trunk/oopsePaper/Mechanics.tex


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

\subsection{\label{integrate}Integrating the Equations of Motion: the
DLM method}

The default method for integrating the equations of motion in {\sc
oopse} is a velocity-Verlet version of the symplectic splitting method
proposed by Dullweber, Leimkuhler and McLachlan
(DLM).\cite{Dullweber1997} 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.\cite{}

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$,\cite{AllenTildesley} leading to
numerical instabilities any time one of the directional atoms or rigid
bodies has an orientation near $\theta=0$ or $\theta=\pi$.  More
modern quaternion-based integration methods have relatively poor
energy conservation.  While quaternions work well for orientational
motion in other ensembles, the microcanonical ensemble has a
constant energy requirement that is quite sensitive to errors in the
equations of motion.  An earlier implementation of {\sc 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
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
feasible, since that the rotation matrix for a single body has nine
elements compared to 3 or 4 elements for Euler angles and quaternions
respectively. System memory has become less costly in recent times,
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,
\begin{equation}
H = \sum_{i} \frac{{\bf v}_i^T m_i {\bf v}_i}{2} + \sum_i
\frac{{\bf j}_i^T \cdot {\bf j}_i}{2 \overleftrightarrow{\mathsf{I}}_i} +
V\left(\left\{{\bf r}\right\}, \left\{\mathsf{A}\right\}\right)
\end{equation}
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$
and $\overleftrightarrow{\mathsf{I}}_i$ are the body-fixed angular
momentum and moment of inertia tensor, respectively.  $\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,
\begin{eqnarray}
\dot{{\bf r}} & = & {\bf v} \\
\dot{{\bf v}} & = & \frac{{\bf f}}{m}
\end{eqnarray}
where ${\bf f}$ is the instantaneous force on the center of mass
of the particle,
\begin{equation}
{\bf f} = \frac{\partial}{\partial
{\bf r}} V(\left\{{\bf r}(t)\right\}, \left\{\mathsf{A}(t)\right\}).
\end{equation}


The equations of motion for the orientational degrees of freedom are
\begin{eqnarray}
\dot{\mathsf{A}} & = & \mathsf{A}
\mbox{ skew}\left(\overleftrightarrow{I}^{-1} \cdot {\bf j}\right) \\
\dot{{\bf j}} & = & {\bf j} \times \left( \overleftrightarrow{I}^{-1}
\cdot {\bf j} \right) - \mbox{ rot}\left(\mathsf{A}^{T} \frac{\partial
V}{\partial \mathsf{A}} \right)
\end{eqnarray}
In these equations of motion, the $\mbox{skew}$ matrix of a vector
${\bf v} = \left( v_1, v_2, v_3 \right)$ is defined:
\begin{equation}
\mbox{skew}\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)
\end{equation}
The $\mbox{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
$\mbox{skew}$ function above:
\begin{equation}
\mbox{rot}\left(\mathsf{A}\right) := \mbox{ skew}^{-1}\left(\mathsf{A}
- \mathsf{A}^{T} \right)
\end{equation}


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 number of extended system integrators for
sampling from other ensembles relevant to chemical physics.  The
integrator can selected with the {\tt ensemble} keyword in the
{\tt .bass} file:

\begin{center}
\begin{tabular}{lll}
{\bf Integrator} & {\bf Ensemble} & {\bf {\tt .bass} line} \\
NVE & microcanonical & {\tt ensemble = ``NVE''; } \\
NVT & canonical & {\tt ensemble = ``NVT''; } \\
NPTi & isobaric-isothermal (with isotropic volume changes) & {\tt
ensemble = ``NPTi'';} \\
NPTf & isobaric-isothermal (with changes to box shape) & {\tt
ensemble = ``NPTf'';} \\
NPTxyz & approximate isobaric-isothermal & {\tt ensemble =
``NPTxyz'';} \\
 &  (with separate barostats on each box dimension) &
\end{tabular}
\end{center}

The relatively well-known Nos\'e-Hoover thermostat is implemented in
{\sc oopse}'s NVT integrator.  This method couples an extra degree of
freedom (the thermostat) to the kinetic energy of the system, and has
been shown to sample the canonical distribution in the system degrees
of freedom while conserving a quantity that is, to within a constant,
the Helmholtz free energy.

NPT algorithms attempt to maintain constant pressure in the system by
coupling the volume of the system to a barostat.  {\sc oopse} contains
three different constant pressure algorithms.  The first two, NPTi and
NPTf have been shown to conserve a quantity that is, to within a
constant, the Gibbs free energy.  The Melchionna modification to the
Hoover barostat is implemented in both NPTi and NPTf.  NPTi allows
only isotropic changes in the simulation box, while box {\it shape}
variations are allowed in NPTf.  The NPTxyz integrator has {\it not}
been shown to sample from the isobaric-isothermal ensemble.  It is
useful, however, in that it maintains orthogonality for the axes of
the simulation box while attempting to equalize pressure along the
three perpendicular directions in the box.

Each of the extended system integrators requires additional keywords
to set target values for the thermodynamic state variables that are
being held constant.  Keywords are also required to set the
characteristic decay times for the dynamics of the extended
variables.

\begin{tabular}{llll}
{\bf variable} & {\bf {\tt .bass} keyword} & {\bf units} & {\bf
default value} \\  
$T_{\mathrm{target}}$ & {\tt targetTemperature = 300;} &  K & none \\
$P_{\mathrm{target}}$ & {\tt targetPressure = 1;} & atm & none \\
$\tau_T$ & {\tt tauThermostat = 1e3;} & fs & none \\
$\tau_B$ & {\tt tauBarostat = 5e3;} & fs  & none \\
         & {\tt resetTime = 200;} & fs & none \\
         & {\tt useInitialExtendedSystemState = ``true'';} & logical &
false
\end{tabular}

Two additional keywords can be used to either clear the extended
system variables periodically ({\tt resetTime}), or to maintain the
state of the extended system variables between simulations ({\tt
useInitialExtendedSystemState}).  More details on these variables
and their use in the integrators follows below.

\subsubsection{\label{sec:noseHooverThermo}Nos\'{e}-Hoover Thermostatting}

The Nos\'e-Hoover equations of motion are given by\cite{Hoover85}
\begin{eqnarray}
\dot{{\bf r}} & = & {\bf v} \\
\dot{{\bf v}} & = & \frac{{\bf f}}{m} - \chi {\bf v} \\
\dot{\mathsf{A}} & = & \\
\dot{{\bf j}} & = & - \chi {\bf j}
\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}^2} \left( \frac{T}{T_{\mathrm{target}}} - 1 \right).
\label{eq:nosehooverext}
\end{equation}

