trunk/oopsePaper/Mechanics.tex


\section{\label{sec:mechanics}Mechanics}

\subsection{\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.


\subsection{\label{sec:extended}Extended Systems for other Ensembles}


{\sc oopse} implements a 


\subsubsection{\label{sec:noseHooverThermo}Nose-Hoover Thermostatting}

To mimic the effects of being in a constant temperature ({\sc nvt})
ensemble, {\sc oopse} uses the Nose-Hoover extended system
approach.\cite{Hoover85} In this method, the equations of motion for
the particle positions and velocities are
\begin{eqnarray}
\dot{{\bf r}} & = & {\bf v} \\
\dot{{\bf v}} & = & \frac{{\bf f}}{m} - \chi {\bf v}
\label{eq:nosehoovereom}
\end{eqnarray}

$\chi$ is an ``extra'' variable included in the extended system, and
it is propagated using the first order equation of motion
\begin{equation}
\dot{\chi} = \frac{1}{\tau_{T}} \left( \frac{T}{T_{target}} - 1 \right)
\label{eq:nosehooverext}
\end{equation}
where $T_{target}$ is the target temperature for the simulation, and
$\tau_{T}$ is a time constant for the thermostat.  

To select the Nose-Hoover {\sc nvt} ensemble, the {\tt ensemble = NVT;} 
command would be used in the simulation's {\sc bass} file.  There is
some subtlety in choosing values for $\tau_{T}$, and it is usually set
to values of a few ps.  Within a {\sc bass} file, $\tau_{T}$ could be
set to 1 ps using the {\tt tauThermostat = 1000; } command.


\subsection{\label{Sec:zcons}Z-Constraint Method}

Based on fluctuatin-dissipation theorem,\bigskip\ force auto-correlation
method was developed to investigate the dynamics of ions inside the ion
channels.\cite{Roux91} Time-dependent friction coefficient can be calculated
from the deviation of the instaneous force from its mean force.

%

\begin{equation}
\xi(z,t)=\langle\delta F(z,t)\delta F(z,0)\rangle/k_{B}T
\end{equation}


where%
\begin{equation}
\delta F(z,t)=F(z,t)-\langle F(z,t)\rangle
\end{equation}


If the time-dependent friction decay rapidly, static friction coefficient can
be approximated by%

\begin{equation}
\xi^{static}(z)=\int_{0}^{\infty}\langle\delta F(z,t)\delta F(z,0)\rangle dt
\end{equation}


Hence, diffusion constant can be estimated by
\begin{equation}
D(z)=\frac{k_{B}T}{\xi^{static}(z)}=\frac{(k_{B}T)^{2}}{\int_{0}^{\infty
}\langle\delta F(z,t)\delta F(z,0)\rangle dt}%
\end{equation}


\bigskip Z-Constraint method, which fixed the z coordinates of the molecules
with respect to the center of the mass of the system, was proposed to obtain
the forces required in force auto-correlation method.\cite{Marrink94} However,
simply resetting the coordinate will move the center of the mass of the whole
system. To avoid this problem,  a new method was used at {\sc oopse}. Instead of
resetting the coordinate, we reset the forces of z-constraint molecules as
well as subtract the total constraint forces from the rest of the system after
force calculation at each time step. 
\begin{verbatim}
$F_{\alpha i}=0$
$V_{\alpha i}=V_{\alpha i}-\frac{\sum\limits_{i}M_{_{\alpha i}}V_{\alpha i}}{\sum\limits_{i}M_{_{\alpha i}}}$
$F_{\alpha i}=F_{\alpha i}-\frac{M_{_{\alpha i}}}{\sum\limits_{\alpha}\sum\limits_{i}M_{_{\alpha i}}}\sum\limits_{\beta}F_{\beta}$
$V_{\alpha i}=V_{\alpha i}-\frac{\sum\limits_{\alpha}\sum\limits_{i}M_{_{\alpha i}}V_{\alpha i}}{\sum\limits_{\alpha}\sum\limits_{i}M_{_{\alpha i}}}$
\end{verbatim}

At the very beginning of the simulation, the molecules may not be at its
constraint position. To move the z-constraint molecule to the specified
position, a simple harmonic potential is used%

\begin{equation}
U(t)=\frac{1}{2}k_{Harmonic}(z(t)-z_{cons})^{2}%
\end{equation}
where $k_{Harmonic}$\bigskip\ is the harmonic force constant, $z(t)$ is
current z coordinate of the center of mass of the z-constraint molecule, and
$z_{cons}$ is the restraint position. Therefore, the harmonic force operated
on the z-constraint molecule at time $t$ can be calculated by%
\begin{equation}
F_{z_{Harmonic}}(t)=-\frac{\partial U(t)}{\partial z(t)}=-k_{Harmonic}%
(z(t)-z_{cons})
\end{equation}
Worthy of mention, other kinds of potential functions can also be used to
drive the z-constraint molecule.
Revision:	899
Committed:	Tue Jan 6 18:53:58 2004 UTC (21 years, 9 months ago) by mmeineke
Content type:	application/x-tex
File size:	8074 byte(s)
Log Message:	reworked the directory structure to reflect the Outline of the paper. All major sections now have their own tex file
#	User	Rev	Content
1	mmeineke	899
2			\section{\label{sec:mechanics}Mechanics}
3
4			\subsection{\label{integrate}Integrating the Equations of Motion: the Symplectic Step Integrator}
5
6			Integration of the equations of motion was carried out using the
7			symplectic splitting method proposed by Dullweber \emph{et
8			al.}.\cite{Dullweber1997} The reason for this integrator selection
9			deals with poor energy conservation of rigid body systems using
10			quaternions. While quaternions work well for orientational motion in
11			alternate ensembles, the microcanonical ensemble has a constant energy
12			requirement that is quite sensitive to errors in the equations of
13			motion. The original implementation of this code utilized quaternions
14			for rotational motion propagation; however, a detailed investigation
15			showed that they resulted in a steady drift in the total energy,
16			something that has been observed by others.\cite{Laird97}
17
18			The key difference in the integration method proposed by Dullweber
19			\emph{et al.} is that the entire rotation matrix is propagated from
20			one time step to the next. In the past, this would not have been as
21			feasible a option, being that the rotation matrix for a single body is
22			nine elements long as opposed to 3 or 4 elements for Euler angles and
23			quaternions respectively. System memory has become much less of an
24			issue in recent times, and this has resulted in substantial benefits
25			in energy conservation. There is still the issue of 5 or 6 additional
26			elements for describing the orientation of each particle, which will
27			increase dump files substantially. Simply translating the rotation
28			matrix into its component Euler angles or quaternions for storage
29			purposes relieves this burden.
30
31			The symplectic splitting method allows for Verlet style integration of
32			both linear and angular motion of rigid bodies. In the integration
33			method, the orientational propagation involves a sequence of matrix
34			evaluations to update the rotation matrix.\cite{Dullweber1997} These
35			matrix rotations end up being more costly computationally than the
36			simpler arithmetic quaternion propagation. With the same time step, a
37			1000 SSD particle simulation shows an average 7\% increase in
38			computation time using the symplectic step method in place of
39			quaternions. This cost is more than justified when comparing the
40			energy conservation of the two methods as illustrated in figure
41			\ref{timestep}.
42
43			\begin{figure}
44			\epsfxsize=6in
45			\epsfbox{timeStep.epsi}
46			\caption{Energy conservation using quaternion based integration versus
47			the symplectic step method proposed by Dullweber \emph{et al.} with
48			increasing time step. For each time step, the dotted line is total
49			energy using the symplectic step integrator, and the solid line comes
50			from the quaternion integrator. The larger time step plots are shifted
51			up from the true energy baseline for clarity.}
52			\label{timestep}
53			\end{figure}
54
55			In figure \ref{timestep}, the resulting energy drift at various time
56			steps for both the symplectic step and quaternion integration schemes
57			is compared. All of the 1000 SSD particle simulations started with the
58			same configuration, and the only difference was the method for
59			handling rotational motion. At time steps of 0.1 and 0.5 fs, both
60			methods for propagating particle rotation conserve energy fairly well,
61			with the quaternion method showing a slight energy drift over time in
62			the 0.5 fs time step simulation. At time steps of 1 and 2 fs, the
63			energy conservation benefits of the symplectic step method are clearly
64			demonstrated. Thus, while maintaining the same degree of energy
65			conservation, one can take considerably longer time steps, leading to
66			an overall reduction in computation time.
67
68			Energy drift in these SSD particle simulations was unnoticeable for
69			time steps up to three femtoseconds. A slight energy drift on the
70			order of 0.012 kcal/mol per nanosecond was observed at a time step of
71			four femtoseconds, and as expected, this drift increases dramatically
72			with increasing time step. To insure accuracy in the constant energy
73			simulations, time steps were set at 2 fs and kept at this value for
74			constant pressure simulations as well.
75
76
77			\subsection{\label{sec:extended}Extended Systems for other Ensembles}
78
79
80			{\sc oopse} implements a
81
82
83			\subsubsection{\label{sec:noseHooverThermo}Nose-Hoover Thermostatting}
84
85			To mimic the effects of being in a constant temperature ({\sc nvt})
86			ensemble, {\sc oopse} uses the Nose-Hoover extended system
87			approach.\cite{Hoover85} In this method, the equations of motion for
88			the particle positions and velocities are
89			\begin{eqnarray}
90			\dot{{\bf r}} & = & {\bf v} \\
91			\dot{{\bf v}} & = & \frac{{\bf f}}{m} - \chi {\bf v}
92			\label{eq:nosehoovereom}
93			\end{eqnarray}
94
95			$\chi$ is an ``extra'' variable included in the extended system, and
96			it is propagated using the first order equation of motion
97			\begin{equation}
98			\dot{\chi} = \frac{1}{\tau_{T}} \left( \frac{T}{T_{target}} - 1 \right)
99			\label{eq:nosehooverext}
100			\end{equation}
101			where $T_{target}$ is the target temperature for the simulation, and
102			$\tau_{T}$ is a time constant for the thermostat.
103
104			To select the Nose-Hoover {\sc nvt} ensemble, the {\tt ensemble = NVT;}
105			command would be used in the simulation's {\sc bass} file. There is
106			some subtlety in choosing values for $\tau_{T}$, and it is usually set
107			to values of a few ps. Within a {\sc bass} file, $\tau_{T}$ could be
108			set to 1 ps using the {\tt tauThermostat = 1000; } command.
109
110
111			\subsection{\label{Sec:zcons}Z-Constraint Method}
112
113			Based on fluctuatin-dissipation theorem,\bigskip\ force auto-correlation
114			method was developed to investigate the dynamics of ions inside the ion
115			channels.\cite{Roux91} Time-dependent friction coefficient can be calculated
116			from the deviation of the instaneous force from its mean force.
117
118			%
119
120			\begin{equation}
121			\xi(z,t)=\langle\delta F(z,t)\delta F(z,0)\rangle/k_{B}T
122			\end{equation}
123
124
125			where%
126			\begin{equation}
127			\delta F(z,t)=F(z,t)-\langle F(z,t)\rangle
128			\end{equation}
129
130
131			If the time-dependent friction decay rapidly, static friction coefficient can
132			be approximated by%
133
134			\begin{equation}
135			\xi^{static}(z)=\int_{0}^{\infty}\langle\delta F(z,t)\delta F(z,0)\rangle dt
136			\end{equation}
137
138
139			Hence, diffusion constant can be estimated by
140			\begin{equation}
141			D(z)=\frac{k_{B}T}{\xi^{static}(z)}=\frac{(k_{B}T)^{2}}{\int_{0}^{\infty
142			}\langle\delta F(z,t)\delta F(z,0)\rangle dt}%
143			\end{equation}
144
145
146			\bigskip Z-Constraint method, which fixed the z coordinates of the molecules
147			with respect to the center of the mass of the system, was proposed to obtain
148			the forces required in force auto-correlation method.\cite{Marrink94} However,
149			simply resetting the coordinate will move the center of the mass of the whole
150			system. To avoid this problem, a new method was used at {\sc oopse}. Instead of
151			resetting the coordinate, we reset the forces of z-constraint molecules as
152			well as subtract the total constraint forces from the rest of the system after
153			force calculation at each time step.
154			\begin{verbatim}
155			$F_{\alpha i}=0$
156			$V_{\alpha i}=V_{\alpha i}-\frac{\sum\limits_{i}M_{_{\alpha i}}V_{\alpha i}}{\sum\limits_{i}M_{_{\alpha i}}}$
157			$F_{\alpha i}=F_{\alpha i}-\frac{M_{_{\alpha i}}}{\sum\limits_{\alpha}\sum\limits_{i}M_{_{\alpha i}}}\sum\limits_{\beta}F_{\beta}$
158			$V_{\alpha i}=V_{\alpha i}-\frac{\sum\limits_{\alpha}\sum\limits_{i}M_{_{\alpha i}}V_{\alpha i}}{\sum\limits_{\alpha}\sum\limits_{i}M_{_{\alpha i}}}$
159			\end{verbatim}
160
161			At the very beginning of the simulation, the molecules may not be at its
162			constraint position. To move the z-constraint molecule to the specified
163			position, a simple harmonic potential is used%
164
165			\begin{equation}
166			U(t)=\frac{1}{2}k_{Harmonic}(z(t)-z_{cons})^{2}%
167			\end{equation}
168			where $k_{Harmonic}$\bigskip\ is the harmonic force constant, $z(t)$ is
169			current z coordinate of the center of mass of the z-constraint molecule, and
170			$z_{cons}$ is the restraint position. Therefore, the harmonic force operated
171			on the z-constraint molecule at time $t$ can be calculated by%
172			\begin{equation}
173			F_{z_{Harmonic}}(t)=-\frac{\partial U(t)}{\partial z(t)}=-k_{Harmonic}%
174			(z(t)-z_{cons})
175			\end{equation}
176			Worthy of mention, other kinds of potential functions can also be used to
177			drive the z-constraint molecule.