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 929 by chrisfen, Mon Jan 12 18:43:56 2004 UTC vs.
Revision 930 by mmeineke, Tue Jan 13 19:24:07 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}
148
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}
147  
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
362   capacity to simulate metallic systems, including some that have
# Line 407 | Line 396 | interactions.
396  
397   \subsection{\label{Sec:pbc}Periodic Boundary Conditions}
398  
399 < \textit{Periodic boundary conditions} are widely used to simulate truly
400 < macroscopic systems with a relatively small number of particles. Simulation
401 < box is replicated throughout space to form an infinite lattice. During the
402 < simulation, when a particle moves in the primary cell, its periodic image
403 < particles in other boxes move in exactly the same direction with exactly the
404 < same orientation.Thus, as a particle leaves the primary cell, one of its
405 < images will enter through the opposite face.If the simulation box is large
406 < enough to avoid "feeling" the symmetric of the periodic lattice,the surface
407 < effect could be ignored. Cubic and parallelepiped are the available periodic
408 < cells. \bigskip In OOPSE, we use the matrix instead of the vector to describe
409 < the property of the simulation box. Therefore, not only the size of the
410 < simulation box could be changed during the simulation, but also the shape of
411 < it. The transformation from box space vector $\overrightarrow{s}$ to its
412 < corresponding real space vector $\overrightarrow{r}$ is defined by
413 < \begin{equation}
414 < \overrightarrow{r}=H\overrightarrow{s}%
415 < \end{equation}
416 <
417 <
418 < where $H=(h_{x},h_{y},h_{z})$ is a transformation matrix made up of the three
419 < box axis vectors. $h_{x},h_{y}$ and $h_{z}$ represent the sides of the
420 < simulation box respectively.
421 <
422 < To find the minimum image, we need to convert the real vector to its
423 < corresponding vector in box space first, \bigskip%
424 < \begin{equation}
425 < \overrightarrow{s}=H^{-1}\overrightarrow{r}%
426 < \end{equation}
427 < And then, each element of $\overrightarrow{s}$ is casted to lie between -0.5
428 < to 0.5,
429 < \begin{equation}
430 < s_{i}^{\prime}=s_{i}-round(s_{i})
431 < \end{equation}
432 < where%
433 <
434 < \begin{equation}
435 < round(x)=\lfloor{x+0.5}\rfloor\text{ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ }if\text{
436 < }x\geqslant0
437 < \end{equation}
438 < %
439 <
440 < \begin{equation}
441 < round(x)=\lceil{x-0.5}\rceil\text{ \ \ \ \ \ \ \ \ \ \ \ \ \ \ }if\text{ }x<0
442 < \end{equation}
443 <
444 <
445 < For example, $round(3.6)=4$,$round(3.1)=3$, $round(-3.6)=-4$, $round(-3.1)=-3$.
446 <
447 < Finally, we could get the minimum image by transforming back to real space,%
448 <
449 < \begin{equation}
450 < \overrightarrow{r^{\prime}}=H^{-1}\overrightarrow{s^{\prime}}%
451 < \end{equation}
399 > \textit{Periodic boundary conditions} are widely used to simulate truly
400 > macroscopic systems with a relatively small number of particles. The
401 > simulation box is replicated throughout space to form an infinite
402 > lattice.  During the simulation, when a particle moves in the primary
403 > cell, its image in other boxes move in exactly the same direction with
404 > exactly the same orientation.Thus, as a particle leaves the primary
405 > cell, one of its images will enter through the opposite face.If the
406 > simulation box is large enough to avoid "feeling" the symmetries of
407 > the periodic lattice, surface effects can be ignored. Cubic,
408 > orthorhombic and parallelepiped are the available periodic cells In
409 > OOPSE. We use a matrix to describe the property of the simulation
410 > box. Therefore, both the size and shape of the simulation box can be
411 > changed during the simulation. The transformation from box space
412 > vector $\mathbf{s}$ to its corresponding real space vector
413 > $\mathbf{r}$ is defined by
414 > \begin{equation}
415 > \mathbf{r}=\underline{\underline{H}}\cdot\mathbf{s}%
416 > \end{equation}
417 >
418 >
419 > where $H=(h_{x},h_{y},h_{z})$ is a transformation matrix made up of
420 > the three box axis vectors. $h_{x},h_{y}$ and $h_{z}$ represent the
421 > three sides of the simulation box respectively.
422 >
423 > To find the minimum image, we convert the real vector to its
424 > corresponding vector in box space first, \bigskip%
425 > \begin{equation}
426 > \mathbf{s}=\underline{\underline{H}}^{-1}\cdot\mathbf{r}%
427 > \end{equation}
428 > And then, each element of $\mathbf{s}$ is wrapped to lie between -0.5 to 0.5,
429 > \begin{equation}
430 > s_{i}^{\prime}=s_{i}-round(s_{i})
431 > \end{equation}
432 > where
433 >
434 > %
435 >
436 > \begin{equation}
437 > round(x)=\left\{
438 > \begin{array}[c]{c}%
439 > \lfloor{x+0.5}\rfloor & \text{if \ }x\geqslant0\\
440 > \lceil{x-0.5}\rceil & \text{otherwise}%
441 > \end{array}
442 > \right.
443 > \end{equation}
444 >
445 >
446 > For example, $round(3.6)=4$,$round(3.1)=3$, $round(-3.6)=-4$,
447 > $round(-3.1)=-3$.
448 >
449 > Finally, we obtain the minimum image coordinates by transforming back
450 > to real space,%
451 >
452 > \begin{equation}
453 > \mathbf{r}^{\prime}=\underline{\underline{H}}^{-1}\cdot\mathbf{s}^{\prime}%
454 > \end{equation}
455 >

Diff Legend

Removed lines
+ Added lines
< Changed lines
> Changed lines