trunk/oopsePaper/oopsePaper.tex

\documentclass[11pt]{article}
\usepackage{amsmath}
\usepackage{amssymb}
\usepackage{endfloat}
\usepackage{listings}
\usepackage{palatino}
\usepackage{graphicx}
\usepackage[ref]{overcite}
\usepackage{setspace}
\usepackage{tabularx}
\pagestyle{plain}
\pagenumbering{arabic}
\oddsidemargin 0.0cm \evensidemargin 0.0cm
\topmargin -21pt \headsep 10pt
\textheight 9.0in \textwidth 6.5in
\brokenpenalty=10000
\renewcommand{\baselinestretch}{1.2}
\renewcommand\citemid{\ } % no comma in optional reference note

\begin{document}
\lstset{language=C,frame=TB,basicstyle=\small,basicstyle=\ttfamily, %
        xleftmargin=0.5in, xrightmargin=0.5in,captionpos=b, %
        abovecaptionskip=0.5cm, belowcaptionskip=0.5cm}
\renewcommand{\lstlistingname}{Scheme}
\title{{\sc oopse}: An Open Source Object-Oriented Parallel Simulation
Engine for Molecular Dynamics}

\author{Matthew A. Meineke, Charles F. Vardeman II, Teng Lin,\\
 Christopher J. Fennell and J. Daniel Gezelter\\
Department of Chemistry and Biochemistry\\
University of Notre Dame\\
Notre Dame, Indiana 46556}

\date{\today}
\maketitle

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

\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 paper 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.  Program design
considerations are presented in Sec.~\ref{oopseSec:design}. And
lastly, Sec.~\ref{oopseSec:conclusion} concludes the chapter.

\section{\label{oopseSec:IOfiles}Concepts \& 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: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.

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:mdlExample} 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 had energy conservation problems when compared to  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}.

\begin{lstlisting}[float,caption={[Defining rigid bodies]A sample definition of a molecule containing 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:miscConcepts}Putting a Script Together}

The actual creation of a {\sc bass} script requires several key components. The first  part of the script needs to be the declaration of all of the molecule prototypes used in the simulation. This is typically done through the inclusion of {\tt .mdl} files. Only the molecules actually present in the simulation need to be declared, however {\sc bass} allows for the declaration of more molecules than are needed. This gives the user the ability to build up a library of commonly used molecules into a single {\tt .mdl} file.

Once all prototypes are declared, the ordering of the rest of the script is less stringent. Typically, the next to follow the molecular prototypes are the component statements. These statements specify which molecules are present within the simulation. The number of components must first be declared before the first component block statement (an example is seen in Sch.~\ref{sch:bassExPrime}).  The component blocks tell {\sc oopse} the number of molecules that will be in the simulation, and the order in which the components blocks are declared sets the ordering of the real atoms within the simulation as well as in the output files. 

The remainder of the script then sets the various simulation parameters for the system of interest. The required set of parameters that must be present in all simulations is given in Table~\ref{table:reqParams}.  The {\tt ensemble} statement is responsible for selecting the integration method used for the calculation of the equations of motion. An in depth discussion of the various methods available in {\sc oopse} can be found in Sec.~\ref{oopseSec:mechanics}. The {\tt forceField} statement is important for the selection of which forces will be used in the course of the simulation. {\sc oopse} supports several force fields, as outlined in Sec.~\ref{oopseSec:empericalEnergy}. The force fields are interchangeable between simulations, with the only requirement being that all atoms needed by the simulation are defined within the selected force field. The time step between force evaluations is set with the {\tt dt} parameter, and {\tt runTime} will set the time length of the simulation. Note, that {\tt runTime} is an absolute time, meaning if the simulation is started at t = 10.0~ns with a {\tt runTime} of 25.0~ns, the simulation will only run for an additional 15.0~ns. The final required parameter, is the {\tt initialConfig} statement. This will set the initial coordinates for the system, as well as the initial time if the {\tt useInitalTime = true;} flag is given. The format of the file specified in {\tt initialConfig}, is given in Sec.~\ref{oopseSec:coordFiles}. Additional parameters are summarized in Table~\ref{table:genParams}.

