  A revised version of the quaternion approach for numerical integration of the
equations of motion for rigid polyatomic molecules is proposed. The modified
approach is based on a formulation of the quaternion dynamics with constraints.
This allows to resolve the rigidity problem rigorously using constraint forces.
It is shown that the procedure for preservation of molecular rigidity can be
realized particularly simply within the Verlet algorithm in velocity form. We
demonstrate that the presented method leads to an improved numerical stability
with respect to the usual quaternion rescaling scheme and it is roughly as good
as the cumbersome atomic-constraint technique.
