| 1 |
\section{\label{integrate}Integrating the Equations of Motion: the Symplectic Step Integrator} |
| 2 |
|
| 3 |
Integration of the equations of motion was carried out using the |
| 4 |
symplectic splitting method proposed by Dullweber \emph{et |
| 5 |
al.}.\cite{Dullweber1997} The reason for this integrator selection |
| 6 |
deals with poor energy conservation of rigid body systems using |
| 7 |
quaternions. While quaternions work well for orientational motion in |
| 8 |
alternate ensembles, the microcanonical ensemble has a constant energy |
| 9 |
requirement that is quite sensitive to errors in the equations of |
| 10 |
motion. The original implementation of this code utilized quaternions |
| 11 |
for rotational motion propagation; however, a detailed investigation |
| 12 |
showed that they resulted in a steady drift in the total energy, |
| 13 |
something that has been observed by others.\cite{Laird97} |
| 14 |
|
| 15 |
The key difference in the integration method proposed by Dullweber |
| 16 |
\emph{et al.} is that the entire rotation matrix is propagated from |
| 17 |
one time step to the next. In the past, this would not have been as |
| 18 |
feasible a option, being that the rotation matrix for a single body is |
| 19 |
nine elements long as opposed to 3 or 4 elements for Euler angles and |
| 20 |
quaternions respectively. System memory has become much less of an |
| 21 |
issue in recent times, and this has resulted in substantial benefits |
| 22 |
in energy conservation. There is still the issue of 5 or 6 additional |
| 23 |
elements for describing the orientation of each particle, which will |
| 24 |
increase dump files substantially. Simply translating the rotation |
| 25 |
matrix into its component Euler angles or quaternions for storage |
| 26 |
purposes relieves this burden. |
| 27 |
|
| 28 |
The symplectic splitting method allows for Verlet style integration of |
| 29 |
both linear and angular motion of rigid bodies. In the integration |
| 30 |
method, the orientational propagation involves a sequence of matrix |
| 31 |
evaluations to update the rotation matrix.\cite{Dullweber1997} These |
| 32 |
matrix rotations end up being more costly computationally than the |
| 33 |
simpler arithmetic quaternion propagation. With the same time step, a |
| 34 |
1000 SSD particle simulation shows an average 7\% increase in |
| 35 |
computation time using the symplectic step method in place of |
| 36 |
quaternions. This cost is more than justified when comparing the |
| 37 |
energy conservation of the two methods as illustrated in figure |
| 38 |
\ref{timestep}. |
| 39 |
|
| 40 |
\begin{figure} |
| 41 |
\epsfxsize=6in |
| 42 |
\epsfbox{timeStep.epsi} |
| 43 |
\caption{Energy conservation using quaternion based integration versus |
| 44 |
the symplectic step method proposed by Dullweber \emph{et al.} with |
| 45 |
increasing time step. For each time step, the dotted line is total |
| 46 |
energy using the symplectic step integrator, and the solid line comes |
| 47 |
from the quaternion integrator. The larger time step plots are shifted |
| 48 |
up from the true energy baseline for clarity.} |
| 49 |
\label{timestep} |
| 50 |
\end{figure} |
| 51 |
|
| 52 |
In figure \ref{timestep}, the resulting energy drift at various time |
| 53 |
steps for both the symplectic step and quaternion integration schemes |
| 54 |
is compared. All of the 1000 SSD particle simulations started with the |
| 55 |
same configuration, and the only difference was the method for |
| 56 |
handling rotational motion. At time steps of 0.1 and 0.5 fs, both |
| 57 |
methods for propagating particle rotation conserve energy fairly well, |
| 58 |
with the quaternion method showing a slight energy drift over time in |
| 59 |
the 0.5 fs time step simulation. At time steps of 1 and 2 fs, the |
| 60 |
energy conservation benefits of the symplectic step method are clearly |
| 61 |
demonstrated. Thus, while maintaining the same degree of energy |
| 62 |
conservation, one can take considerably longer time steps, leading to |
| 63 |
an overall reduction in computation time. |
| 64 |
|
| 65 |
Energy drift in these SSD particle simulations was unnoticeable for |
| 66 |
time steps up to three femtoseconds. A slight energy drift on the |
| 67 |
order of 0.012 kcal/mol per nanosecond was observed at a time step of |
| 68 |
four femtoseconds, and as expected, this drift increases dramatically |
| 69 |
with increasing time step. To insure accuracy in the constant energy |
| 70 |
simulations, time steps were set at 2 fs and kept at this value for |
| 71 |
constant pressure simulations as well. |
