• About OOPSE
• OOPSE News
• Features
• Download
• CVS
• Mailing Lists
• Documentation
• Authors
• Funding
• Results

Integrating the Equations of Motion: the DLM method

The default method for integrating the equations of motion in OOPSE is a velocity-Verlet version of the symplectic splitting method proposed by Dullweber, Leimkuhler and McLachlan (DLM).[42] When there are no directional atoms or rigid bodies present in the simulation, this integrator becomes the standard velocity-Verlet integrator which is known to sample the microcanonical (NVE) ensemble.[43]

Previous integration methods for orientational motion have problems that are avoided in the DLM method. Direct propagation of the Euler angles has a known $1/\sin\theta$ divergence in the equations of motion for $\phi$ and $\psi$ ,[12] leading to numerical instabilities any time one of the directional atoms or rigid bodies has an orientation near $\theta=0$ or $\theta=\pi$ . Quaternion-based integration methods work well for propagating orientational motion; however, energy conservation concerns arise when using the microcanonical (NVE) ensemble. An earlier implementation of OOPSE utilized quaternions for propagation of rotational motion; however, a detailed investigation showed that they resulted in a steady drift in the total energy, something that has been observed by Laird et al.[44]

The key difference in the integration method proposed by Dullweber et al. is that the entire $3 \times 3$ rotation matrix is propagated from one time step to the next. In the past, this would not have been feasible, since the rotation matrix for a single body has nine elements compared with the more memory-efficient methods (using three Euler angles or 4 quaternions). Computer memory has become much less costly in recent years, and this can be translated into substantial benefits in energy conservation.

The basic equations of motion being integrated are derived from the Hamiltonian for conservative systems containing rigid bodies,

$$H = \sum_{i} \left( \frac{1}{2} m_i {\bf v}_i^T \cdot {\bf v}_i +... ...j}_i \right) + V\left(\left\{{\bf r}\right\}, \left\{\mathsf{A}\right\}\right),$$ (4.1)

where ${\bf r}_i$ and ${\bf v}_i$ are the cartesian position vector and velocity of the center of mass of particle $i$ , and ${\bf j}_i$ , $\overleftrightarrow{\mathsf{I}}_i$ are the body-fixed angular momentum and moment of inertia tensor respectively, and the superscript $T$ denotes the transpose of the vector. $\mathsf{A}_i$ is the $3 \times 3$ rotation matrix describing the instantaneous orientation of the particle. $V$ is the potential energy function which may depend on both the positions $\left\{{\bf r}\right\}$ and orientations $\left\{\mathsf{A}\right\}$ of all particles. The equations of motion for the particle centers of mass are derived from Hamilton's equations and are quite simple,

$$\dot{{\bf r}} = {\bf v},$$ (4.2)

$$\dot{{\bf v}} = \frac{{\bf f}}{m},$$ (4.3)

where ${\bf f}$ is the instantaneous force on the center of mass of the particle,

$${\bf f} = - \frac{\partial}{\partial {\bf r}} V(\left\{{\bf r}(t)\right\}, \left\{\mathsf{A}(t)\right\}).$$ (4.4)

The equations of motion for the orientational degrees of freedom are

$$\dot{\mathsf{A}} = \mathsf{A} \cdot \left(\overleftrightarrow{\mathsf{I}}^{-1} \cdot {\bf j}\right),$$ (4.5)

$$\dot{{\bf j}} = {\bf j} \times \left( \overleftrightarrow{\mathsf{I}}^{-1} \cdot {\bf j} \right) - \left(\mathsf{A}^{T} \cdot \frac{\partial V}{\partial \mathsf{A}} \right).$$ (4.6)

In these equations of motion, the skew matrix of a vector ${\bf v} = \left( v_1, v_2, v_3 \right)$ is defined:

$$\left( {\bf v} \right) := \left( \begin{array}{ccc} 0 & v_3 & - v_2 \\ -v_3 & 0 & v_1 \\ v_2 & -v_1 & 0 \end{array} \right).$$ (4.7)

The rot notation refers to the mapping of the $3 \times 3$ rotation matrix to a vector of orientations by first computing the skew-symmetric part $\left(\mathsf{A} - \mathsf{A}^{T}\right)$ and then associating this with a length 3 vector by inverting the skew function above:

