  A new scheme for numerical integration of motion for classical systems
composed of rigid polyatomic molecules is proposed. The scheme is based on a
matrix representation of the rotational degrees of freedom. The equations of
motion are integrated within the Verlet framework in velocity form. It is shown
that, contrary to previous methods, in the approach introduced the rigidity of
molecules can be conserved automatically without any additional
transformations. A comparison of various techniques with respect to numerical
stability is made.
