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.

\begin{equation}
V^{\text{dipole}}_{ij}(\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:dipole}
\end{equation} 

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}

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.

The water model we use to complement this 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}

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

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

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

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

\begin{equation}
V_{\theta_{ijk}} = k_{\theta}( \theta_{ijk} - \theta_0 )^2 \label{eq:bendPot}
\end{equation}

\begin{equation}
V_{\phi_{ijkl}} =  ( k_3 \cos^3 \phi + k_2 \cos^2 \phi + k_1 \cos \phi + k_0)
\label{eq:torsionPot}
\end{equation}


The bonded interactions in the DUFF functional set are limited to the
bend potential and the torsional potential. Bond potentials are not
calculated, instead all bond lengths are fixed to allow for large time
steps to be taken between force evaluations.

The bend functional is of the form:

$k_{\theta}$, the force constant, and $\theta_0$, the equilibrium bend
angle, were taken from the TraPPE forcefield of Siepmann.

The torsion functional has the form:

Here, the authors decided to use a potential in terms of a power
expansion in $\cos \phi$ rather than the typical expansion in
$\phi$. This prevents the need for repeated trigonometric
evaluations. Again, all $k_n$ constants were based on those in TraPPE.

\subsection{Non-Bonded Interactions}
\label{subSec:nonBondedInteractions}

\begin{equation}
V_{\text{LJ}} = \text{internal + external}
\end{equation}


