trunk/oopsePaper/DUFF.tex


\section{\label{sec:DUFF}Dipolar Unified-Atom Force Field}

The \underline{D}ipolar \underline{U}nified-Atom
\underline{F}orce \underline{F}ield (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 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}
\includegraphics[angle=-90,width=\linewidth]{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.


\subsection{\label{subSec:energyFunctions}DUFF Energy Functions}

The total energy of function in 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 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 atoms $i$ and $j$, $\sigma_{ij}$
scales the length of the interaction, and $\epsilon_{ij}$ scales the
energy 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$.
Revision:	740
Committed:	Tue Sep 2 18:40:47 2003 UTC (21 years ago) by mmeineke
Content type:	application/x-tex
File size:	7693 byte(s)
Log Message:	added a charmm citation to the bib file. finished adding all of the equations into DUFF.
#	User	Rev	Content
1	mmeineke	664
2	mmeineke	737	\section{\label{sec:DUFF}Dipolar Unified-Atom Force Field}
3	mmeineke	664
4	mmeineke	737	The \underline{D}ipolar \underline{U}nified-Atom
5			\underline{F}orce \underline{F}ield (DUFF) was developed to
6			simulate lipid bilayers. We needed a model capable of forming
7			bilayers, while still being sufficiently computationally efficient to
8			allow simulations of large systems ($\approx$100's of phospholipids,
9			$\approx$1000's of waters) for long times ($\approx$10's of
10			nanoseconds).
11	mmeineke	710
12	mmeineke	737	With this goal in mind, we have eliminated all point charges. Charge
13			distributions were replaced with dipoles, and charge-neutral
14			distributions were reduced to Lennard-Jones interaction sites. This
15			simplification cuts the length scale of long range interactions from
16			$\frac{1}{r}$ to $\frac{1}{r^3}$, allowing us to avoid the
17			computationally expensive Ewald-Sum. Instead, we can use
18			neighbor-lists and cutoff radii for the dipolar interactions.
19	mmeineke	710
20	mmeineke	737	As an example, lipid head groups in DUFF are represented as point
21			dipole interaction sites.PC and PE Lipid head groups are typically
22			zwitterionic in nature, with charges separated by as much as
23			6~$\mbox{\AA}$. By placing a dipole of 20.6~Debye at the head group
24			center of mass, our model mimics the head group of PC.\cite{Cevc87}
25			Additionally, a Lennard-Jones site is located at the
26			pseudoatom's center of mass. The model is illustrated by the dark grey
27			atom in Fig.~\ref{fig:lipidModel}.
28	mmeineke	713
29	mmeineke	716	\begin{figure}
30			\includegraphics[angle=-90,width=\linewidth]{lipidModel.epsi}
31			\caption{A representation of the lipid model. $\phi$ is the torsion angle, $\theta$ %
32			is the bend angle, $\mu$ is the dipole moment of the head group, and n is the chain length.}
33			\label{fig:lipidModel}
34			\end{figure}
35
36	mmeineke	740	The water model we use to complement the dipoles of the lipids is
37			the soft sticky dipole (SSD) model of Ichiye \emph{et
38			al.}\cite{liu96:new_model} This model is discussed in greater detail
39			in Sec.~\ref{sec:SSD}. In all cases we reduce water to a single
40			Lennard-Jones interaction site. The site also contains a dipole to
41			mimic the partial charges on the hydrogens and the oxygen. However,
42			what makes the SSD model unique is the inclusion of a tetrahedral
43			short range potential to recover the hydrogen bonding of water, an
44			important factor when modeling bilayers, as it has been shown that
45			hydrogen bond network formation is a leading contribution to the
46			entropic driving force towards lipid bilayer formation.\cite{Cevc87}
47
48
49	mmeineke	737	Turning to the tails of the phospholipids, we have used a set of
50			scalable parameters to model the alkyl groups with Lennard-Jones
51			sites. For this, we have used the TraPPE force field of Siepmann
52			\emph{et al}.\cite{Siepmann1998} TraPPE is a unified-atom
53			representation of n-alkanes, which is parametrized against phase
54			equilibria using Gibbs Monte Carlo simulation
55			techniques.\cite{Siepmann1998} One of the advantages of TraPPE is that
56			it generalizes the types of atoms in an alkyl chain to keep the number
57			of pseudoatoms to a minimum; the parameters for an atom such as
58			$\text{CH}_2$ do not change depending on what species are bonded to
59			it.
60	mmeineke	716
61	mmeineke	737	TraPPE also constrains of all bonds to be of fixed length. Typically,
62			bond vibrations are the fastest motions in a molecular dynamic
63			simulation. Small time steps between force evaluations must be used to
64			ensure adequate sampling of the bond potential conservation of
65			energy. By constraining the bond lengths, larger time steps may be
66			used when integrating the equations of motion.
67	mmeineke	716
68	mmeineke	717
69	mmeineke	740	\subsection{\label{subSec:energyFunctions}DUFF Energy Functions}
70	mmeineke	717
71	mmeineke	740	The total energy of function in DUFF is given by the following:
72	mmeineke	737	\begin{equation}
73			V_{\text{Total}} = \sum^{N}_{I=1} V^{I}_{\text{Internal}}
74			+ \sum^{N}_{I=1} \sum^{N}_{J=1} V^{IJ}_{\text{Cross}}
75			\label{eq:totalPotential}
76			\end{equation}
77	mmeineke	740	Where $V^{I}_{\text{Internal}}$ is the internal potential of a molecule:
78	mmeineke	737	\begin{equation}
79			V^{I}_{\text{Internal}} =
80			\sum_{\theta_{ijk} \in I} V_{\text{bend}}(\theta_{ijk})
81			+ \sum_{\phi_{ijkl} \in I} V_{\text{torsion}}(\theta_{ijkl})
82			+ \sum_{i \in I} \sum_{(j>i+4) \in I}
83			\biggl[ V_{\text{LJ}}(r_{ij}) + V_{\text{dipole}}
84			(\mathbf{r}_{ij},\boldsymbol{\Omega}_{i},\boldsymbol{\Omega}_{j})
85			\biggr]
86			\label{eq:internalPotential}
87			\end{equation}
88	mmeineke	740	Here $V_{\text{bend}}$ is the bend potential for all 1, 3 bonded pairs
89			within in the molecule. $V_{\text{torsion}}$ is the torsion The
90			pairwise portions of the internal potential are excluded for pairs
91			that ar closer than three bonds, i.e.~atom pairs farther away than a
92			torsion are included in the pair-wise loop.
93	mmeineke	717
94	mmeineke	740	The cross portion of the total potential, $V^{IJ}_{\text{Cross}}$, is
95			as follows:
96	mmeineke	737	\begin{equation}
97			V^{IJ}_{\text{Cross}} =
98			\sum_{i \in I} \sum_{j \in J}
99			\biggl[ V_{\text{LJ}}(r_{ij}) + V_{\text{dipole}}
100			(\mathbf{r}_{ij},\boldsymbol{\Omega}_{i},\boldsymbol{\Omega}_{j})
101			+ V_{\text{sticky}}
102			(\mathbf{r}_{ij},\boldsymbol{\Omega}_{i},\boldsymbol{\Omega}_{j})
103			\biggr]
104			\label{eq:crossPotentail}
105			\end{equation}
106	mmeineke	740	Where $V_{\text{LJ}}$ is the Lennard Jones potential,
107			$V_{\text{dipole}}$ is the dipole dipole potential, and
108			$V_{\text{sticky}}$ is the sticky potential defined by the SSD model.
109	mmeineke	717
110	mmeineke	740	The bend potential of a molecule is represented by the following function:
111	mmeineke	737	\begin{equation}
112			V_{\theta_{ijk}} = k_{\theta}( \theta_{ijk} - \theta_0 )^2 \label{eq:bendPot}
113			\end{equation}
114	mmeineke	740	Where $\theta_{ijk}$ is the angle defined by atoms $i$, $j$, and $k$
115			(see Fig.~\ref{fig:lipidModel}), and $\theta_0$ is the equilibrium
116			bond angle. $k_{\theta}$ is the force constant which determines the
117			strength of the harmonic bend. The parameters for $k_{\theta}$ and
118			$\theta_0$ are based off of those in TraPPE.\cite{Siepmann1998}
119	mmeineke	717
120	mmeineke	740	The torsion potential and parameters are also taken from TraPPE. It is
121			of the form:
122	mmeineke	698	\begin{equation}
123	mmeineke	740	V_{\text{torsion}}(\phi_{ijkl}) = c_1[1 + \cos \phi]
124			+ c_2[1 + \cos(2\phi)]
125			+ c_3[1 + \cos(3\phi)]
126			\label{eq:origTorsionPot}
127			\end{equation}
128			Here $\phi_{ijkl}$ is the angle defined by four bonded neighbors $i$,
129			$j$, $k$, and $l$ (again, see Fig.~\ref{fig:lipidModel}). However,
130			for computaional efficency, the torsion potentail has been recast
131			after the method of CHARMM\cite{charmm1983} whereby the angle series
132			is converted to a power series of the form:
133			\begin{equation}
134			V_{\text{torsion}}(\phi_{ijkl}) =
135			k_3 \cos^3 \phi + k_2 \cos^2 \phi + k_1 \cos \phi + k_0
136	mmeineke	737	\label{eq:torsionPot}
137	mmeineke	698	\end{equation}
138	mmeineke	740	Where:
139			\begin{align*}
140			k_0 &= c_1 + c_3 \\
141			k_1 &= c_1 - 3c_3 \\
142			k_2 &= 2 c_2 \\
143			k_3 &= 4c_3
144			\end{align*}
145			By recasting the equation to a power series, repeated trigonometric
146			evaluations are avoided during the calculation of the potential.
147	mmeineke	666
148	mmeineke	740	The Lennard-Jones potential is given by:
149			\begin{equation}
150			V_{\text{LJ}}(r_{ij}) =
151			4\epsilon_{ij} \biggl[
152			\biggl(\frac{\sigma_{ij}}{r_{ij}}\biggr)^{12}
153			- \biggl(\frac{\sigma_{ij}}{r_{ij}}\biggr)^{6}
154			\biggr]
155			\label{eq:lennardJonesPot}
156			\end{equation}
157			Where $r_ij$ is the distance between atoms $i$ and $j$, $\sigma_{ij}$
158			scales the length of the interaction, and $\epsilon_{ij}$ scales the
159			energy of the potential.
160	mmeineke	664
161	mmeineke	740	The dipole-dipole potential has the following form:
162	mmeineke	698	\begin{equation}
163	mmeineke	740	V_{\text{dipole}}(\mathbf{r}_{ij},\boldsymbol{\Omega}_{i},
164			\boldsymbol{\Omega}_{j}) = \frac{1}{4\pi\epsilon_{0}} \biggl[
165			\frac{\boldsymbol{\mu}_{i} \cdot \boldsymbol{\mu}_{j}}{r^{3}_{ij}}
166			-
167			\frac{3(\boldsymbol{\mu}_i \cdot \mathbf{r}_{ij}) %
168			(\boldsymbol{\mu}_j \cdot \mathbf{r}_{ij}) }
169			{r^{5}_{ij}} \biggr]
170			\label{eq:dipolePot}
171	mmeineke	698	\end{equation}
172	mmeineke	740	Here $\mathbf{r}_{ij}$ is the vector starting at atom $i$ pointing
173			towards $j$, and $\boldsymbol{\Omega}_i$ and $\boldsymbol{\Omega}_j$
174			are the Euler angles of atom $i$ and $j$
175			respectively. $\boldsymbol{\mu}_i$ is the dipole vector of atom
176			$i$ it takes its orientation from $\boldsymbol{\Omega}_i$.