  We present general arguments and construct a stress tensor operator for
finite lattice spin models. The average value of this operator gives the
Casimir force of the system close to the bulk critical temperature $T_c$. We
verify our arguments via exact results for the force in the two-dimensional
Ising model, $d$-dimensional Gaussian and mean spherical model with $2<d<4$. On
the basis of these exact results and by Monte Carlo simulations for
three-dimensional Ising, XY and Heisenberg models we demonstrate that the
standard deviation of the Casimir force $F_C$ in a slab geometry confining a
critical substance in-between is $k_b T D(T)(A/a^{d-1})^{1/2}$, where $A$ is
the surface area of the plates, $a$ is the lattice spacing and $D(T)$ is a
slowly varying nonuniversal function of the temperature $T$. The numerical
calculations demonstrate that at the critical temperature $T_c$ the force
possesses a Gaussian distribution centered at the mean value of the force
$<F_C>=k_b T_c (d-1)\Delta/(L/a)^{d}$, where $L$ is the distance between the
plates and $\Delta$ is the (universal) Casimir amplitude.
