trunk/mattDisertation/oopse.tex

\chapter{\label{chapt:oopse}OOPSE: AN OPEN SOURCE OBJECT-ORIENTED PARALLEL SIMULATION ENGINE FOR MOLECULAR DYNAMICS}


%% \begin{abstract}
%% We detail the capabilities of a new open-source parallel simulation
%% package ({\sc oopse}) that can perform molecular dynamics simulations
%% on atom types that are missing from other popular packages.  In
%% particular, {\sc oopse} is capable of performing orientational
%% dynamics on dipolar systems, and it can handle simulations of metallic
%% systems using the embedded atom method ({\sc eam}).
%% \end{abstract}

\lstset{language=C,frame=TB,basicstyle=\small,basicstyle=\ttfamily, %
        xleftmargin=0.5in, xrightmargin=0.5in,captionpos=b, %
        abovecaptionskip=0.5cm, belowcaptionskip=0.5cm}

\section{\label{oopseSec:foreword}Foreword}

In this chapter, I present and detail the capabilities of the open
source simulation program {\sc oopse}. It is important to note that a
simulation program of this size and scope would not have been possible
without the collaborative efforts of my colleagues: Charles
F.~Vardeman II, Teng Lin, Christopher J.~Fennell and J.~Daniel
Gezelter. Although my contributions to {\sc oopse} are major,
consideration of my work apart from the others would not give a
complete description to the program's capabilities. As such, all
contributions to {\sc oopse} to date are presented in this chapter.

Charles Vardeman is responsible for the parallelization of the long
range forces in {\sc oopse} (Sec.~\ref{oopseSec:parallelization}) as
well as the inclusion of the embedded-atom potential for transition
metals (Sec.~\ref{oopseSec:eam}). Teng Lin's contributions include
refinement of the periodic boundary conditions
(Sec.~\ref{oopseSec:pbc}), the z-constraint method
(Sec.~\ref{oopseSec:zcons}), refinement of the property analysis
programs (Sec.~\ref{oopseSec:props}), and development in the extended
system integrators (Sec.~\ref{oopseSec:noseHooverThermo}). Christopher
Fennell worked on the symplectic integrator
(Sec.~\ref{oopseSec:integrate}) and the refinement of the {\sc ssd}
water model (Sec.~\ref{oopseSec:SSD}). Daniel Gezelter lent his
talents in the development of the extended system integrators
(Sec.~\ref{oopseSec:noseHooverThermo}) as well as giving general
direction and oversight to the entire project. My responsibilities
covered the creation and specification of {\sc bass}
(Sec.~\ref{oopseSec:IOfiles}), the original development of the single
processor version of {\sc oopse}, contributions to the extended state
integrators (Sec.~\ref{oopseSec:noseHooverThermo}), the implementation
of the Lennard-Jones (Sec.~\ref{sec:LJPot}) and {\sc duff}
(Sec.~\ref{oopseSec:DUFF}) force fields, and initial implementation of
the property analysis (Sec.~\ref{oopseSec:props}) and system
initialization (Sec.~\ref{oopseSec:initCoords}) utility programs. {\sc
oopse}, like many other Molecular Dynamics programs, is a work in
progress, and will continue to be so for many graduate student
lifetimes.

\section{\label{sec:intro}Introduction}

When choosing to simulate a chemical system with molecular dynamics,
there are a variety of options available. For simple systems, one
might consider writing one's own programming code. However, as systems
grow larger and more complex, building and maintaining code for the
simulations becomes a time consuming task. In such cases it is usually
more convenient for a researcher to turn to pre-existing simulation
packages. These packages, such as {\sc amber}\cite{pearlman:1995} and
{\sc charmm}\cite{Brooks83}, provide powerful tools for researchers to
conduct simulations of their systems without spending their time
developing a code base to conduct their research. This then frees them
to perhaps explore experimental analogues to their models. 

Despite their utility, problems with these packages arise when
researchers try to develop techniques or energetic models that the
code was not originally designed to simulate. Examples of techniques
and energetics not commonly implemented include; dipole-dipole
interactions, rigid body dynamics, and metallic potentials. When faced
with these obstacles, a researcher must either develop their own code
or license and extend one of the commercial packages. What we have
elected to do is develop a body of simulation code capable of
implementing the types of models upon which our research is based.

In developing {\sc oopse}, we have adhered to the precepts of Open
Source development, and are releasing our source code with a
permissive license. It is our intent that by doing so, other
researchers might benefit from our work, and add their own
contributions to the package. The license under which {\sc oopse} is
distributed allows any researcher to download and modify the source
code for their own use. In this way further development of {\sc oopse}
is not limited to only the models of interest to ourselves, but also
those of the community of scientists who contribute back to the
project.

We have structured this chapter to first discuss the empirical energy
functions that {\sc oopse } implements in
Sec.~\ref{oopseSec:empiricalEnergy}. Following that is a discussion of
the various input and output files associated with the package
(Sec.~\ref{oopseSec:IOfiles}). Sec.~\ref{oopseSec:mechanics}
elucidates the various Molecular Dynamics algorithms {\sc oopse}
implements in the integration of the Newtonian equations of
motion. Basic analysis of the trajectories obtained from the
simulation is discussed in Sec.~\ref{oopseSec:props}. Program design
considerations are presented in Sec.~\ref{oopseSec:design}. And
lastly, Sec.~\ref{oopseSec:conclusion} concludes the chapter.

\section{\label{oopseSec:empiricalEnergy}The Empirical Energy Functions}

\subsection{\label{oopseSec:atomsMolecules}Atoms, Molecules and Rigid Bodies}

The basic unit of an {\sc oopse} simulation is the atom. The
parameters describing the atom are generalized to make the atom as
flexible a representation as possible. They may represent specific
atoms of an element, or be used for collections of atoms such as
methyl and carbonyl groups. The atoms are also capable of having
directional components associated with them (\emph{e.g.}~permanent
dipoles). Charges, permanent dipoles, and Lennard-Jones parameters for
a given atom type are set in the force field parameter files.

\begin{lstlisting}[float,caption={[Specifier for molecules and atoms] A sample specification of an Ar molecule},label=sch:AtmMole]
molecule{
  name = "Ar";
  nAtoms = 1;
  atom[0]{
    type="Ar";
    position( 0.0, 0.0, 0.0 );
  }
}
\end{lstlisting}


Atoms can be collected into secondary structures such as rigid bodies
or molecules. The molecule is a way for {\sc oopse} to keep track of
the atoms in a simulation in logical manner. Molecular units store the
identities of all the atoms and rigid bodies associated with
themselves, and are responsible for the evaluation of their own
internal interactions (\emph{i.e.}~bonds, bends, and torsions). Scheme
\ref{sch:AtmMole} shows how one creates a molecule in a ``model'' or
\texttt{.mdl} file. The position of the atoms given in the
declaration are relative to the origin of the molecule, and is used
when creating a system containing the molecule.

As stated previously, one of the features that sets {\sc oopse} apart
from most of the current molecular simulation packages is the ability
to handle rigid body dynamics. Rigid bodies are non-spherical
particles or collections of particles that have a constant internal
potential and move collectively.\cite{Goldstein01} They are not
included in most simulation packages because of the algorithmic
complexity involved in propagating orientational degrees of
freedom. Until recently, integrators which propagate orientational
motion have been much worse than those available for translational
motion.

Moving a rigid body involves determination of both the force and
torque applied by the surroundings, which directly affect the
translational and rotational motion in turn. In order to accumulate
the total force on a rigid body, the external forces and torques must
first be calculated for all the internal particles. The total force on
the rigid body is simply the sum of these external forces.
Accumulation of the total torque on the rigid body is more complex
than the force because the torque is applied to the center of mass of
the rigid body. The torque on rigid body $i$ is
\begin{equation}
\boldsymbol{\tau}_i=
        \sum_{a}\biggl[(\mathbf{r}_{ia}-\mathbf{r}_i)\times \mathbf{f}_{ia} 
        + \boldsymbol{\tau}_{ia}\biggr]
\label{eq:torqueAccumulate}
\end{equation}
where $\boldsymbol{\tau}_i$ and $\mathbf{r}_i$ are the torque on and
position of the center of mass respectively, while $\mathbf{f}_{ia}$,
$\mathbf{r}_{ia}$, and $\boldsymbol{\tau}_{ia}$ are the force on,
position of, and torque on the component particles of the rigid body.

The summation of the total torque is done in the body fixed axis of
each rigid body. In order to move between the space fixed and body
fixed coordinate axes, parameters describing the orientation must be
maintained for each rigid body. At a minimum, the rotation matrix
($\mathsf{A}$) can be described by the three Euler angles ($\phi,
\theta,$ and $\psi$), where the elements of $\mathsf{A}$ are composed of
trigonometric operations involving $\phi, \theta,$ and
$\psi$.\cite{Goldstein01} In order to avoid numerical instabilities
inherent in using the Euler angles, the four parameter ``quaternion''
scheme is often used. The elements of $\mathsf{A}$ can be expressed as
arithmetic operations involving the four quaternions ($q_0, q_1, q_2,$
and $q_3$).\cite{allen87:csl} Use of quaternions also leads to
performance enhancements, particularly for very small
systems.\cite{Evans77}

{\sc oopse} utilizes a relatively new scheme that propagates the
entire nine parameter rotation matrix. Further discussion
on this choice can be found in Sec.~\ref{oopseSec:integrate}. An example
definition of a rigid body can be seen in Scheme
\ref{sch:rigidBody}. The positions in the atom definitions are the
placements of the atoms relative to the origin of the rigid body,
which itself has a position relative to the origin of the molecule.

\begin{lstlisting}[float,caption={[Defining rigid bodies]A sample definition of a rigid body},label={sch:rigidBody}]
molecule{
  name = "TIP3P";
  nAtoms = 3;
  atom[0]{
    type = "O_TIP3P";
    position( 0.0, 0.0, -0.06556 );
  }
  atom[1]{
    type = "H_TIP3P";
    position( 0.0, 0.75695, 0.52032 );
  }
  atom[2]{
    type = "H_TIP3P";
    position( 0.0, -0.75695, 0.52032 );
  }

  nRigidBodies = 1;
  rigidBody[0]{
    nMembers = 3;
    members(0, 1, 2);
  }
}
\end{lstlisting}

\subsection{\label{sec:LJPot}The Lennard Jones Force Field}