The instantaneous temperature $T$ is proportional to the total kinetic
energy (both translational and orientational) and is given by
\begin{equation}
T = \frac{2 K}{f k_B}
\end{equation}
Here, $f$ is the total number of degrees of freedom in the system,
\begin{equation}
f = 3 N + 3 N_{\mathrm{orient}} - N_{\mathrm{constraints}}
\end{equation}
and $K$ is the total kinetic energy,
\begin{equation}
K = \sum_{i=1}^{N} \frac{1}{2} m_i {\bf v}_i^T \cdot {\bf v}_i +
\sum_{i=1}^{N_{\mathrm{orient}}} \sum_{\alpha=x,y,z} \frac{{\bf
j}_{i\alpha}^T \cdot {\bf j}_{i\alpha}}{2
\overleftrightarrow{\mathsf{I}}_{i,\alpha \alpha}}
\end{equation}

In eq.(\ref{eq:nosehooverext}), $\tau_T$ is the time constant for
relaxation of the temperature to the target value.  To set values for
$\tau_T$ or $T_{\mathrm{target}}$ in a simulation, one would use the
{\tt tauThermostat} and {\tt targetTemperature} keywords in the {\tt
.bass} file.  The units for {\tt tauThermostat} are fs, and the units
for the {\tt targetTemperature} are degrees K.   The integration of
the equations of motion is carried out in a velocity-Verlet style 2
part algorithm:

{\tt moveA:}
\begin{eqnarray}
T(t) & \leftarrow & \left\{{\bf v}(t)\right\}, \left\{{\bf j}(t)\right\} \\
{\bf v}\left(t + \delta t / 2\right)  & \leftarrow & {\bf
v}(t) + \frac{\delta t}{2} \left( \frac{{\bf f}(t)}{m} - {\bf v}(t)
\chi(t)\right) \\
{\bf r}(t + \delta t) & \leftarrow & {\bf r}(t) + \delta t {\bf
v}\left(t + \delta t / 2 \right) \\
{\bf j}\left(t + \delta t / 2 \right)  & \leftarrow & {\bf
j}(t) + \frac{\delta t}{2} \left( {\bf \tau}^b(t) - {\bf j}(t)
\chi(t) \right) \\
\mathsf{A}(t + \delta t) & \leftarrow & \mathrm{rot}({\bf j}(t +
\delta t / 2) \overleftrightarrow{\mathsf{I}}^{b},
\delta t) \\
\chi\left(t + \delta t / 2 \right) & \leftarrow & \chi(t) +
\frac{\delta t}{2 \tau_T^2} \left( \frac{T(t)}{T_{\mathrm{target}}} - 1
\right) 
\end{eqnarray}

