Rigid-body molecular dynamics simulations typically are performed in a
quaternion representation. The nonseparable form of the Hamiltonian i
n quaternions prevents the use of a standard leapfrog (Verlet) integra
tor, so nonsymplectic Runge-Kutta, multistep, or extrapolation methods
are generally used, This is unfortunate since symplectic methods like
Verlet exhibit superior energy conservation in long-time integrations
. In this article, we describe an alternative method, which we call RS
HAKE (for rotation-SHAKE), in which the entire rotation matrix is evol
ved (using the scheme of McLachlan and Scovel [J. Nonlin. Sci, 16 233
(1995)]) in tandem with the particle positions. We employ a fast appro
ximate Newton solver to preserve the orthogonality of the rotation mat
rix. We test our method on a system of soft-sphere dipoles and compare
with quaternion evolution using a 4th-order predictor-corrector integ
rator, Although the short-time error of the quaternion algorithm is sm
aller for fixed time step than that for RSHAKE, the quaternion scheme
exhibits an energy drift which is not observed in simulations with RSH
AKE, hence a fixed energy tolerance can be achieved by using a larger
time step, The superiority of RSHAKE increases with system size. (C) 1
997 American Institute of Physics.