  We propose new algorithms for numerical integration of the equations of
motion for classical spin systems with fixed spatial site positions. The
algorithms are derived on the basis of a mid-point scheme in conjunction with
the multiple time staging propagation. Contrary to existing predictor-corrector
and decomposition approaches, the algorithms introduced preserve all the
integrals of motion inherent in the basic equations. As is demonstrated for a
lattice ferromagnet model, the present approach appears to be more efficient
even over the recently developed decomposition method.
