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
Laird {\it et al.}\cite{Laird97}      

The key difference in the integration method proposed by Dullweber
\emph{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,
\begin{equation}
H = \sum_{i} \left( \frac{1}{2} m_i {\bf v}_i^T \cdot {\bf v}_i +
\frac{1}{2} {\bf j}_i^T \cdot \overleftrightarrow{\mathsf{I}}_i^{-1} \cdot
{\bf j}_i \right) +
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} \cdot
\mbox{ skew}\left(\overleftrightarrow{\mathsf{I}}^{-1} \cdot {\bf j}\right) \\
\dot{{\bf j}} & = & {\bf j} \times \left( \overleftrightarrow{\mathsf{I}}^{-1}
\cdot {\bf j} \right) - \mbox{ rot}\left(\mathsf{A}^{T} \cdot \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}
Written this way, the $\mbox{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,
\begin{eqnarray}
\dot{j_{x}} & = & \tau^b_x(t)  +
\left(\overleftrightarrow{\mathsf{I}}_{yy} - \overleftrightarrow{\mathsf{I}}_{zz} \right) j_y j_z \\
\dot{j_{y}} & = & \tau^b_y(t) +
\left(\overleftrightarrow{\mathsf{I}}_{zz} - \overleftrightarrow{\mathsf{I}}_{xx} \right) j_z j_x \\
\dot{j_{z}} & = & \tau^b_z(t) +
\left(\overleftrightarrow{\mathsf{I}}_{xx} - \overleftrightarrow{\mathsf{I}}_{yy} \right) j_x j_y 
\end{eqnarray}
which utilize the body-fixed torques, ${\bf \tau}^b$. Torques are
most easily derived in the space-fixed frame, 
\begin{equation}
{\bf \tau}^b(t) = \mathsf{A}(t) \cdot {\bf \tau}^s(t)
\end{equation}
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,
\begin{equation}
{\bf \tau}^s(t) = - \hat{\bf u}(t) \times \left( \frac{\partial}
{\partial \hat{\bf u}} V\left(\left\{ {\bf r}(t) \right\}, \left\{
\mathsf{A}(t) \right\}\right) \right).
\end{equation}
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:
\begin{enumerate}
\item the integrator is area-preserving in phase space (i.e. it is
{\it symplectic}),
\item the integrator is time-{\it reversible}, making it suitable for Hybrid
Monte Carlo applications, and
\item the error for a single time step is of order $O\left(h^3\right)$
for timesteps of length $h$.
\end{enumerate}

The integration of the equations of motion is carried out in a
velocity-Verlet style 2-part algorithm:

{\tt moveA:}
\begin{eqnarray}
{\bf v}\left(t + \delta t / 2\right)  & \leftarrow & {\bf
v}(t) + \frac{\delta t}{2} \left( {\bf f}(t) / m \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} {\bf \tau}^b(t)  \\
\mathsf{A}(t + \delta t) & \leftarrow & \mathrm{rot}\left( \delta t
{\bf j}(t + \delta t / 2) \cdot \overleftrightarrow{\mathsf{I}}^{-1}
\right)
\end{eqnarray}

In this context, the $\mathrm{rot}$ function is the reversible product
of the three body-fixed rotations,
\begin{equation}
\mathrm{rot}({\bf a}) = \mathsf{G}_x(a_x / 2) \cdot
\mathsf{G}_y(a_y / 2) \cdot \mathsf{G}_z(a_z) \cdot \mathsf{G}_y(a_y /
2) \cdot \mathsf{G}_x(a_x /2)
\end{equation}
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$,
\begin{equation}
\mathsf{G}_\alpha( \theta ) = \left\{
\begin{array}{lcl}
\mathsf{A}(t) & \leftarrow & \mathsf{A}(0) \cdot \mathsf{R}_\alpha(\theta)^T \\
{\bf j}(t) & \leftarrow & \mathsf{R}_\alpha(\theta) \cdot {\bf j}(0)
\end{array}
\right.
\end{equation}
$\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
\begin{equation}
\mathsf{R}_x(\theta) \approx \left(
\begin{array}{ccc}
1 & 0 & 0 \\
0 & \frac{1-\theta^2 / 4}{1 + \theta^2 / 4}  & -\frac{\theta}{1+
\theta^2 / 4} \\
0 & \frac{\theta}{1+
\theta^2 / 4} & \frac{1-\theta^2 / 4}{1 + \theta^2 / 4}
\end{array}
\right).
\end{equation}
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

{\tt doForces:}
\begin{eqnarray}
{\bf f}(t + \delta t) & \leftarrow & - \left(\frac{\partial V}{\partial {\bf
r}}\right)_{{\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}

{\sc oopse} automatically updates ${\bf u}$ 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}
{\bf v}\left(t + \delta t \right)  & \leftarrow & {\bf
v}\left(t + \delta t / 2 \right) + \frac{\delta t}{2} \left(
{\bf f}(t + \delta t) / m \right) \\
{\bf j}\left(t + \delta t \right)  & \leftarrow & {\bf
j}\left(t + \delta t / 2 \right) + \frac{\delta t}{2} {\bf
\tau}^b(t + \delta t)  
\end{eqnarray}

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 1000-molecule water simulation
shows an average 7\% increase in computation time using the DLM 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 method proposed by Dullweber \emph{et al.} with increasing time
step. For each time step, the dotted line is total energy using the
DLM 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 DLM and quaternion integration schemes is
compared. All of the 1000 molecule water 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 molecule 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 DLM 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.

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

\begin{center}
\begin{tabular}{llll}
{\bf variable} & {\bf {\tt .bass} keyword} & {\bf units} & {\bf
default value} \\  
$\delta t$ & {\tt dt = 2.0;} & fs & none 
\end{tabular}
\end{center}

\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}} & = & \mathsf{A} \cdot
\mbox{ skew}\left(\overleftrightarrow{\mathsf{I}}^{-1} \cdot {\bf j}\right) \\
\dot{{\bf j}} & = & {\bf j} \times \left( \overleftrightarrow{\mathsf{I}}^{-1}
\cdot {\bf j} \right) - \mbox{ rot}\left(\mathsf{A}^{T} \cdot \frac{\partial
V}{\partial \mathsf{A}} \right) - \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}}}  \frac{1}{2} {\bf j}_i^T \cdot
\overleftrightarrow{\mathsf{I}}_i^{-1} \cdot {\bf j}_i
\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}\left(\delta t *
{\bf j}(t + \delta t / 2) \overleftrightarrow{\mathsf{I}}^{-1} \right) \\
\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}(\delta t * {\bf j}
\overleftrightarrow{\mathsf{I}}^{-1})$ 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 in exactly the same way they are calculated in the {\tt
doForces} portion of the DLM integrator.

