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 signifigant,
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.

{\sc give final breakdown of who wrote which section here.}

\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 convienent 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 analouges 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 do. Examples of uncommonly
implemented techniques and energetics include; dipole-dipole
interactions, rigid body dynamics, and metallic emmbedded
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.

Having written {\sc oopse} we are implementing the concept of Open
Source dcevelopment, and releaseing our source code into the public
domain. 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 emperical energy
functions that {\sc oopse } implements in
Sec.~\ref{oopseSec:empericalEnergy}. Following that is a discusion of
the various input and output files associated with the package
(Sec.~\ref{oopseSec:IOfiles}). In Sec.~\ref{oopseSec:Mechanics}
elucidates the various Molecular Dynamics algorithms {\sc oopse}
mplements 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 as well as the software distribution license is
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 on atoms are not currently supported by {\sc oopse}.

\begin{lstlisting}[float,caption={[Specifier for molecules and atoms] A sample specification of the simple 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 srtructures 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 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}
shws how one creates a molecule in the \texttt{.mdl} files. 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 requirement to
propagate the orientational degrees of freedom. Until recently,
integrators which propagate orientational motion have been lacking.

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 in that it is the torque applied on the center of mass
that dictates rotational motion. The torque on rigid body {\it i} is
\begin{equation}
\boldsymbol{\tau}_i=
        \sum_{a}(\mathbf{r}_{ia}-\mathbf{r}_i)\times \mathbf{f}_{ia} 
        + \boldsymbol{\tau}_{ia},
\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
the 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 internally. Further discussion
on this choice can be found in Sec.~\ref{sec:integrate}. An example
definition of a riged 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 Potential}

The most basic force field implemented in {\sc oopse} is the
Lennard-Jones potential, 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 partial \texttt{.bass} file that
shows a system of 108 Ar particles simulated with 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}]

/* 
 * The Ar molecule is specified 
 * external to the.bass file
 */

#include "argon.mdl" 

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

