ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/trunk/oopsePaper/EmpericalEnergy.tex
(Generate patch)

Comparing trunk/oopsePaper/EmpericalEnergy.tex (file contents):
Revision 806 by mmeineke, Fri Oct 17 05:07:49 2003 UTC vs.
Revision 899 by mmeineke, Tue Jan 6 18:53:58 2004 UTC

# Line 1 | Line 1
1  
2 < \section{The Emperical Energy Functions}
2 > \section{\label{sec:empericalEnergy}The Emperical Energy Functions}
3  
4 < \subsection{Atoms and Molecules}
4 > \subsection{\label{sec:atomsMolecules}Atoms and Molecules}
5  
6 < The basic unit of an OOPSE simulation is the atom. The parameters
6 > The basic unit of an {\sc oopse} simulation is the atom. The parameters
7   describing the atom are generalized to make the atom as flexible a
8   representation as possible. They may represent specific atoms of an
9   element, or be used for collections of atoms such as a methyl
10   group. The atoms are also capable of having a directional component
11   associated with them, often in the form of a dipole. Charges on atoms
12 < are not currently suporrted by OOPSE.
12 > are not currently suporrted by {\sc oopse}.
13  
14   The second most basic building block of a simulation is the
15 < molecule. The molecule is a way for OOPSE to keep track of the atoms
15 > molecule. The molecule is a way for {\sc oopse} to keep track of the atoms
16   in a simulation in logical manner. This particular unit will store the
17   identities of all the atoms associated with itself and is responsible
18   for the evaluation of its own bonded interaction (i.e.~bonds, bends,
19   and torsions).
20 +
21 + \subsection{\label{sec:DUFF}Dipolar Unified-Atom Force Field}
22 +
23 + The \underline{D}ipolar \underline{U}nified-Atom
24 + \underline{F}orce \underline{F}ield ({\sc duff}) was developed to
25 + simulate lipid bilayers. We needed a model capable of forming
26 + bilayers, while still being sufficiently computationally efficient to
27 + allow simulations of large systems ($\approx$100's of phospholipids,
28 + $\approx$1000's of waters) for long times ($\approx$10's of
29 + nanoseconds).
30 +
31 + With this goal in mind, we have eliminated all point charges. Charge
32 + distributions were replaced with dipoles, and charge-neutral
33 + distributions were reduced to Lennard-Jones interaction sites. This
34 + simplification cuts the length scale of long range interactions from
35 + $\frac{1}{r}$ to $\frac{1}{r^3}$, allowing us to avoid the
36 + computationally expensive Ewald-Sum. Instead, we can use
37 + neighbor-lists and cutoff radii for the dipolar interactions.
38 +
39 + As an example, lipid head groups in {\sc duff} are represented as point
40 + dipole interaction sites.PC and PE Lipid head groups are typically
41 + zwitterionic in nature, with charges separated by as much as
42 + 6~$\mbox{\AA}$. By placing a dipole of 20.6~Debye at the head group
43 + center of mass, our model mimics the head group of PC.\cite{Cevc87}
44 + Additionally, a Lennard-Jones site is located at the
45 + pseudoatom's center of mass. The model is illustrated by the dark grey
46 + atom in Fig.~\ref{fig:lipidModel}.
47 +
48 + \begin{figure}
49 + \epsfxsize=6in
50 + \epsfbox{lipidModel.epsi}
51 + \caption{A representation of the lipid model. $\phi$ is the torsion angle, $\theta$ %
52 + is the bend angle, $\mu$ is the dipole moment of the head group, and n is the chain length.}
53 + \label{fig:lipidModel}
54 + \end{figure}
55 +
56 + The water model we use to complement the dipoles of the lipids is
57 + the soft sticky dipole (SSD) model of Ichiye \emph{et
58 + al.}\cite{liu96:new_model} This model is discussed in greater detail
59 + in Sec.~\ref{sec:SSD}. In all cases we reduce water to a single
60 + Lennard-Jones interaction site. The site also contains a dipole to
61 + mimic the partial charges on the hydrogens and the oxygen. However,
62 + what makes the SSD model unique is the inclusion of a tetrahedral
63 + short range potential to recover the hydrogen bonding of water, an
64 + important factor when modeling bilayers, as it has been shown that
65 + hydrogen bond network formation is a leading contribution to the
66 + entropic driving force towards lipid bilayer formation.\cite{Cevc87}
67 +
68 +
69 + Turning to the tails of the phospholipids, we have used a set of
70 + scalable parameters to model the alkyl groups with Lennard-Jones
71 + sites. For this, we have used the TraPPE force field of Siepmann
72 + \emph{et al}.\cite{Siepmann1998} TraPPE is a unified-atom
73 + representation of n-alkanes, which is parametrized against phase
74 + equilibria using Gibbs Monte Carlo simulation
75 + techniques.\cite{Siepmann1998} One of the advantages of TraPPE is that
76 + it generalizes the types of atoms in an alkyl chain to keep the number
77 + of pseudoatoms to a minimum; the parameters for an atom such as
78 + $\text{CH}_2$ do not change depending on what species are bonded to
79 + it.
80 +
81 + TraPPE also constrains of all bonds to be of fixed length. Typically,
82 + bond vibrations are the fastest motions in a molecular dynamic
83 + simulation. Small time steps between force evaluations must be used to
84 + ensure adequate sampling of the bond potential conservation of
85 + energy. By constraining the bond lengths, larger time steps may be
86 + used when integrating the equations of motion.
87 +
88 +
89 + \subsubsection{\label{subSec:energyFunctions}{\sc duff} Energy Functions}
90 +
91 + The total energy of function in {\sc duff} is given by the following:
92 + \begin{equation}
93 + V_{\text{Total}} = \sum^{N}_{I=1} V^{I}_{\text{Internal}}
94 +        + \sum^{N}_{I=1} \sum^{N}_{J=1} V^{IJ}_{\text{Cross}}
95 + \label{eq:totalPotential}
96 + \end{equation}
97 + Where $V^{I}_{\text{Internal}}$ is the internal potential of a molecule:
98 + \begin{equation}
99 + V^{I}_{\text{Internal}} =
100 +        \sum_{\theta_{ijk} \in I} V_{\text{bend}}(\theta_{ijk})
101 +        + \sum_{\phi_{ijkl} \in I} V_{\text{torsion}}(\theta_{ijkl})
102 +        + \sum_{i \in I} \sum_{(j>i+4) \in I}
103 +        \biggl[ V_{\text{LJ}}(r_{ij}) +  V_{\text{dipole}}
104 +        (\mathbf{r}_{ij},\boldsymbol{\Omega}_{i},\boldsymbol{\Omega}_{j})
105 +        \biggr]
106 + \label{eq:internalPotential}
107 + \end{equation}
108 + Here $V_{\text{bend}}$ is the bend potential for all 1, 3 bonded pairs
109 + within in the molecule. $V_{\text{torsion}}$ is the torsion The
110 + pairwise portions of the internal potential are excluded for pairs
111 + that ar closer than three bonds, i.e.~atom pairs farther away than a
112 + torsion are included in the pair-wise loop.
113 +
114 + The cross portion of the total potential, $V^{IJ}_{\text{Cross}}$, is
115 + as follows:
116 + \begin{equation}
117 + V^{IJ}_{\text{Cross}} =
118 +        \sum_{i \in I} \sum_{j \in J}
119 +        \biggl[ V_{\text{LJ}}(r_{ij}) +  V_{\text{dipole}}
120 +        (\mathbf{r}_{ij},\boldsymbol{\Omega}_{i},\boldsymbol{\Omega}_{j})
121 +        + V_{\text{sticky}}
122 +        (\mathbf{r}_{ij},\boldsymbol{\Omega}_{i},\boldsymbol{\Omega}_{j})
123 +        \biggr]
124 + \label{eq:crossPotentail}
125 + \end{equation}
126 + Where $V_{\text{LJ}}$ is the Lennard Jones potential,
127 + $V_{\text{dipole}}$ is the dipole dipole potential, and
128 + $V_{\text{sticky}}$ is the sticky potential defined by the SSD model.
129 +
130 + The bend potential of a molecule is represented by the following function:
131 + \begin{equation}
132 + V_{\theta_{ijk}} = k_{\theta}( \theta_{ijk} - \theta_0 )^2 \label{eq:bendPot}
133 + \end{equation}
134 + Where $\theta_{ijk}$ is the angle defined by atoms $i$, $j$, and $k$
135 + (see Fig.~\ref{fig:lipidModel}), and $\theta_0$ is the equilibrium
136 + bond angle. $k_{\theta}$ is the force constant which determines the
137 + strength of the harmonic bend. The parameters for $k_{\theta}$ and
138 + $\theta_0$ are based off of those in TraPPE.\cite{Siepmann1998}
139 +
140 + The torsion potential and parameters are also taken from TraPPE. It is
141 + of the form:
142 + \begin{equation}
143 + V_{\text{torsion}}(\phi_{ijkl}) = c_1[1 + \cos \phi]
144 +        + c_2[1 + \cos(2\phi)]
145 +        + c_3[1 + \cos(3\phi)]
146 + \label{eq:origTorsionPot}
147 + \end{equation}
148 + Here $\phi_{ijkl}$ is the angle defined by four bonded neighbors $i$,
149 + $j$, $k$, and $l$ (again, see Fig.~\ref{fig:lipidModel}).  However,
150 + for computaional efficency, the torsion potentail has been recast
151 + after the method of CHARMM\cite{charmm1983} whereby the angle series
152 + is converted to a power series of the form:
153 + \begin{equation}
154 + V_{\text{torsion}}(\phi_{ijkl}) =  
155 +        k_3 \cos^3 \phi + k_2 \cos^2 \phi + k_1 \cos \phi + k_0
156 + \label{eq:torsionPot}
157 + \end{equation}
158 + Where:
159 + \begin{align*}
160 + k_0 &= c_1 + c_3 \\
161 + k_1 &= c_1 - 3c_3 \\
162 + k_2 &= 2 c_2 \\
163 + k_3 &= 4c_3
164 + \end{align*}
165 + By recasting the equation to a power series, repeated trigonometric
166 + evaluations are avoided during the calculation of the potential.
167 +
168 + The Lennard-Jones potential is given by:
169 + \begin{equation}
170 + V_{\text{LJ}}(r_{ij}) =
171 +        4\epsilon_{ij} \biggl[
172 +        \biggl(\frac{\sigma_{ij}}{r_{ij}}\biggr)^{12}
173 +        - \biggl(\frac{\sigma_{ij}}{r_{ij}}\biggr)^{6}
174 +        \biggr]
175 + \label{eq:lennardJonesPot}
176 + \end{equation}
177 + Where $r_ij$ is the distance between atoms $i$ and $j$, $\sigma_{ij}$
178 + scales the length of the interaction, and $\epsilon_{ij}$ scales the
179 + energy of the potential.
180 +
181 + The dipole-dipole potential has the following form:
182 + \begin{equation}
183 + V_{\text{dipole}}(\mathbf{r}_{ij},\boldsymbol{\Omega}_{i},
184 +        \boldsymbol{\Omega}_{j}) = \frac{1}{4\pi\epsilon_{0}} \biggl[
185 +        \frac{\boldsymbol{\mu}_{i} \cdot \boldsymbol{\mu}_{j}}{r^{3}_{ij}}
186 +        -
187 +        \frac{3(\boldsymbol{\mu}_i \cdot \mathbf{r}_{ij}) %
188 +                (\boldsymbol{\mu}_j \cdot \mathbf{r}_{ij}) }
189 +                {r^{5}_{ij}} \biggr]
190 + \label{eq:dipolePot}
191 + \end{equation}
192 + Here $\mathbf{r}_{ij}$ is the vector starting at atom $i$ pointing
193 + towards $j$, and $\boldsymbol{\Omega}_i$ and $\boldsymbol{\Omega}_j$
194 + are the Euler angles of atom $i$ and $j$
195 + respectively. $\boldsymbol{\mu}_i$ is the dipole vector of atom
196 + $i$ it takes its orientation from $\boldsymbol{\Omega}_i$.
197 +
198 +
199 + \subsection{\label{sec:SSD}Water Model: SSD and Derivatives}
200 +
201 + In the interest of computational efficiency, the native solvent used
202 + in {\sc oopse} is the Soft Sticky Dipole (SSD) water model. SSD was
203 + developed by Ichiye \emph{et al.} as a modified form of the
204 + hard-sphere water model proposed by Bratko, Blum, and
205 + Luzar.\cite{Bratko85,Bratko95} It consists of a single point dipole
206 + with a Lennard-Jones core and a sticky potential that directs the
207 + particles to assume the proper hydrogen bond orientation in the first
208 + solvation shell. Thus, the interaction between two SSD water molecules
209 + \emph{i} and \emph{j} is given by the potential
210 + \begin{equation}
211 + u_{ij} = u_{ij}^{LJ} (r_{ij})\ + u_{ij}^{dp}
212 + (\mathbf{r}_{ij},\boldsymbol{\Omega}_i,\boldsymbol{\Omega}_j)\ +
213 + u_{ij}^{sp}
214 + (\mathbf{r}_{ij},\boldsymbol{\Omega}_i,\boldsymbol{\Omega}_j),
215 + \end{equation}
216 + where the $\mathbf{r}_{ij}$ is the position vector between molecules
217 + \emph{i} and \emph{j} with magnitude equal to the distance $r_ij$, and
218 + $\boldsymbol{\Omega}_i$ and $\boldsymbol{\Omega}_j$ represent the
219 + orientations of the respective molecules. The Lennard-Jones, dipole,
220 + and sticky parts of the potential are giving by the following
221 + equations,
222 + \begin{equation}
223 + u_{ij}^{LJ}(r_{ij}) = 4\epsilon \left[\left(\frac{\sigma}{r_{ij}}\right)^{12}-\left(\frac{\sigma}{r_{ij}}\right)^{6}\right],
224 + \end{equation}
225 + \begin{equation}
226 + u_{ij}^{dp} = \frac{\boldsymbol{\mu}_i\cdot\boldsymbol{\mu}_j}{r_{ij}^3}-\frac{3(\boldsymbol{\mu}_i\cdot\mathbf{r}_{ij})(\boldsymbol{\mu}_j\cdot\mathbf{r}_{ij})}{r_{ij}^5}\ ,
227 + \end{equation}
228 + \begin{equation}
229 + \begin{split}
230 + u_{ij}^{sp}
231 + (\mathbf{r}_{ij},\boldsymbol{\Omega}_i,\boldsymbol{\Omega}_j)
232 + &=
233 + \frac{\nu_0}{2}[s(r_{ij})w(\mathbf{r}_{ij},\boldsymbol{\Omega}_i,\boldsymbol{\Omega}_j)\\
234 + & \quad \ +
235 + s^\prime(r_{ij})w^\prime(\mathbf{r}_{ij},\boldsymbol{\Omega}_i,\boldsymbol{\Omega}_j)]\ ,
236 + \end{split}
237 + \end{equation}
238 + where $\boldsymbol{\mu}_i$ and $\boldsymbol{\mu}_j$ are the dipole
239 + unit vectors of particles \emph{i} and \emph{j} with magnitude 2.35 D,
240 + $\nu_0$ scales the strength of the overall sticky potential, $s$ and
241 + $s^\prime$ are cubic switching functions. The $w$ and $w^\prime$
242 + functions take the following forms,
243 + \begin{equation}
244 + w(\mathbf{r}_{ij},\boldsymbol{\Omega}_i,\boldsymbol{\Omega}_j)=\sin\theta_{ij}\sin2\theta_{ij}\cos2\phi_{ij},
245 + \end{equation}
246 + \begin{equation}
247 + w^\prime(\mathbf{r}_{ij},\boldsymbol{\Omega}_i,\boldsymbol{\Omega}_j) = (\cos\theta_{ij}-0.6)^2(\cos\theta_{ij}+0.8)^2-w^0,
248 + \end{equation}
249 + where $w^0 = 0.07715$. The $w$ function is the tetrahedral attractive
250 + term that promotes hydrogen bonding orientations within the first
251 + solvation shell, and $w^\prime$ is a dipolar repulsion term that
252 + repels unrealistic dipolar arrangements within the first solvation
253 + shell. A more detailed description of the functional parts and
254 + variables in this potential can be found in other
255 + articles.\cite{liu96:new_model,chandra99:ssd_md}
256 +
257 + Since SSD is a one-site point dipole model, the force calculations are
258 + simplified significantly from a computational standpoint, in that the
259 + number of long-range interactions is dramatically reduced. In the
260 + original Monte Carlo simulations using this model, Ichiye \emph{et
261 + al.} reported a calculation speed up of up to an order of magnitude
262 + over other comparable models while maintaining the structural behavior
263 + of water.\cite{liu96:new_model} In the original molecular dynamics studies of
264 + SSD, it was shown that it actually improves upon the prediction of
265 + water's dynamical properties over TIP3P and SPC/E.\cite{chandra99:ssd_md}
266 +
267 + Recent constant pressure simulations revealed issues in the original
268 + SSD model that led to lower than expected densities at all target
269 + pressures.\cite{Ichiye03,Gezelter04} Reparameterizations of the
270 + original SSD have resulted in improved density behavior, as well as
271 + alterations in the water structure and transport behavior. {\sc oopse} is
272 + easily modified to impliment these new potential parameter sets for
273 + the derivative water models: SSD1, SSD/E, and SSD/RF. All of the
274 + variable parameters are listed in the accompanying BASS file, and
275 + these parameters simply need to be changed to the updated values.
276 +
277 +
278 + \subsection{\label{sec:eam}Embedded Atom Model}
279 +
280 + here there be Monsters

Diff Legend

Removed lines
+ Added lines
< Changed lines
> Changed lines