Once the forces and torques have been obtained at the new time step,
the temperature, velocities, and the extended system variable 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)$, they
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 with 
isotropic box deformations (NPTi)}

To carry out isobaric-isothermal ensemble calculations {\sc oopse}
implements the Melchionna modifications to the Nos\'e-Hoover-Andersen
equations of motion,\cite{Melchionna93}

\begin{eqnarray}
\dot{{\bf r}} & = & {\bf v} + \eta \left( {\bf r} - {\bf R}_0 \right) \\
\dot{{\bf v}} & = & \frac{{\bf f}}{m} - (\eta + \chi) {\bf v} \\
\dot{\mathsf{A}} & = & \mathsf{A} \cdot
\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} \cdot \frac{\partial
V}{\partial \mathsf{A}} \right) - \chi {\bf j} \\
\dot{\chi} & = & \frac{1}{\tau_{T}^2} \left(
\frac{T}{T_{\mathrm{target}}} - 1 \right) \\
\dot{\eta} & = & \frac{1}{\tau_{B}^2 f k_B T_{\mathrm{target}}} V \left( P -
P_{\mathrm{target}} \right) \\
\dot{\mathcal{V}} & = & 3 \mathcal{V} \eta 
\label{eq:melchionna1}
\end{eqnarray}

$\chi$ and $\eta$ are the ``extra'' degrees of freedom in the extended
system.  $\chi$ is a thermostat, and it has the same function as it
does in the Nos\'e-Hoover NVT integrator.  $\eta$ is a barostat which
controls changes to the volume of the simulation box.  ${\bf R}_0$ is
the location of the center of mass for the entire system, and
$\mathcal{V}$ is the volume of the simulation box.  At any time, the
volume can be calculated from the determinant of the matrix which
describes the box shape:
\begin{equation}
\mathcal{V} = \det(\mathsf{H})
\end{equation}