Here $\mathrm{rot}( {\bf j} / {\bf I} )$ is the same symplectic Trotter
factorization of the three rotation operations that was discussed in
the section on the DLM integrator.  Note that this operation modifies
both the rotation matrix $\mathsf{A}$ and the angular momentum ${\bf
j}$.  {\tt moveA} propagates velocities by a half time step, and
positional degrees of freedom by a full time step.  The new positions
(and orientations) are then used to calculate a new set of forces and
torques.

{\tt doForces:}
\begin{eqnarray}
{\bf f}(t + \delta t) & \leftarrow & - \frac{\partial V}{\partial {\bf
r}(t + \delta t)} \\
{\bf \tau}^{s}(t + \delta t) & \leftarrow & {\bf u}(t + \delta t)
\times \frac{\partial V}{\partial {\bf u}} \\
{\bf \tau}^{b}(t + \delta t) & \leftarrow & \mathsf{A}(t + \delta t)
\cdot {\bf \tau}^s(t + \delta t) 
\end{eqnarray}

Here ${\bf u}$ is a unit vector pointing along the principal axis of
the rigid body being propagated, ${\bf \tau}^s$  is the torque in the
space-fixed (laboratory) frame, and ${\bf \tau}^b$ is the torque in
the body-fixed frame.  ${\bf u}$ is automatically calculated when the
rotation matrix $\mathsf{A}$ is calculated in {\tt moveA}.  

Once the forces and torques have been obtained at the new time step,
the velocities can be advanced to the same time value.

{\tt moveB:}
\begin{eqnarray}
T(t + \delta t) & \leftarrow & \left\{{\bf v}(t + \delta t)\right\},
\left\{{\bf j}(t + \delta t)\right\} \\
\chi\left(t + \delta t \right) & \leftarrow & \chi\left(t + \delta t /
2 \right) + \frac{\delta t}{2 \tau_T^2} \left( \frac{T(t+\delta
t)}{T_{\mathrm{target}}} - 1 \right) \\
{\bf v}\left(t + \delta t \right)  & \leftarrow & {\bf
v}\left(t + \delta t / 2 \right) + \frac{\delta t}{2} \left(
\frac{{\bf f}(t + \delta t)}{m} - {\bf v}(t + \delta t)
\chi(t \delta t)\right) \\
{\bf j}\left(t + \delta t \right)  & \leftarrow & {\bf
j}\left(t + \delta t / 2 \right) + \frac{\delta t}{2} \left( {\bf
\tau}^b(t + \delta t) - {\bf j}(t + \delta t)
\chi(t + \delta t) \right) 
\end{eqnarray}

Since ${\bf v}(t + \delta t)$ and ${\bf j}(t + \delta t)$ are required
to caclculate $T(t + \delta t)$ as well as $\chi(t + \delta t)$, the
indirectly depend on their own values at time $t + \delta t$.  {\tt
moveB} is therefore done in an iterative fashion until $\chi(t +
\delta t)$ becomes self-consistent.  The relative tolerance for the
self-consistency check defaults to a value of $\mbox{10}^{-6}$, but
{\sc oopse} will terminate the iteration after 4 loops even if the
consistency check has not been satisfied.

The Nos\'e-Hoover algorithm is known to conserve a Hamiltonian for the
extended system that is, to within a constant, identical to the
Helmholtz free energy,
\begin{equation}
H_{\mathrm{NVT}} = V + K + f k_B T_{\mathrm{target}} \left(
\frac{\tau_{T}^2 \chi^2(t)}{2} + \int_{0}^{t} \chi(t\prime) dt\prime
\right)
\end{equation}
Poor choices of $\delta t$ or $\tau_T$ can result in non-conservation
of $H_{\mathrm{NVT}}$, so the conserved quantity is maintained in the
last column of the {\tt .stat} file to allow checks on the quality of
the integration.

Bond constraints are applied at the end of both the {\tt moveA} and
{\tt moveB} portions of the algorithm.  Details on the constraint
algorithms are given in section \ref{sec:rattle}.

\subsubsection{\label{sec:NPTi}Constant-pressure integration (isotropic box)}


