trunk/oopsePaper/EmpericalEnergy.tex


\section{\label{sec:empericalEnergy}The Emperical 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 a methyl
group. The atoms are also capable of having a directional component
associated with them, often in the form of a dipole. Charges on atoms
are not currently suporrted by {\sc oopse}.

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 in the previously, one of the features that sets 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 many standard simulation packages because of the need to
consider orientational degrees of freedom and include an integrator
that accurately propagates these motions in time.

Moving a rigid body involves determination of both the force and
torque applied by the surroundings, which directly affect the
translation and rotation in turn. In order to accumulate the total
force on a rigid body, the external forces must first be calculated
for all the interal 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 similar to the force in that it is a sum
of the torque applied on each internal particle, mapped onto the
center of mass of the rigid body. 

The application 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 be
maintained for each rigid body. At a minimum, the rotation matrix can
be described and propagated by the three Euler
angles.\cite{Goldstein01} In order to avoid rotational limitations
when using Euler angles, the four parameter ``quaternion'' scheme can
be used instead.\cite{allen87:csl} Use of quaternions also leads to
performance enhancements, particularly for very small
systems.\cite{Evans77} 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{sec:integrate}.

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

The most basic force field implemented in OOPSE is the Lennard-Jones
potential. The Lennard-Jones potential mimics the attractive forces
two charge neutral particles experience when spontaneous dipoles are
induced on each other. This is the predominant intermolecular force in
systems of such as noble gases and simple alkanes.

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
energy well depth of the potential.

Because the potential is calculated between all pairs, the force
evaluation can become computationally expensive for large systems. To
keep the pair evaluation to a manegable number, OOPSE employs the use
of a cut-off radius.\cite{allen87:csl} The cutoff radius is set to be
$2.5\sigma_{ii}$, where $\sigma_{ii}$ is the largest self self 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 equation 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 \underline{D}ipolar \underline{U}nified-Atom
\underline{F}orce \underline{F}ield ({\sc duff}) was developed to
simulate lipid bilayers. We needed 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, we have eliminated all point charges. Charge
distributions were replaced with dipoles, and charge-neutral
distributions 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 and cutoff radii for the dipolar interactions.

As an example, lipid head groups in {\sc duff} are represented as point
dipole interaction sites.PC and PE Lipid head groups are typically
zwitterionic in nature, with charges separated by as much as
6~$\mbox{\AA}$. By placing a dipole of 20.6~Debye at the head group
center of mass, our model mimics the head group of PC.\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}.

\begin{figure}
\epsfxsize=6in
\epsfbox{lipidModel.epsi}
\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}

