trunk/oopsePaper/EmpericalEnergy.tex


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

\subsection{\label{sec: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}[caption={[Specifier for molecules and atoms] An example specifying 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}

The second most basic building block of a simulation is the
molecule. The molecule is a way for {\sc oopse} to keep track of the
atoms in a simulation in logical manner. This particular unit will
store the identities of all the atoms associated with itself and is
responsible for the evaluation of its own bonded interaction
(i.e.~bonds, bends, and torsions).

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

\subsection{\label{sec:LJPot}The Lennard Jones Potential}

The most basic force field implemented in {\sc oopse} is the Lennard-Jones
potential. 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 particle $i$ and $j$,
$\sigma_{ij}$ scales the length of the interaction, and
$\epsilon_{ij}$ scales the well depth of the potential.

Because this potential is calculated between all pairs, the force
evaluation can become computationally expensive for large systems. To
keep the pair evaluation 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 in the system. 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 entire potential. This causes
the potential to go to zero 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{sec:DUFF}Dipolar Unified-Atom Force Field}

The Dipolar Unified-atom Force Field ({\sc duff}) was developed to
simulate lipid bilayers. The systems require a model capable of forming
bilayers, while still being sufficiently computationally efficient to
allow simulations of large systems ($\approx$100's of phospholipids,
$\approx$1000's of waters) 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 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 out
reparameterization of the soft sticky dipole (SSD) model of Ichiye
\emph{et al.}\cite{liu96:new_model}

\begin{figure}
\epsfxsize=\linewidth
\epsfbox{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{fig:lipidModel}
\end{figure}

Turning to the tails of the phospholipids, we have used a set of
scalable parameters to model the alkyl groups with Lennard-Jones
sites. For this, we have used 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 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 of 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 conservation of
energy. By constraining the bond lengths, larger time steps may be
used when integrating the equations of motion.


\subsubsection{\label{subSec:energyFunctions}{\sc duff} Energy Functions}

The total energy of function in {\sc duff} is given by the following:
\begin{equation}
V_{\text{Total}} = \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 a molecule:
\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, 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_{\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}), and $\theta_0$ is the equilibrium
bond angle. $k_{\theta}$ is the force constant which determines the
strength of the harmonic bend. The parameters for $k_{\theta}$ and
$\theta_0$ are based off of those in TraPPE.\cite{Siepmann1998}

The torsion potential and parameters are also taken 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_{ijkl}$ 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} whereby the angle series is
converted to a power series of the form:
\begin{equation}
V_{\text{torsion}}(\phi_{ijkl}) =  
        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 equation to a power series, repeated trigonometric
evaluations are avoided during the calculation of the potential.


The cross portion of the total potential, $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. 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 Euler angles of atom $i$ and $j$ respectively. $|\mu_i|$ is
the magnitude of the dipole moment of atom $i$ and
$\boldsymbol{\hat{u}}_i$ is the standard unit orientation vector of
$\boldsymbol{\Omega}_i$.


\subsection{\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}.

!!!Place a {\sc BASS} scheme showing SSD parameters around here!!!

\subsection{\label{sec: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{Sec: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 "feeling" the symmetries of
the periodic lattice, surface effects can be ignored. Cubic,
orthorhombic and parallelepiped are the available periodic cells In
{\sc 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{\underline{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, we convert the real vector to its
corresponding vector in box space first, \bigskip%
\begin{equation}
\mathbf{s}=\underline{\underline{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\geqslant0\\
\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 by transforming back
to real space,%

\begin{equation}
\mathbf{r}^{\prime}=\underline{\underline{H}}^{-1}\cdot\mathbf{s}^{\prime}%
\end{equation}

Revision:	967
Committed:	Tue Jan 20 14:17:21 2004 UTC (21 years, 3 months ago) by chrisfen
Content type:	application/x-tex
File size:	21394 byte(s)
Log Message:	Corrected .bib file entries