The NPTi integrator requires an instantaneous pressure. This quantity
is calculated via the pressure tensor,
\begin{equation}
\overleftrightarrow{\mathsf{P}}(t) = \frac{1}{\mathcal{V}(t)} \left(
\sum_{i=1}^{N} m_i {\bf v}_i(t) \otimes {\bf v}_i(t) \right) +
\overleftrightarrow{\mathsf{W}}(t)
\end{equation}
The kinetic contribution to the pressure tensor utilizes the {\it
outer} product of the velocities denoted by the $\otimes$ symbol.  The
stress tensor is calculated from another outer product of the
inter-atomic separation vectors (${\bf r}_{ij} = {\bf r}_j - {\bf
r}_i$) with the forces between the same two atoms,
\begin{equation}
\overleftrightarrow{\mathsf{W}}(t) = \sum_{i} \sum_{j>i} {\bf r}_{ij}(t)
\otimes {\bf f}_{ij}(t)
\end{equation}
The instantaneous pressure is then simply obtained from the trace of
the Pressure tensor,
\begin{equation}
P(t) = \frac{1}{3} \mathrm{Tr} \left( \overleftrightarrow{\mathsf{P}}(t)
\right)
\end{equation}

In eq.(\ref{eq:melchionna1}), $\tau_B$ is the time constant for
relaxation of the pressure to the target value.  To set values for
$\tau_B$ or $P_{\mathrm{target}}$ in a simulation, one would use the
{\tt tauBarostat} and {\tt targetPressure} keywords in the {\tt .bass}
file.  The units for {\tt tauBarostat} are fs, and the units for the
{\tt targetPressure} are atmospheres.  Like in the NVT integrator, 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\} \\
P(t) & \leftarrow & \left\{{\bf r}(t)\right\}, \left\{{\bf v}(t)\right\}, \left\{{\bf f}(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)
\left(\chi(t) + \eta(t) \right) \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}\left(\delta t *
{\bf j}(t + \delta t / 2) \overleftrightarrow{\mathsf{I}}^{-1} \right) \\
\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) \\
\eta(t + \delta t / 2) & \leftarrow & \eta(t) + \frac{\delta t \mathcal{V}(t)}{2 N k_B
T(t) \tau_B^2} \left( P(t) - P_{\mathrm{target}} \right) \\ 
{\bf r}(t + \delta t) & \leftarrow & {\bf r}(t) + \delta t \left\{ {\bf
v}\left(t + \delta t / 2 \right) + \eta(t + \delta t / 2)\left[ {\bf
r}(t + \delta t) - {\bf R}_0 \right] \right\} \\
\mathsf{H}(t + \delta t) & \leftarrow & e^{-\delta t \eta(t + \delta t
/ 2)} \mathsf{H}(t)
\end{eqnarray}

Most of these equations are identical to their counterparts in the NVT
integrator, but the propagation of positions to time $t + \delta t$
depends on the positions at the same time.  {\sc oopse} carries out
this step iteratively (with a limit of 5 passes through the iterative
loop).  Also, the simulation box $\mathsf{H}$ is scaled uniformly for
one full time step by an exponential factor that depends on the value
of $\eta$ at time $t +
\delta t / 2$.  Reshaping the box uniformly also scales the volume of
the box by
\begin{equation}
\mathcal{V}(t + \delta t) \leftarrow e^{ - 3 \delta t \eta(t + \delta t /2)}
\mathcal{V}(t)
\end{equation}