The water model we use to complement the dipoles of the lipids is
the soft sticky dipole (SSD) model of Ichiye \emph{et
al.}\cite{liu96:new_model} This model is discussed in greater detail
in Sec.~\ref{sec:SSD}. In all cases we reduce water to a single
Lennard-Jones interaction site. The site also contains a dipole to
mimic the partial charges on the hydrogens and the oxygen. However,
what makes the SSD model unique is the inclusion of a tetrahedral
short range potential to recover the hydrogen bonding of water, an
important factor when modeling bilayers, as it has been shown that
hydrogen bond network formation is a leading contribution to the
entropic driving force towards lipid bilayer formation.\cite{Cevc87}


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^{N}_{J=1} 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}}(\theta_{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 in the molecule. $V_{\text{torsion}}$ is the torsion The
pairwise portions of the internal potential are excluded for pairs
that ar closer than three bonds, i.e.~atom pairs farther away than a
torsion are included in the pair-wise loop. 

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.

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_{ijkl}) = 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}).  However,
for computaional efficency, the torsion potentail 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 dipole-dipole potential has the following form:
\begin{equation}
V_{\text{dipole}}(\mathbf{r}_{ij},\boldsymbol{\Omega}_{i},
        \boldsymbol{\Omega}_{j}) = \frac{1}{4\pi\epsilon_{0}} \biggl[
        \frac{\boldsymbol{\mu}_{i} \cdot \boldsymbol{\mu}_{j}}{r^{3}_{ij}}
        -
        \frac{3(\boldsymbol{\mu}_i \cdot \mathbf{r}_{ij}) %
                (\boldsymbol{\mu}_j \cdot \mathbf{r}_{ij}) }
                {r^{5}_{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. $\boldsymbol{\mu}_i$ is the dipole vector of atom
$i$ it takes its orientation from $\boldsymbol{\Omega}_i$.


\subsection{\label{sec:SSD}Water Model: SSD and Derivatives}

In the interest of computational efficiency, the native solvent used
in {\sc oopse} is the Soft Sticky Dipole (SSD) water model. SSD was
developed by Ichiye \emph{et al.} 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}
u_{ij} = u_{ij}^{LJ} (r_{ij})\ + u_{ij}^{dp}
(\mathbf{r}_{ij},\boldsymbol{\Omega}_i,\boldsymbol{\Omega}_j)\ +
u_{ij}^{sp}
(\mathbf{r}_{ij},\boldsymbol{\Omega}_i,\boldsymbol{\Omega}_j),
\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, dipole,
and sticky parts of the potential are giving by the following
equations,
\begin{equation}
u_{ij}^{LJ}(r_{ij}) = 4\epsilon \left[\left(\frac{\sigma}{r_{ij}}\right)^{12}-\left(\frac{\sigma}{r_{ij}}\right)^{6}\right],
\end{equation}
\begin{equation}
u_{ij}^{dp} = \frac{\boldsymbol{\mu}_i\cdot\boldsymbol{\mu}_j}{r_{ij}^3}-\frac{3(\boldsymbol{\mu}_i\cdot\mathbf{r}_{ij})(\boldsymbol{\mu}_j\cdot\mathbf{r}_{ij})}{r_{ij}^5}\ ,
\end{equation}
\begin{equation}
\begin{split}
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)\\
& \quad \ +
s^\prime(r_{ij})w^\prime(\mathbf{r}_{ij},\boldsymbol{\Omega}_i,\boldsymbol{\Omega}_j)]\ ,
\end{split}
\end{equation}
where $\boldsymbol{\mu}_i$ and $\boldsymbol{\mu}_j$ are the dipole
unit vectors of particles \emph{i} and \emph{j} with magnitude 2.35 D,
$\nu_0$ scales the strength of the overall sticky potential, $s$ and
$s^\prime$ are cubic switching functions. The $w$ and $w^\prime$
functions take the following forms,
\begin{equation}
w(\mathbf{r}_{ij},\boldsymbol{\Omega}_i,\boldsymbol{\Omega}_j)=\sin\theta_{ij}\sin2\theta_{ij}\cos2\phi_{ij},
\end{equation}
\begin{equation}
w^\prime(\mathbf{r}_{ij},\boldsymbol{\Omega}_i,\boldsymbol{\Omega}_j) = (\cos\theta_{ij}-0.6)^2(\cos\theta_{ij}+0.8)^2-w^0,
\end{equation}
where $w^0 = 0.07715$. The $w$ function is the tetrahedral attractive
term that promotes hydrogen bonding orientations within the first
solvation shell, and $w^\prime$ is a dipolar repulsion term that
repels unrealistic dipolar arrangements within the first solvation
shell. A more detailed description of the functional parts and
variables in this potential can be found in other
articles.\cite{liu96:new_model,chandra99:ssd_md}

Since SSD is a one-site point dipole model, the force calculations are
simplified significantly from a computational standpoint, in that the
number of long-range interactions is dramatically reduced. In the
original Monte Carlo simulations using this model, Ichiye \emph{et
al.} reported a calculation speed up of up to an order of magnitude
over other comparable models while maintaining the structural behavior
of water.\cite{liu96:new_model} In the original molecular dynamics studies of
SSD, it was shown that it actually improves upon the prediction of
water's dynamical properties over TIP3P and SPC/E.\cite{chandra99:ssd_md}

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} Reparameterizations of the
original SSD have resulted in improved density behavior, as well as
alterations in the water structure and transport behavior. {\sc oopse} is
easily modified to impliment these new potential parameter sets for
the derivative water models: SSD1, SSD/E, and SSD/RF. All of the
variable parameters are listed in the accompanying BASS file, and
these parameters simply need to be changed to the updated values.


\subsection{\label{sec:eam}Embedded Atom Model}

here there be Monsters

