trunk/oopsePaper/integrator.tex

\section{\label{integrate}Integrating the Equations of Motion: the Symplectic Step Integrator}

Integration of the equations of motion was carried out using the
symplectic splitting method proposed by Dullweber \emph{et
al.}.\cite{Dullweber1997} The reason for this integrator selection
deals with poor energy conservation of rigid body systems using
quaternions. While quaternions work well for orientational motion in
alternate ensembles, the microcanonical ensemble has a constant energy
requirement that is quite sensitive to errors in the equations of
motion. The original implementation of this code utilized quaternions
for rotational motion propagation; however, a detailed investigation
showed that they resulted in a steady drift in the total energy,
something that has been observed by others.\cite{Laird97}

The key difference in the integration method proposed by Dullweber
\emph{et al.} is that the entire rotation matrix is propagated from
one time step to the next. In the past, this would not have been as
feasible a option, being that the rotation matrix for a single body is
nine elements long as opposed to 3 or 4 elements for Euler angles and
quaternions respectively. System memory has become much less of an
issue in recent times, and this has resulted in substantial benefits
in energy conservation. There is still the issue of 5 or 6 additional
elements for describing the orientation of each particle, which will
increase dump files substantially. Simply translating the rotation
matrix into its component Euler angles or quaternions for storage
purposes relieves this burden.

The symplectic splitting method allows for Verlet style integration of
both linear and angular motion of rigid bodies. In the integration
method, the orientational propagation involves a sequence of matrix
evaluations to update the rotation matrix.\cite{Dullweber1997} These
matrix rotations end up being more costly computationally than the
simpler arithmetic quaternion propagation. With the same time step, a
1000 SSD particle simulation shows an average 7\% increase in
computation time using the symplectic step method in place of
quaternions. This cost is more than justified when comparing the
energy conservation of the two methods as illustrated in figure
\ref{timestep}.

\begin{figure}
\epsfxsize=6in
\epsfbox{timeStep.epsi}
\caption{Energy conservation using quaternion based integration versus 
the symplectic step method proposed by Dullweber \emph{et al.} with
increasing time step. For each time step, the dotted line is total
energy using the symplectic step integrator, and the solid line comes
from the quaternion integrator. The larger time step plots are shifted
up from the true energy baseline for clarity.}
\label{timestep}
\end{figure}

In figure \ref{timestep}, the resulting energy drift at various time
steps for both the symplectic step and quaternion integration schemes
is compared. All of the 1000 SSD particle simulations started with the
same configuration, and the only difference was the method for
handling rotational motion. At time steps of 0.1 and 0.5 fs, both
methods for propagating particle rotation conserve energy fairly well,
with the quaternion method showing a slight energy drift over time in
the 0.5 fs time step simulation. At time steps of 1 and 2 fs, the
energy conservation benefits of the symplectic step method are clearly
demonstrated. Thus, while maintaining the same degree of energy
conservation, one can take considerably longer time steps, leading to
an overall reduction in computation time.

Energy drift in these SSD particle simulations was unnoticeable for
time steps up to three femtoseconds. A slight energy drift on the
order of 0.012 kcal/mol per nanosecond was observed at a time step of
four femtoseconds, and as expected, this drift increases dramatically
with increasing time step. To insure accuracy in the constant energy
simulations, time steps were set at 2 fs and kept at this value for
constant pressure simulations as well.
Revision:	818
Committed:	Fri Oct 24 21:27:59 2003 UTC (21 years, 7 months ago) by gezelter
Content type:	application/x-tex
File size:	4010 byte(s)
Log Message:	Formatting changes to get out of RevTex nonsense
#	User	Rev	Content
1	chrisfen	775	\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	gezelter	818	\epsfxsize=6in
42			\epsfbox{timeStep.epsi}
43	chrisfen	775	\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.