trunk/oopsePaper/DUFF.tex


\section{\label{sec:DUFF}The DUFF Force Field}

The DUFF (\underline{D}ipolar \underline{U}nified-atom
\underline{F}orce \underline{F}ield) force field was developed to
simulate lipid bilayer formation and equilibrium dynamics. We needed a
model capable of forming bilayers, while still being sufficiently
computationally efficient allowing simulations of large systems
(\~100's of phospholipids, \~1000's of waters) for long times (\~10's
of nanoseconds).

With this goal in mind, we decided to eliminate all charged
interactions within the force field. Charge distributions were
replaced with dipolar entities, 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}$ (Eq.~\ref{eq:dipole} vs.~ Eq.~\ref{eq:coloumb}),
allowing us to avoid the computationally expensive Ewald-Sum. Instead,
we can use neighbor-lists and cutoff radii for the dipolar
interactions.

\begin{align}
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} \\
V^{\text{ch}}_{ij}(\mathbf{r}_{ij}) &= \frac{q_{i}q_{j}}%
        {4\pi\epsilon_{0} r_{ij}} \label{eq:coloumb}
\end{align} 

Applying this standard to the lipid model, we decided to represent the
lipid model as a point dipole interaction site. Lipid head groups are
typically zwitterionic in nature, with sometimes full integer charges
separated by only 5 to 6~$\mbox{\AA}$. By placing a dipole of
20.6~Debye at the head groups center of mass, our model mimics the
dipole of DMPC.\cite{Cevc87} Then, to account for the steric hindrance
of the head group, a Lennard-Jones interaction site is also located at
the pseudoatom's center of mass. The model is illustrated 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 tail chains of the phospholipids, we needed a set of
scalable parameters to model the alkyl groups as Lennard-Jones
interaction sites. For this, we used the TraPPE force field of
Siepmann \emph{et al}.\cite{Siepmann1998} The force field is a
unified-atom representation of n-alkanes. It is parametrized against
phase equilibria using Gibbs Monte Carlo simulation techniques. One of
the advantages of TraPPE is that is generalizes the types of atoms in
an alkyl chain to keep the number of pseudoatoms to a minimum; the
$\mbox{CH}_2$ in propane is the same as the central and offset
$\mbox{CH}_2$'s in pentane, meaning the pseudoatom type does not
change according to the atom's environment.

Another advantage of using TraPPE is the constraining of all bonds to
be of fixed length. Typically, bond vibrations are the motions in a
molecular dynamic simulation. This necessitates a small time step
between force evaluations be used to ensure adequate sampling of the
bond potential. Failure to do so will result in loss of energy
conservation within the microcanonical ensemble. By constraining this
degree of freedom, time steps larger than were previously allowable
are able to be used when integrating the equations of motion.

After developing the model for the phospholipids, we needed a model
for water that would complement our lipid. For this we turned to 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}. The basic idea of the model is to 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}

BREAK 

END OF CURRENT REVISIONS 

BREAK


The main energy function in OOPSE is DUFF (the Dipolar Unified-atom
Force Field). DUFF is a collection of parameters taken from Seipmann
 and  The total energy of interaction is given by
Eq.~\ref{eq:totalPotential}:
\begin{equation}
V_{\text{Total}} =
        \overbrace{V_{\theta} + V_{\phi}}^{\text{bonded}} +
        \underbrace{V_{\text{LJ}} + V_{\text{Dp}} + %
        V_{\text{SSD}}}_{\text{non-bonded}} \label{eq:totalPotential}
\end{equation}