The {\tt doForces} step for the NPTi integrator is exactly the same as
in both the DLM and NVT integrators.  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\} \\
P(t + \delta t) & \leftarrow & \left\{{\bf r}(t + \delta t)\right\},
\left\{{\bf v}(t + \delta t)\right\}, \left\{{\bf f}(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) \\
\eta(t + \delta t) & \leftarrow & \eta(t + \delta t / 2) +
\frac{\delta t \mathcal{V}(t + \delta t)}{2 N k_B T(t + \delta t) \tau_B^2}
\left( P(t + \delta t) - P_{\mathrm{target}} 
\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) + \eta(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}

Once again, since ${\bf v}(t + \delta t)$ and ${\bf j}(t + \delta t)$
are required to caclculate $T(t + \delta t)$, $P(t + \delta t)$, $\chi(t +
\delta t)$, and $\eta(t + \delta t)$, they 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)$ and $\eta(t + \delta
t)$ become 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 Melchionna modification of the Nos\'e-Hoover-Andersen algorithm is
known to conserve a Hamiltonian for the extended system that is, to
within a constant, identical to the Gibbs free energy,
\begin{equation}
H_{\mathrm{NPTi}} = 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) + P_{\mathrm{target}} \mathcal{V}(t).
\end{equation}
Poor choices of $\delta t$, $\tau_T$, or $\tau_B$ can result in
non-conservation of $H_{\mathrm{NPTi}}$, so the conserved quantity is
maintained in the last column of the {\tt .stat} file to allow checks
on the quality of the integration.  It is also known that this
algorithm samples the equilibrium distribution for the enthalpy
(including contributions for the thermostat and barostat), 
\begin{equation}
H_{\mathrm{NPTi}} = V + K + \frac{f k_B T_{\mathrm{target}}}{2} \left(
\chi^2 \tau_T^2 + \eta^2 \tau_B^2 \right) +  P_{\mathrm{target}}
\mathcal{V}(t). 
\end{equation}

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:NPTf}Constant-pressure integration with a
flexible box (NPTf)} 

There is a relatively simple generalization of the
Nos\'e-Hoover-Andersen method to include changes in the simulation box
{\it shape} as well as in the volume of the box.  This method utilizes
the full $3 \times 3$ pressure tensor and introduces a tensor of
extended variables ($\overleftrightarrow{\eta}$) to control changes to
the box shape.  The equations of motion for this method are
\begin{eqnarray}
\dot{{\bf r}} & = & {\bf v} + \overleftrightarrow{\eta} \cdot \left( {\bf r} - {\bf R}_0 \right) \\
\dot{{\bf v}} & = & \frac{{\bf f}}{m} - (\overleftrightarrow{\eta} +
\chi \mathsf{1}) {\bf v} \\
\dot{\mathsf{A}} & = & \mathsf{A} \cdot
\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} \cdot \frac{\partial
V}{\partial \mathsf{A}} \right) - \chi {\bf j} \\
\dot{\chi} & = & \frac{1}{\tau_{T}^2} \left(
\frac{T}{T_{\mathrm{target}}} - 1 \right) \\
\dot{\overleftrightarrow{eta}} & = & \frac{1}{\tau_{B}^2 f k_B
T_{\mathrm{target}}} V \left( \overleftrightarrow{\mathsf{P}} - P_{\mathrm{target}}\mathsf{1} \right) \\
\dot{\mathsf{H}} & = &  \overleftrightarrow{\eta} \cdot \mathsf{H} 
\label{eq:melchionna2}
\end{eqnarray}

Here, $\mathsf{1}$ is the unit matrix and $\overleftrightarrow{\mathsf{P}}$
is the pressure tensor.  Again, the volume, $\mathcal{V} = \det
\mathsf{H}$. 

The propagation of the equations of motion is nearly identical to the
NPTi integration:

