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 bilaers, 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
seperated 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 henderanc
of the head group, a Lennard-Jones interaction site is also oacted at
the psuedoatom'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.
%( $ \mbox{CH_3} $ %-$\mathbf{\mbox{CH_2}}$-$\mbox{CH_3}$ is the same as

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 neccesitates 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.

The main energy function in OOPSE is DUFF (the Dipolar Unified-atom
Force Field). DUFF is a collection of parameters taken from Seipmann
 and Ichiye \emph{et
al.}\cite{liu96:new_model} 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 trigonemtric
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:	716
Committed:	Sun Aug 24 04:00:44 2003 UTC (21 years ago) by mmeineke
Content type:	application/x-tex
File size:	5116 byte(s)
Log Message:	added the lipid figure, and did some 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			model capable of forming bilaers, while still being sufficiently
8			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			seperated by only 5 to 6~$\mbox{\AA}$. By placing a dipole of
39			20.6~Debye at the head groups center of mass, our model mimics the
40			dipole of DMPC.\cite{Cevc87} Then, to account for the steric henderanc
41			of the head group, a Lennard-Jones interaction site is also oacted at
42			the psuedoatom's center of mass. The model is illustrated in
43			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			an alkyl chain to keep the number of pseudoatoms to a minimum.
60			%( $ \mbox{CH_3} $ %-$\mathbf{\mbox{CH_2}}$-$\mbox{CH_3}$ is the same as
61
62			Another advantage of using TraPPE is the constraining of all bonds to
63			be of fixed length. Typically, bond vibrations are the motions in a
64			molecular dynamic simulation. This neccesitates a small time step
65			between force evaluations be used to ensure adequate sampling of the
66			bond potential. Failure to do so will result in loss of energy
67			conservation within the microcanonical ensemble. By constraining this
68			degree of freedom, time steps larger than were previously allowable
69			are able to be used when integrating the equations of motion.
70
71			The main energy function in OOPSE is DUFF (the Dipolar Unified-atom
72			Force Field). DUFF is a collection of parameters taken from Seipmann
73			and Ichiye \emph{et
74	mmeineke	664	al.}\cite{liu96:new_model} The total energy of interaction is given by
75	mmeineke	666	Eq.~\ref{eq:totalPotential}:
76	mmeineke	698	\begin{equation}
77			V_{\text{Total}} =
78			\overbrace{V_{\theta} + V_{\phi}}^{\text{bonded}} +
79			\underbrace{V_{\text{LJ}} + V_{\text{Dp}} + %
80			V_{\text{SSD}}}_{\text{non-bonded}} \label{eq:totalPotential}
81			\end{equation}
82	mmeineke	666
83	mmeineke	698	\subsection{Bonded Interactions}
84			\label{subSec:bondedInteractions}
85	mmeineke	664
86	mmeineke	698	The bonded interactions in the DUFF functional set are limited to the
87			bend potential and the torsional potential. Bond potentials are not
88			calculated, instead all bond lengths are fixed to allow for large time
89			steps to be taken between force evaluations.
90	mmeineke	666
91	mmeineke	698	The bend functional is of the form:
92			\begin{equation}
93			V_{\theta} = \sum k_{\theta}( \theta - \theta_0 )^2 \label{eq:bendPot}
94			\end{equation}
95			$k_{\theta}$, the force constant, and $\theta_0$, the equilibrium bend
96			angle, were taken from the TraPPE forcefield of Siepmann.
97
98			The torsion functional has the form:
99			\begin{equation}
100	mmeineke	709	V_{\phi} = \sum ( k_3 \cos^3 \phi + k_2 \cos^2 \phi + k_1 \cos \phi + k_0)
101	mmeineke	698	\label{eq:torsionPot}
102			\end{equation}
103			Here, the authors decided to use a potential in terms of a power
104			expansion in $\cos \phi$ rather than the typical expansion in
105			$\phi$. This prevents the need for repeated trigonemtric
106			evaluations. Again, all $k_n$ constants were based on those in TraPPE.
107
108			\subsection{Non-Bonded Interactions}
109			\label{subSec:nonBondedInteractions}
110
111			\begin{equation}
112			V_{\text{LJ}} = \text{internal + external}
113			\end{equation}
114
115