\subsection{Bonded Interactions}
\label{subSec:bondedInteractions}

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:
\begin{equation}
V_{\theta} = \sum k_{\theta}( \theta - \theta_0 )^2 \label{eq:bendPot}
\end{equation}
$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:
\begin{equation}
V_{\phi} =  \sum ( k_3 \cos^3 \phi + k_2 \cos^2 \phi + k_1 \cos \phi + k_0)
\label{eq:torsionPot}
\end{equation}
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:	717
Committed:	Mon Aug 25 20:25:12 2003 UTC (21 years ago) by mmeineke
Content type:	application/x-tex
File size:	6034 byte(s)
Log Message:	did some more work on DUFF.
#	User	Rev	Content
1	mmeineke	664
2	mmeineke	713	\section{\label{sec:DUFF}The DUFF Force Field}
3	mmeineke	664
4	mmeineke	713	The DUFF (\underline{D}ipolar \underline{U}nified-atom
5			\underline{F}orce \underline{F}ield) force field was developed to
6			simulate lipid bilayer formation and equilibrium dynamics. We needed a
7	mmeineke	717	model capable of forming bilayers, while still being sufficiently
8	mmeineke	713	computationally efficient allowing simulations of large systems
9			(\~100's of phospholipids, \~1000's of waters) for long times (\~10's
10			of nanoseconds).
11	mmeineke	710
12	mmeineke	713	With this goal in mind, we decided to eliminate all charged
13			interactions within the force field. Charge distributions were
14			replaced with dipolar entities, and charge neutral distributions were
15			reduced to Lennard-Jones interaction sites. This simplification cuts
16			the length scale of long range interactions from $\frac{1}{r}$ to
17	mmeineke	716	$\frac{1}{r^3}$ (Eq.~\ref{eq:dipole} vs.~ Eq.~\ref{eq:coloumb}),
18			allowing us to avoid the computationally expensive Ewald-Sum. Instead,
19			we can use neighbor-lists and cutoff radii for the dipolar
20			interactions.
21	mmeineke	710
22	mmeineke	713	\begin{align}
23			V^{\text{dipole}}_{ij}(\mathbf{r}_{ij},\boldsymbol{\Omega}_{i},
24			\boldsymbol{\Omega}_{j}) &=
25			\frac{1}{4\pi\epsilon_{0}} \biggl[
26			\frac{\boldsymbol{\mu}_{i} \cdot \boldsymbol{\mu}_{j}}{r^{3}_{ij}}
27			-
28			\frac{3(\boldsymbol{\mu}_i \cdot \mathbf{r}_{ij}) %
29			(\boldsymbol{\mu}_j \cdot \mathbf{r}_{ij}) }
30			{r^{5}_{ij}} \biggr]\label{eq:dipole} \\
31			V^{\text{ch}}_{ij}(\mathbf{r}_{ij}) &= \frac{q_{i}q_{j}}%
32			{4\pi\epsilon_{0} r_{ij}} \label{eq:coloumb}
33			\end{align}
34
35	mmeineke	716	Applying this standard to the lipid model, we decided to represent the
36			lipid model as a point dipole interaction site. Lipid head groups are
37			typically zwitterionic in nature, with sometimes full integer charges
38	mmeineke	717	separated by only 5 to 6~$\mbox{\AA}$. By placing a dipole of
39	mmeineke	716	20.6~Debye at the head groups center of mass, our model mimics the
40	mmeineke	717	dipole of DMPC.\cite{Cevc87} Then, to account for the steric hindrance
41			of the head group, a Lennard-Jones interaction site is also located at
42			the pseudoatom's center of mass. The model is illustrated in
43	mmeineke	716	Fig.~\ref{fig:lipidModel}.
44	mmeineke	713
45	mmeineke	716	\begin{figure}
46			\includegraphics[angle=-90,width=\linewidth]{lipidModel.epsi}
47			\caption{A representation of the lipid model. $\phi$ is the torsion angle, $\theta$ %
48			is the bend angle, $\mu$ is the dipole moment of the head group, and n is the chain length.}
49			\label{fig:lipidModel}
50			\end{figure}
51
52			Turning to the tail chains of the phospholipids, we needed a set of
53			scalable parameters to model the alkyl groups as Lennard-Jones
54			interaction sites. For this, we used the TraPPE force field of
55			Siepmann \emph{et al}.\cite{Siepmann1998} The force field is a
56			unified-atom representation of n-alkanes. It is parametrized against
57			phase equilibria using Gibbs Monte Carlo simulation techniques. One of
58			the advantages of TraPPE is that is generalizes the types of atoms in
59	mmeineke	717	an alkyl chain to keep the number of pseudoatoms to a minimum; the
60			$\mbox{CH}_2$ in propane is the same as the central and offset
61			$\mbox{CH}_2$'s in pentane, meaning the pseudoatom type does not
62			change according to the atom's environment.
63	mmeineke	716
64			Another advantage of using TraPPE is the constraining of all bonds to
65			be of fixed length. Typically, bond vibrations are the motions in a
66	mmeineke	717	molecular dynamic simulation. This necessitates a small time step
67	mmeineke	716	between force evaluations be used to ensure adequate sampling of the
68			bond potential. Failure to do so will result in loss of energy
69			conservation within the microcanonical ensemble. By constraining this
70			degree of freedom, time steps larger than were previously allowable
71			are able to be used when integrating the equations of motion.
72
73	mmeineke	717	After developing the model for the phospholipids, we needed a model
74			for water that would complement our lipid. For this we turned to the
75			soft sticky dipole (SSD) model of Ichiye \emph{et
76			al.}\cite{liu96:new_model} This model is discussed in greater detail
77			in Sec.~\ref{sec:SSD}. The basic idea of the model is to reduce water
78			to a single Lennard-Jones interaction site. The site also contains a
79			dipole to mimic the partial charges on the hydrogens and the
80			oxygen. However, what makes the SSD model unique is the inclusion of a
81			tetrahedral short range potential to recover the hydrogen bonding of
82			water, an important factor when modeling bilayers, as it has been
83			shown that hydrogen bond network formation is a leading contribution
84			to the entropic driving force towards lipid bilayer
85			formation.\cite{Cevc87}
86
87			BREAK
88
89			END OF CURRENT REVISIONS
90
91			BREAK
92
93
94
95
96
97	mmeineke	716	The main energy function in OOPSE is DUFF (the Dipolar Unified-atom
98			Force Field). DUFF is a collection of parameters taken from Seipmann
99	mmeineke	717	and The total energy of interaction is given by
100	mmeineke	666	Eq.~\ref{eq:totalPotential}:
101	mmeineke	698	\begin{equation}
102			V_{\text{Total}} =
103			\overbrace{V_{\theta} + V_{\phi}}^{\text{bonded}} +
104			\underbrace{V_{\text{LJ}} + V_{\text{Dp}} + %
105			V_{\text{SSD}}}_{\text{non-bonded}} \label{eq:totalPotential}
106			\end{equation}
107	mmeineke	666
108	mmeineke	698	\subsection{Bonded Interactions}
109			\label{subSec:bondedInteractions}
110	mmeineke	664
111	mmeineke	698	The bonded interactions in the DUFF functional set are limited to the
112			bend potential and the torsional potential. Bond potentials are not
113			calculated, instead all bond lengths are fixed to allow for large time
114			steps to be taken between force evaluations.
115	mmeineke	666
116	mmeineke	698	The bend functional is of the form:
117			\begin{equation}
118			V_{\theta} = \sum k_{\theta}( \theta - \theta_0 )^2 \label{eq:bendPot}
119			\end{equation}
120			$k_{\theta}$, the force constant, and $\theta_0$, the equilibrium bend
121			angle, were taken from the TraPPE forcefield of Siepmann.
122
123			The torsion functional has the form:
124			\begin{equation}
125	mmeineke	709	V_{\phi} = \sum ( k_3 \cos^3 \phi + k_2 \cos^2 \phi + k_1 \cos \phi + k_0)
126	mmeineke	698	\label{eq:torsionPot}
127			\end{equation}
128			Here, the authors decided to use a potential in terms of a power
129			expansion in $\cos \phi$ rather than the typical expansion in
130	mmeineke	717	$\phi$. This prevents the need for repeated trigonometric
131	mmeineke	698	evaluations. Again, all $k_n$ constants were based on those in TraPPE.
132
133			\subsection{Non-Bonded Interactions}
134			\label{subSec:nonBondedInteractions}
135
136			\begin{equation}
137			V_{\text{LJ}} = \text{internal + external}
138			\end{equation}
139
140