\subsection{\label{Sec:pbc}Periodic Boundary Conditions} 
 
\textit{Periodic boundary conditions} are widely used to simulate truly 
macroscopic systems with a relatively small number of particles. Simulation 
box is replicated throughout space to form an infinite lattice. During the 
simulation, when a particle moves in the primary cell, its periodic image 
particles 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 symmetric of the periodic lattice,the surface 
effect could be ignored. Cubic and parallelepiped are the available periodic 
cells. \bigskip In OOPSE, we use the matrix instead of the vector to describe 
the property of the simulation box. Therefore, not only the size of the 
simulation box could be changed during the simulation, but also the shape of 
it. The transformation from box space vector $\overrightarrow{s}$ to its 
corresponding real space vector $\overrightarrow{r}$ is defined by 
\begin{equation} 
\overrightarrow{r}=H\overrightarrow{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 sides of the 
simulation box respectively. 
 
To find the minimum image, we need to convert the real vector to its 
corresponding vector in box space first, \bigskip% 
\begin{equation} 
\overrightarrow{s}=H^{-1}\overrightarrow{r}% 
\end{equation} 
And then, each element of $\overrightarrow{s}$ is casted to lie between -0.5 
to 0.5, 
\begin{equation} 
s_{i}^{\prime}=s_{i}-round(s_{i}) 
\end{equation} 
where% 
 
\begin{equation} 
round(x)=\lfloor{x+0.5}\rfloor\text{ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ }if\text{ 
}x\geqslant0 
\end{equation} 
% 
 
\begin{equation} 
round(x)=\lceil{x-0.5}\rceil\text{ \ \ \ \ \ \ \ \ \ \ \ \ \ \ }if\text{ }x<0 
\end{equation} 
 
 
For example, $round(3.6)=4$,$round(3.1)=3$, $round(-3.6)=-4$, $round(-3.1)=-3$. 
 
Finally, we could get the minimum image by transforming back to real space,% 
 
\begin{equation} 
\overrightarrow{r^{\prime}}=H^{-1}\overrightarrow{s^{\prime}}% 
\end{equation} 
Revision:	915
Committed:	Fri Jan 9 20:25:50 2004 UTC (21 years, 4 months ago) by mmeineke
Content type:	application/x-tex
File size:	18489 byte(s)
Log Message:	added a lennard jones section to the Emperical Energy section
#	User	Rev	Content
1	mmeineke	806
2	mmeineke	899	\section{\label{sec:empericalEnergy}The Emperical Energy Functions}
3	mmeineke	806
4	mmeineke	915	\subsection{\label{sec:atomsMolecules}Atoms, Molecules and Rigid Bodies}
5	mmeineke	806
6	gezelter	818	The basic unit of an {\sc oopse} simulation is the atom. The parameters
7	mmeineke	806	describing the atom are generalized to make the atom as flexible a
8			representation as possible. They may represent specific atoms of an
9			element, or be used for collections of atoms such as a methyl
10			group. The atoms are also capable of having a directional component
11			associated with them, often in the form of a dipole. Charges on atoms
12	gezelter	818	are not currently suporrted by {\sc oopse}.
13	mmeineke	806
14			The second most basic building block of a simulation is the
15	gezelter	818	molecule. The molecule is a way for {\sc oopse} to keep track of the atoms
16	mmeineke	806	in a simulation in logical manner. This particular unit will store the
17			identities of all the atoms associated with itself and is responsible
18			for the evaluation of its own bonded interaction (i.e.~bonds, bends,
19			and torsions).
20	mmeineke	899
21	mmeineke	915	As stated in the previously, one of the features that sets OOPSE apart
22			from most of the current molecular simulation packages is the ability
23			to handle rigid body dynamics. Rigid bodies are non-spherical
24			particles or collections of particles that have a constant internal
25			potential and move collectively.\cite{Goldstein01} They are not
26			included in many standard simulation packages because of the need to
27			consider orientational degrees of freedom and include an integrator
28			that accurately propagates these motions in time.
29
30			Moving a rigid body involves determination of both the force and
31			torque applied by the surroundings, which directly affect the
32			translation and rotation in turn. In order to accumulate the total
33			force on a rigid body, the external forces must first be calculated
34			for all the interal particles. The total force on the rigid body is
35			simply the sum of these external forces. Accumulation of the total
36			torque on the rigid body is similar to the force in that it is a sum
37			of the torque applied on each internal particle, mapped onto the
38			center of mass of the rigid body.
39
40			The application of the total torque is done in the body fixed axis of
41			the rigid body. In order to move between the space fixed and body
42			fixed coordinate axes, parameters describing the orientation be
43			maintained for each rigid body. At a minimum, the rotation matrix can
44			be described and propagated by the three Euler
45			angles.\cite{Goldstein01} In order to avoid rotational limitations
46			when using Euler angles, the four parameter ``quaternion'' scheme can
47			be used instead.\cite{allen87:csl} Use of quaternions also leads to
48			performance enhancements, particularly for very small
49			systems.\cite{Evans77} OOPSE utilizes a relatively new scheme that
50			propagates the entire nine parameter rotation matrix. Further
51			discussion on this choice can be found in Sec.~\ref{sec:integrate}.
52
53			\subsection{\label{sec:LJPot}The Lennard Jones Potential}
54
55			The most basic force field implemented in OOPSE is the Lennard-Jones
56			potential. The Lennard-Jones potential mimics the attractive forces
57			two charge neutral particles experience when spontaneous dipoles are
58			induced on each other. This is the predominant intermolecular force in
59			systems of such as noble gases and simple alkanes.
60
61			The Lennard-Jones potential is given by:
62			\begin{equation}
63			V_{\text{LJ}}(r_{ij}) =
64			4\epsilon_{ij} \biggl[
65			\biggl(\frac{\sigma_{ij}}{r_{ij}}\biggr)^{12}
66			- \biggl(\frac{\sigma_{ij}}{r_{ij}}\biggr)^{6}
67			\biggr]
68			\label{eq:lennardJonesPot}
69			\end{equation}
70			Where $r_ij$ is the distance between particle $i$ and $j$, $\sigma_{ij}$
71			scales the length of the interaction, and $\epsilon_{ij}$ scales the
72			energy well depth of the potential.
73
74			Because the potential is calculated between all pairs, the force
75			evaluation can become computationally expensive for large systems. To
76			keep the pair evaluation to a manegable number, OOPSE employs the use
77			of a cut-off radius.\cite{allen87:csl} The cutoff radius is set to be
78			$2.5\sigma_{ii}$, where $\sigma_{ii}$ is the largest self self length
79			parameter in the system. Truncating the calculation at
80			$r_{\text{cut}}$ introduces a discontinuity into the potential
81			energy. To offset this discontinuity, the energy value at
82			$r_{\text{cut}}$ is subtracted from the entire potential. This causes
83			the equation to go to zero at the cut-off radius.
84
85			Interactions between dissimilar particles requires the generation of
86			cross term parameters for $\sigma$ and $\epsilon$. These are
87			calculated through the Lorentz-Berthelot mixing
88			rules:\cite{allen87:csl}
89			\begin{equation}
90			\sigma_{ij} = \frac{1}{2}[\sigma_{ii} + \sigma_{jj}]
91			\label{eq:sigmaMix}
92			\end{equation}
93			and
94			\begin{equation}
95			\epsilon_{ij} = \sqrt{\epsilon_{ii} \epsilon_{jj}}
96			\label{eq:epsilonMix}
97			\end{equation}
98
99
100	mmeineke	899	\subsection{\label{sec:DUFF}Dipolar Unified-Atom Force Field}
101
102			The \underline{D}ipolar \underline{U}nified-Atom
103			\underline{F}orce \underline{F}ield ({\sc duff}) was developed to
104			simulate lipid bilayers. We needed a model capable of forming
105			bilayers, while still being sufficiently computationally efficient to
106			allow simulations of large systems ($\approx$100's of phospholipids,
107			$\approx$1000's of waters) for long times ($\approx$10's of
108			nanoseconds).
109
110			With this goal in mind, we have eliminated all point charges. Charge
111			distributions were replaced with dipoles, and charge-neutral
112			distributions were reduced to Lennard-Jones interaction sites. This
113			simplification cuts the length scale of long range interactions from
114			$\frac{1}{r}$ to $\frac{1}{r^3}$, allowing us to avoid the
115			computationally expensive Ewald-Sum. Instead, we can use
116			neighbor-lists and cutoff radii for the dipolar interactions.
117
118			As an example, lipid head groups in {\sc duff} are represented as point
119			dipole interaction sites.PC and PE Lipid head groups are typically
120			zwitterionic in nature, with charges separated by as much as
121			6~$\mbox{\AA}$. By placing a dipole of 20.6~Debye at the head group
122			center of mass, our model mimics the head group of PC.\cite{Cevc87}
123			Additionally, a Lennard-Jones site is located at the
124			pseudoatom's center of mass. The model is illustrated by the dark grey
125			atom in Fig.~\ref{fig:lipidModel}.
126
127			\begin{figure}
128			\epsfxsize=6in
129			\epsfbox{lipidModel.epsi}
130			\caption{A representation of the lipid model. $\phi$ is the torsion angle, $\theta$ %
131			is the bend angle, $\mu$ is the dipole moment of the head group, and n is the chain length.}
132			\label{fig:lipidModel}
133			\end{figure}
134
135			The water model we use to complement the dipoles of the lipids is
136			the soft sticky dipole (SSD) model of Ichiye \emph{et
137			al.}\cite{liu96:new_model} This model is discussed in greater detail
138			in Sec.~\ref{sec:SSD}. In all cases we reduce water to a single
139			Lennard-Jones interaction site. The site also contains a dipole to
140			mimic the partial charges on the hydrogens and the oxygen. However,
141			what makes the SSD model unique is the inclusion of a tetrahedral
142			short range potential to recover the hydrogen bonding of water, an
143			important factor when modeling bilayers, as it has been shown that
144			hydrogen bond network formation is a leading contribution to the
145			entropic driving force towards lipid bilayer formation.\cite{Cevc87}
146
147
148			Turning to the tails of the phospholipids, we have used a set of
149			scalable parameters to model the alkyl groups with Lennard-Jones
150			sites. For this, we have used the TraPPE force field of Siepmann
151			\emph{et al}.\cite{Siepmann1998} TraPPE is a unified-atom
152			representation of n-alkanes, which is parametrized against phase
153			equilibria using Gibbs Monte Carlo simulation
154			techniques.\cite{Siepmann1998} One of the advantages of TraPPE is that
155			it generalizes the types of atoms in an alkyl chain to keep the number
156			of pseudoatoms to a minimum; the parameters for an atom such as
157			$\text{CH}_2$ do not change depending on what species are bonded to
158			it.
159
160			TraPPE also constrains of all bonds to be of fixed length. Typically,
161			bond vibrations are the fastest motions in a molecular dynamic
162			simulation. Small time steps between force evaluations must be used to
163			ensure adequate sampling of the bond potential conservation of
164			energy. By constraining the bond lengths, larger time steps may be
165			used when integrating the equations of motion.
166
167
168			\subsubsection{\label{subSec:energyFunctions}{\sc duff} Energy Functions}
169
170			The total energy of function in {\sc duff} is given by the following:
171			\begin{equation}
172			V_{\text{Total}} = \sum^{N}_{I=1} V^{I}_{\text{Internal}}
173			+ \sum^{N}_{I=1} \sum^{N}_{J=1} V^{IJ}_{\text{Cross}}
174			\label{eq:totalPotential}
175			\end{equation}
176			Where $V^{I}_{\text{Internal}}$ is the internal potential of a molecule:
177			\begin{equation}
178			V^{I}_{\text{Internal}} =
179			\sum_{\theta_{ijk} \in I} V_{\text{bend}}(\theta_{ijk})
180			+ \sum_{\phi_{ijkl} \in I} V_{\text{torsion}}(\theta_{ijkl})
181			+ \sum_{i \in I} \sum_{(j>i+4) \in I}
182			\biggl[ V_{\text{LJ}}(r_{ij}) + V_{\text{dipole}}
183			(\mathbf{r}_{ij},\boldsymbol{\Omega}_{i},\boldsymbol{\Omega}_{j})
184			\biggr]
185			\label{eq:internalPotential}
186			\end{equation}
187			Here $V_{\text{bend}}$ is the bend potential for all 1, 3 bonded pairs
188			within in the molecule. $V_{\text{torsion}}$ is the torsion The
189			pairwise portions of the internal potential are excluded for pairs
190			that ar closer than three bonds, i.e.~atom pairs farther away than a
191			torsion are included in the pair-wise loop.
192
193			The cross portion of the total potential, $V^{IJ}_{\text{Cross}}$, is
194			as follows:
195			\begin{equation}
196			V^{IJ}_{\text{Cross}} =
197			\sum_{i \in I} \sum_{j \in J}
198			\biggl[ V_{\text{LJ}}(r_{ij}) + V_{\text{dipole}}
199			(\mathbf{r}_{ij},\boldsymbol{\Omega}_{i},\boldsymbol{\Omega}_{j})
200			+ V_{\text{sticky}}
201			(\mathbf{r}_{ij},\boldsymbol{\Omega}_{i},\boldsymbol{\Omega}_{j})
202			\biggr]
203			\label{eq:crossPotentail}
204			\end{equation}
205			Where $V_{\text{LJ}}$ is the Lennard Jones potential,
206			$V_{\text{dipole}}$ is the dipole dipole potential, and
207			$V_{\text{sticky}}$ is the sticky potential defined by the SSD model.
208
209			The bend potential of a molecule is represented by the following function:
210			\begin{equation}
211			V_{\theta_{ijk}} = k_{\theta}( \theta_{ijk} - \theta_0 )^2 \label{eq:bendPot}
212			\end{equation}
213			Where $\theta_{ijk}$ is the angle defined by atoms $i$, $j$, and $k$
214			(see Fig.~\ref{fig:lipidModel}), and $\theta_0$ is the equilibrium
215			bond angle. $k_{\theta}$ is the force constant which determines the
216			strength of the harmonic bend. The parameters for $k_{\theta}$ and
217			$\theta_0$ are based off of those in TraPPE.\cite{Siepmann1998}
218
219			The torsion potential and parameters are also taken from TraPPE. It is
220			of the form:
221			\begin{equation}
222			V_{\text{torsion}}(\phi_{ijkl}) = c_1[1 + \cos \phi]
223			+ c_2[1 + \cos(2\phi)]
224			+ c_3[1 + \cos(3\phi)]
225			\label{eq:origTorsionPot}
226			\end{equation}
227			Here $\phi_{ijkl}$ is the angle defined by four bonded neighbors $i$,
228			$j$, $k$, and $l$ (again, see Fig.~\ref{fig:lipidModel}). However,
229			for computaional efficency, the torsion potentail has been recast
230			after the method of CHARMM\cite{charmm1983} whereby the angle series
231			is converted to a power series of the form:
232			\begin{equation}
233			V_{\text{torsion}}(\phi_{ijkl}) =
234			k_3 \cos^3 \phi + k_2 \cos^2 \phi + k_1 \cos \phi + k_0
235			\label{eq:torsionPot}
236			\end{equation}
237			Where:
238			\begin{align*}
239			k_0 &= c_1 + c_3 \\
240			k_1 &= c_1 - 3c_3 \\
241			k_2 &= 2 c_2 \\
242			k_3 &= 4c_3
243			\end{align*}
244			By recasting the equation to a power series, repeated trigonometric
245			evaluations are avoided during the calculation of the potential.
246
247
248	mmeineke	915
249	mmeineke	899	The dipole-dipole potential has the following form:
250			\begin{equation}
251			V_{\text{dipole}}(\mathbf{r}_{ij},\boldsymbol{\Omega}_{i},
252			\boldsymbol{\Omega}_{j}) = \frac{1}{4\pi\epsilon_{0}} \biggl[
253			\frac{\boldsymbol{\mu}_{i} \cdot \boldsymbol{\mu}_{j}}{r^{3}_{ij}}
254			-
255			\frac{3(\boldsymbol{\mu}_i \cdot \mathbf{r}_{ij}) %
256			(\boldsymbol{\mu}_j \cdot \mathbf{r}_{ij}) }
257			{r^{5}_{ij}} \biggr]
258			\label{eq:dipolePot}
259			\end{equation}
260			Here $\mathbf{r}_{ij}$ is the vector starting at atom $i$ pointing
261			towards $j$, and $\boldsymbol{\Omega}_i$ and $\boldsymbol{\Omega}_j$
262			are the Euler angles of atom $i$ and $j$
263			respectively. $\boldsymbol{\mu}_i$ is the dipole vector of atom
264			$i$ it takes its orientation from $\boldsymbol{\Omega}_i$.
265
266
267			\subsection{\label{sec:SSD}Water Model: SSD and Derivatives}
268
269			In the interest of computational efficiency, the native solvent used
270			in {\sc oopse} is the Soft Sticky Dipole (SSD) water model. SSD was
271			developed by Ichiye \emph{et al.} as a modified form of the
272			hard-sphere water model proposed by Bratko, Blum, and
273			Luzar.\cite{Bratko85,Bratko95} It consists of a single point dipole
274			with a Lennard-Jones core and a sticky potential that directs the
275			particles to assume the proper hydrogen bond orientation in the first
276			solvation shell. Thus, the interaction between two SSD water molecules
277			\emph{i} and \emph{j} is given by the potential
278			\begin{equation}
279			u_{ij} = u_{ij}^{LJ} (r_{ij})\ + u_{ij}^{dp}
280			(\mathbf{r}_{ij},\boldsymbol{\Omega}_i,\boldsymbol{\Omega}_j)\ +
281			u_{ij}^{sp}
282			(\mathbf{r}_{ij},\boldsymbol{\Omega}_i,\boldsymbol{\Omega}_j),
283			\end{equation}
284			where the $\mathbf{r}_{ij}$ is the position vector between molecules
285			\emph{i} and \emph{j} with magnitude equal to the distance $r_ij$, and
286			$\boldsymbol{\Omega}_i$ and $\boldsymbol{\Omega}_j$ represent the
287			orientations of the respective molecules. The Lennard-Jones, dipole,
288			and sticky parts of the potential are giving by the following
289			equations,
290			\begin{equation}
291			u_{ij}^{LJ}(r_{ij}) = 4\epsilon \left[\left(\frac{\sigma}{r_{ij}}\right)^{12}-\left(\frac{\sigma}{r_{ij}}\right)^{6}\right],
292			\end{equation}
293			\begin{equation}
294			u_{ij}^{dp} = \frac{\boldsymbol{\mu}_i\cdot\boldsymbol{\mu}_j}{r_{ij}^3}-\frac{3(\boldsymbol{\mu}_i\cdot\mathbf{r}_{ij})(\boldsymbol{\mu}_j\cdot\mathbf{r}_{ij})}{r_{ij}^5}\ ,
295			\end{equation}
296			\begin{equation}
297			\begin{split}
298			u_{ij}^{sp}
299			(\mathbf{r}_{ij},\boldsymbol{\Omega}_i,\boldsymbol{\Omega}_j)
300			&=
301			\frac{\nu_0}{2}[s(r_{ij})w(\mathbf{r}_{ij},\boldsymbol{\Omega}_i,\boldsymbol{\Omega}_j)\\
302			& \quad \ +
303			s^\prime(r_{ij})w^\prime(\mathbf{r}_{ij},\boldsymbol{\Omega}_i,\boldsymbol{\Omega}_j)]\ ,
304			\end{split}
305			\end{equation}
306			where $\boldsymbol{\mu}_i$ and $\boldsymbol{\mu}_j$ are the dipole
307			unit vectors of particles \emph{i} and \emph{j} with magnitude 2.35 D,
308			$\nu_0$ scales the strength of the overall sticky potential, $s$ and
309			$s^\prime$ are cubic switching functions. The $w$ and $w^\prime$
310			functions take the following forms,
311			\begin{equation}
312			w(\mathbf{r}_{ij},\boldsymbol{\Omega}_i,\boldsymbol{\Omega}_j)=\sin\theta_{ij}\sin2\theta_{ij}\cos2\phi_{ij},
313			\end{equation}
314			\begin{equation}
315			w^\prime(\mathbf{r}_{ij},\boldsymbol{\Omega}_i,\boldsymbol{\Omega}_j) = (\cos\theta_{ij}-0.6)^2(\cos\theta_{ij}+0.8)^2-w^0,
316			\end{equation}
317			where $w^0 = 0.07715$. The $w$ function is the tetrahedral attractive
318			term that promotes hydrogen bonding orientations within the first
319			solvation shell, and $w^\prime$ is a dipolar repulsion term that
320			repels unrealistic dipolar arrangements within the first solvation
321			shell. A more detailed description of the functional parts and
322			variables in this potential can be found in other
323			articles.\cite{liu96:new_model,chandra99:ssd_md}
324
325			Since SSD is a one-site point dipole model, the force calculations are
326			simplified significantly from a computational standpoint, in that the
327			number of long-range interactions is dramatically reduced. In the
328			original Monte Carlo simulations using this model, Ichiye \emph{et
329			al.} reported a calculation speed up of up to an order of magnitude
330			over other comparable models while maintaining the structural behavior
331			of water.\cite{liu96:new_model} In the original molecular dynamics studies of
332			SSD, it was shown that it actually improves upon the prediction of
333			water's dynamical properties over TIP3P and SPC/E.\cite{chandra99:ssd_md}
334
335			Recent constant pressure simulations revealed issues in the original
336			SSD model that led to lower than expected densities at all target
337			pressures.\cite{Ichiye03,Gezelter04} Reparameterizations of the
338			original SSD have resulted in improved density behavior, as well as
339			alterations in the water structure and transport behavior. {\sc oopse} is
340			easily modified to impliment these new potential parameter sets for
341			the derivative water models: SSD1, SSD/E, and SSD/RF. All of the
342			variable parameters are listed in the accompanying BASS file, and
343			these parameters simply need to be changed to the updated values.
344
345
346			\subsection{\label{sec:eam}Embedded Atom Model}
347
348			here there be Monsters
349	mmeineke	915
350			\subsection{\label{Sec:pbc}Periodic Boundary Conditions}
351
352			\textit{Periodic boundary conditions} are widely used to simulate truly
353			macroscopic systems with a relatively small number of particles. Simulation
354			box is replicated throughout space to form an infinite lattice. During the
355			simulation, when a particle moves in the primary cell, its periodic image
356			particles in other boxes move in exactly the same direction with exactly the
357			same orientation.Thus, as a particle leaves the primary cell, one of its
358			images will enter through the opposite face.If the simulation box is large
359			enough to avoid "feeling" the symmetric of the periodic lattice,the surface
360			effect could be ignored. Cubic and parallelepiped are the available periodic
361			cells. \bigskip In OOPSE, we use the matrix instead of the vector to describe
362			the property of the simulation box. Therefore, not only the size of the
363			simulation box could be changed during the simulation, but also the shape of
364			it. The transformation from box space vector $\overrightarrow{s}$ to its
365			corresponding real space vector $\overrightarrow{r}$ is defined by
366			\begin{equation}
367			\overrightarrow{r}=H\overrightarrow{s}%
368			\end{equation}
369
370
371			where $H=(h_{x},h_{y},h_{z})$ is a transformation matrix made up of the three
372			box axis vectors. $h_{x},h_{y}$ and $h_{z}$ represent the sides of the
373			simulation box respectively.
374
375			To find the minimum image, we need to convert the real vector to its
376			corresponding vector in box space first, \bigskip%
377			\begin{equation}
378			\overrightarrow{s}=H^{-1}\overrightarrow{r}%
379			\end{equation}
380			And then, each element of $\overrightarrow{s}$ is casted to lie between -0.5
381			to 0.5,
382			\begin{equation}
383			s_{i}^{\prime}=s_{i}-round(s_{i})
384			\end{equation}
385			where%
386
387			\begin{equation}
388			round(x)=\lfloor{x+0.5}\rfloor\text{ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ }if\text{
389			}x\geqslant0
390			\end{equation}
391			%
392
393			\begin{equation}
394			round(x)=\lceil{x-0.5}\rceil\text{ \ \ \ \ \ \ \ \ \ \ \ \ \ \ }if\text{ }x<0
395			\end{equation}
396
397
398			For example, $round(3.6)=4$,$round(3.1)=3$, $round(-3.6)=-4$, $round(-3.1)=-3$.
399
400			Finally, we could get the minimum image by transforming back to real space,%
401
402			\begin{equation}
403			\overrightarrow{r^{\prime}}=H^{-1}\overrightarrow{s^{\prime}}%
404			\end{equation}