Integration of the equations of motion

The velocity Verlet algorithm [8] was initially introduced to improve the numerical stability of the leap frog scheme [15]. The velocity Verlet algorithm has subsequently been derived in a systematic way by means of a time-reversible partitioning of the Louville operator due to Tuckerman, et al. [9] and is now widely used in simulations. It is an example of a second order symplectic integrator. It has the form

$$x(\delta t) = x(0) + \dot x(0) \delta t + {(\delta t^2) \over 2} a(0)$$ (9)

$$\dot x(\delta t) = \dot x(0) + {\delta t\over 2}[ a(0) + a (\delta t)]$$

where a(0) is the acceleration term evaluated using x(0).

While the systematic derivation for translational degrees of freedom does not apply to rotation of a rigid body, one can still propose a velocity Verlet like algorithm for the quaternions. We adopt here the scheme proposed by Omelyan [10]. The conditions on the quaternions, $Q_\alpha Q_\alpha =1$ and $Q_\alpha \dot Q_\alpha =0$, are incorporated into the coefficient $\Lambda$ of a constraint force with the form $f_\alpha = -2\Lambda Q_\alpha$ so that the integrator for $Q_\alpha$ takes the form

$$Q_\alpha (\delta t) = Q_\alpha (0) + \dot Q_\alpha (0) \del... ...er 2} \ddot Q_\alpha(0) + f_\alpha(0) {(\delta t)^2 \over 2}.$$ (10)

The condition $Q_\alpha (\delta t)Q_\alpha (\delta t) = 1$ leads to an explicit expression for the coefficient $\Lambda$, namely

$$(\delta t)^2 \Lambda = 1-s_1 (\delta t)^2/2 -\sqrt{1 -s_1 (\delta t)^2 -s_2 (\delta t)^3 -(s_3 -s_1^2)(\delta t)^4 /4},$$ (11)

where the s_i terms are sums; $s_1 = \dot Q_\alpha (0)\dot Q_\alpha (0)$, $s_2 = \dot Q_\alpha (0)\ddot Q_\alpha (0)$, and $s_3 = \ddot Q_\alpha (0)\ddot Q_\alpha (0)$. For small $\delta t$, $\Lambda \rightarrow s_2 \delta t /2$.

The updated values for $\dot Q_\alpha(\delta t)$ and $\ddot Q_\alpha(\delta t)$ are obtained using Eqs. (4) and (8) with values for $\omega_{px}$, $\omega_{py}$, and $\omega_{pz}$ at $\delta t$ obtained by solving the Euler equations, Eqs. ( 3). Since $\dot \omega _{p\alpha}$ is proportional to $\omega _{p\beta} \omega _{p\gamma}$, it is necessary to iterate the second member of Eq. (3) in order to obtain a self-consistent result. Here we suggest a scheme that converges rapidly.

First determine the $\omega$-independent part of $\dot \omega _{p\alpha} (\delta t)$which involves just the torques, Eq. (3), and call it $T_\alpha (\delta t)$. A zeroth estimate for $\omega_{p\alpha} (\delta t)$ is then

$$\omega _{p\alpha} ^{(0)} (\delta t) = \omega _{p\alpha} (0) ... ...t \over 2} [\dot \omega_{p\alpha} (0) + T_\alpha (\delta t)].$$ (12)

This estimate for $\omega_{p\alpha} (\delta t)$ is then used to estimate the $\omega$-dependent part of the right-hand side of Eq. (3), say $g_\alpha ^{(0)}[\omega^{(0)}(\delta t)]$. The first estimate for $\omega_{p\alpha} (\delta t)$ is then

$$\omega_{p\alpha}^{(1)} (\delta t) = \omega _{p\alpha}^{(0)}(\delta t) + {\delta t\over 2} g_\alpha ^{(0)}.$$ (13)

Now use $\omega _{p\alpha}^{(1)} (\delta t)$ to construct $g_\alpha ^{(1)}[\omega ^{(1)}(\delta t)]$and then to generate the second estimate for $\omega_{p\alpha} (\delta t)$,

$$\omega _{p\alpha}^{(2)} (\delta t) = \omega _{p\alpha} ^{(0)}(\delta t) + {\delta t\over 2} g_\alpha ^{(1)}.$$ (14)

This process can be continued until the desired level of convergence has been reached. We find that three iterations is sufficient for the examples discussed in the next section. Eqs. (10) through (14) constitute the "constraint force algorithm."

A related algorithm for integrating the equations of motion for quaternions that is patterned after the original Verlet algorithm has been described by Svanberg. [16] His "mid-step implicit algorithm" is similar to a velocity Verlet algorithm that iterates $\dot Q_{\alpha}$and imposes the $Q_\alpha Q_\alpha =1$ and $Q_\alpha \dot Q_\alpha =0$conditions by scaling. Omelyan [10] has shown that a velocity Verlet algorithm with scaling for quaternions is inferior to the version described above. This is illustrated by the example discussed next.