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 968 by tim, Tue Jan 20 16:49:22 2004 UTC vs.
Revision 973 by mmeineke, Wed Jan 21 19:00:23 2004 UTC

# Line 11 | Line 11 | dipoles). Charges on atoms are not currently supported
11   directional components associated with them (\emph{e.g.}~permanent
12   dipoles). Charges on atoms are not currently supported by {\sc oopse}.
13  
14 < \begin{lstlisting}[caption={[Specifier for molecules and atoms] An example specifying the simple Ar molecule},label=sch:AtmMole]
14 > \begin{lstlisting}[caption={[Specifier for molecules and atoms] A sample specification of the simple Ar molecule},label=sch:AtmMole]
15   molecule{
16    name = "Ar";
17    nAtoms = 1;
18    atom[0]{
19 <     type="Ar";
20 <     position( 0.0, 0.0, 0.0 );
19 >    type="Ar";
20 >    position( 0.0, 0.0, 0.0 );
21    }
22   }
23   \end{lstlisting}
24  
25 < The second most basic building block of a simulation is the
26 < molecule. The molecule is a way for {\sc oopse} to keep track of the
27 < atoms in a simulation in logical manner. This particular unit will
28 < store the identities of all the atoms associated with itself and is
29 < responsible for the evaluation of its own bonded interaction
30 < (i.e.~bonds, bends, and torsions).
25 > Atoms can be collected into secondary srtructures such as rigid bodies
26 > or molecules. The molecule is a way for {\sc oopse} to keep track of
27 > the atoms in a simulation in logical manner. Molecular units store the
28 > identities of all the atoms associated with themselves, and are
29 > responsible for the evaluation of their own internal interactions
30 > (\emph{i.e.}~bonds, bends, and torsions). Scheme \ref{sch:AtmMole}
31 > shws how one creates a molecule in the \texttt{.mdl} files. The
32 > position of the atoms given in the declaration are relative to the
33 > origin of the molecule, and is used when creating a system containing
34 > the molecule.
35  
36   As stated previously, one of the features that sets {\sc oopse} apart
37   from most of the current molecular simulation packages is the ability
# Line 75 | Line 79 | entire nine parameter rotation matrix internally. Furt
79  
80   {\sc oopse} utilizes a relatively new scheme that propagates the
81   entire nine parameter rotation matrix internally. Further discussion
82 < on this choice can be found in Sec.~\ref{sec:integrate}.
82 > on this choice can be found in Sec.~\ref{sec:integrate}. An example
83 > definition of a riged body can be seen in Scheme
84 > \ref{sch:rigidBody}. The positions in the atom definitions are the
85 > placements of the atoms relative to the origin of the rigid body,
86 > which itself has a position relative to the origin of the molecule.
87  
88 + \begin{lstlisting}[caption={[Defining rigid bodies]A sample definition of a rigid body},label={sch:rigidBody}]
89 + molecule{
90 +  name = "TIP3P_water";
91 +  nRigidBodies = 1;
92 +  rigidBody[0]{
93 +    nAtoms = 3;
94 +    atom[0]{
95 +      type = "O_TIP3P";
96 +      position( 0.0, 0.0, -0.06556 );    
97 +    }                                    
98 +    atom[1]{
99 +      type = "H_TIP3P";
100 +      position( 0.0, 0.75695, 0.52032 );
101 +    }
102 +    atom[2]{
103 +      type = "H_TIP3P";
104 +      position( 0.0, -0.75695, 0.52032 );
105 +    }
106 +    position( 0.0, 0.0, 0.0 );
107 +    orientation( 0.0, 0.0, 1.0 );
108 +  }
109 + }
110 + \end{lstlisting}
111 +
112   \subsection{\label{sec:LJPot}The Lennard Jones Potential}
113  
114 < The most basic force field implemented in {\sc oopse} is the Lennard-Jones
115 < potential. The Lennard-Jones potential. Which mimics the Van der Waals
116 < interaction at long distances, and uses an empirical repulsion at
117 < short distances. The Lennard-Jones potential is given by:
114 > The most basic force field implemented in {\sc oopse} is the
115 > Lennard-Jones potential, which mimics the van der Waals interaction at
116 > long distances, and uses an empirical repulsion at short
117 > distances. The Lennard-Jones potential is given by:
118   \begin{equation}
119   V_{\text{LJ}}(r_{ij}) =
120          4\epsilon_{ij} \biggl[
# Line 91 | Line 123 | V_{\text{LJ}}(r_{ij}) =
123          \biggr]
124   \label{eq:lennardJonesPot}
125   \end{equation}
126 < Where $r_{ij}$ is the distance between particle $i$ and $j$,
126 > Where $r_{ij}$ is the distance between particles $i$ and $j$,
127   $\sigma_{ij}$ scales the length of the interaction, and
128 < $\epsilon_{ij}$ scales the well depth of the potential.
128 > $\epsilon_{ij}$ scales the well depth of the potential. Scheme
129 > \ref{sch:LJFF} gives and example partial \texttt{.bass} file that
130 > shows a system of 108 Ar particles simulated with the Lennard-Jones
131 > force field.
132  
133 + \begin{lstlisting}[caption={[Invocation of the Lennard-Jones force field] A sample system using the Lennard-Jones force field.},label={sch:LJFF}]
134 +
135 + /*
136 + * The Ar molecule is specified
137 + * external to the.bass file
138 + */
139 +
140 + #include "argon.mdl"
141 +
142 + nComponents = 1;
143 + component{
144 +  type = "Ar";
145 +  nMol = 108;
146 + }
147 +
148 + /*
149 + * The initial configuration is generated
150 + * before the simulation is invoked.
151 + */
152 +
153 + initialConfig = "./argon.init";
154 +
155 + forceField = "LJ";
156 + \end{lstlisting}
157 +
158   Because this potential is calculated between all pairs, the force
159   evaluation can become computationally expensive for large systems. To
160 < keep the pair evaluation to a manageable number, {\sc oopse} employs a
161 < cut-off radius.\cite{allen87:csl} The cutoff radius is set to be
162 < $2.5\sigma_{ii}$, where $\sigma_{ii}$ is the largest Lennard-Jones length
163 < parameter in the system. Truncating the calculation at
164 < $r_{\text{cut}}$ introduces a discontinuity into the potential
160 > keep the pair evaluations to a manageable number, {\sc oopse} employs
161 > a cut-off radius.\cite{allen87:csl} The cutoff radius is set to be
162 > $2.5\sigma_{ii}$, where $\sigma_{ii}$ is the largest Lennard-Jones
163 > length parameter present in the simulation. Truncating the calculation
164 > at $r_{\text{cut}}$ introduces a discontinuity into the potential
165   energy. To offset this discontinuity, the energy value at
166 < $r_{\text{cut}}$ is subtracted from the entire potential. This causes
167 < the potential to go to zero at the cut-off radius.
166 > $r_{\text{cut}}$ is subtracted from the potential. This causes the
167 > potential to go to zero smoothly at the cut-off radius.
168  
169   Interactions between dissimilar particles requires the generation of
170   cross term parameters for $\sigma$ and $\epsilon$. These are
# Line 121 | Line 181 | and
181   \end{equation}
182  
183  
184 +
185   \subsection{\label{sec:DUFF}Dipolar Unified-Atom Force Field}
186  
187 < The Dipolar Unified-atom Force Field ({\sc duff}) was developed to
188 < simulate lipid bilayers. The systems require a model capable of forming
189 < bilayers, while still being sufficiently computationally efficient to
190 < allow simulations of large systems ($\approx$100's of phospholipids,
191 < $\approx$1000's of waters) for long times ($\approx$10's of
192 < nanoseconds).
187 > The dipolar unified-atom force field ({\sc duff}) was developed to
188 > simulate lipid bilayers. The simulations require a model capable of
189 > forming bilayers, while still being sufficiently computationally
190 > efficient to allow large systems ($\approx$100's of phospholipids,
191 > $\approx$1000's of waters) to be simulated for long times
192 > ($\approx$10's of nanoseconds).
193  
194 < With this goal in mind, {\sc duff} has no point charges. Charge
195 < neutral distributions were replaced with dipoles, while most atoms and
196 < groups of atoms were reduced to Lennard-Jones interaction sites. This
197 < simplification cuts the length scale of long range interactions from
198 < $\frac{1}{r}$ to $\frac{1}{r^3}$, allowing us to avoid the
199 < computationally expensive Ewald sum. Instead, we can use
194 > With this goal in mind, {\sc duff} has no point
195 > charges. Charge-neutral distributions were replaced with dipoles,
196 > while most atoms and groups of atoms were reduced to Lennard-Jones
197 > interaction sites. This simplification cuts the length scale of long
198 > range interactions from $\frac{1}{r}$ to $\frac{1}{r^3}$, allowing us
199 > to avoid the computationally expensive Ewald sum. Instead, we can use
200   neighbor-lists, reaction field, and cutoff radii for the dipolar
201   interactions.
202  
203   As an example, lipid head-groups in {\sc duff} are represented as
204   point dipole interaction sites. By placing a dipole of 20.6~Debye at
205   the head group center of mass, our model mimics the head group of
206 < phosphatidylcholine.\cite{Cevc87} Additionally, a Lennard-Jones site
207 < is located at the pseudoatom's center of mass. The model is
208 < 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
206 > phosphatidylcholine.\cite{Cevc87} Additionally, a large Lennard-Jones
207 > site is located at the pseudoatom's center of mass. The model is
208 > illustrated by the dark grey atom in Fig.~\ref{fig:lipidModel}. The
209 > water model we use to complement the dipoles of the lipids is our
210   reparameterization of the soft sticky dipole (SSD) model of Ichiye
211   \emph{et al.}\cite{liu96:new_model}
212  
# Line 157 | Line 219 | is the chain length.}
219   \label{fig:lipidModel}
220   \end{figure}
221  
222 < Turning to the tails of the phospholipids, we have used a set of
223 < scalable parameters to model the alkyl groups with Lennard-Jones
224 < sites. For this, we have used the TraPPE force field of Siepmann
222 > We have used a set of scalable parameters to model the alkyl groups
223 > with Lennard-Jones sites. For this, we have borrowed parameters from
224 > the TraPPE force field of Siepmann
225   \emph{et al}.\cite{Siepmann1998} TraPPE is a unified-atom
226   representation of n-alkanes, which is parametrized against phase
227 < equilibria using Gibbs Monte Carlo simulation
227 > equilibria using Gibbs ensemble Monte Carlo simulation
228   techniques.\cite{Siepmann1998} One of the advantages of TraPPE is that
229   it generalizes the types of atoms in an alkyl chain to keep the number
230   of pseudoatoms to a minimum; the parameters for an atom such as
231   $\text{CH}_2$ do not change depending on what species are bonded to
232   it.
233  
234 < TraPPE also constrains of all bonds to be of fixed length. Typically,
234 > TraPPE also constrains all bonds to be of fixed length. Typically,
235   bond vibrations are the fastest motions in a molecular dynamic
236   simulation. Small time steps between force evaluations must be used to
237 < ensure adequate sampling of the bond potential conservation of
238 < energy. By constraining the bond lengths, larger time steps may be
239 < used when integrating the equations of motion.
237 > ensure adequate sampling of the bond potential to ensure conservation
238 > of energy. By constraining the bond lengths, larger time steps may be
239 > used when integrating the equations of motion. A simulation using {\sc
240 > duff} is illustrated in Scheme \ref{sch:DUFF}.
241  
242 + \begin{lstlisting}[caption={[Invocation of {\sc duff}]Sample \texttt{.bass} file showing a simulation utilizing {\sc duff}},label={sch:DUFF}]
243  
244 + #include "water.mdl"
245 + #include "lipid.mdl"
246 +
247 + nComponents = 2;
248 + component{
249 +  type = "simpleLipid_16";
250 +  nMol = 60;
251 + }
252 +
253 + component{
254 +  type = "SSD_water";
255 +  nMol = 1936;
256 + }
257 +
258 + initialConfig = "bilayer.init";
259 +
260 + forceField = "DUFF";
261 +
262 + \end{lstlisting}
263 +
264   \subsubsection{\label{subSec:energyFunctions}{\sc duff} Energy Functions}
265  
266 < The total energy of function in {\sc duff} is given by the following:
266 > The total potential energy function in {\sc duff} is
267   \begin{equation}
268 < V_{\text{Total}} = \sum^{N}_{I=1} V^{I}_{\text{Internal}}
268 > V = \sum^{N}_{I=1} V^{I}_{\text{Internal}}
269          + \sum^{N}_{I=1} \sum_{J>I} V^{IJ}_{\text{Cross}}
270   \label{eq:totalPotential}
271   \end{equation}
272 < Where $V^{I}_{\text{Internal}}$ is the internal potential of a molecule:
272 > Where $V^{I}_{\text{Internal}}$ is the internal potential of molecule $I$:
273   \begin{equation}
274   V^{I}_{\text{Internal}} =
275          \sum_{\theta_{ijk} \in I} V_{\text{bend}}(\theta_{ijk})
# Line 197 | Line 281 | Here $V_{\text{bend}}$ is the bend potential for all 1
281   \label{eq:internalPotential}
282   \end{equation}
283   Here $V_{\text{bend}}$ is the bend potential for all 1, 3 bonded pairs
284 < within the molecule, and $V_{\text{torsion}}$ is the torsion potential
284 > within the molecule $I$, and $V_{\text{torsion}}$ is the torsion potential
285   for all 1, 4 bonded pairs. The pairwise portions of the internal
286   potential are excluded for pairs that are closer than three bonds,
287   i.e.~atom pairs farther away than a torsion are included in the
# Line 206 | Line 290 | The bend potential of a molecule is represented by the
290  
291   The bend potential of a molecule is represented by the following function:
292   \begin{equation}
293 < V_{\theta_{ijk}} = k_{\theta}( \theta_{ijk} - \theta_0 )^2 \label{eq:bendPot}
293 > V_{\text{bend}}(\theta_{ijk}) = k_{\theta}( \theta_{ijk} - \theta_0 )^2 \label{eq:bendPot}
294   \end{equation}
295   Where $\theta_{ijk}$ is the angle defined by atoms $i$, $j$, and $k$
296 < (see Fig.~\ref{fig:lipidModel}), and $\theta_0$ is the equilibrium
297 < bond angle. $k_{\theta}$ is the force constant which determines the
296 > (see Fig.~\ref{fig:lipidModel}), $\theta_0$ is the equilibrium
297 > bond angle, and $k_{\theta}$ is the force constant which determines the
298   strength of the harmonic bend. The parameters for $k_{\theta}$ and
299 < $\theta_0$ are based off of those in TraPPE.\cite{Siepmann1998}
299 > $\theta_0$ are borrowed from those in TraPPE.\cite{Siepmann1998}
300  
301 < The torsion potential and parameters are also taken from TraPPE. It is
301 > The torsion potential and parameters are also borrowed from TraPPE. It is
302   of the form:
303   \begin{equation}
304   V_{\text{torsion}}(\phi) = c_1[1 + \cos \phi]
# Line 222 | Line 306 | V_{\text{torsion}}(\phi) = c_1[1 + \cos \phi]
306          + c_3[1 + \cos(3\phi)]
307   \label{eq:origTorsionPot}
308   \end{equation}
309 < Here $\phi_{ijkl}$ is the angle defined by four bonded neighbors $i$,
309 > Here $\phi$ is the angle defined by four bonded neighbors $i$,
310   $j$, $k$, and $l$ (again, see Fig.~\ref{fig:lipidModel}). For
311   computational efficiency, the torsion potential has been recast after
312 < the method of CHARMM\cite{charmm1983} whereby the angle series is
312 > the method of CHARMM,\cite{charmm1983} in which the angle series is
313   converted to a power series of the form:
314   \begin{equation}
315 < V_{\text{torsion}}(\phi_{ijkl}) =  
315 > V_{\text{torsion}}(\phi) =  
316          k_3 \cos^3 \phi + k_2 \cos^2 \phi + k_1 \cos \phi + k_0
317   \label{eq:torsionPot}
318   \end{equation}
# Line 239 | Line 323 | k_3 &= 4c_3
323   k_2 &= 2 c_2 \\
324   k_3 &= 4c_3
325   \end{align*}
326 < By recasting the equation to a power series, repeated trigonometric
327 < evaluations are avoided during the calculation of the potential.
326 > By recasting the potential as a power series, repeated trigonometric
327 > evaluations are avoided during the calculation of the potential energy.
328  
329  
330 < The cross portion of the total potential, $V^{IJ}_{\text{Cross}}$, is
330 > The cross potential between molecules $I$ and $J$, $V^{IJ}_{\text{Cross}}$, is
331   as follows:
332   \begin{equation}
333   V^{IJ}_{\text{Cross}} =
# Line 257 | Line 341 | $V_{\text{dipole}}$ is the dipole dipole potential, an
341   \end{equation}
342   Where $V_{\text{LJ}}$ is the Lennard Jones potential,
343   $V_{\text{dipole}}$ is the dipole dipole potential, and
344 < $V_{\text{sticky}}$ is the sticky potential defined by the SSD
345 < model. Note that not all atom types include all interactions.
344 > $V_{\text{sticky}}$ is the sticky potential defined by the SSD model
345 > (Sec.~\ref{sec:SSD}). Note that not all atom types include all
346 > interactions.
347  
348   The dipole-dipole potential has the following form:
349   \begin{equation}
# Line 273 | Line 358 | towards $j$, and $\boldsymbol{\Omega}_i$ and $\boldsym
358   \end{equation}
359   Here $\mathbf{r}_{ij}$ is the vector starting at atom $i$ pointing
360   towards $j$, and $\boldsymbol{\Omega}_i$ and $\boldsymbol{\Omega}_j$
361 < are the Euler angles of atom $i$ and $j$ respectively. $|\mu_i|$ is
362 < the magnitude of the dipole moment of atom $i$ and
363 < $\boldsymbol{\hat{u}}_i$ is the standard unit orientation vector of
364 < $\boldsymbol{\Omega}_i$.
361 > are the orientational degrees of freedom for atoms $i$ and $j$
362 > respectively. $|\mu_i|$ is the magnitude of the dipole moment of atom
363 > $i$, $\boldsymbol{\hat{u}}_i$ is the standard unit orientation
364 > vector of $\boldsymbol{\Omega}_i$, and $\boldsymbol{\hat{r}}_{ij}$ is
365 > the unit vector pointing along $\mathbf{r}_{ij}$.
366  
367  
368 < \subsection{\label{sec:SSD}The {\sc duff} Water Models: SSD/E and SSD/RF}
368 > \subsubsection{\label{sec:SSD}The {\sc duff} Water Models: SSD/E and SSD/RF}
369  
370   In the interest of computational efficiency, the default solvent used
371   by {\sc oopse} is the extended Soft Sticky Dipole (SSD/E) water
# Line 368 | Line 454 | density corrected SSD models can be found in reference
454   density corrected SSD models can be found in reference
455   \ref{Gezelter04}.
456  
457 < !!!Place a {\sc BASS} scheme showing SSD parameters around here!!!
457 > \begin{lstlisting}[caption={[A simulation of {\sc ssd} water]An example file showing a simulation including {\sc ssd} water.},label={sch:ssd}]
458  
459 + #include "water.mdl"
460 +
461 + nComponents = 1;
462 + component{
463 +  type = "SSD_water";
464 +  nMol = 864;
465 + }
466 +
467 + initialConfig = "liquidWater.init";
468 +
469 + forceField = "DUFF";
470 +
471 + /*
472 + * The reactionField flag toggles reaction
473 + * field corrections.
474 + */
475 +
476 + reactionField = false; // defaults to false
477 + dielectric = 80.0; // dielectric for reaction field
478 +
479 + /*
480 + * The following two flags set the cutoff
481 + * radius for the electrostatic forces
482 + * as well as the skin thickness of the switching
483 + * function.
484 + */
485 +
486 + electrostaticCutoffRadius  = 9.2;
487 + electrostaticSkinThickness = 1.38;
488 +
489 + \end{lstlisting}
490 +
491 +
492   \subsection{\label{sec:eam}Embedded Atom Method}
493  
494   Several other molecular dynamics packages\cite{dynamo86} exist which have the
# Line 412 | Line 531 | interactions. Foiles et al. fit EAM potentials for fcc
531   \subsection{\label{Sec:pbc}Periodic Boundary Conditions}
532  
533   \newcommand{\roundme}{\operatorname{round}}
534 <
534 >
535   \textit{Periodic boundary conditions} are widely used to simulate truly
536   macroscopic systems with a relatively small number of particles. The
537 < simulation box is replicated throughout space to form an infinite
538 < lattice.  During the simulation, when a particle moves in the primary
539 < cell, its images in other boxes move in exactly the same direction with
540 < exactly the same orientation. So, as a particle leaves the primary
541 < cell, one of its images will enter through the opposite face.If the
542 < simulation box is large enough to avoid \textquotedblleft feeling\textquotedblright\ the symmetries of
543 < the periodic lattice, surface effects can be ignored. Cubic,
544 < orthorhombic and parallelepiped are the available periodic cells in
545 < {\sc oopse}. We use a matrix to describe the property of the simulation
546 < box. Both the size and shape of the simulation box can be
547 < changed during the simulation. The transformation from box space
548 < vector $\mathbf{s}$ to its corresponding real space vector
430 < $\mathbf{r}$ is defined by
537 > simulation box is replicated throughout space to form an infinite lattice.
538 > During the simulation, when a particle moves in the primary cell, its image in
539 > other boxes move in exactly the same direction with exactly the same
540 > orientation.Thus, as a particle leaves the primary cell, one of its images
541 > will enter through the opposite face.If the simulation box is large enough to
542 > avoid \textquotedblleft feeling\textquotedblright\ the symmetries of the
543 > periodic lattice, surface effects can be ignored. Cubic, orthorhombic and
544 > parallelepiped are the available periodic cells In OOPSE. We use a matrix to
545 > describe the property of the simulation box. Therefore, both the size and
546 > shape of the simulation box can be changed during the simulation. The
547 > transformation from box space vector $\mathbf{s}$ to its corresponding real
548 > space vector $\mathbf{r}$ is defined by
549   \begin{equation}
550   \mathbf{r}=\underline{\mathbf{H}}\cdot\mathbf{s}%
551   \end{equation}
552  
553  
554 < where $H=(h_{x},h_{y},h_{z})$ is a transformation matrix made up of
555 < the three box axis vectors. $h_{x},h_{y}$ and $h_{z}$ represent the
556 < three sides of the simulation box respectively.
554 > where $H=(h_{x},h_{y},h_{z})$ is a transformation matrix made up of the three
555 > box axis vectors. $h_{x},h_{y}$ and $h_{z}$ represent the three sides of the
556 > simulation box respectively.
557  
558 < To find the minimum image of a vector $\mathbf{r}$, we convert the real vector to its
559 < corresponding vector in box space first, \bigskip%
558 > To find the minimum image of a vector $\mathbf{r}$, we convert the real vector
559 > to its corresponding vector in box space first, \bigskip%
560   \begin{equation}
561   \mathbf{s}=\underline{\mathbf{H}}^{-1}\cdot\mathbf{r}%
562   \end{equation}
# Line 452 | Line 570 | where
570  
571   \begin{equation}
572   \roundme(x)=\left\{
573 < \begin{array}{cc}
573 > \begin{array}{cc}%
574   \lfloor{x+0.5}\rfloor & \text{if \ }x\geqslant0\\
575   \lceil{x-0.5}\rceil & \text{otherwise}%
576   \end{array}
577   \right.
578   \end{equation}
461 – For example, $\roundme(3.6)=4$, $\roundme(3.1)=3$, $\roundme(-3.6)=-4$,
462 – $\roundme(-3.1)=-3$.
579  
464 – Finally, we obtain the minimum image coordinates $\mathbf{r}^{\prime}$ by transforming back
465 – to real space,%
580  
581 + For example, $\roundme(3.6)=4$,$\roundme(3.1)=3$, $\roundme(-3.6)=-4$, $\roundme(-3.1)=-3$.
582 +
583 + Finally, we obtain the minimum image coordinates $\mathbf{r}^{\prime}$ by
584 + transforming back to real space,%
585 +
586   \begin{equation}
587   \mathbf{r}^{\prime}=\underline{\mathbf{H}}^{-1}\cdot\mathbf{s}^{\prime}%
588   \end{equation}

Diff Legend

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