Revision:	737
Committed:	Mon Sep 1 19:50:07 2003 UTC (21 years ago) by mmeineke
Content type:	application/x-tex
File size:	5857 byte(s)
Log Message:	added some text to DUFF, changed the preamble from two column to preprint format.
#	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	\begin{equation}
21	mmeineke	713	V^{\text{dipole}}_{ij}(\mathbf{r}_{ij},\boldsymbol{\Omega}_{i},
22	mmeineke	737	\boldsymbol{\Omega}_{j}) =
23	mmeineke	713	\frac{1}{4\pi\epsilon_{0}} \biggl[
24			\frac{\boldsymbol{\mu}_{i} \cdot \boldsymbol{\mu}_{j}}{r^{3}_{ij}}
25			-
26			\frac{3(\boldsymbol{\mu}_i \cdot \mathbf{r}_{ij}) %
27			(\boldsymbol{\mu}_j \cdot \mathbf{r}_{ij}) }
28	mmeineke	737	{r^{5}_{ij}} \biggr]\label{eq:dipole}
29			\end{equation}
30	mmeineke	713
31	mmeineke	737	As an example, lipid head groups in DUFF are represented as point
32			dipole interaction sites.PC and PE Lipid head groups are typically
33			zwitterionic in nature, with charges separated by as much as
34			6~$\mbox{\AA}$. By placing a dipole of 20.6~Debye at the head group
35			center of mass, our model mimics the head group of PC.\cite{Cevc87}
36			Additionally, a Lennard-Jones site is located at the
37			pseudoatom's center of mass. The model is illustrated by the dark grey
38			atom in Fig.~\ref{fig:lipidModel}.
39	mmeineke	713
40	mmeineke	716	\begin{figure}
41			\includegraphics[angle=-90,width=\linewidth]{lipidModel.epsi}
42			\caption{A representation of the lipid model. $\phi$ is the torsion angle, $\theta$ %
43			is the bend angle, $\mu$ is the dipole moment of the head group, and n is the chain length.}
44			\label{fig:lipidModel}
45			\end{figure}
46
47	mmeineke	737	Turning to the tails of the phospholipids, we have used a set of
48			scalable parameters to model the alkyl groups with Lennard-Jones
49			sites. For this, we have used the TraPPE force field of Siepmann
50			\emph{et al}.\cite{Siepmann1998} TraPPE is a unified-atom
51			representation of n-alkanes, which is parametrized against phase
52			equilibria using Gibbs Monte Carlo simulation
53			techniques.\cite{Siepmann1998} One of the advantages of TraPPE is that
54			it generalizes the types of atoms in an alkyl chain to keep the number
55			of pseudoatoms to a minimum; the parameters for an atom such as
56			$\text{CH}_2$ do not change depending on what species are bonded to
57			it.
58	mmeineke	716
59	mmeineke	737	TraPPE also constrains of all bonds to be of fixed length. Typically,
60			bond vibrations are the fastest motions in a molecular dynamic
61			simulation. Small time steps between force evaluations must be used to
62			ensure adequate sampling of the bond potential conservation of
63			energy. By constraining the bond lengths, larger time steps may be
64			used when integrating the equations of motion.
65	mmeineke	716
66	mmeineke	737	The water model we use to complement this the dipoles of the lipids is
67			the soft sticky dipole (SSD) model of Ichiye \emph{et
68	mmeineke	717	al.}\cite{liu96:new_model} This model is discussed in greater detail
69	mmeineke	737	in Sec.~\ref{sec:SSD}. In all cases we reduce water to a single
70			Lennard-Jones interaction site. The site also contains a dipole to
71			mimic the partial charges on the hydrogens and the oxygen. However,
72			what makes the SSD model unique is the inclusion of a tetrahedral
73			short range potential to recover the hydrogen bonding of water, an
74			important factor when modeling bilayers, as it has been shown that
75			hydrogen bond network formation is a leading contribution to the
76			entropic driving force towards lipid bilayer formation.\cite{Cevc87}
77	mmeineke	717
78	mmeineke	737	\subsection{\label{subSec:energyFunctions}Energy Functions}
79	mmeineke	717
80	mmeineke	737	\begin{equation}
81			V_{\text{Total}} = \sum^{N}_{I=1} V^{I}_{\text{Internal}}
82			+ \sum^{N}_{I=1} \sum^{N}_{J=1} V^{IJ}_{\text{Cross}}
83			\label{eq:totalPotential}
84			\end{equation}
85	mmeineke	717
86	mmeineke	737	\begin{equation}
87			V^{I}_{\text{Internal}} =
88			\sum_{\theta_{ijk} \in I} V_{\text{bend}}(\theta_{ijk})
89			+ \sum_{\phi_{ijkl} \in I} V_{\text{torsion}}(\theta_{ijkl})
90			+ \sum_{i \in I} \sum_{(j>i+4) \in I}
91			\biggl[ V_{\text{LJ}}(r_{ij}) + V_{\text{dipole}}
92			(\mathbf{r}_{ij},\boldsymbol{\Omega}_{i},\boldsymbol{\Omega}_{j})
93			\biggr]
94			\label{eq:internalPotential}
95			\end{equation}
96	mmeineke	717
97	mmeineke	737	\begin{equation}
98			V^{IJ}_{\text{Cross}} =
99			\sum_{i \in I} \sum_{j \in J}
100			\biggl[ V_{\text{LJ}}(r_{ij}) + V_{\text{dipole}}
101			(\mathbf{r}_{ij},\boldsymbol{\Omega}_{i},\boldsymbol{\Omega}_{j})
102			+ V_{\text{sticky}}
103			(\mathbf{r}_{ij},\boldsymbol{\Omega}_{i},\boldsymbol{\Omega}_{j})
104			\biggr]
105			\label{eq:crossPotentail}
106			\end{equation}
107	mmeineke	717
108	mmeineke	737	\begin{equation}
109			V_{\theta_{ijk}} = k_{\theta}( \theta_{ijk} - \theta_0 )^2 \label{eq:bendPot}
110			\end{equation}
111	mmeineke	717
112	mmeineke	698	\begin{equation}
113	mmeineke	737	V_{\phi_{ijkl}} = ( k_3 \cos^3 \phi + k_2 \cos^2 \phi + k_1 \cos \phi + k_0)
114			\label{eq:torsionPot}
115	mmeineke	698	\end{equation}
116	mmeineke	666
117	mmeineke	664
118	mmeineke	698	The bonded interactions in the DUFF functional set are limited to the
119			bend potential and the torsional potential. Bond potentials are not
120			calculated, instead all bond lengths are fixed to allow for large time
121			steps to be taken between force evaluations.
122	mmeineke	666
123	mmeineke	698	The bend functional is of the form:
124	mmeineke	737
125	mmeineke	698	$k_{\theta}$, the force constant, and $\theta_0$, the equilibrium bend
126			angle, were taken from the TraPPE forcefield of Siepmann.
127
128			The torsion functional has the form:
129	mmeineke	737
130	mmeineke	698	Here, the authors decided to use a potential in terms of a power
131			expansion in $\cos \phi$ rather than the typical expansion in
132	mmeineke	717	$\phi$. This prevents the need for repeated trigonometric
133	mmeineke	698	evaluations. Again, all $k_n$ constants were based on those in TraPPE.
134
135			\subsection{Non-Bonded Interactions}
136			\label{subSec:nonBondedInteractions}
137
138			\begin{equation}
139			V_{\text{LJ}} = \text{internal + external}
140			\end{equation}
141
142