$$\left(\mathsf{A}\right) := ^{-1}\left(\mathsf{A} - \mathsf{A}^{T} \right).$$ (4.8)

Written this way, the rot operation creates a set of conjugate angle coordinates to the body-fixed angular momenta represented by ${\bf j}$ . This equation of motion for angular momenta is equivalent to the more familiar body-fixed forms,

$$\dot{j_{x}} = \tau^b_x(t) - \left(\overleftrightarrow{\mathsf{I}}_{yy}^{-1} - \overleftrightarrow{\mathsf{I}}_{zz}^{-1} \right) j_y j_z,$$ (4.9)

$$\dot{j_{y}} = \tau^b_y(t) - \left(\overleftrightarrow{\mathsf{I}}_{zz}^{-1} - \overleftrightarrow{\mathsf{I}}_{xx}^{-1} \right) j_z j_x,$$ (4.10)

$$\dot{j_{z}} = \tau^b_z(t) - \left(\overleftrightarrow{\mathsf{I}}_{xx}^{-1} - \overleftrightarrow{\mathsf{I}}_{yy}^{-1} \right) j_x j_y,$$ (4.11)

which utilize the body-fixed torques, ${\bf\tau}^b$ . Torques are most easily derived in the space-fixed frame,

$${\bf\tau}^b(t) = \mathsf{A}(t) \cdot {\bf\tau}^s(t),$$ (4.12)

where the torques are either derived from the forces on the constituent atoms of the rigid body, or for directional atoms, directly from derivatives of the potential energy,

