A symplectic method for rigid-body molecular simulation

Abstract

Rigid-body molecular dynamics simulations typically are performed in a quaternion representation. The nonseparable form of the Hamiltonian in quaternions prevents the use of a standard leapfrog ͑Verlet͒ integrator, 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 RSHAKE ͑for rotation-SHAKE͒, in which the entire rotation matrix is evolved ͑using the scheme of McLachlan and Scovel ͓J. Nonlin. Sci. 16 233 ͑1995͔͒͒ in tandem with the particle positions. We employ a fast approximate Newton solver to preserve the orthogonality of the rotation matrix. We test our method on a system of soft-sphere dipoles and compare with quaternion evolution using a 4th-order predictor-corrector integrator. Although the short-time error of the quaternion algorithm is smaller for fixed time step than that for RSHAKE, the quaternion scheme exhibits an energy drift which is not observed in simulations with RSHAKE, hence a fixed energy tolerance can be achieved by using a larger time step. The superiority of RSHAKE increases with system size.

Implementation of Particle-scale Rotation in the 3-D Lattice Solid Model

“…In recent years, some more advanced and complicated algorithms have been proposed which focus on conserving the symplectic structure, time reversible properties and energy (OMELYAN, 1998a, b;DULLWEBER et al, 1997;KOL et al, 1997).…”

Symplectic quaternion scheme for biophysical molecular dynamics

“…However, existing methods designed to treat sets of rigid bodies are either ͑a͒ solved iteratively and are, hence, not reversible 6,7 and cannot be used in hybrid Monte Carlo ͑HMC͒ calculations; ͑b͒ are not symplectic and, hence, not stable at long times; 8,9 or ͑c͒ introduce many extra parameters that must themselves be constrained and, hence, go beyond the four parameters required to define a nonsingular mapping of rigid body rotations. 10,11 Finally, using modern methods, phase space volume preserving but nonsymplectic integrators have been developed 12 but theoretical tools to formally assess the utility of this and similar approaches have been lacking. ͑Note, a symplectic integrator both possesses a time-step-dependent energy, which is invariant along the trajectory produced by the integrator, and is phase space volume preserving.…”

A new algorithm to model the dynamics of 3-D bonded rigid bodies with rotations

“…Different numerical algorithms are employed to discover the optimal approach to integrate the rotational equations, such as high-order Gear methods [1], leapfroglike algorithms [10] and more recently symplectic splitting method [15]. However, previous studies have not addressed the third equally important issue which is the calculation of interactions (or potential) between bonded bodies caused by the relative rotation between bodies.…”

“…For example, it is found that performing a straightforward parameterization of the orientational degrees of freedom of a 3-D rigid body, using a Euler angle representation, is not numerically efficient because of a singularity inherent in this description [8,9]. To avoid this singularity, orientations are typically expressed in terms of a unit quaternion [8,9,14,21], rotation matrices [15,22] and vectors [3,4,19].…”

“…In recently years, more advanced and accurate algorithms have been proposed which focus on conserving the symplectic structure, time reversible and energy [4,7,15,20,21].…”