/*
 * The initial configuration is generated
 * before the simulation is invoked.
 */

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 is set to be
$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. To offset this discontinuity, 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.

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 ($\approx$100's of phospholipids,
$\approx$1000's of waters) to be simulated for long times
($\approx$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}$, allowing us
to avoid the computationally expensive Ewald sum. Instead, we can use
neighbor-lists, reaction field, and cutoff radii for the dipolar
interactions.

As an example, lipid head-groups in {\sc duff} are represented as
point dipole interaction sites. By placing a dipole of 20.6~Debye at
the head group center of mass, our model mimics the head group of
phosphatidylcholine.\cite{Cevc87} Additionally, a large Lennard-Jones
site is located at the pseudoatom's center of mass. The model is
illustrated by the dark grey atom in Fig.~\ref{fig: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 an 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 sampling of the bond potential to ensure conservation
of energy. 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}_{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{fig: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}
Here $\phi$ is the angle defined by four bonded neighbors $i$,
$j$, $k$, and $l$ (again, see Fig.~\ref{fig:lipidModel}). For
computational efficiency, the torsion potential has been recast after
the method of CHARMM,\cite{charmm1983} 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{sec: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}
        -
        \frac{3(\boldsymbol{\hat{u}}_i \cdot \mathbf{r}_{ij}) %
                (\boldsymbol{\hat{u}}_j \cdot \mathbf{r}_{ij}) }
                {r^{2}_{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}$.


\subsubsection{\label{sec: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 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 maintains reasonable agreement with the Soper
diffraction data for the structural features of liquid
water.\cite{Soper86,liu96:new_model} Additionally, the dynamical properties
exhibited by SSD 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 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
{\sc 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
\ref{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 reactionField flag toggles reaction 
 * field corrections.
 */

reactionField = false; // defaults to false
dielectric = 80.0; // dielectric for reaction field

/*
 * 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}

Several other molecular dynamics packages\cite{dynamo86} exist 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 a
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 done
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}S

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{FDB86}. 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 truly
macroscopic systems 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 boxes 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 \textquotedblleft feeling\textquotedblright\ the symmetries of the
periodic lattice, surface effects can be ignored. Cubic, orthorhombic and
parallelepiped are the available periodic cells In OOPSE. We use a matrix to
describe the property of the simulation box. Therefore, both the size and
shape of the simulation box can be changed during the simulation. The
transformation from box space vector $\mathbf{s}$ to its corresponding real
space vector $\mathbf{r}$ is defined by
\begin{equation}
\mathbf{r}=\underline{\mathbf{H}}\cdot\mathbf{s}%
\end{equation}


where $H=(h_{x},h_{y},h_{z})$ is a transformation matrix made up of the three
box axis vectors. $h_{x},h_{y}$ and $h_{z}$ represent the three sides of the
simulation box respectively.

To find the minimum image of a vector $\mathbf{r}$, we convert the real vector
to its corresponding vector in box space first, \bigskip%
\begin{equation}
\mathbf{s}=\underline{\mathbf{H}}^{-1}\cdot\mathbf{r}%
\end{equation}
And then, each element of $\mathbf{s}$ is wrapped to lie between -0.5 to 0.5,
\begin{equation}
s_{i}^{\prime}=s_{i}-\roundme(s_{i})
\end{equation}
where

%

\begin{equation}
\roundme(x)=\left\{
\begin{array}{cc}%
\lfloor{x+0.5}\rfloor & \text{if \ }x\geqslant 0 \\
\lceil{x-0.5}\rceil & \text{otherwise}%
\end{array}
\right.
\end{equation}


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}=\underline{\mathbf{H}}^{-1}\cdot\mathbf{s}^{\prime}%
\end{equation}


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

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

Every {\sc oopse} simuation begins with a {\sc bass} file. {\sc bass}
(\underline{B}izarre \underline{A}tom \underline{S}imulation
\underline{S}yntax) 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 are 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
Fig.~\ref{fig:bassExample}.

\begin{figure}
\centering
\framebox[\linewidth]{\rule{0cm}{0.75\linewidth}I'm a {\sc bass} file!}
\caption{Here is an example \texttt{.bass} file}
\label{fig:bassExample}
\end{figure}

Within the \texttt{.bass} file it is neccassary 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
inconvient to duplicate the simulation at the begining 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.

\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
App.~\ref{appCoordFormat}. 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.

There are three major files used by {\sc oopse} written in the coordinate
format, they are as follows: the initialization file, the simulation
trajectory file, and the final coordinates of the simulation. The
initialization file is neccassary for {\sc oopse} to start the simulation
with the proper coordinates. It is typically denoted with the
extension \texttt{.init}. 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. The
trajectory file is given the extension \texttt{.dump}. The final
coordinate file is the end of run or \texttt{.eor} file. The
\texttt{.eor} file stores the final configuration of teh 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 or
recover the previous simulation.

\subsection{Generation of Initial Coordinates}

As was stated in Sec.~\ref{subSec: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} is {\sc bass}
aware, and will recognize arguments and parameters in the
\texttt{.bass} file that would otherwise be ignored by the
simulation. The program itself is under contiunual development, and is
offered here as a helper tool only.

\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. The file allows the
user to observe the system variables as a function od 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{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 this integrator selection
deals with poor energy conservation of rigid body systems using
quaternions. 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 this code 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.} is that the entire rotation matrix is propagated from
one time step to the next. In the past, this would not have been as
feasible a option, being that the rotation matrix for a single body is
nine elements long as opposed to 3 or 4 elements for Euler angles and
quaternions respectively. System memory has become much less of an
issue in recent times, and this has resulted in substantial benefits
in energy conservation. There is still the issue of 5 or 6 additional
elements for describing the orientation of each particle, which will
increase dump files substantially. Simply translating the rotation
matrix into its component Euler angles or quaternions for storage
purposes relieves this burden.

The symplectic splitting 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 end up being more costly computationally than the
simpler arithmetic quaternion propagation. With the same time step, a
1000 SSD particle simulation shows an average 7\% increase in
computation time using the symplectic step method in place of
quaternions. This cost is more than justified when comparing the
energy conservation of the two methods as illustrated in figure
\ref{timestep}.

\begin{figure}
\centering
\includegraphics[width=\linewidth]{timeStep.eps}
\caption{Energy conservation using quaternion based integration versus 
the symplectic step method proposed by Dullweber \emph{et al.} with
increasing time step. For each time step, the dotted line is total
energy using the symplectic step integrator, and the solid line comes
from the quaternion integrator. The larger time step plots are shifted
up from the true energy baseline for clarity.}
\label{timestep}
\end{figure}

In figure \ref{timestep}, the resulting energy drift at various time
steps for both the symplectic step and quaternion integration schemes
is compared. All of the 1000 SSD particle simulations started with the
same configuration, and the only difference was the method for
handling rotational motion. At time steps of 0.1 and 0.5 fs, both
methods for propagating particle rotation conserve energy fairly well,
with the quaternion method showing a slight energy drift over time in
the 0.5 fs time step simulation. At time steps of 1 and 2 fs, the
energy conservation benefits of the symplectic step method are clearly
demonstrated. Thus, while maintaining the same degree of energy
conservation, one can take considerably longer time steps, leading to
an overall reduction in computation time.

Energy drift in these SSD particle simulations was unnoticeable for
time steps up to three femtoseconds. A slight energy drift on the
order of 0.012 kcal/mol per nanosecond was observed at a time step of
four femtoseconds, and as expected, this drift increases dramatically
with increasing time step. To insure accuracy in the constant energy
simulations, time steps were set at 2 fs and kept at this value for
constant pressure simulations as well.


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


{\sc oopse} implements a 


\subsubsection{\label{sec: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{Sec:zcons}Z-Constraint Method}

Based on fluctuatin-dissipation theorem,\bigskip\ force auto-correlation
method was developed to investigate the dynamics of ions inside the ion
channels.\cite{Roux91} Time-dependent friction coefficient can be calculated
from the deviation of the instaneous force from its mean force.

%

\begin{equation}
\xi(z,t)=\langle\delta F(z,t)\delta F(z,0)\rangle/k_{B}T
\end{equation}
where%
\begin{equation}
\delta F(z,t)=F(z,t)-\langle F(z,t)\rangle
\end{equation}


If the time-dependent friction decay rapidly, static friction coefficient can
be approximated by%

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


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


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

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

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

\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{fig: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.  

\begin{figure} 
\centering
\includegraphics[width=\linewidth]{dynamicPropsMem.eps}
\caption{This diagram illustrates the dynamic memory allocation 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{fig:dynamicPropsMemory}
\end{figure}

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

\subsection{\label{sec:architecture} 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 parseing engine and {\texttt libmdtools} is the
software library developed around the simulation engine.

\subsection{\label{sec:parallelization} Parallelization of OOPSE}

Although processor power is doubling roughly every 18 months according
to the famous Moore's Law\cite{moore}, 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{openSource}Open Source and Distribution License}


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

\begin{itemize}
        
\item Restate capabilities

\item recap major structure / design choices

        \begin{itemize}
        
        \item parallel
        \item symplectic integration
        \item languages

        \end{itemize}

\item How well does it meet the primary goal

\end{itemize}
Revision:	1051
Committed:	Fri Feb 13 21:24:07 2004 UTC (20 years, 5 months ago) by mmeineke
Content type:	application/x-tex
File size:	51028 byte(s)
Log Message:	added an introduction to oopse, and started smoothing out the wrinkles.