$${\bf\tau}^s(t) = - \hat{\bf u}(t) \times \left( \frac{\partial} {... ...eft(\left\{ {\bf r}(t) \right\}, \left\{ \mathsf{A}(t) \right\}\right) \right).$$ (4.13)

Here $\hat{\bf u}$ is a unit vector pointing along the principal axis of the particle in the space-fixed frame.

The DLM method uses a Trotter factorization of the orientational propagator. This has three effects:

1. the integrator is area-preserving in phase space (i.e. it is symplectic),
2. the integrator is time-reversible, making it suitable for Hybrid Monte Carlo applications, and
3. the error for a single time step is of order $\mathcal{O}\left(h^4\right)$ for timesteps of length $h$ .

The integration of the equations of motion is carried out in a velocity-Verlet style 2-part algorithm, where $h= \delta t$ :

moveA:

$${\bf v}\left(t + h / 2\right) \leftarrow {\bf v}(t) + \frac{h}{2} \left( {\bf f}(t) / m \right),$$

$${\bf r}(t + h) \leftarrow {\bf r}(t) + h {\bf v}\left(t + h / 2 \right),$$

$${\bf j}\left(t + h / 2 \right) \leftarrow {\bf j}(t) + \frac{h}{2} {\bf\tau}^b(t),$$

$$\mathsf{A}(t + h) \leftarrow \mathrm{rotate}\left( h {\bf j} (t + h / 2) \cdot \overleftrightarrow{\mathsf{I}}^{-1} \right).$$

In this context, the $\mathrm{rotate}$ function is the reversible product of the three body-fixed rotations,

$$\mathrm{rotate}({\bf a}) = \mathsf{G}_x(a_x / 2) \cdot \mathsf{G}... ...\cdot \mathsf{G}_z(a_z) \cdot \mathsf{G}_y(a_y / 2) \cdot \mathsf{G}_x(a_x /2),$$ (4.14)

where each rotational propagator, $\mathsf{G}_\alpha(\theta)$ , rotates both the rotation matrix ( $\mathsf{A}$ ) and the body-fixed angular momentum (${\bf j}$ ) by an angle $\theta$ around body-fixed axis $\alpha$ ,

$$\mathsf{G}_\alpha( \theta ) = \left\{ \begin{array}{lcl} \mathsf{... ... & \leftarrow & \mathsf{R}_\alpha(\theta) \cdot {\bf j}(0). \end{array} \right.$$ (4.15)

$\mathsf{R}_\alpha$ is a quadratic approximation to the single-axis rotation matrix. For example, in the small-angle limit, the rotation matrix around the body-fixed x-axis can be approximated as

$$\mathsf{R}_x(\theta) \approx \left( \begin{array}{ccc} 1 & 0 & 0 ... ...1+ \theta^2 / 4} & \frac{1-\theta^2 / 4}{1 + \theta^2 / 4} \end{array} \right).$$ (4.16)

All other rotations follow in a straightforward manner.

After the first part of the propagation, the forces and body-fixed torques are calculated at the new positions and orientations

doForces:

$${\bf f}(t + h) \leftarrow - \left(\frac{\partial V}{\partial {\bf r}}\right)_{{\bf r}(t + h)},$$

$${\bf\tau}^{s}(t + h) \leftarrow {\bf u}(t + h) \times \frac{\partial V}{\partial {\bf u}},$$

$${\bf\tau}^{b}(t + h) \leftarrow \mathsf{A}(t + h) \cdot {\bf\tau}^s(t + h).$$

OOPSE automatically updates ${\bf u}$ when the rotation matrix $\mathsf{A}$ is calculated in moveA. Once the forces and torques have been obtained at the new time step, the velocities can be advanced to the same time value.

moveB:

$${\bf v}\left(t + h \right) \leftarrow {\bf v}\left(t + h / 2 \right) + \frac{h}{2} \left( {\bf f}(t + h) / m \right),$$

$${\bf j}\left(t + h \right) \leftarrow {\bf j}\left(t + h / 2 \right) + \frac{h}{2} {\bf\tau}^b(t + h) .$$

The matrix rotations used in the DLM method end up being more costly computationally than the simpler arithmetic quaternion propagation. With the same time step, a 1024-molecule water simulation incurs an average 12% increase in computation time using the DLM method in place of quaternions. This cost is more than justified when comparing the energy conservation achieved by the two methods. Figure 4.1 provides a comparative analysis of the DLM method versus the traditional quaternion scheme.

Figure 4.1: Analysis of the energy conservation of the DLM and quaternion integration methods. $\delta \mathrm{E}_1$ is the linear drift in energy over time and $\delta \mathrm{E}_0$ is the standard deviation of energy fluctuations around this drift. All simulations were of a 1024-molecule simulation of SSD water at 298 K starting from the same initial configuration. Note that the DLM method provides more than an order of magnitude improvement in both the energy drift and the size of the energy fluctuations when compared with the quaternion method at any given time step. At time steps larger than 4 fs, the quaternion scheme resulted in rapidly rising energies which eventually lead to simulation failure. Using the DLM method, time steps up to 8 fs can be taken before this behavior is evident.

In Fig. 4.1, $\delta$   E$_1$ is a measure of the linear energy drift in units of kcal mol$^{-1}$ per particle over a nanosecond of simulation time, and $\delta$   E$_0$ is the standard deviation of the energy fluctuations in units of kcal mol$^{-1}$ per particle. In the top plot, it is apparent that the energy drift is reduced by a significant amount (2 to 3 orders of magnitude improvement at all tested time steps) by chosing the DLM method over the simple non-symplectic quaternion integration method. In addition to this improvement in energy drift, the fluctuations in the total energy are also dampened by 1 to 2 orders of magnitude by utilizing the DLM method.

Although the DLM method is more computationally expensive than the traditional quaternion scheme for propagating a single time step, consideration of the computational cost for a long simulation with a particular level of energy conservation is in order. A plot of energy drift versus computational cost was generated (Fig. 4.2). This figure provides an estimate of the CPU time required under the two integration schemes for 1 nanosecond of simulation time for the model 1024-molecule system. By chosing a desired energy drift value it is possible to determine the CPU time required for both methods. If a $\delta$   E$_1$ of $1 \times 10^{-3}$   kcal mol$^{-1}$ per particle is desired, a nanosecond of simulation time will require 19 hours of CPU time with the DLM integrator, while the quaternion scheme will require 154 hours of CPU time. This demonstrates the computational advantage of the integration scheme utilized in OOPSE.

Figure 4.2: Energy drift as a function of required simulation run time. $\delta \mathrm{E}_1$ is the linear drift in energy over time. Simulations were performed on a single 2.5 GHz Pentium 4 processor. Simulation time comparisons can be made by tracing horizontally from one curve to the other. For example, a simulation that takes 24 hours using the DLM method will take roughly 210 hours using the simple quaternion method if the same degree of energy conservation is desired.

There is only one specific keyword relevant to the default integrator, and that is the time step for integrating the equations of motion.

variable	Meta-data keyword	units	default value
$h$	dt = 2.0;	fs	none

Copyright © 2006 - University of Notre Dame

Updated on January 16, 2006