{\tt moveA:}
\begin{eqnarray}
T(t) & \leftarrow & \left\{{\bf v}(t)\right\}, \left\{{\bf j}(t)\right\} \\
\overleftrightarrow{\mathsf{P}}(t) & \leftarrow & \left\{{\bf r}(t)\right\}, \left\{{\bf v}(t)\right\}, \left\{{\bf f}(t)\right\} \\
{\bf v}\left(t + \delta t / 2\right)  & \leftarrow & {\bf
v}(t) + \frac{\delta t}{2} \left( \frac{{\bf f}(t)}{m} - 
\left(\chi(t)\mathsf{1} + \overleftrightarrow{\eta}(t) \right) \cdot
{\bf v}(t) \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}\left(\delta t *
{\bf j}(t + \delta t / 2) \overleftrightarrow{\mathsf{I}}^{-1} \right) \\
\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) \\
\overleftrightarrow{\eta}(t + \delta t / 2) & \leftarrow & \overleftrightarrow{\eta}(t) + \frac{\delta t \mathcal{V}(t)}{2 N k_B
T(t) \tau_B^2} \left( \overleftrightarrow{\mathsf{P}}(t) - P_{\mathrm{target}}\mathsf{1} \right) \\ 
{\bf r}(t + \delta t) & \leftarrow & {\bf r}(t) + \delta t \left\{ {\bf
v}\left(t + \delta t / 2 \right) + \overleftrightarrow{\eta}(t +
\delta t / 2) \cdot \left[ {\bf
r}(t + \delta t) - {\bf R}_0 \right] \right\} \\
\mathsf{H}(t + \delta t) & \leftarrow & \mathsf{H}(t) \cdot e^{-\delta t
\overleftrightarrow{\eta}(t + \delta t / 2)} 
\end{eqnarray}
{\sc oopse} uses a power series expansion truncated at second order
for the exponential operation which scales the simulation box.

The {\tt moveB} portion of the algorithm is largely unchanged from the
NPTi integrator:

{\tt moveB:}
\begin{eqnarray}
T(t + \delta t) & \leftarrow & \left\{{\bf v}(t + \delta t)\right\},
\left\{{\bf j}(t + \delta t)\right\} \\
\overleftrightarrow{\mathsf{P}}(t + \delta t) & \leftarrow & \left\{{\bf r}(t + \delta t)\right\},
\left\{{\bf v}(t + \delta t)\right\}, \left\{{\bf f}(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) \\
\overleftrightarrow{\eta}(t + \delta t) & \leftarrow & \overleftrightarrow{\eta}(t + \delta t / 2) +
\frac{\delta t \mathcal{V}(t + \delta t)}{2 N k_B T(t + \delta t) \tau_B^2}
\left( \overleftrightarrow{P}(t + \delta t) - P_{\mathrm{target}}\mathsf{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} - 
(\chi(t + \delta t)\mathsf{1} + \overleftrightarrow{\eta}(t + \delta
t)) \right) \cdot {\bf v}(t + \delta t) \\
{\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}

The iterative schemes for both {\tt moveA} and {\tt moveB} are
identical to those described for the NPTi integrator.

The NPTf integrator is known to conserve the following Hamiltonian:
\begin{equation}
H_{\mathrm{NPTf}} = 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) + P_{\mathrm{target}} \mathcal{V}(t) + \frac{f k_B
T_{\mathrm{target}}}{2}
\mathrm{Tr}\left[\overleftrightarrow{\eta}(t)\right]^2 \tau_B^2.
\end{equation}

This integrator must be used with care, particularly in liquid
simulations.  Liquids have very small restoring forces in the
off-diagonal directions, and the simulation box can very quickly form
elongated and sheared geometries which become smaller than the
electrostatic or Lennard-Jones cutoff radii.  It finds most use in
simulating crystals or liquid crystals which assume non-orthorhombic
geometries.

\subsubsection{\label{nptxyz}Constant pressure in 3 axes (NPTxyz)}

There is one additional extended system integrator which is somewhat
simpler than the NPTf method described above.  In this case, the three
axes have independent barostats which each attempt to preserve the
target pressure along the box walls perpendicular to that particular
axis.  The lengths of the box axes are allowed to fluctuate
independently, but the angle between the box axes does not change.
The equations of motion are identical to those described above, but
only the {\it diagonal} elements of $\overleftrightarrow{\eta}$ are
computed.  The off-diagonal elements are set to zero (even when the
pressure tensor has non-zero off-diagonal elements).

It should be noted that the NPTxyz integrator is {\it not} known to
preserve any Hamiltonian of interest to the chemical physics
community.  The integrator is extremely useful, however, in generating
initial conditions for other integration methods.  It {\it is} suitable
for use with liquid simulations, or in cases where there is
orientational anisotropy in the system (i.e. in lipid bilayer
simulations).

\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:	1082
Committed:	Wed Mar 3 23:16:18 2004 UTC (21 years, 4 months ago) by gezelter
Content type:	application/x-tex
File size:	37228 byte(s)
Log Message:	NPTf complete, NPTxyz complete