\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:	1076
Committed:	Mon Mar 1 22:29:28 2004 UTC (21 years, 8 months ago) by gezelter
Content type:	application/x-tex
File size:	18401 byte(s)
Log Message:	more changes for integrators
#	User	Rev	Content
1	mmeineke	899
2			\section{\label{sec:mechanics}Mechanics}
3
4	gezelter	1069	\subsection{\label{integrate}Integrating the Equations of Motion: the
5			DLM method}
6	mmeineke	899
7	gezelter	1069	The default method for integrating the equations of motion in {\sc
8	gezelter	1076	oopse} is a velocity-Verlet version of the symplectic splitting method
9			proposed by Dullweber, Leimkuhler and McLachlan
10			(DLM).\cite{Dullweber1997} When there are no directional atoms or
11			rigid bodies present in the simulation, this integrator becomes the
12			standard velocity-Verlet integrator which is known to sample the
13			microcanonical (NVE) ensemble.\cite{}
14	mmeineke	899
15	gezelter	1069	Previous integration methods for orientational motion have problems
16			that are avoided in the DLM method. Direct propagation of the Euler
17			angles has a known $1/\sin\theta$ divergence in the equations of
18			motion for $\phi$ and $\psi$,\cite{AllenTildesley} leading to
19			numerical instabilities any time one of the directional atoms or rigid
20			bodies has an orientation near $\theta=0$ or $\theta=\pi$. More
21			modern quaternion-based integration methods have relatively poor
22			energy conservation. While quaternions work well for orientational
23	gezelter	1076	motion in other ensembles, the microcanonical ensemble has a
24	gezelter	1069	constant energy requirement that is quite sensitive to errors in the
25			equations of motion. An earlier implementation of {\sc oopse}
26			utilized quaternions for propagation of rotational motion; however, a
27			detailed investigation showed that they resulted in a steady drift in
28			the total energy, something that has been observed by
29			others.\cite{Laird97}
30
31	mmeineke	899	The key difference in the integration method proposed by Dullweber
32			\emph{et al.} is that the entire rotation matrix is propagated from
33	gezelter	1076	one time step to the next. In the past, this would not have been
34			feasible, since that the rotation matrix for a single body has nine
35			elements compared to 3 or 4 elements for Euler angles and quaternions
36			respectively. System memory has become less costly in recent times,
37			and this can be translated into substantial benefits in energy
38			conservation.
39	mmeineke	899
40	gezelter	1069	The basic equations of motion being integrated are derived from the
41			Hamiltonian for conservative systems containing rigid bodies,
42			\begin{equation}
43	gezelter	1076	H = \sum_{i} \frac{{\bf v}_i^T m_i {\bf v}_i}{2} + \sum_i
44			\frac{{\bf j}_i^T \cdot {\bf j}_i}{2 \overleftrightarrow{\mathsf{I}}_i} +
45			V\left(\left\{{\bf r}\right\}, \left\{\mathsf{A}\right\}\right)
46	gezelter	1069	\end{equation}
47	gezelter	1076	where ${\bf r}_i$ and ${\bf v}_i$ are the cartesian position vector
48			and velocity of the center of mass of particle $i$, and ${\bf j}_i$
49			and $\overleftrightarrow{\mathsf{I}}_i$ are the body-fixed angular
50			momentum and moment of inertia tensor, respectively. $\mathsf{A}_i$
51			is the $3 \times 3$ rotation matrix describing the instantaneous
52			orientation of the particle. $V$ is the potential energy function
53			which may depend on both the positions $\left\{{\bf r}\right\}$ and
54			orientations $\left\{\mathsf{A}\right\}$ of all particles. The
55			equations of motion for the particle centers of mass are derived from
56			Hamilton's equations and are quite simple,
57	gezelter	1069	\begin{eqnarray}
58	gezelter	1076	\dot{{\bf r}} & = & {\bf v} \\
59			\dot{{\bf v}} & = & \frac{{\bf f}}{m}
60	gezelter	1069	\end{eqnarray}
61	gezelter	1076	where ${\bf f}$ is the instantaneous force on the center of mass
62			of the particle,
63	gezelter	1069	\begin{equation}
64	gezelter	1076	{\bf f} = \frac{\partial}{\partial
65			{\bf r}} V(\left\{{\bf r}(t)\right\}, \left\{\mathsf{A}(t)\right\}).
66	gezelter	1069	\end{equation}
67
68	gezelter	1076
69			The equations of motion for the orientational degrees of freedom are
70	gezelter	1069	\begin{eqnarray}
71	gezelter	1076	\dot{\mathsf{A}} & = & \mathsf{A}
72			\mbox{ skew}\left(\overleftrightarrow{I}^{-1} \cdot {\bf j}\right) \\
73			\dot{{\bf j}} & = & {\bf j} \times \left( \overleftrightarrow{I}^{-1}
74			\cdot {\bf j} \right) - \mbox{ rot}\left(\mathsf{A}^{T} \frac{\partial
75			V}{\partial \mathsf{A}} \right)
76	gezelter	1069	\end{eqnarray}
77	gezelter	1076	In these equations of motion, the $\mbox{skew}$ matrix of a vector
78			${\bf v} = \left( v_1, v_2, v_3 \right)$ is defined:
79			\begin{equation}
80			\mbox{skew}\left( {\bf v} \right) := \left(
81			\begin{array}{ccc}
82			0 & v_3 & - v_2 \\
83			-v_3 & 0 & v_1 \\
84			v_2 & -v_1 & 0
85			\end{array}
86			\right)
87			\end{equation}
88			The $\mbox{rot}$ notation refers to the mapping of the $3 \times 3$
89			rotation matrix to a vector of orientations by first computing the
90			skew-symmetric part $\left(\mathsf{A} - \mathsf{A}^{T}\right)$ and
91			then associating this with a length 3 vector by inverting the
92			$\mbox{skew}$ function above:
93			\begin{equation}
94			\mbox{rot}\left(\mathsf{A}\right) := \mbox{ skew}^{-1}\left(\mathsf{A}
95			- \mathsf{A}^{T} \right)
96			\end{equation}
97	gezelter	1069
98
99
100			In the integration method, the
101			orientational propagation involves a sequence of matrix evaluations to
102			update the rotation matrix.\cite{Dullweber1997} These matrix rotations
103			end up being more costly computationally than the simpler arithmetic
104			quaternion propagation. With the same time step, a 1000 SSD particle
105			simulation shows an average 7\% increase in computation time using the
106			symplectic step method in place of quaternions. This cost is more than
107			justified when comparing the energy conservation of the two methods as
108			illustrated in figure
109	mmeineke	899	\ref{timestep}.
110
111			\begin{figure}
112			\epsfxsize=6in
113			\epsfbox{timeStep.epsi}
114			\caption{Energy conservation using quaternion based integration versus
115			the symplectic step method proposed by Dullweber \emph{et al.} with
116			increasing time step. For each time step, the dotted line is total
117			energy using the symplectic step integrator, and the solid line comes
118			from the quaternion integrator. The larger time step plots are shifted
119			up from the true energy baseline for clarity.}
120			\label{timestep}
121			\end{figure}
122
123			In figure \ref{timestep}, the resulting energy drift at various time
124			steps for both the symplectic step and quaternion integration schemes
125			is compared. All of the 1000 SSD particle simulations started with the
126			same configuration, and the only difference was the method for
127			handling rotational motion. At time steps of 0.1 and 0.5 fs, both
128			methods for propagating particle rotation conserve energy fairly well,
129			with the quaternion method showing a slight energy drift over time in
130			the 0.5 fs time step simulation. At time steps of 1 and 2 fs, the
131			energy conservation benefits of the symplectic step method are clearly
132			demonstrated. Thus, while maintaining the same degree of energy
133			conservation, one can take considerably longer time steps, leading to
134			an overall reduction in computation time.
135
136			Energy drift in these SSD particle simulations was unnoticeable for
137			time steps up to three femtoseconds. A slight energy drift on the
138			order of 0.012 kcal/mol per nanosecond was observed at a time step of
139			four femtoseconds, and as expected, this drift increases dramatically
140			with increasing time step. To insure accuracy in the constant energy
141			simulations, time steps were set at 2 fs and kept at this value for
142			constant pressure simulations as well.
143
144
145			\subsection{\label{sec:extended}Extended Systems for other Ensembles}
146
147	gezelter	1076	{\sc oopse} implements a number of extended system integrators for
148			sampling from other ensembles relevant to chemical physics. The
149			integrator can selected with the {\tt ensemble} keyword in the
150			{\tt .bass} file:
151	mmeineke	899
152	gezelter	1076	\begin{center}
153			\begin{tabular}{lll}
154			{\bf Integrator} & {\bf Ensemble} & {\bf {\tt .bass} line} \\
155			NVE & microcanonical & {\tt ensemble = ``NVE''; } \\
156			NVT & canonical & {\tt ensemble = ``NVT''; } \\
157			NPTi & isobaric-isothermal (with isotropic volume changes) & {\tt
158			ensemble = ``NPTi'';} \\
159			NPTf & isobaric-isothermal (with changes to box shape) & {\tt
160			ensemble = ``NPTf'';} \\
161			NPTxyz & approximate isobaric-isothermal & {\tt ensemble =
162			``NPTxyz'';} \\
163			& (with separate barostats on each box dimension) &
164			\end{tabular}
165			\end{center}
166	gezelter	1069
167	gezelter	1076	The relatively well-known Nos\'e-Hoover thermostat is implemented in
168			{\sc oopse}'s NVT integrator. This method couples an extra degree of
169			freedom (the thermostat) to the kinetic energy of the system, and has
170			been shown to sample the canonical distribution in the system degrees
171			of freedom while conserving a quantity that is, to within a constant,
172			the Helmholtz free energy.
173	gezelter	1069
174	gezelter	1076	NPT algorithms attempt to maintain constant pressure in the system by
175			coupling the volume of the system to a barostat. {\sc oopse} contains
176			three different constant pressure algorithms. The first two, NPTi and
177			NPTf have been shown to conserve a quantity that is, to within a
178			constant, the Gibbs free energy. The Melchionna modification to the
179			Hoover barostat is implemented in both NPTi and NPTf. NPTi allows
180			only isotropic changes in the simulation box, while box {\it shape}
181			variations are allowed in NPTf. The NPTxyz integrator has {\it not}
182			been shown to sample from the isobaric-isothermal ensemble. It is
183			useful, however, in that it maintains orthogonality for the axes of
184			the simulation box while attempting to equalize pressure along the
185			three perpendicular directions in the box.
186	mmeineke	899
187	gezelter	1076	Each of the extended system integrators requires additional keywords
188			to set target values for the thermodynamic state variables that are
189			being held constant. Keywords are also required to set the
190			characteristic decay times for the dynamics of the extended
191			variables.
192	mmeineke	899
193	gezelter	1076	\begin{tabular}{llll}
194			{\bf variable} & {\bf {\tt .bass} keyword} & {\bf units} & {\bf
195			default value} \\
196			$T_{\mathrm{target}}$ & {\tt targetTemperature = 300;} & K & none \\
197			$P_{\mathrm{target}}$ & {\tt targetPressure = 1;} & atm & none \\
198			$\tau_T$ & {\tt tauThermostat = 1e3;} & fs & none \\
199			$\tau_B$ & {\tt tauBarostat = 5e3;} & fs & none \\
200			& {\tt resetTime = 200;} & fs & none \\
201			& {\tt useInitialExtendedSystemState = ``true'';} & logical &
202			false
203			\end{tabular}
204	mmeineke	899
205	gezelter	1076	Two additional keywords can be used to either clear the extended
206			system variables periodically ({\tt resetTime}), or to maintain the
207			state of the extended system variables between simulations ({\tt
208			useInitialExtendedSystemState}). More details on these variables
209			and their use in the integrators follows below.
210
211			\subsubsection{\label{sec:noseHooverThermo}Nos\'{e}-Hoover Thermostatting}
212
213			The Nos\'e-Hoover equations of motion are given by\cite{Hoover85}
214	mmeineke	899	\begin{eqnarray}
215			\dot{{\bf r}} & = & {\bf v} \\
216	gezelter	1076	\dot{{\bf v}} & = & \frac{{\bf f}}{m} - \chi {\bf v} \\
217			\dot{\mathsf{A}} & = & \\
218			\dot{{\bf j}} & = & - \chi {\bf j}
219	mmeineke	899	\label{eq:nosehoovereom}
220			\end{eqnarray}
221
222			$\chi$ is an ``extra'' variable included in the extended system, and
223			it is propagated using the first order equation of motion
224			\begin{equation}
225	gezelter	1076	\dot{\chi} = \frac{1}{\tau_{T}^2} \left( \frac{T}{T_{\mathrm{target}}} - 1 \right).
226	mmeineke	899	\label{eq:nosehooverext}
227			\end{equation}
228
229	gezelter	1076	The instantaneous temperature $T$ is proportional to the total kinetic
230			energy (both translational and orientational) and is given by
231			\begin{equation}
232			T = \frac{2 K}{f k_B}
233			\end{equation}
234			Here, $f$ is the total number of degrees of freedom in the system,
235			\begin{equation}
236			f = 3 N + 3 N_{\mathrm{orient}} - N_{\mathrm{constraints}}
237			\end{equation}
238			and $K$ is the total kinetic energy,
239			\begin{equation}
240			K = \sum_{i=1}^{N} \frac{1}{2} m_i {\bf v}_i^T \cdot {\bf v}_i +
241			\sum_{i=1}^{N_{\mathrm{orient}}} \sum_{\alpha=x,y,z} \frac{{\bf
242			j}_{i\alpha}^T \cdot {\bf j}_{i\alpha}}{2
243			\overleftrightarrow{\mathsf{I}}_{i,\alpha \alpha}}
244			\end{equation}
245	mmeineke	899
246	gezelter	1076	In eq.(\ref{eq:nosehooverext}), $\tau_T$ is the time constant for
247			relaxation of the temperature to the target value. To set values for
248			$\tau_T$ or $T_{\mathrm{target}}$ in a simulation, one would use the
249			{\tt tauThermostat} and {\tt targetTemperature} keywords in the {\tt
250			.bass} file. The units for {\tt tauThermostat} are fs, and the units
251			for the {\tt targetTemperature} are degrees K. The integration of
252			the equations of motion is carried out in a velocity-Verlet style 2
253			part algorithm:
254	mmeineke	899
255	gezelter	1076	{\tt moveA:}
256			\begin{eqnarray}
257			T(t) & \leftarrow & \left\{{\bf v}(t)\right\}, \left\{{\bf j}(t)\right\} \\
258			{\bf v}\left(t + \delta t / 2\right) & \leftarrow & {\bf
259			v}(t) + \frac{\delta t}{2} \left( \frac{{\bf f}(t)}{m} - {\bf v}(t)
260			\chi(t)\right) \\
261			{\bf r}(t + \delta t) & \leftarrow & {\bf r}(t) + \delta t {\bf
262			v}\left(t + \delta t / 2 \right) \\
263			{\bf j}\left(t + \delta t / 2 \right) & \leftarrow & {\bf
264			j}(t) + \frac{\delta t}{2} \left( {\bf \tau}^b(t) - {\bf j}(t)
265			\chi(t) \right) \\
266			\mathsf{A}(t + \delta t) & \leftarrow & \mathrm{rot}({\bf j}(t +
267			\delta t / 2) \overleftrightarrow{\mathsf{I}}^{b},
268			\delta t) \\
269			\chi\left(t + \delta t / 2 \right) & \leftarrow & \chi(t) +
270			\frac{\delta t}{2 \tau_T^2} \left( \frac{T(t)}{T_{\mathrm{target}}} - 1
271			\right)
272			\end{eqnarray}
273
274			Here $\mathrm{rot}( {\bf j} / {\bf I} )$ is the same symplectic Trotter
275			factorization of the three rotation operations that was discussed in
276			the section on the DLM integrator. Note that this operation modifies
277			both the rotation matrix $\mathsf{A}$ and the angular momentum ${\bf
278			j}$. {\tt moveA} propagates velocities by a half time step, and
279			positional degrees of freedom by a full time step. The new positions
280			(and orientations) are then used to calculate a new set of forces and
281			torques.
282
283			{\tt doForces:}
284			\begin{eqnarray}
285			{\bf f}(t + \delta t) & \leftarrow & - \frac{\partial V}{\partial {\bf
286			r}(t + \delta t)} \\
287			{\bf \tau}^{s}(t + \delta t) & \leftarrow & {\bf u}(t + \delta t)
288			\times \frac{\partial V}{\partial {\bf u}} \\
289			{\bf \tau}^{b}(t + \delta t) & \leftarrow & \mathsf{A}(t + \delta t)
290			\cdot {\bf \tau}^s(t + \delta t)
291			\end{eqnarray}
292
293			Here ${\bf u}$ is a unit vector pointing along the principal axis of
294			the rigid body being propagated, ${\bf \tau}^s$ is the torque in the
295			space-fixed (laboratory) frame, and ${\bf \tau}^b$ is the torque in
296			the body-fixed frame. ${\bf u}$ is automatically calculated when the
297			rotation matrix $\mathsf{A}$ is calculated in {\tt moveA}.
298
299			Once the forces and torques have been obtained at the new time step,
300			the velocities can be advanced to the same time value.
301
302			{\tt moveB:}
303			\begin{eqnarray}
304			T(t + \delta t) & \leftarrow & \left\{{\bf v}(t + \delta t)\right\},
305			\left\{{\bf j}(t + \delta t)\right\} \\
306			\chi\left(t + \delta t \right) & \leftarrow & \chi\left(t + \delta t /
307			2 \right) + \frac{\delta t}{2 \tau_T^2} \left( \frac{T(t+\delta
308			t)}{T_{\mathrm{target}}} - 1 \right) \\
309			{\bf v}\left(t + \delta t \right) & \leftarrow & {\bf
310			v}\left(t + \delta t / 2 \right) + \frac{\delta t}{2} \left(
311			\frac{{\bf f}(t + \delta t)}{m} - {\bf v}(t + \delta t)
312			\chi(t \delta t)\right) \\
313			{\bf j}\left(t + \delta t \right) & \leftarrow & {\bf
314			j}\left(t + \delta t / 2 \right) + \frac{\delta t}{2} \left( {\bf
315			\tau}^b(t + \delta t) - {\bf j}(t + \delta t)
316			\chi(t + \delta t) \right)
317			\end{eqnarray}
318
319			Since ${\bf v}(t + \delta t)$ and ${\bf j}(t + \delta t)$ are required
320			to caclculate $T(t + \delta t)$ as well as $\chi(t + \delta t)$, the
321			indirectly depend on their own values at time $t + \delta t$. {\tt
322			moveB} is therefore done in an iterative fashion until $\chi(t +
323			\delta t)$ becomes self-consistent. The relative tolerance for the
324			self-consistency check defaults to a value of $\mbox{10}^{-6}$, but
325			{\sc oopse} will terminate the iteration after 4 loops even if the
326			consistency check has not been satisfied.
327
328			The Nos\'e-Hoover algorithm is known to conserve a Hamiltonian for the
329			extended system that is, to within a constant, identical to the
330			Helmholtz free energy,
331			\begin{equation}
332			H_{\mathrm{NVT}} = V + K + f k_B T_{\mathrm{target}} \left(
333			\frac{\tau_{T}^2 \chi^2(t)}{2} + \int_{0}^{t} \chi(t\prime) dt\prime
334			\right)
335			\end{equation}
336			Poor choices of $\delta t$ or $\tau_T$ can result in non-conservation
337			of $H_{\mathrm{NVT}}$, so the conserved quantity is maintained in the
338			last column of the {\tt .stat} file to allow checks on the quality of
339			the integration.
340
341			Bond constraints are applied at the end of both the {\tt moveA} and
342			{\tt moveB} portions of the algorithm. Details on the constraint
343			algorithms are given in section \ref{sec:rattle}.
344
345			\subsubsection{\label{sec:NPTi}Constant-pressure integration (isotropic box)}
346
347
348	mmeineke	899	\subsection{\label{Sec:zcons}Z-Constraint Method}
349
350			Based on fluctuatin-dissipation theorem,\bigskip\ force auto-correlation
351			method was developed to investigate the dynamics of ions inside the ion
352			channels.\cite{Roux91} Time-dependent friction coefficient can be calculated
353			from the deviation of the instaneous force from its mean force.
354
355			%
356
357			\begin{equation}
358			\xi(z,t)=\langle\delta F(z,t)\delta F(z,0)\rangle/k_{B}T
359			\end{equation}
360
361
362			where%
363			\begin{equation}
364			\delta F(z,t)=F(z,t)-\langle F(z,t)\rangle
365			\end{equation}
366
367
368			If the time-dependent friction decay rapidly, static friction coefficient can
369			be approximated by%
370
371			\begin{equation}
372			\xi^{static}(z)=\int_{0}^{\infty}\langle\delta F(z,t)\delta F(z,0)\rangle dt
373			\end{equation}
374
375
376			Hence, diffusion constant can be estimated by
377			\begin{equation}
378			D(z)=\frac{k_{B}T}{\xi^{static}(z)}=\frac{(k_{B}T)^{2}}{\int_{0}^{\infty
379			}\langle\delta F(z,t)\delta F(z,0)\rangle dt}%
380			\end{equation}
381
382
383			\bigskip Z-Constraint method, which fixed the z coordinates of the molecules
384			with respect to the center of the mass of the system, was proposed to obtain
385			the forces required in force auto-correlation method.\cite{Marrink94} However,
386			simply resetting the coordinate will move the center of the mass of the whole
387			system. To avoid this problem, a new method was used at {\sc oopse}. Instead of
388			resetting the coordinate, we reset the forces of z-constraint molecules as
389			well as subtract the total constraint forces from the rest of the system after
390			force calculation at each time step.
391			\begin{verbatim}
392			$F_{\alpha i}=0$
393			$V_{\alpha i}=V_{\alpha i}-\frac{\sum\limits_{i}M_{_{\alpha i}}V_{\alpha i}}{\sum\limits_{i}M_{_{\alpha i}}}$
394			$F_{\alpha i}=F_{\alpha i}-\frac{M_{_{\alpha i}}}{\sum\limits_{\alpha}\sum\limits_{i}M_{_{\alpha i}}}\sum\limits_{\beta}F_{\beta}$
395			$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}}}$
396			\end{verbatim}
397
398			At the very beginning of the simulation, the molecules may not be at its
399			constraint position. To move the z-constraint molecule to the specified
400			position, a simple harmonic potential is used%
401
402			\begin{equation}
403			U(t)=\frac{1}{2}k_{Harmonic}(z(t)-z_{cons})^{2}%
404			\end{equation}
405			where $k_{Harmonic}$\bigskip\ is the harmonic force constant, $z(t)$ is
406			current z coordinate of the center of mass of the z-constraint molecule, and
407			$z_{cons}$ is the restraint position. Therefore, the harmonic force operated
408			on the z-constraint molecule at time $t$ can be calculated by%
409			\begin{equation}
410			F_{z_{Harmonic}}(t)=-\frac{\partial U(t)}{\partial z(t)}=-k_{Harmonic}%
411			(z(t)-z_{cons})
412			\end{equation}
413			Worthy of mention, other kinds of potential functions can also be used to
414			drive the z-constraint molecule.