The most basic force field implemented in {\sc oopse} is the
Lennard-Jones force field, which mimics the van der Waals interaction at
long distances, and uses an empirical repulsion at short
distances. The Lennard-Jones potential is given by:
\begin{equation}
V_{\text{LJ}}(r_{ij}) = 
        4\epsilon_{ij} \biggl[
        \biggl(\frac{\sigma_{ij}}{r_{ij}}\biggr)^{12}
        - \biggl(\frac{\sigma_{ij}}{r_{ij}}\biggr)^{6}
        \biggr]
\label{eq:lennardJonesPot}
\end{equation}
Where $r_{ij}$ is the distance between particles $i$ and $j$,
$\sigma_{ij}$ scales the length of the interaction, and
$\epsilon_{ij}$ scales the well depth of the potential. Scheme
\ref{sch:LJFF} gives an example \texttt{.bass} file that
sets up a system of 108 Ar particles to be simulated using the
Lennard-Jones force field.

\begin{lstlisting}[float,caption={[Invocation of the Lennard-Jones force field] A sample system using the Lennard-Jones force field.},label={sch:LJFF}]

#include "argon.mdl" 

nComponents = 1;
component{
  type = "Ar";
  nMol = 108;
}

initialConfig = "./argon.init";

forceField = "LJ";
\end{lstlisting}

Because this potential is calculated between all pairs, the force
evaluation can become computationally expensive for large systems. To
keep the pair evaluations to a manageable number, {\sc oopse} employs
a cut-off radius.\cite{allen87:csl} The cutoff radius can either be
specified in the \texttt{.bass} file, or left as its default value of
$2.5\sigma_{ii}$, where $\sigma_{ii}$ is the largest Lennard-Jones
length parameter present in the simulation. Truncating the calculation
at $r_{\text{cut}}$ introduces a discontinuity into the potential
energy and the force. To offset this discontinuity in the potential,
the energy value at $r_{\text{cut}}$ is subtracted from the
potential. This causes the potential to go to zero smoothly at the
cut-off radius, and preserves conservation of energy in integrating
the equations of motion. There still remains a discontinuity in the derivative (the forces), however, this does not significantly affect the dynamics.

Interactions between dissimilar particles requires the generation of
cross term parameters for $\sigma$ and $\epsilon$. These are
calculated through the Lorentz-Berthelot mixing
rules:\cite{allen87:csl}
\begin{equation}
\sigma_{ij} = \frac{1}{2}[\sigma_{ii} + \sigma_{jj}]
\label{eq:sigmaMix}
\end{equation}
and
\begin{equation}
\epsilon_{ij} = \sqrt{\epsilon_{ii} \epsilon_{jj}}
\label{eq:epsilonMix}
\end{equation}

\subsection{\label{oopseSec:DUFF}Dipolar Unified-Atom Force Field}

