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 package {\sc oopse}. It is important to note, that a
simulation package 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 package'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 uncommonly
implemented techniques and energetics include; dipole-dipole
interactions, rigid body dynamics, and metallic embedded
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 package 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
(\textbf{A}) can be described by the three Euler angles ($\phi,
\theta,$ and $\psi$), where the elements of \textbf{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 \textbf{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_water";
  nRigidBodies = 1;
  rigidBody[0]{ 
    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 );
    }
    position( 0.0, 0.0, 0.0 );
    orientation( 0.0, 0.0, 1.0 );
  }
}
\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 and 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.

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
center of mass, our model mimics the charge separation found in common
phospholipids 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]{lipidModel.eps}
\caption{A representation of the lipid model. $\phi$ is the torsion angle, $\theta$ %
is the bend angle, $\mu$ is the dipole moment of the head group, and n
is the chain length.}
\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}]Sample \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 pairs that are closer than three bonds,
i.e.~atom pairs farther away than a torsion are included in the
pair-wise loop.


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. 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]An example 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 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 et al. fit EAM potentials for fcc metals Cu, Ag, Au, Ni, Pd, Pt and alloys of these metals\cite{FBD86}. These potential fits are in the DYNAMO 86 format and are included with {\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, $\mathbf{H}$, to describe the shape and
size of the simulation box. $\mathbf{H}$ is defined:
\begin{equation}
\mathbf{H} = ( \mathbf{h}_x, \mathbf{h}_y, \mathbf{h}_z )
\end{equation}
Where $\mathbf{h}_j$ is the column vector of the $j$th axis of the
box.  During the course of the simulation both the size and shape of
the box can be changed to allow volume fluctations 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} &= \mathbf{H}^{-1} \mathbf{r} \\
\mathbf{r} &= \mathbf{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 on the 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}=\mathbf{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 interatomic 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 enesemble
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"

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";
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 $\mathbf{H}$ column vectors. 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 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. 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 a program called
\texttt{sysBuilder} to aid in the creation of the \texttt{.init}
file. \texttt{sysBuilder} uses {\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 Symplectic Step Integrator}

Integration of the equations of motion was carried out using the
symplectic splitting method proposed by Dullweber \emph{et
al.}.\cite{Dullweber1997} The reason for the selection of this
integrator, is the poor energy conservation of rigid body systems
using quaternion dynamics. While quaternions work well for
orientational motion in alternate ensembles, the microcanonical
ensemble has a constant energy requirement that is quite sensitive to
errors in the equations of motion. The original implementation of {\sc
oopse} utilized quaternions for rotational motion propagation;
however, a detailed investigation showed that they resulted in a
steady drift in the total energy, something that has been observed by
others.\cite{Laird97}

The key difference in the integration method proposed by Dullweber
\emph{et al}.~({\sc dlm}) is that the entire rotation matrix is propagated from
one time step to the next. In the past, this would not have been a
feasible option, since the rotation matrix for a single body is nine
elements long as opposed to three or four elements for Euler angles
and quaternions respectively. System memory has become much less of an
issue in recent times, and the {\sc dlm} method has used memory in
exchange for substantial benefits in energy conservation.

The {\sc dlm} method allows for Verlet style integration of both
linear and angular motion of rigid bodies. In the integration method,
the orientational propagation involves a sequence of matrix
evaluations to update the rotation matrix.\cite{Dullweber1997} These
matrix rotations are more costly computationally than the simpler
arithmetic quaternion propagation. With the same time step, a 1000 SSD
particle simulation shows an average 7\% increase in computation time
using the {\sc 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 {\sc dlm} dynamics]{Energy conservation using quaternion based integration versus 
the {\sc dlm} method with
increasing time step. For each time step, the dotted line is total
energy using the {\sc 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 {\sc dlm} and quaternion integration schemes
is compared. All of the 1000 SSD particle simulations started with the
same configuration, and the only difference was the method for
handling rotational motion. At time steps of 0.1 and 0.5 fs, both
methods for propagating particle rotation conserve energy fairly well,
with the quaternion method showing a slight energy drift over time in
the 0.5 fs time step simulation. At time steps of 1 and 2 fs, the
energy conservation benefits of the {\sc 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.

Energy drift in these SSD particle simulations was unnoticeable for
time steps up to three femtoseconds. A slight energy drift on the
order of 0.012 kcal/mol per nanosecond was observed at a time step of
four femtoseconds, and as expected, this drift increases dramatically
with increasing time step.


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


{\sc oopse} implements a 


\subsection{\label{oopseSec:noseHooverThermo}Nose-Hoover Thermostatting}

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

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

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

\subsection{\label{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 on the equations of
motion.\cite{fowles99:lagrange} 

Consider a system described by qoordinates $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 neccassary 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 statisitics 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 enesemble
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 perpindicular 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,
\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}


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

Based on fluctuation-dissipation theorem, a force auto-correlation
method was developed to investigate the dynamics of ions inside the ion
channels.\cite{Roux91} Time-dependent friction coefficient can be calculated
from the deviation of the 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 decay rapidly, static friction coefficient can
be approximated by%

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


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


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

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

\begin{equation}
U(t)=\frac{1}{2}k_{Harmonic}(z(t)-z_{cons})^{2}%
\end{equation}
where $k_{Harmonic}$\bigskip\ is the harmonic force constant, $z(t)$ is
current z coordinate of the center of mass of the z-constraint molecule, and
$z_{cons}$ is the restraint position. Therefore, the harmonic force operated
on the z-constraint molecule at time $t$ can be calculated by%
\begin{equation}
F_{z_{Harmonic}}(t)=-\frac{\partial U(t)}{\partial z(t)}=-k_{Harmonic}%
(z(t)-z_{cons})
\end{equation}
Worthy of mention, other kinds of potential functions can also be used to
drive the z-constraint molecule.

\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 the following
pair correlations between species A and B:
\begin{itemize}
        \item $g_{\text{AB}}(r)$: Eq.~\ref{eq:gofr}
        \item $g_{\text{AB}}(r, \cos \theta)$: Eq.~\ref{eq:gofrCosTheta}
        \item $g_{\text{AB}}(r, \cos \omega)$: Eq.~\ref{eq:gofrCosOmega}
        \item $g_{\text{AB}}(x, y, z)$: Eq.~\ref{eq:gofrXYZ}
        \item $\langle \cos \omega \rangle_{\text{AB}}(r)$: 
                Eq.~\ref{eq:cosOmegaOfR}
\end{itemize}

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.

All static properties are calculated on a frame by frame basis. The
trajectory is read a single frame at a time, and the appropriate
calculations are done on each frame. Once one frame is finished, the
next frame is read in, and a running average of the property being
calculated is accumulated in each frame. The program allows for the
user to specify more than one property be calculated in single run,
preventing the need to read a file multiple times.

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

The dynamic properties of a trajectory are calculated with the program
\texttt{dynamicProps}. The program will calculate 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. Eq.~\ref{eq:velCorr} and Eq.~\ref{eq:angularVelCorr} are the
velocity and angular velocity correlation functions respectively. The
latter is only applicable to directional species in the simulation.

The \texttt{dynamicProps} program handles he file in a manner different from
\texttt{staticProps}. As the properties calculated by this program are time
dependent, multiple frames must be read in simultaneously by the
program. For small trajectories this is no problem, and the entire
trajectory is read into memory. However, for long trajectories of
large systems, the files can be quite large. In order to accommodate
large files, \texttt{dynamicProps} adopts a scheme whereby two blocks of memory
are allocated to read in several frames each.

In this two block scheme, the correlation functions are first
calculated within each memory block, then the cross correlations
between the frames contained within the two blocks are
calculated. Once completed, the memory blocks are incremented, and the
process is repeated. A diagram illustrating the process is shown in
Fig.~\ref{oopseFig:dynamicPropsMemory}. As was the case with
\texttt{staticProps}, multiple properties may be calculated in a
single run to avoid multiple reads on the same file.


\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 then 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 month by month, 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 method's have been developed including atomic, spatial
and force decomposition methods.

Algorithmically simplest of the three method's 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
$O(N/P)$ for atomic decomposition. Unfortunately all processors must
communicate positions and forces with all other processors leading
communication to scale as an unfavorable $O(N)$ 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 $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 $(N/P)^{2/3}$ in three
dimensions.

Force decomposition assigns particles to processors based on a block
decomposition of the force matrix. Processors are split into a
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 then the spatial method but still scales computationally
as $O(N/P)$ and scales as $(N/\sqrt{p})$ in communication
cost. Plimpton also found that force decompositions scales more
favorably then spatial decomposition up to 10,000 atoms and favorably
competes with spatial methods for up to 100,000 atoms.

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

\subsection{\label{openSource}Open Source and Distribution License}

\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 MPI. It also implements several integration
ensembles 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
development package. This allows researchers to not only benefit from
{\sc oopse}, but also contribute to {\sc oopse}'s development as
well.Documentation and source code for {\sc oopse} can be downloaded
from \texttt{http://www.openscience.org/oopse/}.

Revision:	1071
Committed:	Thu Feb 26 21:55:43 2004 UTC (21 years, 8 months ago) by mmeineke
Content type:	application/x-tex
File size:	66434 byte(s)
Log Message:	added citations to the bib, and started on a honking huge section of RATTLE in the oopse chapter.