\begin{table}
\caption{The Global Keywords: Required Parameters}
\label{table:reqParams}
\begin{center}
% Note when adding or removing columns, the \hsize numbers must add up to the total number
% of columns.
\begin{tabularx}{\linewidth}%
  {>{\setlength{\hsize}{1.00\hsize}}X%
  >{\setlength{\hsize}{0.4\hsize}}X%
  >{\setlength{\hsize}{1.2\hsize}}X%
  >{\setlength{\hsize}{1.4\hsize}}X}

{\bf keyword} & {\bf units} & {\bf use} & {\bf remarks} \\ \hline

{\tt forceField} & string & Sets the force field. & Possible force fields are "DUFF", "LJ", and "EAM". \\
{\tt ensemble} & string & Sets the ensemble. & Possible ensembles are "NVE", "NVT", "NPTi", "NPTf", and "NPTxyz".\\ 
{\tt dt} & fs & Sets the time step. & Selection of {\tt dt} should be small enough to sample the fastest motion of the simulation. \\
{\tt nComponents} & integer & Sets the number of components. & Needs to appear before the first {\tt Component} block. \\
{\tt initialConfig} & string & Sets the file containing the initial configuration. & Can point to any file containing the configuration in the correct order. \\
{\tt runTime} & fs & Sets the time at which the simulation should end. & This is an absolute time, and will end the simulation when the current time meets or exceeds the {\tt runTime}. \\


\end{tabularx}
\end{center}
\end{table}

\begin{table}
\caption{The Global Keywords: General Parameters}
\label{table:genParams}
\begin{center}
% Note when adding or removing columns, the \hsize numbers must add up to the total number
% of columns.
\begin{tabularx}{\linewidth}%
  {>{\setlength{\hsize}{1.00\hsize}}X%
  >{\setlength{\hsize}{0.4\hsize}}X%
  >{\setlength{\hsize}{1.2\hsize}}X%
  >{\setlength{\hsize}{1.4\hsize}}X}

{\bf keyword} & {\bf units} & {\bf use} & {\bf remarks} \\ \hline

{\tt finalConfig} & string & Option to set the name of the final output file. & Useful when stringing simulations together. Defaults to the {\tt .bass} file with an {\tt .eor} extension. \\
{\tt useInitialTime} & logical & Sets whether the initial time is taken from the {\tt .init} file. & Useful when recovering a simulation from a crashed processor. Default is false. \\
{\tt sampleTime} & fs & Sets the frequency at which the {\tt .dump} file is written. & Default sets the frequency to the {\tt runTime}. \\
{\tt statusTime} & fs & Sets the frequency at which the {\tt .stat} file is written. & Defaults sets the frequency to the {\tt sampleTime}. \\
{\tt LJrcut} & $\mbox{\AA}$ & Manually sets the Lennard-Jones cutoff. & Defaults to $2.5\sigma_L$, where $\sigma_L$ is the largest LJ $\sigma$ in the simulation. \\
{\tt electrostaticCutoffRadius}& & & \\
      & $\mbox{\AA}$ & Manually sets the cutoff used by the electrostatic potentials. & Defaults to $15\mbox{\AA}$ \\
{\tt electrostaticSkinThickness} & & & \\
     & $\mbox{\AA}$  & Manually sets the skin thickness for the electrostatic switching function. & Defaults to 5~\% of the {\tt electrostaticSkinThickness}. \\
{\tt useReactionField} & logical & Turns the reaction field correction on/off. & Default is "false". \\
{\tt dielectric} & unitless & Sets the dielectric constant for reaction field. & If {\tt useReactionField} is true, then {\tt dielectric} must be set. \\
{\tt usePeriodicBoundaryConditions} & & & \\
        & logical & Turns periodic boundary conditions on/off. & Default is "true". \\
{\tt seed } & integer & Sets the seed value for the random number generator. & The seed needs to be at least 9 digits long. The default is to take the seed from the CPU clock.