The dipolar unified-atom force field ({\sc duff}) was developed to
simulate lipid bilayers. The simulations require a model capable of
forming bilayers, while still being sufficiently computationally
efficient to allow large systems ($\sim$100's of phospholipids,
$\sim$1000's of waters) to be simulated for long times
($\sim$10's of nanoseconds).

With this goal in mind, {\sc duff} has no point
charges. Charge-neutral distributions were replaced with dipoles,
while most atoms and groups of atoms were reduced to Lennard-Jones
interaction sites. This simplification cuts the length scale of long
range interactions from $\frac{1}{r}$ to $\frac{1}{r^3}$, and allows
us to avoid the computationally expensive Ewald sum. Instead, we can
use neighbor-lists and cutoff radii for the dipolar interactions, or
include a reaction field to mimic larger range interactions.

As an example, lipid head-groups in {\sc duff} are represented as
point dipole interaction sites. By placing a dipole at the head
group's center of mass, our model mimics the charge separation found
in common phospholipid head groups such as
phosphatidylcholine.\cite{Cevc87} Additionally, a large Lennard-Jones
site is located at the pseudoatom's center of mass. The model is
illustrated by the red atom in Fig.~\ref{oopseFig:lipidModel}. The
water model we use to complement the dipoles of the lipids is our
reparameterization of the soft sticky dipole (SSD) model of Ichiye
\emph{et al.}\cite{liu96:new_model}

\begin{figure}
\centering
\includegraphics[width=\linewidth]{twoChainFig.eps}
\caption[A representation of a lipid model in {\sc duff}]{A representation of the lipid model. $\phi$ is the torsion angle, $\theta$ %
is the bend angle, and $\mu$ is the dipole moment of the head group.}
\label{oopseFig:lipidModel}
\end{figure}

We have used a set of scalable parameters to model the alkyl groups
with Lennard-Jones sites. For this, we have borrowed parameters from
the TraPPE force field of Siepmann
\emph{et al}.\cite{Siepmann1998} TraPPE is a unified-atom
representation of n-alkanes, which is parametrized against phase
equilibria using Gibbs ensemble Monte Carlo simulation
techniques.\cite{Siepmann1998} One of the advantages of TraPPE is that
it generalizes the types of atoms in an alkyl chain to keep the number
of pseudoatoms to a minimum; the parameters for a unified atom such as
$\text{CH}_2$ do not change depending on what species are bonded to
it.

TraPPE also constrains all bonds to be of fixed length. Typically,
bond vibrations are the fastest motions in a molecular dynamic
simulation. Small time steps between force evaluations must be used to
ensure adequate energy conservation in the bond degrees of freedom. By
constraining the bond lengths, larger time steps may be used when
integrating the equations of motion. A simulation using {\sc duff} is
illustrated in Scheme \ref{sch:DUFF}.

\begin{lstlisting}[float,caption={[Invocation of {\sc duff}]A portion of a \texttt{.bass} file showing a simulation utilizing {\sc duff}},label={sch:DUFF}]

#include "water.mdl"
#include "lipid.mdl"

nComponents = 2;
component{
  type = "simpleLipid_16";
  nMol = 60;
}

component{
  type = "SSD_water";
  nMol = 1936;
}

initialConfig = "bilayer.init";

forceField = "DUFF";

\end{lstlisting}

\subsection{\label{oopseSec:energyFunctions}{\sc duff} Energy Functions}

The total potential energy function in {\sc duff} is
\begin{equation}
V = \sum^{N}_{I=1} V^{I}_{\text{Internal}}
        + \sum^{N-1}_{I=1} \sum_{J>I} V^{IJ}_{\text{Cross}}
\label{eq:totalPotential}
\end{equation}
Where $V^{I}_{\text{Internal}}$ is the internal potential of molecule $I$:
\begin{equation}
 V^{I}_{\text{Internal}} = 
        \sum_{\theta_{ijk} \in I} V_{\text{bend}}(\theta_{ijk})
        + \sum_{\phi_{ijkl} \in I} V_{\text{torsion}}(\phi_{ijkl})
        + \sum_{i \in I} \sum_{(j>i+4) \in I} 
        \biggl[ V_{\text{LJ}}(r_{ij}) +  V_{\text{dipole}}
        (\mathbf{r}_{ij},\boldsymbol{\Omega}_{i},\boldsymbol{\Omega}_{j})
        \biggr]
\label{eq:internalPotential}
\end{equation}
Here $V_{\text{bend}}$ is the bend potential for all 1, 3 bonded pairs
within the molecule $I$, and $V_{\text{torsion}}$ is the torsion potential
for all 1, 4 bonded pairs. The pairwise portions of the internal
potential are excluded for atom pairs that are involved in the same bond, bend, or torsion. All other atom pairs within the molecule are subject to the LJ pair potential.


The bend potential of a molecule is represented by the following function:
\begin{equation}
V_{\text{bend}}(\theta_{ijk}) = k_{\theta}( \theta_{ijk} - \theta_0 )^2 \label{eq:bendPot}
\end{equation}
Where $\theta_{ijk}$ is the angle defined by atoms $i$, $j$, and $k$
(see Fig.~\ref{oopseFig:lipidModel}), $\theta_0$ is the equilibrium
bond angle, and $k_{\theta}$ is the force constant which determines the
strength of the harmonic bend. The parameters for $k_{\theta}$ and
$\theta_0$ are borrowed from those in TraPPE.\cite{Siepmann1998}

The torsion potential and parameters are also borrowed from TraPPE. It is
of the form:
\begin{equation}
V_{\text{torsion}}(\phi) = c_1[1 + \cos \phi] 
        + c_2[1 + \cos(2\phi)] 
        + c_3[1 + \cos(3\phi)]
\label{eq:origTorsionPot}
\end{equation}
Where:
\begin{equation}
\cos\phi = (\hat{\mathbf{r}}_{ij} \times \hat{\mathbf{r}}_{jk}) \cdot
        (\hat{\mathbf{r}}_{jk} \times \hat{\mathbf{r}}_{kl})
\label{eq:torsPhi}
\end{equation}
Here, $\hat{\mathbf{r}}_{\alpha\beta}$ are the set of unit bond
vectors between atoms $i$, $j$, $k$, and $l$. For computational
efficiency, the torsion potential has been recast after the method of
{\sc charmm},\cite{Brooks83} in which the angle series is converted to
a power series of the form:
\begin{equation}
V_{\text{torsion}}(\phi) =  
        k_3 \cos^3 \phi + k_2 \cos^2 \phi + k_1 \cos \phi + k_0
\label{eq:torsionPot}
\end{equation}
Where:
\begin{align*}
k_0 &= c_1 + c_3 \\
k_1 &= c_1 - 3c_3 \\
k_2 &= 2 c_2 \\
k_3 &= 4c_3
\end{align*}
By recasting the potential as a power series, repeated trigonometric
evaluations are avoided during the calculation of the potential energy.


The cross potential between molecules $I$ and $J$, $V^{IJ}_{\text{Cross}}$, is
as follows:
\begin{equation}
V^{IJ}_{\text{Cross}} = 
        \sum_{i \in I} \sum_{j \in J}
        \biggl[ V_{\text{LJ}}(r_{ij}) +  V_{\text{dipole}}
        (\mathbf{r}_{ij},\boldsymbol{\Omega}_{i},\boldsymbol{\Omega}_{j})
        + V_{\text{sticky}}
        (\mathbf{r}_{ij},\boldsymbol{\Omega}_{i},\boldsymbol{\Omega}_{j})
        \biggr]
\label{eq:crossPotentail}
\end{equation}
Where $V_{\text{LJ}}$ is the Lennard Jones potential,
$V_{\text{dipole}}$ is the dipole dipole potential, and
$V_{\text{sticky}}$ is the sticky potential defined by the SSD model
(Sec.~\ref{oopseSec:SSD}). Note that not all atom types include all
interactions.

The dipole-dipole potential has the following form:
\begin{equation}
V_{\text{dipole}}(\mathbf{r}_{ij},\boldsymbol{\Omega}_{i},
        \boldsymbol{\Omega}_{j}) = \frac{|\mu_i||\mu_j|}{4\pi\epsilon_{0}r_{ij}^{3}} \biggl[
        \boldsymbol{\hat{u}}_{i} \cdot \boldsymbol{\hat{u}}_{j}
        -
        3(\boldsymbol{\hat{u}}_i \cdot \hat{\mathbf{r}}_{ij}) %
                (\boldsymbol{\hat{u}}_j \cdot \hat{\mathbf{r}}_{ij}) \biggr]
\label{eq:dipolePot}
\end{equation}
Here $\mathbf{r}_{ij}$ is the vector starting at atom $i$ pointing
towards $j$, and $\boldsymbol{\Omega}_i$ and $\boldsymbol{\Omega}_j$
are the orientational degrees of freedom for atoms $i$ and $j$
respectively. $|\mu_i|$ is the magnitude of the dipole moment of atom
$i$, $\boldsymbol{\hat{u}}_i$ is the standard unit orientation vector
of $\boldsymbol{\Omega}_i$, and $\boldsymbol{\hat{r}}_{ij}$ is the
unit vector pointing along $\mathbf{r}_{ij}$
($\boldsymbol{\hat{r}}_{ij}=\mathbf{r}_{ij}/|\mathbf{r}_{ij}|$).

To improve computational efficiency of the dipole-dipole interactions,
{\sc oopse} employs an electrostatic cutoff radius. This parameter can
be set in the \texttt{.bass} file, and controls the length scale over
which dipole interactions are felt. To compensate for the
discontinuity in the potential and the forces at the cutoff radius, we
have implemented a switching function to smoothly scale the
dipole-dipole interaction at the cutoff.
\begin{equation}
S(r_{ij}) = 
        \begin{cases}
        1 & \text{if $r_{ij} \le r_t$},\\
        \frac{(r_{\text{cut}} + 2r_{ij} - 3r_t)(r_{\text{cut}} - r_{ij})^2}
        {(r_{\text{cut}} - r_t)^2} 
        & \text{if $r_t < r_{ij} \le r_{\text{cut}}$}, \\
        0 & \text{if $r_{ij} > r_{\text{cut}}$.}
        \end{cases}
\label{eq:dipoleSwitching}
\end{equation}
Here $S(r_{ij})$ scales the potential at a given $r_{ij}$, and $r_t$
is the taper radius some given thickness less than the electrostatic
cutoff. The switching thickness can be set in the \texttt{.bass} file.

\subsection{\label{oopseSec:SSD}The {\sc duff} Water Models: SSD/E and SSD/RF}

In the interest of computational efficiency, the default solvent used
by {\sc oopse} is the extended Soft Sticky Dipole (SSD/E) water
model.\cite{Gezelter04} The original SSD was developed by Ichiye
\emph{et al.}\cite{liu96:new_model} as a modified form of the hard-sphere 
water model proposed by Bratko, Blum, and
Luzar.\cite{Bratko85,Bratko95} It consists of a single point dipole
with a Lennard-Jones core and a sticky potential that directs the
particles to assume the proper hydrogen bond orientation in the first
solvation shell. Thus, the interaction between two SSD water molecules
\emph{i} and \emph{j} is given by the potential
\begin{equation}
V_{ij} = 
        V_{ij}^{LJ} (r_{ij})\ + V_{ij}^{dp}
        (\mathbf{r}_{ij},\boldsymbol{\Omega}_i,\boldsymbol{\Omega}_j)\ +
        V_{ij}^{sp}
        (\mathbf{r}_{ij},\boldsymbol{\Omega}_i,\boldsymbol{\Omega}_j),
\label{eq:ssdPot}
\end{equation}
where the $\mathbf{r}_{ij}$ is the position vector between molecules
\emph{i} and \emph{j} with magnitude equal to the distance $r_{ij}$, and
$\boldsymbol{\Omega}_i$ and $\boldsymbol{\Omega}_j$ represent the
orientations of the respective molecules. The Lennard-Jones and dipole
parts of the potential are given by equations \ref{eq:lennardJonesPot}
and \ref{eq:dipolePot} respectively. The sticky part is described by
the following,
\begin{equation}
u_{ij}^{sp}(\mathbf{r}_{ij},\boldsymbol{\Omega}_i,\boldsymbol{\Omega}_j)=
        \frac{\nu_0}{2}[s(r_{ij})w(\mathbf{r}_{ij},
        \boldsymbol{\Omega}_i,\boldsymbol{\Omega}_j) +
        s^\prime(r_{ij})w^\prime(\mathbf{r}_{ij},
        \boldsymbol{\Omega}_i,\boldsymbol{\Omega}_j)]\ ,
\label{eq:stickyPot}
\end{equation}
where $\nu_0$ is a strength parameter for the sticky potential, and
$s$ and $s^\prime$ are cubic switching functions which turn off the
sticky interaction beyond the first solvation shell. The $w$ function
can be thought of as an attractive potential with tetrahedral
geometry:
\begin{equation}
w({\bf r}_{ij},{\bf \Omega}_i,{\bf \Omega}_j)=
        \sin\theta_{ij}\sin2\theta_{ij}\cos2\phi_{ij},
\label{eq:stickyW}
\end{equation}
while the $w^\prime$ function counters the normal aligned and
anti-aligned structures favored by point dipoles:
\begin{equation}
w^\prime({\bf r}_{ij},{\bf \Omega}_i,{\bf \Omega}_j)=
        (\cos\theta_{ij}-0.6)^2(\cos\theta_{ij}+0.8)^2-w^0,
\label{eq:stickyWprime}
\end{equation}
It should be noted that $w$ is proportional to the sum of the $Y_3^2$
and $Y_3^{-2}$ spherical harmonics (a linear combination which
enhances the tetrahedral geometry for hydrogen bonded structures),
while $w^\prime$ is a purely empirical function.  A more detailed
description of the functional parts and variables in this potential
can be found in the original SSD
articles.\cite{liu96:new_model,liu96:monte_carlo,chandra99:ssd_md,Ichiye03}

Since SSD/E is a single-point {\it dipolar} model, the force
calculations are simplified significantly relative to the standard
{\it charged} multi-point models. In the original Monte Carlo
simulations using this model, Ichiye {\it et al.} reported that using
SSD decreased computer time by a factor of 6-7 compared to other
models.\cite{liu96:new_model} What is most impressive is that these savings
did not come at the expense of accurate depiction of the liquid state
properties.  Indeed, SSD/E maintains reasonable agreement with the Head-Gordon
diffraction data for the structural features of liquid
water.\cite{hura00,liu96:new_model} Additionally, the dynamical properties
exhibited by SSD/E agree with experiment better than those of more
computationally expensive models (like TIP3P and
SPC/E).\cite{chandra99:ssd_md} The combination of speed and accurate depiction
of solvent properties makes SSD/E a very attractive model for the
simulation of large scale biochemical simulations.

Recent constant pressure simulations revealed issues in the original
SSD model that led to lower than expected densities at all target
pressures.\cite{Ichiye03,Gezelter04} The default model in {\sc oopse}
is therefore SSD/E, a density corrected derivative of SSD that
exhibits improved liquid structure and transport behavior. If the use
of a reaction field long-range interaction correction is desired, it
is recommended that the parameters be modified to those of the SSD/RF
model (an SSD variant  parameterized for reaction field). Solvent parameters can be easily modified in an accompanying
\texttt{.bass} file as illustrated in the scheme below. A table of the
parameter values and the drawbacks and benefits of the different
density corrected SSD models can be found in
reference~\cite{Gezelter04}.

\begin{lstlisting}[float,caption={[A simulation of {\sc ssd} water]A portion of a \texttt{.bass} file showing a simulation including {\sc ssd} water.},label={sch:ssd}]

#include "water.mdl"

nComponents = 1;
component{
  type = "SSD_water";
  nMol = 864;
}

initialConfig = "liquidWater.init";

forceField = "DUFF";

/*
 * The following two flags set the cutoff 
 * radius for the electrostatic forces 
 * as well as the skin thickness of the switching
 * function.
 */

electrostaticCutoffRadius  = 9.2; 
electrostaticSkinThickness = 1.38;

\end{lstlisting}


\subsection{\label{oopseSec:eam}Embedded Atom Method}

There are Molecular Dynamics packages which have the
capacity to simulate metallic systems, including some that have
parallel computational abilities\cite{plimpton93}. Potentials that
describe bonding transition metal
systems\cite{Finnis84,Ercolessi88,Chen90,Qi99,Ercolessi02} have an
attractive interaction which models  ``Embedding''
a positively charged metal ion in the electron density due to the
free valance ``sea'' of electrons created by the surrounding atoms in
the system. A mostly-repulsive pairwise part of the potential
describes the interaction of the positively charged metal core ions
with one another. A particular potential description called the
Embedded Atom Method\cite{Daw84,FBD86,johnson89,Lu97}({\sc eam}) that has
particularly wide adoption has been selected for inclusion in {\sc oopse}. A
good review of {\sc eam} and other metallic potential formulations was written
by Voter.\cite{voter}

The {\sc eam} potential has the form:
\begin{eqnarray}
V & = & \sum_{i} F_{i}\left[\rho_{i}\right] + \sum_{i} \sum_{j \neq i}
\phi_{ij}({\bf r}_{ij})  \\
\rho_{i}  & = & \sum_{j \neq i} f_{j}({\bf r}_{ij})
\end{eqnarray}
where $F_{i} $ is the embedding function that equates the energy
required to embed a positively-charged core ion $i$ into a linear
superposition of spherically averaged atomic electron densities given
by $\rho_{i}$.  $\phi_{ij}$ is a primarily repulsive pairwise
interaction between atoms $i$ and $j$. In the original formulation of
{\sc eam}\cite{Daw84}, $\phi_{ij}$ was an entirely repulsive term,
however in later refinements to {\sc eam} have shown that non-uniqueness
between $F$ and $\phi$ allow for more general forms for
$\phi$.\cite{Daw89} There is a cutoff distance, $r_{cut}$, which
limits the summations in the {\sc eam} equation to the few dozen atoms
surrounding atom $i$ for both the density $\rho$ and pairwise $\phi$
interactions. Foiles \emph{et al}.~fit {\sc eam} potentials for the fcc
metals Cu, Ag, Au, Ni, Pd, Pt and alloys of these metals.\cite{FBD86}
These fits are included in {\sc oopse}.

\subsection{\label{oopseSec:pbc}Periodic Boundary Conditions} 

\newcommand{\roundme}{\operatorname{round}}

\textit{Periodic boundary conditions} are widely used to simulate bulk properties with a relatively small number of particles. The
simulation box is replicated throughout space to form an infinite
lattice.  During the simulation, when a particle moves in the primary
cell, its image in other cells move in exactly the same direction with
exactly the same orientation. Thus, as a particle leaves the primary
cell, one of its images will enter through the opposite face. If the
simulation box is large enough to avoid ``feeling'' the symmetries of
the periodic lattice, surface effects can be ignored. The available
periodic cells in OOPSE are cubic, orthorhombic and parallelepiped. We
use a $3 \times 3$ matrix, $\mathsf{H}$, to describe the shape and
size of the simulation box. $\mathsf{H}$ is defined:
\begin{equation}
\mathsf{H} = ( \mathbf{h}_x, \mathbf{h}_y, \mathbf{h}_z )
\end{equation}
Where $\mathbf{h}_{\alpha}$ is the column vector of the $\alpha$ axis of the
box.  During the course of the simulation both the size and shape of
the box can be changed to allow volume fluctuations when constraining
the pressure.

A real space vector, $\mathbf{r}$ can be transformed in to a box space
vector, $\mathbf{s}$, and back through the following transformations:
\begin{align}
\mathbf{s} &= \mathsf{H}^{-1} \mathbf{r} \\
\mathbf{r} &= \mathsf{H} \mathbf{s}
\end{align}
The vector $\mathbf{s}$ is now a vector expressed as the number of box
lengths in the $\mathbf{h}_x$, $\mathbf{h}_y$, and $\mathbf{h}_z$
directions. To find the minimum image of a vector $\mathbf{r}$, we
first convert it to its corresponding vector in box space, and then,
cast each element to lie in the range $[-0.5,0.5]$:
\begin{equation}
s_{i}^{\prime}=s_{i}-\roundme(s_{i})
\end{equation}
Where $s_i$ is the $i$th element of $\mathbf{s}$, and
$\roundme(s_i)$ is given by
\begin{equation}
\roundme(x) =
        \begin{cases}
        \lfloor x+0.5 \rfloor & \text{if $x \ge 0$} \\
        \lceil x-0.5 \rceil & \text{if $x < 0$ }
        \end{cases}
\end{equation}
Here $\lfloor x \rfloor$ is the floor operator, and gives the largest
integer value that is not greater than $x$, and $\lceil x \rceil$ is
the ceiling operator, and gives the smallest integer that is not less
than $x$.  For example, $\roundme(3.6)=4$, $\roundme(3.1)=3$,
$\roundme(-3.6)=-4$, $\roundme(-3.1)=-3$.

Finally, we obtain the minimum image coordinates $\mathbf{r}^{\prime}$ by
transforming back to real space,
\begin{equation}
\mathbf{r}^{\prime}=\mathsf{H}^{-1}\mathbf{s}^{\prime}%
\end{equation}
In this way, particles are allowed to diffuse freely in $\mathbf{r}$,
but their minimum images, $\mathbf{r}^{\prime}$ are used to compute
the inter-atomic forces.


\section{\label{oopseSec:IOfiles}Input and Output Files}

\subsection{{\sc bass} and Model Files}

Every {\sc oopse} simulation begins with a Bizarre Atom Simulation
Syntax ({\sc bass}) file. {\sc bass} is a script syntax that is parsed
by {\sc oopse} at runtime. The {\sc bass} file allows for the user to
completely describe the system they wish to simulate, as well as tailor
{\sc oopse}'s behavior during the simulation. {\sc bass} files are
denoted with the extension
\texttt{.bass}, an example file is shown in
Scheme~\ref{sch:bassExample}.

\begin{lstlisting}[float,caption={[An example of a complete {\sc bass} file] An example showing a complete {\sc bass} file.},label={sch:bassExample}]

molecule{
  name = "Ar";
  nAtoms = 1;
  atom[0]{
    type="Ar";
    position( 0.0, 0.0, 0.0 );
  }
}

nComponents = 1;
component{
  type = "Ar";
  nMol = 108;
}

initialConfig = "./argon.init";

forceField = "LJ";
ensemble = "NVE"; // specify the simulation ensemble
dt = 1.0;         // the time step for integration
runTime = 1e3;    // the total simulation run time
sampleTime = 100; // trajectory file frequency
statusTime = 50;  // statistics file frequency

\end{lstlisting}

Within the \texttt{.bass} file it is necessary to provide a complete
description of the molecule before it is actually placed in the
simulation. The {\sc bass} syntax was originally developed with this
goal in mind, and allows for the specification of all the atoms in a
molecular prototype, as well as any bonds, bends, or torsions. These
descriptions can become lengthy for complex molecules, and it would be
inconvenient to duplicate the simulation at the beginning of each {\sc
bass} script. Addressing this issue {\sc bass} allows for the
inclusion of model files at the top of a \texttt{.bass} file. These
model files, denoted with the \texttt{.mdl} extension, allow the user
to describe a molecular prototype once, then simply include it into
each simulation containing that molecule. Returning to the example in
Scheme~\ref{sch:bassExample}, the \texttt{.mdl} file's contents would
be Scheme~\ref{sch:mdlExample}, and the new \texttt{.bass} file would
become Scheme~\ref{sch:bassExPrime}.

\begin{lstlisting}[float,caption={An example \texttt{.mdl} file.},label={sch:mdlExample}]

molecule{
  name = "Ar";
  nAtoms = 1;
  atom[0]{
    type="Ar";
    position( 0.0, 0.0, 0.0 );
  }
}

\end{lstlisting}

\begin{lstlisting}[float,caption={Revised {\sc bass} example.},label={sch:bassExPrime}]

#include "argon.mdl"

nComponents = 1;
component{
  type = "Ar";
  nMol = 108;
}

initialConfig = "./argon.init";

forceField = "LJ";
ensemble = "NVE";
dt = 1.0;
runTime = 1e3;
sampleTime = 100;
statusTime = 50; 

\end{lstlisting}

\subsection{\label{oopseSec:coordFiles}Coordinate Files}

The standard format for storage of a systems coordinates is a modified
xyz-file syntax, the exact details of which can be seen in
Scheme~\ref{sch:dumpFormat}. As all bonding and molecular information
is stored in the \texttt{.bass} and \texttt{.mdl} files, the
coordinate files are simply the complete set of coordinates for each
atom at a given simulation time. One important note, although the
simulation propagates the complete rotation matrix, directional
entities are written out using quanternions, to save space in the
output files.

\begin{lstlisting}[float,caption={[The format of the coordinate files]Shows the format of the coordinate files. The fist line is the number of atoms. The second line begins with the time stamp followed by the three $\mathsf{H}$ column vectors. It is important to note, that for extended system ensembles, additional information pertinent to the integrators may be stored on this line as well. The next lines are the atomic coordinates for all atoms in the system. First is the name followed by position, velocity, quanternions, and lastly angular velocities.},label=sch:dumpFormat]

nAtoms
time; Hxx Hyx Hzx; Hxy Hyy Hzy; Hxz Hyz Hzz;
Name1 x y z vx vy vz q0 q1 q2 q3 jx jy jz
Name2 x y z vx vy vz q0 q1 q2 q3 jx jy jz
etc...

\end{lstlisting}


There are three major files used by {\sc oopse} written in the
coordinate format, they are as follows: the initialization file
(\texttt{.init}), the simulation trajectory file (\texttt{.dump}), and
the final coordinates of the simulation. The initialization file is
necessary for {\sc oopse} to start the simulation with the proper
coordinates, and is generated before the simulation run. The
trajectory file is created at the beginning of the simulation, and is
used to store snapshots of the simulation at regular intervals. The
first frame is a duplication of the
\texttt{.init} file, and each subsequent frame is appended to the file
at an interval specified in the \texttt{.bass} file with the
\texttt{sampleTime} flag. The final coordinate file is the end of run file. The
\texttt{.eor} file stores the final configuration of the system for a
given simulation. The file is updated at the same time as the
\texttt{.dump} file, however, it only contains the most recent
frame. In this way, an \texttt{.eor} file may be used as the
initialization file to a second simulation in order to continue a
simulation or recover one from a processor that has crashed during the
course of the run.

\subsection{\label{oopseSec:initCoords}Generation of Initial Coordinates}

As was stated in Sec.~\ref{oopseSec:coordFiles}, an initialization
file is needed to provide the starting coordinates for a
simulation. The {\sc oopse} package provides several system building
programs to aid in the creation of the \texttt{.init}
file. The programs use {\sc bass}, and will recognize
arguments and parameters in the \texttt{.bass} file that would
otherwise be ignored by the simulation.

\subsection{The Statistics File}

The last output file generated by {\sc oopse} is the statistics
file. This file records such statistical quantities as the
instantaneous temperature, volume, pressure, etc. It is written out
with the frequency specified in the \texttt{.bass} file with the
\texttt{statusTime} keyword. The file allows the user to observe the
system variables as a function of simulation time while the simulation
is in progress. One useful function the statistics file serves is to
monitor the conserved quantity of a given simulation ensemble, this
allows the user to observe the stability of the integrator. The
statistics file is denoted with the \texttt{.stat} file extension.

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

\subsection{\label{oopseSec: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{Frenkel1996}

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{allen87:csl} 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$,
$\overleftrightarrow{\mathsf{I}}_i$ are the body-fixed angular
momentum and moment of inertia tensor respectively, and the
superscript $T$ denotes the transpose of the vector.  $\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 $\mathcal{O}\left(h^4\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, where $h= \delta t$:

{\tt moveA:}
\begin{align*}
{\bf v}\left(t + h / 2\right)  &\leftarrow  {\bf v}(t) 
        + \frac{h}{2} \left( {\bf f}(t) / m \right) \\
%
{\bf r}(t + h) &\leftarrow {\bf r}(t) 
        + h  {\bf v}\left(t + h / 2 \right) \\
%
{\bf j}\left(t + h / 2 \right)  &\leftarrow {\bf j}(t) 
        + \frac{h}{2} {\bf \tau}^b(t)  \\
%
\mathsf{A}(t + h) &\leftarrow \mathrm{rotate}\left( h {\bf j}
        (t + h / 2) \cdot \overleftrightarrow{\mathsf{I}}^{-1} \right)
\end{align*}

In this context, the $\mathrm{rotate}$ function is the reversible product
of the three body-fixed rotations,
\begin{equation}
\mathrm{rotate}({\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{align*}
{\bf f}(t + h) &\leftarrow  
        - \left(\frac{\partial V}{\partial {\bf r}}\right)_{{\bf r}(t + h)} \\
%
{\bf \tau}^{s}(t + h) &\leftarrow {\bf u}(t + h)
        \times \frac{\partial V}{\partial {\bf u}} \\
%
{\bf \tau}^{b}(t + h) &\leftarrow \mathsf{A}(t + h)
        \cdot {\bf \tau}^s(t + h) 
\end{align*}

{\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{align*}
{\bf v}\left(t + h \right)  &\leftarrow  {\bf v}\left(t + h / 2 \right) 
        + \frac{h}{2} \left( {\bf f}(t + h) / m \right) \\
%
{\bf j}\left(t + h \right)  &\leftarrow {\bf j}\left(t + h / 2 \right) 
        + \frac{h}{2} {\bf \tau}^b(t + h)  
\end{align*}

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
Fig.~\ref{timestep}.

\begin{figure}
\centering
\includegraphics[width=\linewidth]{timeStep.eps}
\caption[Energy conservation for quaternion versus DLM dynamics]{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 Fig.~\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} \\  
$h$ & {\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 & {\tt ensemble = NPTi;} \\
  &  (with isotropic volume changes) & \\
NPTf & isobaric-isothermal & {\tt ensemble = NPTf;} \\
  & (with changes to box shape) & \\
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\cite{Hoover85} 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.\cite{melchionna93}

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.\cite{melchionna93} 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{center}
\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 &
true
\end{tabular}
\end{center}

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.

\subsection{\label{oopseSec: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{align*}
T(t) &\leftarrow \left\{{\bf v}(t)\right\}, \left\{{\bf j}(t)\right\} \\
%
{\bf v}\left(t + h / 2\right)  &\leftarrow {\bf v}(t) 
        + \frac{h}{2} \left( \frac{{\bf f}(t)}{m} - {\bf v}(t)
        \chi(t)\right) \\
%
{\bf r}(t + h) &\leftarrow {\bf r}(t) 
        + h {\bf v}\left(t + h / 2 \right) \\
%
{\bf j}\left(t + h / 2 \right)  &\leftarrow {\bf j}(t) 
        + \frac{h}{2} \left( {\bf \tau}^b(t) - {\bf j}(t)
        \chi(t) \right) \\
%
\mathsf{A}(t + h) &\leftarrow \mathrm{rotate}
        \left(h * {\bf j}(t + h / 2) 
        \overleftrightarrow{\mathsf{I}}^{-1} \right) \\
%
\chi\left(t + h / 2 \right) &\leftarrow \chi(t) 
        + \frac{h}{2 \tau_T^2} \left( \frac{T(t)}
        {T_{\mathrm{target}}} - 1 \right) 
\end{align*}

Here $\mathrm{rotate}(h * {\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{align*}
T(t + h) &\leftarrow \left\{{\bf v}(t + h)\right\}, 
        \left\{{\bf j}(t + h)\right\} \\
%
\chi\left(t + h \right) &\leftarrow \chi\left(t + h /
        2 \right) + \frac{h}{2 \tau_T^2} \left( \frac{T(t+h)}
        {T_{\mathrm{target}}} - 1 \right) \\
%
{\bf v}\left(t + h \right)  &\leftarrow {\bf v}\left(t 
        + h / 2 \right) + \frac{h}{2} \left(
        \frac{{\bf f}(t + h)}{m} - {\bf v}(t + h)
        \chi(t h)\right) \\
%
{\bf j}\left(t + h \right) &\leftarrow {\bf j}\left(t
        + h / 2 \right) + \frac{h}{2} 
        \left( {\bf \tau}^b(t + h) - {\bf j}(t + h) 
        \chi(t + h) \right) 
\end{align*}

Since ${\bf v}(t + h)$ and ${\bf j}(t + h)$ are required to caclculate
$T(t + h)$ as well as $\chi(t + h)$, they indirectly depend on their
own values at time $t + h$.  {\tt moveB} is therefore done in an
iterative fashion until $\chi(t + h)$ 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,\cite{melchionna93}
\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 $h$ 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{oopseSec:rattle}.

\subsection{\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{align*}
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\} \\
%
{\bf v}\left(t + h / 2\right)  &\leftarrow {\bf v}(t) 
        + \frac{h}{2} \left( \frac{{\bf f}(t)}{m} - {\bf v}(t)
        \left(\chi(t) + \eta(t) \right) \right) \\
%
{\bf j}\left(t + h / 2 \right)  &\leftarrow {\bf j}(t) 
        + \frac{h}{2} \left( {\bf \tau}^b(t) - {\bf j}(t)
        \chi(t) \right) \\
%
\mathsf{A}(t + h) &\leftarrow \mathrm{rotate}\left(h *
        {\bf j}(t + h / 2) \overleftrightarrow{\mathsf{I}}^{-1} 
        \right) \\
%
\chi\left(t + h / 2 \right) &\leftarrow \chi(t) +
        \frac{h}{2 \tau_T^2} \left( \frac{T(t)}{T_{\mathrm{target}}} - 1
        \right) \\
%
\eta(t + h / 2) &\leftarrow \eta(t) + \frac{h 
        \mathcal{V}(t)}{2 N k_B T(t) \tau_B^2} \left( P(t) 
        - P_{\mathrm{target}} \right) \\ 
%
{\bf r}(t + h) &\leftarrow {\bf r}(t) + h 
        \left\{ {\bf v}\left(t + h / 2 \right) 
        + \eta(t + h / 2)\left[ {\bf r}(t + h) 
        - {\bf R}_0 \right] \right\} \\
%
\mathsf{H}(t + h) &\leftarrow e^{-h \eta(t + h / 2)} 
        \mathsf{H}(t)
\end{align*}

Most of these equations are identical to their counterparts in the NVT
integrator, but the propagation of positions to time $t + h$
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 +
h / 2$.  Reshaping the box uniformly also scales the volume of
the box by
\begin{equation}
\mathcal{V}(t + h) \leftarrow e^{ - 3 h \eta(t + h /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{align*}
T(t + h) &\leftarrow \left\{{\bf v}(t + h)\right\},
        \left\{{\bf j}(t + h)\right\} \\
%
P(t + h) &\leftarrow  \left\{{\bf r}(t + h)\right\},
        \left\{{\bf v}(t + h)\right\} \\
%
\chi\left(t + h \right) &\leftarrow \chi\left(t + h /
        2 \right) + \frac{h}{2 \tau_T^2} \left( \frac{T(t+h)}
        {T_{\mathrm{target}}} - 1 \right) \\
%
\eta(t + h) &\leftarrow \eta(t + h / 2) +
        \frac{h \mathcal{V}(t + h)}{2 N k_B T(t + h) 
        \tau_B^2} \left( P(t + h) - P_{\mathrm{target}} \right) \\ 
%
{\bf v}\left(t + h \right)  &\leftarrow {\bf v}\left(t 
        + h / 2 \right) + \frac{h}{2} \left(
        \frac{{\bf f}(t + h)}{m} - {\bf v}(t + h)
        (\chi(t + h) + \eta(t + h)) \right) \\
%
{\bf j}\left(t + h \right)  &\leftarrow {\bf j}\left(t 
        + h / 2 \right) + \frac{h}{2} \left( {\bf
        \tau}^b(t + h) - {\bf j}(t + h)
        \chi(t + h) \right) 
\end{align*}

Once again, since ${\bf v}(t + h)$ and ${\bf j}(t + h)$ are required
to caclculate $T(t + h)$, $P(t + h)$, $\chi(t + h)$, and $\eta(t +
h)$, they indirectly depend on their own values at time $t + h$.  {\tt
moveB} is therefore done in an iterative fashion until $\chi(t + h)$
and $\eta(t + h)$ 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{oopseSec:rattle}.

\subsection{\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 \cdot \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{align*}
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\} \\
%
{\bf v}\left(t + h / 2\right)  &\leftarrow {\bf v}(t) 
        + \frac{h}{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 + h / 2 \right)  &\leftarrow {\bf j}(t) 
        + \frac{h}{2} \left( {\bf \tau}^b(t) - {\bf j}(t)
        \chi(t) \right) \\
%
\mathsf{A}(t + h) &\leftarrow \mathrm{rotate}\left(h *
        {\bf j}(t + h / 2) \overleftrightarrow{\mathsf{I}}^{-1} 
        \right) \\
%
\chi\left(t + h / 2 \right) &\leftarrow \chi(t) +
        \frac{h}{2 \tau_T^2} \left( \frac{T(t)}{T_{\mathrm{target}}} 
        - 1 \right) \\
%
\overleftrightarrow{\eta}(t + h / 2) &\leftarrow 
        \overleftrightarrow{\eta}(t) + \frac{h \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 + h) &\leftarrow {\bf r}(t) + h \left\{ {\bf v}
        \left(t + h / 2 \right) + \overleftrightarrow{\eta}(t +
        h / 2) \cdot \left[ {\bf r}(t + h) 
        - {\bf R}_0 \right] \right\} \\
%
\mathsf{H}(t + h) &\leftarrow \mathsf{H}(t) \cdot e^{-h
        \overleftrightarrow{\eta}(t + h / 2)} 
\end{align*}
{\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{align*}
T(t + h) &\leftarrow \left\{{\bf v}(t + h)\right\},
        \left\{{\bf j}(t + h)\right\} \\
%
\overleftrightarrow{\mathsf{P}}(t + h) &\leftarrow \left\{{\bf r}
        (t + h)\right\}, \left\{{\bf v}(t 
        + h)\right\}, \left\{{\bf f}(t + h)\right\} \\
%
\chi\left(t + h \right) &\leftarrow \chi\left(t + h /
        2 \right) + \frac{h}{2 \tau_T^2} \left( \frac{T(t+
        h)}{T_{\mathrm{target}}} - 1 \right) \\
%
\overleftrightarrow{\eta}(t + h) &\leftarrow 
        \overleftrightarrow{\eta}(t + h / 2) +
        \frac{h \mathcal{V}(t + h)}{2 N k_B T(t + h) 
        \tau_B^2} \left( \overleftrightarrow{P}(t + h) 
        - P_{\mathrm{target}}\mathsf{1} \right) \\ 
%
{\bf v}\left(t + h \right)  &\leftarrow {\bf v}\left(t 
        + h / 2 \right) + \frac{h}{2} \left(
        \frac{{\bf f}(t + h)}{m} - 
        (\chi(t + h)\mathsf{1} + \overleftrightarrow{\eta}(t 
        + h)) \right) \cdot {\bf v}(t + h) \\
%
{\bf j}\left(t + h \right)  &\leftarrow {\bf j}\left(t 
        + h / 2 \right) + \frac{h}{2} \left( {\bf \tau}^b(t 
        + h) - {\bf j}(t + h) \chi(t + h) \right) 
\end{align*}

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.  The NPTf integrator
finds most use in simulating crystals or liquid crystals which assume
non-orthorhombic geometries.

\subsection{\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{oopseSec:rattle}The {\sc rattle} Method for Bond 
        Constraints}

In order to satisfy the constraints of fixed bond lengths within {\sc
oopse}, we have implemented the {\sc rattle} algorithm of
Andersen.\cite{andersen83} The algorithm is a velocity verlet
formulation of the {\sc shake} method\cite{ryckaert77} of iteratively
solving the Lagrange multipliers of constraint. The system of Lagrange
multipliers allows one to reformulate the equations of motion with
explicit constraint forces.\cite{fowles99:lagrange} 

Consider a system described by coordinates $q_1$ and $q_2$ subject to an
equation of constraint:
\begin{equation}
\sigma(q_1, q_2,t) = 0
\label{oopseEq:lm1}
\end{equation}
The Lagrange formulation of the equations of motion can be written:
\begin{equation}
\delta\int_{t_1}^{t_2}L\, dt = 
        \int_{t_1}^{t_2} \sum_i \biggl [ \frac{\partial L}{\partial q_i}
        - \frac{d}{dt}\biggl(\frac{\partial L}{\partial \dot{q}_i} 
        \biggr ) \biggr] \delta q_i \, dt = 0
\label{oopseEq:lm2}
\end{equation}
Here, $\delta q_i$ is not independent for each $q$, as $q_1$ and $q_2$
are linked by $\sigma$. However, $\sigma$ is fixed at any given
instant of time, giving:
\begin{align}
\delta\sigma &= \biggl( \frac{\partial\sigma}{\partial q_1} \delta q_1 %
        + \frac{\partial\sigma}{\partial q_2} \delta q_2 \biggr) = 0 \\
%
\frac{\partial\sigma}{\partial q_1} \delta q_1 &= %
        - \frac{\partial\sigma}{\partial q_2} \delta q_2 \\
%
\delta q_2 &= - \biggl(\frac{\partial\sigma}{\partial q_1} \bigg / %
        \frac{\partial\sigma}{\partial q_2} \biggr) \delta q_1
\end{align}
Substituted back into Eq.~\ref{oopseEq:lm2},
\begin{equation}
\int_{t_1}^{t_2}\biggl [ \biggl(\frac{\partial L}{\partial q_1}
        - \frac{d}{dt}\,\frac{\partial L}{\partial \dot{q}_1} 
        \biggr)
        - \biggl( \frac{\partial L}{\partial q_1}
        - \frac{d}{dt}\,\frac{\partial L}{\partial \dot{q}_1} 
        \biggr) \biggl(\frac{\partial\sigma}{\partial q_1} \bigg / %
        \frac{\partial\sigma}{\partial q_2} \biggr)\biggr] \delta q_1 \, dt = 0
\label{oopseEq:lm3}
\end{equation}
Leading to,
\begin{equation}
\frac{\biggl(\frac{\partial L}{\partial q_1}
        - \frac{d}{dt}\,\frac{\partial L}{\partial \dot{q}_1} 
        \biggr)}{\frac{\partial\sigma}{\partial q_1}} = 
\frac{\biggl(\frac{\partial L}{\partial q_2}
        - \frac{d}{dt}\,\frac{\partial L}{\partial \dot{q}_2} 
        \biggr)}{\frac{\partial\sigma}{\partial q_2}}
\label{oopseEq:lm4}
\end{equation}
This relation can only be statisfied, if both are equal to a single
function $-\lambda(t)$,
\begin{align}
\frac{\biggl(\frac{\partial L}{\partial q_1}
        - \frac{d}{dt}\,\frac{\partial L}{\partial \dot{q}_1} 
        \biggr)}{\frac{\partial\sigma}{\partial q_1}} &= -\lambda(t) \\
%
\frac{\partial L}{\partial q_1}
        - \frac{d}{dt}\,\frac{\partial L}{\partial \dot{q}_1} &= 
         -\lambda(t)\,\frac{\partial\sigma}{\partial q_1} \\
%
\frac{\partial L}{\partial q_1}
        - \frac{d}{dt}\,\frac{\partial L}{\partial \dot{q}_1} 
         + \mathcal{G}_i &= 0
\end{align}
Where $\mathcal{G}_i$, the force of constraint on $i$, is:
\begin{equation}
\mathcal{G}_i = \lambda(t)\,\frac{\partial\sigma}{\partial q_1}
\label{oopseEq:lm5}
\end{equation}

In a simulation, this would involve the solution of a set of $(m + n)$
number of equations. Where $m$ is the number of constraints, and $n$
is the number of constrained coordinates. In practice, this is not
done, as the matrix inversion necessary to solve the system of
equations would be very time consuming to solve. Additionally, the
numerical error in the solution of the set of $\lambda$'s would be
compounded by the error inherent in propagating by the Velocity Verlet
algorithm ($\Delta t^4$). The Verlet propagation error is negligible
in an unconstrained system, as one is interested in the statistics of
the run, and not that the run be numerically exact to the ``true''
integration. This relates back to the ergodic hypothesis that a time
integral of a valid trajectory will still give the correct ensemble
average. However, in the case of constraints, if the equations of
motion leave the ``true'' trajectory, they are departing from the
constrained surface. The method that is used, is to iteratively solve
for $\lambda(t)$ at each time step.

In {\sc rattle} the equations of motion are modified subject to the
following two constraints:
\begin{align}
\sigma_{ij}[\mathbf{r}(t)] \equiv 
        [ \mathbf{r}_i(t) - \mathbf{r}_j(t)]^2  - d_{ij}^2 &= 0 %
        \label{oopseEq:c1} \\
%
[\mathbf{\dot{r}}_i(t) - \mathbf{\dot{r}}_j(t)] \cdot 
        [\mathbf{r}_i(t) - \mathbf{r}_j(t)] &= 0 \label{oopseEq:c2}
\end{align}
Eq.~\ref{oopseEq:c1} is the set of bond constraints, where $d_{ij}$ is
the constrained distance between atom $i$ and
$j$. Eq.~\ref{oopseEq:c2} constrains the velocities of $i$ and $j$ to
be perpendicular to the bond vector, so that the bond can neither grow
nor shrink. The constrained dynamics equations become:
\begin{equation}
m_i \mathbf{\ddot{r}}_i = \mathbf{F}_i + \mathbf{\mathcal{G}}_i
\label{oopseEq:r1}
\end{equation}
Where,$\mathbf{\mathcal{G}}_i$ are the forces of constraint on $i$,
and are defined:
\begin{equation}
\mathbf{\mathcal{G}}_i = - \sum_j \lambda_{ij}(t)\,\nabla \sigma_{ij}
\label{oopseEq:r2}
\end{equation}

In Velocity Verlet, if $\Delta t = h$, the propagation can be written:
\begin{align}
\mathbf{r}_i(t+h) &=
        \mathbf{r}_i(t) + h\mathbf{\dot{r}}(t) + 
        \frac{h^2}{2m_i}\,\Bigl[ \mathbf{F}_i(t) + 
        \mathbf{\mathcal{G}}_{Ri}(t) \Bigr] \label{oopseEq:vv1} \\
%
\mathbf{\dot{r}}_i(t+h) &= 
        \mathbf{\dot{r}}_i(t) + \frac{h}{2m_i}
        \Bigl[ \mathbf{F}_i(t) + \mathbf{\mathcal{G}}_{Ri}(t) + 
        \mathbf{F}_i(t+h) + \mathbf{\mathcal{G}}_{Vi}(t+h) \Bigr] %
        \label{oopseEq:vv2}
\end{align}
Where:
\begin{align}
\mathbf{\mathcal{G}}_{Ri}(t) &= 
        -2 \sum_j \lambda_{Rij}(t) \mathbf{r}_{ij}(t) \\
%
\mathbf{\mathcal{G}}_{Vi}(t+h) &=
        -2 \sum_j \lambda_{Vij}(t+h) \mathbf{r}(t+h)
\end{align}
Next, define:
\begin{align}
g_{ij} &= h \lambda_{Rij}(t) \\
k_{ij} &= h \lambda_{Vij}(t+h) \\
\mathbf{q}_i &= \mathbf{\dot{r}}_i(t) + \frac{h}{2m_i} \mathbf{F}_i(t)
        - \frac{1}{m_i}\sum_j g_{ij}\mathbf{r}_{ij}(t)
\end{align}
Using these definitions, Eq.~\ref{oopseEq:vv1} and \ref{oopseEq:vv2}
can be rewritten as,
\begin{align}
\mathbf{r}_i(t+h) &= \mathbf{r}_i(t) + h \mathbf{q}_i \\
%
\mathbf{\dot{r}}(t+h) &= \mathbf{q}_i + \frac{h}{2m_i}\mathbf{F}_i(t+h)
        -\frac{1}{m_i}\sum_j k_{ij} \mathbf{r}_{ij}(t+h)
\end{align}

To integrate the equations of motion, the {\sc rattle} algorithm first
solves for $\mathbf{r}(t+h)$. Let,
\begin{equation}
\mathbf{q}_i = \mathbf{\dot{r}}(t) + \frac{h}{2m_i}\mathbf{F}_i(t)
\end{equation}
Here $\mathbf{q}_i$ corresponds to an initial unconstrained move. Next
pick a constraint $j$, and let,
\begin{equation}
\mathbf{s} = \mathbf{r}_i(t) + h\mathbf{q}_i(t) 
        - \mathbf{r}_j(t) + h\mathbf{q}_j(t)
\label{oopseEq:ra1}
\end{equation}
If
\begin{equation}
\Big| |\mathbf{s}|^2 - d_{ij}^2 \Big| > \text{tolerance},
\end{equation}
then the constraint is unsatisfied, and corrections are made to the
positions. First we define a test corrected configuration as,
\begin{align}
\mathbf{r}_i^T(t+h) = \mathbf{r}_i(t) + h\biggl[\mathbf{q}_i - 
        g_{ij}\,\frac{\mathbf{r}_{ij}(t)}{m_i} \biggr] \\
%
\mathbf{r}_j^T(t+h) = \mathbf{r}_j(t) + h\biggl[\mathbf{q}_j +
        g_{ij}\,\frac{\mathbf{r}_{ij}(t)}{m_j} \biggr] 
\end{align}
And we chose $g_{ij}$ such that, $|\mathbf{r}_i^T - \mathbf{r}_j^T|^2
= d_{ij}^2$. Solving the quadratic for $g_{ij}$ we obtain the
approximation,
\begin{equation}
g_{ij} = \frac{(s^2 - d^2)}{2h[\mathbf{s}\cdot\mathbf{r}_{ij}(t)]
        (\frac{1}{m_i} + \frac{1}{m_j})}
\end{equation}
Although not an exact solution for $g_{ij}$, as this is an iterative
scheme overall, the eventual solution will converge. With a trial
$g_{ij}$, the new $\mathbf{q}$'s become,
\begin{align}
\mathbf{q}_i &= \mathbf{q}^{\text{old}}_i - g_{ij}\,
        \frac{\mathbf{r}_{ij}(t)}{m_i} \\
%
\mathbf{q}_j &= \mathbf{q}^{\text{old}}_j + g_{ij}\,
        \frac{\mathbf{r}_{ij}(t)}{m_j} 
\end{align}
The whole algorithm is then repeated from Eq.~\ref{oopseEq:ra1} until
all constraints are satisfied.

The second step of {\sc rattle}, is to then update the velocities. The
step starts with,
\begin{equation}
\mathbf{\dot{r}}_i(t+h) = \mathbf{q}_i + \frac{h}{2m_i}\mathbf{F}_i(t+h)
\end{equation}
Next we pick a constraint $j$, and calculate the dot product $\ell$.
\begin{equation}
\ell = \mathbf{r}_{ij}(t+h) \cdot \mathbf{\dot{r}}_{ij}(t+h)
\label{oopseEq:rv1}
\end{equation}
Here if constraint Eq.~\ref{oopseEq:c2} holds, $\ell$ should be
zero. Therefore if $\ell$ is greater than some tolerance, then
corrections are made to the $i$ and $j$ velocities.
\begin{align}
\mathbf{\dot{r}}_i^T &= \mathbf{\dot{r}}_i(t+h) - k_{ij}
        \frac{\mathbf{\dot{r}}_{ij}(t+h)}{m_i} \\
%
\mathbf{\dot{r}}_j^T &= \mathbf{\dot{r}}_j(t+h) + k_{ij}
        \frac{\mathbf{\dot{r}}_{ij}(t+h)}{m_j}
\end{align}
Like in the previous step, we select a value for $k_{ij}$ such that
$\ell$ is zero.
\begin{equation}
k_{ij} = \frac{\ell}{d^2_{ij}(\frac{1}{m_i} + \frac{1}{m_j})}
\end{equation}
The test velocities, $\mathbf{\dot{r}}^T_i$ and
$\mathbf{\dot{r}}^T_j$, then replace their respective velocities, and
the algorithm is iterated from Eq.~\ref{oopseEq:rv1} until all
constraints are satisfied.


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

Based on the fluctuation-dissipation theorem, a force auto-correlation
method was developed by Roux and Karplus to investigate the dynamics
of ions inside ion channels.\cite{Roux91} The time-dependent friction
coefficient can be calculated from the deviation of the instantaneous
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 decays rapidly, the static friction
coefficient can be approximated by
\begin{equation}
\xi_{\text{static}}(z)=\int_{0}^{\infty}\langle\delta F(z,t)\delta F(z,0)\rangle dt
\end{equation}
Allowing diffusion constant to then be calculated through the
Einstein relation:\cite{Marrink94}
\begin{equation}
D(z)=\frac{k_{B}T}{\xi_{\text{static}}(z)}=\frac{(k_{B}T)^{2}}{\int_{0}^{\infty
}\langle\delta F(z,t)\delta F(z,0)\rangle dt}%
\end{equation}

The Z-Constraint method, which fixes the z coordinates of the
molecules with respect to the center of the mass of the system, has
been a method suggested to obtain the forces required for the force
auto-correlation calculation.\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 in {\sc oopse}. Instead of
resetting the coordinate, we reset the forces of z-constrained
molecules as well as subtract the total constraint forces from the
rest of the system after the force calculation at each time step.

After the force calculation, define $G_\alpha$ as
\begin{equation}
G_{\alpha} = \sum_i F_{\alpha i}
\label{oopseEq:zc1}
\end{equation}
Where $F_{\alpha i}$ is the force in the z direction of atom $i$ in
z-constrained molecule $\alpha$. The forces of the z constrained
molecule are then set to:
\begin{equation}
F_{\alpha i} = F_{\alpha i} - 
        \frac{m_{\alpha i} G_{\alpha}}{\sum_i m_{\alpha i}}
\end{equation}
Here, $m_{\alpha i}$ is the mass of atom $i$ in the z-constrained
molecule. Having rescaled the forces, the velocities must also be
rescaled to subtract out any center of mass velocity in the z
direction.
\begin{equation}
v_{\alpha i} = v_{\alpha i} -
        \frac{\sum_i m_{\alpha i} v_{\alpha i}}{\sum_i m_{\alpha i}}
\end{equation}
Where $v_{\alpha i}$ is the velocity of atom $i$ in the z direction.
Lastly, all of the accumulated z constrained forces must be subtracted
from the system to keep the system center of mass from drifting.
\begin{equation}
F_{\beta i} = F_{\beta i} - \frac{m_{\beta i} \sum_{\alpha} G_{\alpha}}
        {\sum_{\beta}\sum_i m_{\beta i}}
\end{equation}
Where $\beta$ are all of the unconstrained molecules in the
system. Similarly, the velocities of the unconstrained molecules must
also be scaled.
\begin{equation}
v_{\beta i} = v_{\beta i} + \sum_{\alpha}
        \frac{\sum_i m_{\alpha i} v_{\alpha i}}{\sum_i m_{\alpha i}}
\end{equation}

At the very beginning of the simulation, the molecules may not be at their
constrained positions. To move a z-constrained molecule to its specified
position, a simple harmonic potential is used
\begin{equation}
U(t)=\frac{1}{2}k_{\text{Harmonic}}(z(t)-z_{\text{cons}})^{2}%
\end{equation}
where $k_{\text{Harmonic}}$ is the harmonic force constant, $z(t)$ is the
current $z$ coordinate of the center of mass of the constrained molecule, and
$z_{\text{cons}}$ is the constrained position. The harmonic force operating
on the z-constrained molecule at time $t$ can be calculated by
\begin{equation}
F_{z_{\text{Harmonic}}}(t)=-\frac{\partial U(t)}{\partial z(t)}=
        -k_{\text{Harmonic}}(z(t)-z_{\text{cons}})
\end{equation}

\section{\label{oopseSec:props}Trajectory Analysis}

\subsection{\label{oopseSec:staticProps}Static Property Analysis}

The static properties of the trajectories are analyzed with the
program \texttt{staticProps}. The code is capable of calculating a
number of pair correlations between species A and B. Some of which
only apply to directional entities. The summary of pair correlations
can be found in Table~\ref{oopseTb:gofrs}

\begin{table}
\caption[The list of pair correlations in \texttt{staticProps}]{THE DIFFERENT PAIR CORRELATIONS IN \texttt{staticProps}}
\label{oopseTb:gofrs}
\begin{center}
\begin{tabular}{|l|c|c|}
\hline
Name      & Equation & Directional Atom \\ \hline
$g_{\text{AB}}(r)$              & Eq.~\ref{eq:gofr}         & neither \\ \hline
$g_{\text{AB}}(r, \cos \theta)$ & Eq.~\ref{eq:gofrCosTheta} & A \\ \hline
$g_{\text{AB}}(r, \cos \omega)$ & Eq.~\ref{eq:gofrCosOmega} & both \\ \hline
$g_{\text{AB}}(x, y, z)$        & Eq.~\ref{eq:gofrXYZ}      & neither \\ \hline
$\langle \cos \omega \rangle_{\text{AB}}(r)$ & Eq.~\ref{eq:cosOmegaOfR} &%
        both \\ \hline
\end{tabular}
\begin{minipage}{\linewidth}
\centering
\vspace{2mm}
The third column specifies which atom, if any, need be a directional entity.
\end{minipage}
\end{center}
\end{table}

The first pair correlation, $g_{\text{AB}}(r)$, is defined as follows:
\begin{equation}
g_{\text{AB}}(r) = \frac{V}{N_{\text{A}}N_{\text{B}}}\langle %%
        \sum_{i \in \text{A}} \sum_{j \in \text{B}} %%
        \delta( r - |\mathbf{r}_{ij}|) \rangle \label{eq:gofr}
\end{equation}
Where $\mathbf{r}_{ij}$ is the vector
\begin{equation*}
\mathbf{r}_{ij} = \mathbf{r}_j - \mathbf{r}_i \notag
\end{equation*}
and $\frac{V}{N_{\text{A}}N_{\text{B}}}$ normalizes the average over
the expected pair density at a given $r$.

The next two pair correlations, $g_{\text{AB}}(r, \cos \theta)$ and
$g_{\text{AB}}(r, \cos \omega)$, are similar in that they are both two
dimensional histograms. Both use $r$ for the primary axis then a
$\cos$ for the secondary axis ($\cos \theta$ for
Eq.~\ref{eq:gofrCosTheta} and $\cos \omega$ for
Eq.~\ref{eq:gofrCosOmega}). This allows for the investigator to
correlate alignment on directional entities. $g_{\text{AB}}(r, \cos
\theta)$ is defined as follows:
\begin{equation}
g_{\text{AB}}(r, \cos \theta) = \frac{V}{N_{\text{A}}N_{\text{B}}}\langle  
\sum_{i \in \text{A}} \sum_{j \in \text{B}}  
\delta( \cos \theta - \cos \theta_{ij}) 
\delta( r - |\mathbf{r}_{ij}|) \rangle
\label{eq:gofrCosTheta}
\end{equation}
Where
\begin{equation*}
\cos \theta_{ij} = \mathbf{\hat{i}} \cdot \mathbf{\hat{r}}_{ij}
\end{equation*}
Here $\mathbf{\hat{i}}$ is the unit directional vector of species $i$
and $\mathbf{\hat{r}}_{ij}$ is the unit vector associated with vector
$\mathbf{r}_{ij}$.

The second two dimensional histogram is of the form:
\begin{equation}
g_{\text{AB}}(r, \cos \omega) = 
        \frac{V}{N_{\text{A}}N_{\text{B}}}\langle 
        \sum_{i \in \text{A}} \sum_{j \in \text{B}} 
        \delta( \cos \omega - \cos \omega_{ij})
        \delta( r - |\mathbf{r}_{ij}|) \rangle \label{eq:gofrCosOmega}
\end{equation}
Here
\begin{equation*}
\cos \omega_{ij} = \mathbf{\hat{i}} \cdot \mathbf{\hat{j}}
\end{equation*}
Again, $\mathbf{\hat{i}}$ and $\mathbf{\hat{j}}$ are the unit
directional vectors of species $i$ and $j$.

The static analysis code is also cable of calculating a three
dimensional pair correlation of the form:
\begin{equation}\label{eq:gofrXYZ}
g_{\text{AB}}(x, y, z) = 
        \frac{V}{N_{\text{A}}N_{\text{B}}}\langle 
        \sum_{i \in \text{A}} \sum_{j \in \text{B}} 
        \delta( x - x_{ij})
        \delta( y - y_{ij})
        \delta( z - z_{ij}) \rangle
\end{equation}
Where $x_{ij}$, $y_{ij}$, and $z_{ij}$ are the $x$, $y$, and $z$
components respectively of vector $\mathbf{r}_{ij}$.

The final pair correlation is similar to
Eq.~\ref{eq:gofrCosOmega}. $\langle \cos \omega
\rangle_{\text{AB}}(r)$ is calculated in the following way:
\begin{equation}\label{eq:cosOmegaOfR}
\langle \cos \omega \rangle_{\text{AB}}(r)  = 
        \langle \sum_{i \in \text{A}} \sum_{j \in \text{B}}
        (\cos \omega_{ij}) \delta( r - |\mathbf{r}_{ij}|) \rangle
\end{equation}
Here $\cos \omega_{ij}$ is defined in the same way as in
Eq.~\ref{eq:gofrCosOmega}. This equation is a single dimensional pair
correlation that gives the average correlation of two directional
entities as a function of their distance from each other.

\subsection{\label{dynamicProps}Dynamic Property Analysis}

The dynamic properties of a trajectory are calculated with the program
\texttt{dynamicProps}. The program calculates the following properties:
\begin{gather}
\langle | \mathbf{r}(t) - \mathbf{r}(0) |^2 \rangle \label{eq:rms}\\
\langle \mathbf{v}(t) \cdot \mathbf{v}(0) \rangle \label{eq:velCorr} \\
\langle \mathbf{j}(t) \cdot \mathbf{j}(0) \rangle \label{eq:angularVelCorr}
\end{gather}

Eq.~\ref{eq:rms} is the root mean square displacement function. Which
allows one to observe the average displacement of an atom as a
function of time. The quantity is useful when calculating diffusion
coefficients because of the Einstein Relation, which is valid at long
times.\cite{allen87:csl}
\begin{equation}
2tD = \langle | \mathbf{r}(t) - \mathbf{r}(0) |^2 \rangle
\label{oopseEq:einstein}
\end{equation}

Eq.~\ref{eq:velCorr} and \ref{eq:angularVelCorr} are the translational
velocity and angular velocity correlation functions respectively. The
latter is only applicable to directional species in the
simulation. The velocity autocorrelation functions are useful when
determining vibrational information about the system of interest.

\section{\label{oopseSec:design}Program Design}

\subsection{\label{sec:architecture} {\sc oopse} Architecture}

The core of OOPSE is divided into two main object libraries:
\texttt{libBASS} and \texttt{libmdtools}. \texttt{libBASS} is the
library developed around the parsing engine and \texttt{libmdtools}
is the software library developed around the simulation engine. These
two libraries are designed to encompass all the basic functions and
tools that {\sc oopse} provides. Utility programs, such as the
property analyzers, need only link against the software libraries to
gain access to parsing, force evaluation, and input / output
routines.

Contained in \texttt{libBASS} are all the routines associated with
reading and parsing the \texttt{.bass} input files. Given a
\texttt{.bass} file, \texttt{libBASS} will open it and any associated
\texttt{.mdl} files; then create structures in memory that are
templates of all the molecules specified in the input files. In
addition, any simulation parameters set in the \texttt{.bass} file
will be placed in a structure for later query by the controlling
program.

Located in \texttt{libmdtools} are all other routines necessary to a
Molecular Dynamics simulation. The library uses the main data
structures returned by \texttt{libBASS} to initialize the various
parts of the simulation: the atom structures and positions, the force
field, the integrator, \emph{et cetera}. After initialization, the
library can be used to perform a variety of tasks: integrate a
Molecular Dynamics trajectory, query phase space information from a
specific frame of a completed trajectory, or even recalculate force or
energetic information about specific frames from a completed
trajectory.

With these core libraries in place, several programs have been
developed to utilize the routines provided by \texttt{libBASS} and
\texttt{libmdtools}. The main program of the package is \texttt{oopse}
and the corresponding parallel version \texttt{oopse\_MPI}. These two
programs will take the \texttt{.bass} file, and create (and integrate)
the simulation specified in the script. The two analysis programs
\texttt{staticProps} and \texttt{dynamicProps} utilize the core
libraries to initialize and read in trajectories from previously
completed simulations, in addition to the ability to use functionality
from \texttt{libmdtools} to recalculate forces and energies at key
frames in the trajectories. Lastly, the family of system building
programs (Sec.~\ref{oopseSec:initCoords}) also use the libraries to
store and output the system configurations they create.

\subsection{\label{oopseSec:parallelization} Parallelization of {\sc oopse}}

Although processor power is continually growing roughly following
Moore's Law, it is still unreasonable to simulate systems of more then
a 1000 atoms on a single processor. To facilitate study of larger
system sizes or smaller systems on long time scales in a reasonable
period of time, parallel methods were developed allowing multiple
CPU's to share the simulation workload. Three general categories of
parallel decomposition methods have been developed including atomic,
spatial and force decomposition methods.

Algorithmically simplest of the three methods is atomic decomposition
where N particles in a simulation are split among P processors for the
duration of the simulation. Computational cost scales as an optimal
$\mathcal{O}(N/P)$ for atomic decomposition. Unfortunately all
processors must communicate positions and forces with all other
processors at every force evaluation, leading communication costs to
scale as an unfavorable $\mathcal{O}(N)$, \emph{independent of the
number of processors}. This communication bottleneck led to the
development of spatial and force decomposition methods in which
communication among processors scales much more favorably. Spatial or
domain decomposition divides the physical spatial domain into 3D boxes
in which each processor is responsible for calculation of forces and
positions of particles located in its box. Particles are reassigned to
different processors as they move through simulation space. To
calculate forces on a given particle, a processor must know the
positions of particles within some cutoff radius located on nearby
processors instead of the positions of particles on all
processors. Both communication between processors and computation
scale as $\mathcal{O}(N/P)$ in the spatial method. However, spatial
decomposition adds algorithmic complexity to the simulation code and
is not very efficient for small N since the overall communication
scales as the surface to volume ratio $\mathcal{O}(N/P)^{2/3}$ in
three dimensions.

The parallelization method used in {\sc oopse} is the force
decomposition method.  Force decomposition assigns particles to
processors based on a block decomposition of the force
matrix. Processors are split into an optimally square grid forming row
and column processor groups. Forces are calculated on particles in a
given row by particles located in that processors column
assignment. Force decomposition is less complex to implement than the
spatial method but still scales computationally as $\mathcal{O}(N/P)$
and scales as $\mathcal{O}(N/\sqrt{P})$ in communication
cost. Plimpton has also found that force decompositions scale more
favorably than spatial decompositions for systems up to 10,000 atoms
and favorably compete with spatial methods up to 100,000
atoms.\cite{plimpton95}

\subsection{\label{oopseSec:memAlloc}Memory Issues in Trajectory Analysis}

For large simulations, the trajectory files can sometimes reach sizes
in excess of several gigabytes. In order to effectively analyze that
amount of data, two memory management schemes have been devised for
\texttt{staticProps} and for \texttt{dynamicProps}. The first scheme,
developed for \texttt{staticProps}, is the simplest. As each frame's
statistics are calculated independent of each other, memory is
allocated for each frame, then freed once correlation calculations are
complete for the snapshot. To prevent multiple passes through a
potentially large file, \texttt{staticProps} is capable of calculating
all requested correlations per frame with only a single pair loop in
each frame and a single read of the file.

The second, more advanced memory scheme, is used by
\texttt{dynamicProps}. Here, the program must have multiple frames in
memory to calculate time dependent correlations. In order to prevent a
situation where the program runs out of memory due to large
trajectories, the user is able to specify that the trajectory be read
in blocks. The number of frames in each block is specified by the
user, and upon reading a block of the trajectory,
\texttt{dynamicProps} will calculate all of the time correlation frame
pairs within the block. After in-block correlations are complete, a
second block of the trajectory is read, and the cross correlations are
calculated between the two blocks. This second block is then freed and
then incremented and the process repeated until the end of the
trajectory. Once the end is reached, the first block is freed then
incremented, and the again the internal time correlations are
calculated. The algorithm with the second block is then repeated with
the new origin block, until all frame pairs have been correlated in
time. This process is illustrated in
Fig.~\ref{oopseFig:dynamicPropsMemory}.

\begin{figure} 
\centering
\includegraphics[width=\linewidth]{dynamicPropsMem.eps}
\caption[A representation of the block correlations in \texttt{dynamicProps}]{This diagram illustrates the memory management used by \texttt{dynamicProps}, which follows the scheme: $\sum^{N_{\text{memory blocks}}}_{i=1}[ \operatorname{self}(i) + \sum^{N_{\text{memory blocks}}}_{j>i} \operatorname{cross}(i,j)]$. The shaded region represents the self correlation of the memory block, and the open blocks are read one at a time and the cross correlations between blocks are calculated.}
\label{oopseFig:dynamicPropsMemory}
\end{figure}

\section{\label{oopseSec:conclusion}Conclusion}

We have presented the design and implementation of our open source
simulation package {\sc oopse}. The package offers novel capabilities
to the field of Molecular Dynamics simulation packages in the form of
dipolar force fields, and symplectic integration of rigid body
dynamics. It is capable of scaling across multiple processors through
the use of force based decomposition using MPI. It also implements
several advanced integrators allowing the end user control over
temperature and pressure. In addition, it is capable of integrating
constrained dynamics through both the {\sc rattle} algorithm and the
z-constraint method.

These features are all brought together in a single open-source
program. This allows researchers to not only benefit from
{\sc oopse}, but also contribute to {\sc oopse}'s development as
well.

Revision:	1112
Committed:	Thu Apr 15 02:45:10 2004 UTC (21 years ago) by mmeineke
Content type:	application/x-tex
File size:	99402 byte(s)
Log Message:	some more updates.