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 ({\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.


\subsection{\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 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:	818
Committed:	Fri Oct 24 21:27:59 2003 UTC (22 years ago) by gezelter
Content type:	application/x-tex
File size:	7696 byte(s)
Log Message:	Formatting changes to get out of RevTex nonsense
#	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	gezelter	818	\underline{F}orce \underline{F}ield ({\sc duff}) was developed to
6	mmeineke	737	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	gezelter	818	As an example, lipid head groups in {\sc duff} are represented as point
21	mmeineke	737	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	gezelter	818	\epsfxsize=6in
31			\epsfbox{lipidModel.epsi}
32	mmeineke	716	\caption{A representation of the lipid model. $\phi$ is the torsion angle, $\theta$ %
33			is the bend angle, $\mu$ is the dipole moment of the head group, and n is the chain length.}
34			\label{fig:lipidModel}
35			\end{figure}
36
37	mmeineke	740	The water model we use to complement the dipoles of the lipids is
38			the soft sticky dipole (SSD) model of Ichiye \emph{et
39			al.}\cite{liu96:new_model} This model is discussed in greater detail
40			in Sec.~\ref{sec:SSD}. In all cases we reduce water to a single
41			Lennard-Jones interaction site. The site also contains a dipole to
42			mimic the partial charges on the hydrogens and the oxygen. However,
43			what makes the SSD model unique is the inclusion of a tetrahedral
44			short range potential to recover the hydrogen bonding of water, an
45			important factor when modeling bilayers, as it has been shown that
46			hydrogen bond network formation is a leading contribution to the
47			entropic driving force towards lipid bilayer formation.\cite{Cevc87}
48
49
50	mmeineke	737	Turning to the tails of the phospholipids, we have used a set of
51			scalable parameters to model the alkyl groups with Lennard-Jones
52			sites. For this, we have used the TraPPE force field of Siepmann
53			\emph{et al}.\cite{Siepmann1998} TraPPE is a unified-atom
54			representation of n-alkanes, which is parametrized against phase
55			equilibria using Gibbs Monte Carlo simulation
56			techniques.\cite{Siepmann1998} One of the advantages of TraPPE is that
57			it generalizes the types of atoms in an alkyl chain to keep the number
58			of pseudoatoms to a minimum; the parameters for an atom such as
59			$\text{CH}_2$ do not change depending on what species are bonded to
60			it.
61	mmeineke	716
62	mmeineke	737	TraPPE also constrains of all bonds to be of fixed length. Typically,
63			bond vibrations are the fastest motions in a molecular dynamic
64			simulation. Small time steps between force evaluations must be used to
65			ensure adequate sampling of the bond potential conservation of
66			energy. By constraining the bond lengths, larger time steps may be
67			used when integrating the equations of motion.
68	mmeineke	716
69	mmeineke	717
70	gezelter	818	\subsection{\label{subSec:energyFunctions}{\sc duff} Energy Functions}
71	mmeineke	717
72	gezelter	818	The total energy of function in {\sc duff} is given by the following:
73	mmeineke	737	\begin{equation}
74			V_{\text{Total}} = \sum^{N}_{I=1} V^{I}_{\text{Internal}}
75			+ \sum^{N}_{I=1} \sum^{N}_{J=1} V^{IJ}_{\text{Cross}}
76			\label{eq:totalPotential}
77			\end{equation}
78	mmeineke	740	Where $V^{I}_{\text{Internal}}$ is the internal potential of a molecule:
79	mmeineke	737	\begin{equation}
80			V^{I}_{\text{Internal}} =
81			\sum_{\theta_{ijk} \in I} V_{\text{bend}}(\theta_{ijk})
82			+ \sum_{\phi_{ijkl} \in I} V_{\text{torsion}}(\theta_{ijkl})
83			+ \sum_{i \in I} \sum_{(j>i+4) \in I}
84			\biggl[ V_{\text{LJ}}(r_{ij}) + V_{\text{dipole}}
85			(\mathbf{r}_{ij},\boldsymbol{\Omega}_{i},\boldsymbol{\Omega}_{j})
86			\biggr]
87			\label{eq:internalPotential}
88			\end{equation}
89	mmeineke	740	Here $V_{\text{bend}}$ is the bend potential for all 1, 3 bonded pairs
90			within in the molecule. $V_{\text{torsion}}$ is the torsion The
91			pairwise portions of the internal potential are excluded for pairs
92			that ar closer than three bonds, i.e.~atom pairs farther away than a
93			torsion are included in the pair-wise loop.
94	mmeineke	717
95	mmeineke	740	The cross portion of the total potential, $V^{IJ}_{\text{Cross}}$, is
96			as follows:
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	740	Where $V_{\text{LJ}}$ is the Lennard Jones potential,
108			$V_{\text{dipole}}$ is the dipole dipole potential, and
109			$V_{\text{sticky}}$ is the sticky potential defined by the SSD model.
110	mmeineke	717
111	mmeineke	740	The bend potential of a molecule is represented by the following function:
112	mmeineke	737	\begin{equation}
113			V_{\theta_{ijk}} = k_{\theta}( \theta_{ijk} - \theta_0 )^2 \label{eq:bendPot}
114			\end{equation}
115	mmeineke	740	Where $\theta_{ijk}$ is the angle defined by atoms $i$, $j$, and $k$
116			(see Fig.~\ref{fig:lipidModel}), and $\theta_0$ is the equilibrium
117			bond angle. $k_{\theta}$ is the force constant which determines the
118			strength of the harmonic bend. The parameters for $k_{\theta}$ and
119			$\theta_0$ are based off of those in TraPPE.\cite{Siepmann1998}
120	mmeineke	717
121	mmeineke	740	The torsion potential and parameters are also taken from TraPPE. It is
122			of the form:
123	mmeineke	698	\begin{equation}
124	mmeineke	740	V_{\text{torsion}}(\phi_{ijkl}) = c_1[1 + \cos \phi]
125			+ c_2[1 + \cos(2\phi)]
126			+ c_3[1 + \cos(3\phi)]
127			\label{eq:origTorsionPot}
128			\end{equation}
129			Here $\phi_{ijkl}$ is the angle defined by four bonded neighbors $i$,
130			$j$, $k$, and $l$ (again, see Fig.~\ref{fig:lipidModel}). However,
131			for computaional efficency, the torsion potentail has been recast
132			after the method of CHARMM\cite{charmm1983} whereby the angle series
133			is converted to a power series of the form:
134			\begin{equation}
135			V_{\text{torsion}}(\phi_{ijkl}) =
136			k_3 \cos^3 \phi + k_2 \cos^2 \phi + k_1 \cos \phi + k_0
137	mmeineke	737	\label{eq:torsionPot}
138	mmeineke	698	\end{equation}
139	mmeineke	740	Where:
140			\begin{align*}
141			k_0 &= c_1 + c_3 \\
142			k_1 &= c_1 - 3c_3 \\
143			k_2 &= 2 c_2 \\
144			k_3 &= 4c_3
145			\end{align*}
146			By recasting the equation to a power series, repeated trigonometric
147			evaluations are avoided during the calculation of the potential.
148	mmeineke	666
149	mmeineke	740	The Lennard-Jones potential is given by:
150			\begin{equation}
151			V_{\text{LJ}}(r_{ij}) =
152			4\epsilon_{ij} \biggl[
153			\biggl(\frac{\sigma_{ij}}{r_{ij}}\biggr)^{12}
154			- \biggl(\frac{\sigma_{ij}}{r_{ij}}\biggr)^{6}
155			\biggr]
156			\label{eq:lennardJonesPot}
157			\end{equation}
158			Where $r_ij$ is the distance between atoms $i$ and $j$, $\sigma_{ij}$
159			scales the length of the interaction, and $\epsilon_{ij}$ scales the
160			energy of the potential.
161	mmeineke	664
162	mmeineke	740	The dipole-dipole potential has the following form:
163	mmeineke	698	\begin{equation}
164	mmeineke	740	V_{\text{dipole}}(\mathbf{r}_{ij},\boldsymbol{\Omega}_{i},
165			\boldsymbol{\Omega}_{j}) = \frac{1}{4\pi\epsilon_{0}} \biggl[
166			\frac{\boldsymbol{\mu}_{i} \cdot \boldsymbol{\mu}_{j}}{r^{3}_{ij}}
167			-
168			\frac{3(\boldsymbol{\mu}_i \cdot \mathbf{r}_{ij}) %
169			(\boldsymbol{\mu}_j \cdot \mathbf{r}_{ij}) }
170			{r^{5}_{ij}} \biggr]
171			\label{eq:dipolePot}
172	mmeineke	698	\end{equation}
173	mmeineke	740	Here $\mathbf{r}_{ij}$ is the vector starting at atom $i$ pointing
174			towards $j$, and $\boldsymbol{\Omega}_i$ and $\boldsymbol{\Omega}_j$
175			are the Euler angles of atom $i$ and $j$
176			respectively. $\boldsymbol{\mu}_i$ is the dipole vector of atom
177			$i$ it takes its orientation from $\boldsymbol{\Omega}_i$.