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 925 by chrisfen, Mon Jan 12 18:43:56 2004 UTC vs.
Revision 933 by chuckv, Tue Jan 13 20:03:21 2004 UTC

# Line 9 | Line 9 | are not currently suporrted by {\sc oopse}.
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 {\sc oopse}.
12 > are not currently suported by {\sc oopse}.
13  
14   The second most basic building block of a simulation is the
15   molecule. The molecule is a way for {\sc oopse} to keep track of the
# Line 68 | Line 68 | potential. The Lennard-Jones potential mimics the attr
68   \subsection{\label{sec:LJPot}The Lennard Jones Potential}
69  
70   The most basic force field implemented in OOPSE is the Lennard-Jones
71 < potential. The Lennard-Jones potential mimics the attractive forces
72 < two charge neutral particles experience when spontaneous dipoles are
73 < induced on each other. This is the predominant intermolecular force in
74 < systems of such as noble gases and simple alkanes.
75 <
76 < The Lennard-Jones potential is given by:
71 > potential. The Lennard-Jones potential. Which mimics the Van der Waals
72 > interaction at long distances, and uses an emperical repulsion at
73 > short distances. The Lennard-Jones potential is given by:
74   \begin{equation}
75   V_{\text{LJ}}(r_{ij}) =
76          4\epsilon_{ij} \biggl[
# Line 82 | Line 79 | Where $r_ij$ is the distance between particle $i$ and
79          \biggr]
80   \label{eq:lennardJonesPot}
81   \end{equation}
82 < Where $r_ij$ is the distance between particle $i$ and $j$, $\sigma_{ij}$
83 < scales the length of the interaction, and $\epsilon_{ij}$ scales the
84 < energy well depth of the potential.
82 > Where $r_{ij}$ is the distance between particle $i$ and $j$,
83 > $\sigma_{ij}$ scales the length of the interaction, and
84 > $\epsilon_{ij}$ scales the well depth of the potential.
85  
86 < Because the potential is calculated between all pairs, the force
86 > Because this potential is calculated between all pairs, the force
87   evaluation can become computationally expensive for large systems. To
88 < keep the pair evaluation to a manegable number, OOPSE employs the use
89 < of a cut-off radius.\cite{allen87:csl} The cutoff radius is set to be
90 < $2.5\sigma_{ii}$, where $\sigma_{ii}$ is the largest self self length
88 > keep the pair evaluation to a manegable number, OOPSE employs a
89 > cut-off radius.\cite{allen87:csl} The cutoff radius is set to be
90 > $2.5\sigma_{ii}$, where $\sigma_{ii}$ is the largest Lennard-Jones length
91   parameter in the system. Truncating the calculation at
92   $r_{\text{cut}}$ introduces a discontinuity into the potential
93   energy. To offset this discontinuity, the energy value at
94   $r_{\text{cut}}$ is subtracted from the entire potential. This causes
95 < the equation to go to zero at the cut-off radius.
95 > the potential to go to zero at the cut-off radius.
96  
97   Interactions between dissimilar particles requires the generation of
98   cross term parameters for $\sigma$ and $\epsilon$. These are
# Line 114 | Line 111 | The \underline{D}ipolar \underline{U}nified-Atom
111  
112   \subsection{\label{sec:DUFF}Dipolar Unified-Atom Force Field}
113  
114 < The \underline{D}ipolar \underline{U}nified-Atom
115 < \underline{F}orce \underline{F}ield ({\sc duff}) was developed to
119 < simulate lipid bilayers. We needed a model capable of forming
114 > The Dipolar Unified-atom Force Field ({\sc duff}) was developed to
115 > simulate lipid bilayers. The systems require a model capable of forming
116   bilayers, while still being sufficiently computationally efficient to
117   allow simulations of large systems ($\approx$100's of phospholipids,
118   $\approx$1000's of waters) for long times ($\approx$10's of
119   nanoseconds).
120  
121 < With this goal in mind, we have eliminated all point charges. Charge
122 < distributions were replaced with dipoles, and charge-neutral
123 < distributions were reduced to Lennard-Jones interaction sites. This
121 > With this goal in mind, {\sc duff} has no point charges. Charge
122 > neutral distributions were replaced with dipoles, while most atoms and
123 > groups of atoms were reduced to Lennard-Jones interaction sites. This
124   simplification cuts the length scale of long range interactions from
125   $\frac{1}{r}$ to $\frac{1}{r^3}$, allowing us to avoid the
126 < computationally expensive Ewald-Sum. Instead, we can use
127 < neighbor-lists and cutoff radii for the dipolar interactions.
126 > computationally expensive Ewald sum. Instead, we can use
127 > neighbor-lists, reaction field, and cutoff radii for the dipolar
128 > interactions.
129  
130 < As an example, lipid head groups in {\sc duff} are represented as point
131 < dipole interaction sites.PC and PE Lipid head groups are typically
132 < zwitterionic in nature, with charges separated by as much as
133 < 6~$\mbox{\AA}$. By placing a dipole of 20.6~Debye at the head group
134 < center of mass, our model mimics the head group of PC.\cite{Cevc87}
135 < Additionally, a Lennard-Jones site is located at the
136 < pseudoatom's center of mass. The model is illustrated by the dark grey
137 < atom in Fig.~\ref{fig:lipidModel}.
130 > As an example, lipid head-groups in {\sc duff} are represented as
131 > point dipole interaction sites. By placing a dipole of 20.6~Debye at
132 > the head group center of mass, our model mimics the head group of
133 > phosphatidylcholine.\cite{Cevc87} Additionally, a Lennard-Jones site
134 > is located at the pseudoatom's center of mass. The model is
135 > illustrated by the dark grey atom in Fig.~\ref{fig:lipidModel}. The water model we use to complement the dipoles of the lipids is out
136 > repaarameterization of the soft sticky dipole (SSD) model of Ichiye
137 > \emph{et al.}\cite{liu96:new_model}
138  
139   \begin{figure}
140 + \epsfxsize=\linewidth
141   \epsfbox{lipidModel.eps}
142   \caption{A representation of the lipid model. $\phi$ is the torsion angle, $\theta$ %
143 < is the bend angle, $\mu$ is the dipole moment of the head group, and n is the chain length.}
143 > is the bend angle, $\mu$ is the dipole moment of the head group, and n
144 > is the chain length.}
145   \label{fig:lipidModel}
146   \end{figure}
147  
149 – The water model we use to complement the dipoles of the lipids is
150 – the soft sticky dipole (SSD) model of Ichiye \emph{et
151 – al.}\cite{liu96:new_model} This model is discussed in greater detail
152 – in Sec.~\ref{sec:SSD}. In all cases we reduce water to a single
153 – Lennard-Jones interaction site. The site also contains a dipole to
154 – mimic the partial charges on the hydrogens and the oxygen. However,
155 – what makes the SSD model unique is the inclusion of a tetrahedral
156 – short range potential to recover the hydrogen bonding of water, an
157 – important factor when modeling bilayers, as it has been shown that
158 – hydrogen bond network formation is a leading contribution to the
159 – entropic driving force towards lipid bilayer formation.\cite{Cevc87}
160 –
161 –
148   Turning to the tails of the phospholipids, we have used a set of
149   scalable parameters to model the alkyl groups with Lennard-Jones
150   sites. For this, we have used the TraPPE force field of Siepmann
# Line 184 | Line 170 | V_{\text{Total}} = \sum^{N}_{I=1} V^{I}_{\text{Interna
170   The total energy of function in {\sc duff} is given by the following:
171   \begin{equation}
172   V_{\text{Total}} = \sum^{N}_{I=1} V^{I}_{\text{Internal}}
173 <        + \sum^{N}_{I=1} \sum^{N}_{J=1} V^{IJ}_{\text{Cross}}
173 >        + \sum^{N}_{I=1} \sum_{J>I} V^{IJ}_{\text{Cross}}
174   \label{eq:totalPotential}
175   \end{equation}
176   Where $V^{I}_{\text{Internal}}$ is the internal potential of a molecule:
177   \begin{equation}
178   V^{I}_{\text{Internal}} =
179          \sum_{\theta_{ijk} \in I} V_{\text{bend}}(\theta_{ijk})
180 <        + \sum_{\phi_{ijkl} \in I} V_{\text{torsion}}(\theta_{ijkl})
180 >        + \sum_{\phi_{ijkl} \in I} V_{\text{torsion}}(\phi_{ijkl})
181          + \sum_{i \in I} \sum_{(j>i+4) \in I}
182          \biggl[ V_{\text{LJ}}(r_{ij}) +  V_{\text{dipole}}
183          (\mathbf{r}_{ij},\boldsymbol{\Omega}_{i},\boldsymbol{\Omega}_{j})
# Line 199 | Line 185 | within in the molecule. $V_{\text{torsion}}$ is the to
185   \label{eq:internalPotential}
186   \end{equation}
187   Here $V_{\text{bend}}$ is the bend potential for all 1, 3 bonded pairs
188 < within in the molecule. $V_{\text{torsion}}$ is the torsion The
189 < pairwise portions of the internal potential are excluded for pairs
190 < that ar closer than three bonds, i.e.~atom pairs farther away than a
191 < torsion are included in the pair-wise loop.
188 > within the molecule, and $V_{\text{torsion}}$ is the torsion potential
189 > for all 1, 4 bonded pairs. The pairwise portions of the internal
190 > potential are excluded for pairs that are closer than three bonds,
191 > i.e.~atom pairs farther away than a torsion are included in the
192 > pair-wise loop.
193  
207 – The cross portion of the total potential, $V^{IJ}_{\text{Cross}}$, is
208 – as follows:
209 – \begin{equation}
210 – V^{IJ}_{\text{Cross}} =
211 –        \sum_{i \in I} \sum_{j \in J}
212 –        \biggl[ V_{\text{LJ}}(r_{ij}) +  V_{\text{dipole}}
213 –        (\mathbf{r}_{ij},\boldsymbol{\Omega}_{i},\boldsymbol{\Omega}_{j})
214 –        + V_{\text{sticky}}
215 –        (\mathbf{r}_{ij},\boldsymbol{\Omega}_{i},\boldsymbol{\Omega}_{j})
216 –        \biggr]
217 – \label{eq:crossPotentail}
218 – \end{equation}
219 – Where $V_{\text{LJ}}$ is the Lennard Jones potential,
220 – $V_{\text{dipole}}$ is the dipole dipole potential, and
221 – $V_{\text{sticky}}$ is the sticky potential defined by the SSD model.
194  
195   The bend potential of a molecule is represented by the following function:
196   \begin{equation}
# Line 233 | Line 205 | V_{\text{torsion}}(\phi_{ijkl}) = c_1[1 + \cos \phi]
205   The torsion potential and parameters are also taken from TraPPE. It is
206   of the form:
207   \begin{equation}
208 < V_{\text{torsion}}(\phi_{ijkl}) = c_1[1 + \cos \phi]
208 > V_{\text{torsion}}(\phi) = c_1[1 + \cos \phi]
209          + c_2[1 + \cos(2\phi)]
210          + c_3[1 + \cos(3\phi)]
211   \label{eq:origTorsionPot}
212   \end{equation}
213   Here $\phi_{ijkl}$ is the angle defined by four bonded neighbors $i$,
214 < $j$, $k$, and $l$ (again, see Fig.~\ref{fig:lipidModel}).  However,
215 < for computaional efficency, the torsion potentail has been recast
216 < after the method of CHARMM\cite{charmm1983} whereby the angle series
217 < is converted to a power series of the form:
214 > $j$, $k$, and $l$ (again, see Fig.~\ref{fig:lipidModel}). For
215 > computaional efficency, the torsion potential has been recast after
216 > the method of CHARMM\cite{charmm1983} whereby the angle series is
217 > converted to a power series of the form:
218   \begin{equation}
219   V_{\text{torsion}}(\phi_{ijkl}) =  
220          k_3 \cos^3 \phi + k_2 \cos^2 \phi + k_1 \cos \phi + k_0
# Line 259 | Line 231 | evaluations are avoided during the calculation of the
231   evaluations are avoided during the calculation of the potential.
232  
233  
234 + The cross portion of the total potential, $V^{IJ}_{\text{Cross}}$, is
235 + as follows:
236 + \begin{equation}
237 + V^{IJ}_{\text{Cross}} =
238 +        \sum_{i \in I} \sum_{j \in J}
239 +        \biggl[ V_{\text{LJ}}(r_{ij}) +  V_{\text{dipole}}
240 +        (\mathbf{r}_{ij},\boldsymbol{\Omega}_{i},\boldsymbol{\Omega}_{j})
241 +        + V_{\text{sticky}}
242 +        (\mathbf{r}_{ij},\boldsymbol{\Omega}_{i},\boldsymbol{\Omega}_{j})
243 +        \biggr]
244 + \label{eq:crossPotentail}
245 + \end{equation}
246 + Where $V_{\text{LJ}}$ is the Lennard Jones potential,
247 + $V_{\text{dipole}}$ is the dipole dipole potential, and
248 + $V_{\text{sticky}}$ is the sticky potential defined by the SSD
249 + model. Note that not all atom types include all interactions.
250  
251   The dipole-dipole potential has the following form:
252   \begin{equation}
253   V_{\text{dipole}}(\mathbf{r}_{ij},\boldsymbol{\Omega}_{i},
254 <        \boldsymbol{\Omega}_{j}) = \frac{1}{4\pi\epsilon_{0}} \biggl[
255 <        \frac{\boldsymbol{\mu}_{i} \cdot \boldsymbol{\mu}_{j}}{r^{3}_{ij}}
254 >        \boldsymbol{\Omega}_{j}) = \frac{|\mu_i||\mu_j|}{4\pi\epsilon_{0}r_{ij}^{3}} \biggl[
255 >        \boldsymbol{\hat{u}}_{i} \cdot \boldsymbol{\hat{u}}_{j}
256          -
257 <        \frac{3(\boldsymbol{\mu}_i \cdot \mathbf{r}_{ij}) %
258 <                (\boldsymbol{\mu}_j \cdot \mathbf{r}_{ij}) }
259 <                {r^{5}_{ij}} \biggr]
257 >        \frac{3(\boldsymbol{\hat{u}}_i \cdot \mathbf{r}_{ij}) %
258 >                (\boldsymbol{\hat{u}}_j \cdot \mathbf{r}_{ij}) }
259 >                {r^{2}_{ij}} \biggr]
260   \label{eq:dipolePot}
261   \end{equation}
262   Here $\mathbf{r}_{ij}$ is the vector starting at atom $i$ pointing
263   towards $j$, and $\boldsymbol{\Omega}_i$ and $\boldsymbol{\Omega}_j$
264 < are the Euler angles of atom $i$ and $j$
265 < respectively. $\boldsymbol{\mu}_i$ is the dipole vector of atom
266 < $i$ it takes its orientation from $\boldsymbol{\Omega}_i$.
264 > are the Euler angles of atom $i$ and $j$ respectively. $|\mu_i|$ is
265 > the magnitude of the dipole moment of atom $i$ and
266 > $\boldsymbol{\hat{u}}_i$ is the standard unit orientation vector of
267 > $\boldsymbol{\Omega}_i$.
268  
269  
270   \subsection{\label{sec:SSD}The {\sc DUFF} Water Models: SSD/E and SSD/RF}
# Line 367 | Line 356 | models can be found in reference \ref{Gezelter04}.
356  
357   !!!Place a {\sc BASS} scheme showing SSD parameters around here!!!
358  
359 < \subsection{\label{sec:eam}Embedded Atom Model}
359 > \subsection{\label{sec:eam}Embedded Atom Method}
360  
361 < Several molecular dynamics codes\cite{dynamo86} exist which have the
361 > Several other molecular dynamics packages\cite{dynamo86} exist which have the
362   capacity to simulate metallic systems, including some that have
363   parallel computational abilities\cite{plimpton93}. Potentials that
364   describe bonding transition metal
365   systems\cite{Finnis84,Ercolessi88,Chen90,Qi99,Ercolessi02} have a
366 < attractive interaction which models the stabilization of ``Embedding''
367 < a positively charged metal ion in an electron density created by the
366 > attractive interaction which models  ``Embedding''
367 > a positively charged metal ion in the electron density due to the
368   free valance ``sea'' of electrons created by the surrounding atoms in
369   the system. A mostly repulsive pairwise part of the potential
370   describes the interaction of the positively charged metal core ions
371   with one another. A particular potential description called the
372 < Embedded Atom Method\cite{Daw84,FBD86,johnson89,Lu97}(EAM) that has
372 > Embedded Atom Method\cite{Daw84,FBD86,johnson89,Lu97}({\sc eam}) that has
373   particularly wide adoption has been selected for inclusion in OOPSE. A
374 < good review of EAM and other metallic potential formulations was done
374 > good review of {\sc eam} and other metallic potential formulations was done
375   by Voter.\cite{voter}
376  
377   The {\sc eam} potential has the form:
# Line 390 | Line 379 | V & = & \sum_{i} F_{i}\left[\rho_{i}\right] + \sum_{i}
379   V & = & \sum_{i} F_{i}\left[\rho_{i}\right] + \sum_{i} \sum_{j \neq i}
380   \phi_{ij}({\bf r}_{ij})  \\
381   \rho_{i}  & = & \sum_{j \neq i} f_{j}({\bf r}_{ij})
382 < \end{eqnarray}
382 > \end{eqnarray}S
383  
384 < where $\phi_{ij}$ is a primarily repulsive pairwise interaction
385 < between atoms $i$ and $j$.In the origional formulation of
397 < EAM\cite{Daw84}, $\phi_{ij}$ was an entirely repulsive term, however
398 < in later refinements to EAM have shown that nonuniqueness between $F$
399 < and $\phi$ allow for more general forms for $\phi$.\cite{Daw89} The
400 < embedding function $F_{i}$ is the energy required to embedded an
401 < positively-charged core ion $i$ into a linear supeposition of
384 > where $F_{i} $ is the embedding function that equates the energy required to embed a
385 > positively-charged core ion $i$ into a linear superposition of
386   sperically averaged atomic electron densities given by
387 < $\rho_{i}$. There is a cutoff distance, $r_{cut}$, which limits the
387 > $\rho_{i}$.  $\phi_{ij}$ is a primarily repulsive pairwise interaction
388 > between atoms $i$ and $j$. In the original formulation of
389 > {\sc eam} cite{Daw84}, $\phi_{ij}$ was an entirely repulsive term, however
390 > in later refinements to EAM have shown that non-uniqueness between $F$
391 > and $\phi$ allow for more general forms for $\phi$.\cite{Daw89}
392 > There is a cutoff distance, $r_{cut}$, which limits the
393   summations in the {\sc eam} equation to the few dozen atoms
394   surrounding atom $i$ for both the density $\rho$ and pairwise $\phi$
395 < interactions.
395 > interactions. Foiles et al. fit EAM potentials for fcc metals Cu, Ag, Au, Ni, Pd, Pt and alloys of these metals\cite{FDB86}. These potential fits are in the DYNAMO 86 format and are included with {\sc oopse}.
396  
397 +
398   \subsection{\label{Sec:pbc}Periodic Boundary Conditions}
399  
400 < \textit{Periodic boundary conditions} are widely used to simulate truly
401 < macroscopic systems with a relatively small number of particles. Simulation
402 < box is replicated throughout space to form an infinite lattice. During the
403 < simulation, when a particle moves in the primary cell, its periodic image
404 < particles in other boxes move in exactly the same direction with exactly the
405 < same orientation.Thus, as a particle leaves the primary cell, one of its
406 < images will enter through the opposite face.If the simulation box is large
407 < enough to avoid "feeling" the symmetric of the periodic lattice,the surface
408 < effect could be ignored. Cubic and parallelepiped are the available periodic
409 < cells. \bigskip In OOPSE, we use the matrix instead of the vector to describe
410 < the property of the simulation box. Therefore, not only the size of the
411 < simulation box could be changed during the simulation, but also the shape of
412 < it. The transformation from box space vector $\overrightarrow{s}$ to its
413 < corresponding real space vector $\overrightarrow{r}$ is defined by
414 < \begin{equation}
415 < \overrightarrow{r}=H\overrightarrow{s}%
416 < \end{equation}
417 <
418 <
419 < where $H=(h_{x},h_{y},h_{z})$ is a transformation matrix made up of the three
420 < box axis vectors. $h_{x},h_{y}$ and $h_{z}$ represent the sides of the
421 < simulation box respectively.
422 <
423 < To find the minimum image, we need to convert the real vector to its
424 < corresponding vector in box space first, \bigskip%
425 < \begin{equation}
426 < \overrightarrow{s}=H^{-1}\overrightarrow{r}%
427 < \end{equation}
428 < And then, each element of $\overrightarrow{s}$ is casted to lie between -0.5
429 < to 0.5,
430 < \begin{equation}
431 < s_{i}^{\prime}=s_{i}-round(s_{i})
432 < \end{equation}
433 < where%
434 <
435 < \begin{equation}
436 < round(x)=\lfloor{x+0.5}\rfloor\text{ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ }if\text{
437 < }x\geqslant0
438 < \end{equation}
439 < %
440 <
441 < \begin{equation}
442 < round(x)=\lceil{x-0.5}\rceil\text{ \ \ \ \ \ \ \ \ \ \ \ \ \ \ }if\text{ }x<0
443 < \end{equation}
444 <
445 <
446 < For example, $round(3.6)=4$,$round(3.1)=3$, $round(-3.6)=-4$, $round(-3.1)=-3$.
447 <
448 < Finally, we could get the minimum image by transforming back to real space,%
449 <
450 < \begin{equation}
451 < \overrightarrow{r^{\prime}}=H^{-1}\overrightarrow{s^{\prime}}%
452 < \end{equation}
400 > \textit{Periodic boundary conditions} are widely used to simulate truly
401 > macroscopic systems with a relatively small number of particles. The
402 > simulation box is replicated throughout space to form an infinite
403 > lattice.  During the simulation, when a particle moves in the primary
404 > cell, its image in other boxes move in exactly the same direction with
405 > exactly the same orientation.Thus, as a particle leaves the primary
406 > cell, one of its images will enter through the opposite face.If the
407 > simulation box is large enough to avoid "feeling" the symmetries of
408 > the periodic lattice, surface effects can be ignored. Cubic,
409 > orthorhombic and parallelepiped are the available periodic cells In
410 > OOPSE. We use a matrix to describe the property of the simulation
411 > box. Therefore, both the size and shape of the simulation box can be
412 > changed during the simulation. The transformation from box space
413 > vector $\mathbf{s}$ to its corresponding real space vector
414 > $\mathbf{r}$ is defined by
415 > \begin{equation}
416 > \mathbf{r}=\underline{\underline{H}}\cdot\mathbf{s}%
417 > \end{equation}
418 >
419 >
420 > where $H=(h_{x},h_{y},h_{z})$ is a transformation matrix made up of
421 > the three box axis vectors. $h_{x},h_{y}$ and $h_{z}$ represent the
422 > three sides of the simulation box respectively.
423 >
424 > To find the minimum image, we convert the real vector to its
425 > corresponding vector in box space first, \bigskip%
426 > \begin{equation}
427 > \mathbf{s}=\underline{\underline{H}}^{-1}\cdot\mathbf{r}%
428 > \end{equation}
429 > And then, each element of $\mathbf{s}$ is wrapped to lie between -0.5 to 0.5,
430 > \begin{equation}
431 > s_{i}^{\prime}=s_{i}-round(s_{i})
432 > \end{equation}
433 > where
434 >
435 > %
436 >
437 > \begin{equation}
438 > round(x)=\left\{
439 > \begin{array}[c]{c}%
440 > \lfloor{x+0.5}\rfloor & \text{if \ }x\geqslant0\\
441 > \lceil{x-0.5}\rceil & \text{otherwise}%
442 > \end{array}
443 > \right.
444 > \end{equation}
445 >
446 >
447 > For example, $round(3.6)=4$,$round(3.1)=3$, $round(-3.6)=-4$,
448 > $round(-3.1)=-3$.
449 >
450 > Finally, we obtain the minimum image coordinates by transforming back
451 > to real space,%
452 >
453 > \begin{equation}
454 > \mathbf{r}^{\prime}=\underline{\underline{H}}^{-1}\cdot\mathbf{s}^{\prime}%
455 > \end{equation}
456 >

Diff Legend

– Removed lines
+ Added lines
< Changed lines
> Changed lines