\end{tabularx}
\end{center}
\end{table}


\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, body fixed angular momentum.},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 (\texttt{.eor}). 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.  Several helper programs are provided with {\sc oopse} to illustrate possible build routes. However, as each simulation is different, system creation is left to the end user. The {\tt .init} file must list the atoms in the correct order or {\sc oopse} will give an atom mismatch error. 

The correct ordering of the atoms relies on the ordering of atoms and molecules within the model and {\sc bass} scripts. {\sc oopse} expects the order to comply with the following guidelines:
\begin{enumerate}
\item All of the molecules of the first declared component are given before proceeding to the molecules of the second component, and so on for all declared components.
\item The ordering of the atoms for each molecule follows the order declared in the molecule's declaration within the model file.
\end{enumerate}
An example is given in Scheme~\ref{sch:initEx1} resulting in the {\tt .init} file shown in Scheme~\ref{sch:initEx2}.

\begin{lstlisting}[float,caption={This scheme illustrates the declaration of the $\text{I}_2$ molecule and the HCl molecule. The two molecules are then included into a simulation.}, label=sch:initEx1]

molecule{
  name = "I2";
  nAtoms = 2;
  atom[0]{
    type = "I";
  }
  atom[1]{
    type = "I";
  }
  nBonds = 1;
  bond[0]{
    members( 0, 1);
  }
}

molecule{
  name = "HCl"
  nAtoms = 2;
  atom[0]{
    type = "H";
  }
  atom[1]{
    type = "Cl";
  }
  nBonds = 1;
  bond[0]{
    members( 0, 1);
  }
}

nComponents = 2;
component{
  type = "HCl";
  nMol = 4;
}
component{
  type = "I2";
  nMol = 1;
}

initialConfig = "mixture.init";

\end{lstlisting}

\begin{lstlisting}[float,caption={This is the contents of the {\tt mixture.init} file matching the declarations in Scheme~\ref{sch:initEx1}. Note that even though $\text{I}_2$ is declared before HCl, the {\tt .init} file follows the order in which the components were included.},label=sch:initEx2]

10
0.0;  10.0  0.0  0.0;  0.0  10.0  0.0;  0.0  0.0  10.0;
H  ...
Cl ...
H  ...
Cl ...
H  ...
Cl ...
H  ...
Cl ...
I  ...
I  ...

\end{lstlisting}


\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:empiricalEnergy}The Empirical Energy Functions}

\
\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}$, removing the need for the computationally expensive Ewald sum. Instead, we Verlet neighbor-lists and cutoff radii are used for the dipolar interactions, or a reaction field is added to mimic longer 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\cite{fennell04} of the soft sticky dipole (SSD) model of Ichiye
\emph{et al.}\cite{liu96:new_model}

\begin{figure}
\centering
\includegraphics[width=\linewidth]{twoChainFig.pdf}
\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 and {\sc duff} also constrain 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}

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

\subsubsection{\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{fennell04} 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,fennell04} 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{fennell04}.

\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}

{\sc oopse} implements a potential that
describes bonding transition metal
systems\cite{Finnis84,Ercolessi88,Chen90,Qi99,Ercolessi02} and has 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: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.pdf}
\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{sec:constraints}Constraint Methods}

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

\subsubsection{\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: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}

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


\newpage
\section{Acknowledgments}
The authors would like to thank the Notre Dame BoB computer cluster where much of this project was tested. Additionally, the authors would like to acknowledge their funding from {\LARGE FIX ME}.

\bibliographystyle{achemso}
\bibliography{oopsePaper}

\end{document}
Revision:	1168
Committed:	Wed May 12 19:32:47 2004 UTC (20 years, 1 month ago) by mmeineke
Content type:	application/x-tex
File size:	89020 byte(s)
Log Message:	added more revisions to the paper