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 899 by mmeineke, Tue Jan 6 18:53:58 2004 UTC vs.
Revision 973 by mmeineke, Wed Jan 21 19:00:23 2004 UTC

# Line 1 | Line 1
1  
2 < \section{\label{sec:empericalEnergy}The Emperical Energy Functions}
2 > \section{\label{sec:empiricalEnergy}The Empirical Energy Functions}
3  
4 < \subsection{\label{sec:atomsMolecules}Atoms and Molecules}
4 > \subsection{\label{sec:atomsMolecules}Atoms, Molecules and Rigid Bodies}
5  
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 {\sc oopse}.
6 > The basic unit of an {\sc oopse} simulation is the atom. The
7 > parameters describing the atom are generalized to make the atom as
8 > flexible a representation as possible. They may represent specific
9 > atoms of an element, or be used for collections of atoms such as
10 > methyl and carbonyl groups. The atoms are also capable of having
11 > directional components associated with them (\emph{e.g.}~permanent
12 > dipoles). Charges on atoms are not currently supported 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 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).
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 );
21 >  }
22 > }
23 > \end{lstlisting}
24  
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
38 + to handle rigid body dynamics. Rigid bodies are non-spherical
39 + particles or collections of particles that have a constant internal
40 + potential and move collectively.\cite{Goldstein01} They are not
41 + included in most simulation packages because of the requirement to
42 + propagate the orientational degrees of freedom. Until recently,
43 + integrators which propagate orientational motion have been lacking.
44 +
45 + Moving a rigid body involves determination of both the force and
46 + torque applied by the surroundings, which directly affect the
47 + translational and rotational motion in turn. In order to accumulate
48 + the total force on a rigid body, the external forces and torques must
49 + first be calculated for all the internal particles. The total force on
50 + the rigid body is simply the sum of these external forces.
51 + Accumulation of the total torque on the rigid body is more complex
52 + than the force in that it is the torque applied on the center of mass
53 + that dictates rotational motion. The torque on rigid body {\it i} is
54 + \begin{equation}
55 + \boldsymbol{\tau}_i=
56 +        \sum_{a}(\mathbf{r}_{ia}-\mathbf{r}_i)\times \mathbf{f}_{ia}
57 +        + \boldsymbol{\tau}_{ia},
58 + \label{eq:torqueAccumulate}
59 + \end{equation}
60 + where $\boldsymbol{\tau}_i$ and $\mathbf{r}_i$ are the torque on and
61 + position of the center of mass respectively, while $\mathbf{f}_{ia}$,
62 + $\mathbf{r}_{ia}$, and $\boldsymbol{\tau}_{ia}$ are the force on,
63 + position of, and torque on the component particles of the rigid body.
64 +
65 + The summation of the total torque is done in the body fixed axis of
66 + the rigid body. In order to move between the space fixed and body
67 + fixed coordinate axes, parameters describing the orientation must be
68 + maintained for each rigid body. At a minimum, the rotation matrix
69 + (\textbf{A}) can be described by the three Euler angles ($\phi,
70 + \theta,$ and $\psi$), where the elements of \textbf{A} are composed of
71 + trigonometric operations involving $\phi, \theta,$ and
72 + $\psi$.\cite{Goldstein01} In order to avoid numerical instabilities
73 + inherent in using the Euler angles, the four parameter ``quaternion''
74 + scheme is often used. The elements of \textbf{A} can be expressed as
75 + arithmetic operations involving the four quaternions ($q_0, q_1, q_2,$
76 + and $q_3$).\cite{allen87:csl} Use of quaternions also leads to
77 + performance enhancements, particularly for very small
78 + systems.\cite{Evans77}
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}. 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
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[
121 +        \biggl(\frac{\sigma_{ij}}{r_{ij}}\biggr)^{12}
122 +        - \biggl(\frac{\sigma_{ij}}{r_{ij}}\biggr)^{6}
123 +        \biggr]
124 + \label{eq:lennardJonesPot}
125 + \end{equation}
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. 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 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 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
171 + calculated through the Lorentz-Berthelot mixing
172 + rules:\cite{allen87:csl}
173 + \begin{equation}
174 + \sigma_{ij} = \frac{1}{2}[\sigma_{ii} + \sigma_{jj}]
175 + \label{eq:sigmaMix}
176 + \end{equation}
177 + and
178 + \begin{equation}
179 + \epsilon_{ij} = \sqrt{\epsilon_{ii} \epsilon_{jj}}
180 + \label{eq:epsilonMix}
181 + \end{equation}
182 +
183 +
184 +
185   \subsection{\label{sec:DUFF}Dipolar Unified-Atom Force Field}
186  
187 < The \underline{D}ipolar \underline{U}nified-Atom
188 < \underline{F}orce \underline{F}ield ({\sc duff}) was developed to
189 < simulate lipid bilayers. We needed a model capable of forming
190 < bilayers, while still being sufficiently computationally efficient to
191 < allow simulations of large systems ($\approx$100's of phospholipids,
192 < $\approx$1000's of waters) for long times ($\approx$10's of
29 < 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, we have eliminated all point charges. Charge
195 < distributions were replaced with dipoles, and charge-neutral
196 < distributions 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
200 < neighbor-lists and cutoff radii for the dipolar interactions.
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 point
204 < dipole interaction sites.PC and PE Lipid head groups are typically
205 < zwitterionic in nature, with charges separated by as much as
206 < 6~$\mbox{\AA}$. By placing a dipole of 20.6~Debye at the head group
207 < center of mass, our model mimics the head group of PC.\cite{Cevc87}
208 < Additionally, a Lennard-Jones site is located at the
209 < pseudoatom's center of mass. The model is illustrated by the dark grey
210 < atom in Fig.~\ref{fig:lipidModel}.
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 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  
213   \begin{figure}
214 < \epsfxsize=6in
215 < \epsfbox{lipidModel.epsi}
214 > \epsfxsize=\linewidth
215 > \epsfbox{lipidModel.eps}
216   \caption{A representation of the lipid model. $\phi$ is the torsion angle, $\theta$ %
217 < is the bend angle, $\mu$ is the dipole moment of the head group, and n is the chain length.}
217 > is the bend angle, $\mu$ is the dipole moment of the head group, and n
218 > is the chain length.}
219   \label{fig:lipidModel}
220   \end{figure}
221  
222 < The water model we use to complement the dipoles of the lipids is
223 < the soft sticky dipole (SSD) model of Ichiye \emph{et
224 < 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
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}}
269 <        + \sum^{N}_{I=1} \sum^{N}_{J=1} V^{IJ}_{\text{Cross}}
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})
276 <        + \sum_{\phi_{ijkl} \in I} V_{\text{torsion}}(\theta_{ijkl})
276 >        + \sum_{\phi_{ijkl} \in I} V_{\text{torsion}}(\phi_{ijkl})
277          + \sum_{i \in I} \sum_{(j>i+4) \in I}
278          \biggl[ V_{\text{LJ}}(r_{ij}) +  V_{\text{dipole}}
279          (\mathbf{r}_{ij},\boldsymbol{\Omega}_{i},\boldsymbol{\Omega}_{j})
# Line 106 | 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 in the molecule. $V_{\text{torsion}}$ is the torsion The
285 < pairwise portions of the internal potential are excluded for pairs
286 < that ar closer than three bonds, i.e.~atom pairs farther away than a
287 < torsion are included in the pair-wise loop.
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
288 > pair-wise loop.
289  
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.
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_{ijkl}) = c_1[1 + \cos \phi]
304 > V_{\text{torsion}}(\phi) = c_1[1 + \cos \phi]
305          + c_2[1 + \cos(2\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$,
310 < $j$, $k$, and $l$ (again, see Fig.~\ref{fig:lipidModel}).  However,
311 < for computaional efficency, the torsion potentail has been recast
312 < after the method of CHARMM\cite{charmm1983} whereby the angle series
313 < is converted to a power series of the form:
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} 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 162 | 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 < The Lennard-Jones potential is given by:
329 >
330 > The cross potential between molecules $I$ and $J$, $V^{IJ}_{\text{Cross}}$, is
331 > as follows:
332   \begin{equation}
333 < V_{\text{LJ}}(r_{ij}) =
334 <        4\epsilon_{ij} \biggl[
335 <        \biggl(\frac{\sigma_{ij}}{r_{ij}}\biggr)^{12}
336 <        - \biggl(\frac{\sigma_{ij}}{r_{ij}}\biggr)^{6}
333 > V^{IJ}_{\text{Cross}} =
334 >        \sum_{i \in I} \sum_{j \in J}
335 >        \biggl[ V_{\text{LJ}}(r_{ij}) +  V_{\text{dipole}}
336 >        (\mathbf{r}_{ij},\boldsymbol{\Omega}_{i},\boldsymbol{\Omega}_{j})
337 >        + V_{\text{sticky}}
338 >        (\mathbf{r}_{ij},\boldsymbol{\Omega}_{i},\boldsymbol{\Omega}_{j})
339          \biggr]
340 < \label{eq:lennardJonesPot}
340 > \label{eq:crossPotentail}
341   \end{equation}
342 < Where $r_ij$ is the distance between atoms $i$ and $j$, $\sigma_{ij}$
343 < scales the length of the interaction, and $\epsilon_{ij}$ scales the
344 < energy of the potential.
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 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}
350   V_{\text{dipole}}(\mathbf{r}_{ij},\boldsymbol{\Omega}_{i},
351 <        \boldsymbol{\Omega}_{j}) = \frac{1}{4\pi\epsilon_{0}} \biggl[
352 <        \frac{\boldsymbol{\mu}_{i} \cdot \boldsymbol{\mu}_{j}}{r^{3}_{ij}}
351 >        \boldsymbol{\Omega}_{j}) = \frac{|\mu_i||\mu_j|}{4\pi\epsilon_{0}r_{ij}^{3}} \biggl[
352 >        \boldsymbol{\hat{u}}_{i} \cdot \boldsymbol{\hat{u}}_{j}
353          -
354 <        \frac{3(\boldsymbol{\mu}_i \cdot \mathbf{r}_{ij}) %
355 <                (\boldsymbol{\mu}_j \cdot \mathbf{r}_{ij}) }
356 <                {r^{5}_{ij}} \biggr]
354 >        \frac{3(\boldsymbol{\hat{u}}_i \cdot \mathbf{r}_{ij}) %
355 >                (\boldsymbol{\hat{u}}_j \cdot \mathbf{r}_{ij}) }
356 >                {r^{2}_{ij}} \biggr]
357   \label{eq:dipolePot}
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$
362 < respectively. $\boldsymbol{\mu}_i$ is the dipole vector of atom
363 < $i$ it takes its orientation from $\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}Water Model: SSD and Derivatives}
368 > \subsubsection{\label{sec:SSD}The {\sc duff} Water Models: SSD/E and SSD/RF}
369  
370 < In the interest of computational efficiency, the native solvent used
371 < in {\sc oopse} is the Soft Sticky Dipole (SSD) water model. SSD was
372 < developed by Ichiye \emph{et al.} as a modified form of the
373 < hard-sphere water model proposed by Bratko, Blum, and
370 > In the interest of computational efficiency, the default solvent used
371 > by {\sc oopse} is the extended Soft Sticky Dipole (SSD/E) water
372 > model.\cite{Gezelter04} The original SSD was developed by Ichiye
373 > \emph{et al.}\cite{liu96:new_model} as a modified form of the hard-sphere
374 > water model proposed by Bratko, Blum, and
375   Luzar.\cite{Bratko85,Bratko95} It consists of a single point dipole
376   with a Lennard-Jones core and a sticky potential that directs the
377   particles to assume the proper hydrogen bond orientation in the first
378   solvation shell. Thus, the interaction between two SSD water molecules
379   \emph{i} and \emph{j} is given by the potential
380   \begin{equation}
381 < u_{ij} = u_{ij}^{LJ} (r_{ij})\ + u_{ij}^{dp}
382 < (\mathbf{r}_{ij},\boldsymbol{\Omega}_i,\boldsymbol{\Omega}_j)\ +
383 < u_{ij}^{sp}
384 < (\mathbf{r}_{ij},\boldsymbol{\Omega}_i,\boldsymbol{\Omega}_j),
381 > V_{ij} =
382 >        V_{ij}^{LJ} (r_{ij})\ + V_{ij}^{dp}
383 >        (\mathbf{r}_{ij},\boldsymbol{\Omega}_i,\boldsymbol{\Omega}_j)\ +
384 >        V_{ij}^{sp}
385 >        (\mathbf{r}_{ij},\boldsymbol{\Omega}_i,\boldsymbol{\Omega}_j),
386 > \label{eq:ssdPot}
387   \end{equation}
388   where the $\mathbf{r}_{ij}$ is the position vector between molecules
389 < \emph{i} and \emph{j} with magnitude equal to the distance $r_ij$, and
389 > \emph{i} and \emph{j} with magnitude equal to the distance $r_{ij}$, and
390   $\boldsymbol{\Omega}_i$ and $\boldsymbol{\Omega}_j$ represent the
391 < orientations of the respective molecules. The Lennard-Jones, dipole,
392 < and sticky parts of the potential are giving by the following
393 < equations,
391 > orientations of the respective molecules. The Lennard-Jones and dipole
392 > parts of the potential are given by equations \ref{eq:lennardJonesPot}
393 > and \ref{eq:dipolePot} respectively. The sticky part is described by
394 > the following,
395   \begin{equation}
396 < u_{ij}^{LJ}(r_{ij}) = 4\epsilon \left[\left(\frac{\sigma}{r_{ij}}\right)^{12}-\left(\frac{\sigma}{r_{ij}}\right)^{6}\right],
396 > u_{ij}^{sp}(\mathbf{r}_{ij},\boldsymbol{\Omega}_i,\boldsymbol{\Omega}_j)=
397 >        \frac{\nu_0}{2}[s(r_{ij})w(\mathbf{r}_{ij},
398 >        \boldsymbol{\Omega}_i,\boldsymbol{\Omega}_j) +
399 >        s^\prime(r_{ij})w^\prime(\mathbf{r}_{ij},
400 >        \boldsymbol{\Omega}_i,\boldsymbol{\Omega}_j)]\ ,
401 > \label{eq:stickyPot}
402   \end{equation}
403 + where $\nu_0$ is a strength parameter for the sticky potential, and
404 + $s$ and $s^\prime$ are cubic switching functions which turn off the
405 + sticky interaction beyond the first solvation shell. The $w$ function
406 + can be thought of as an attractive potential with tetrahedral
407 + geometry:
408   \begin{equation}
409 < 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}\ ,
409 > w({\bf r}_{ij},{\bf \Omega}_i,{\bf \Omega}_j)=
410 >        \sin\theta_{ij}\sin2\theta_{ij}\cos2\phi_{ij},
411 > \label{eq:stickyW}
412   \end{equation}
413 + while the $w^\prime$ function counters the normal aligned and
414 + anti-aligned structures favored by point dipoles:
415   \begin{equation}
416 < \begin{split}
417 < u_{ij}^{sp}
418 < (\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}
416 > w^\prime({\bf r}_{ij},{\bf \Omega}_i,{\bf \Omega}_j)=
417 >        (\cos\theta_{ij}-0.6)^2(\cos\theta_{ij}+0.8)^2-w^0,
418 > \label{eq:stickyWprime}
419   \end{equation}
420 < where $\boldsymbol{\mu}_i$ and $\boldsymbol{\mu}_j$ are the dipole
421 < unit vectors of particles \emph{i} and \emph{j} with magnitude 2.35 D,
422 < $\nu_0$ scales the strength of the overall sticky potential, $s$ and
423 < $s^\prime$ are cubic switching functions. The $w$ and $w^\prime$
424 < functions take the following forms,
420 > It should be noted that $w$ is proportional to the sum of the $Y_3^2$
421 > and $Y_3^{-2}$ spherical harmonics (a linear combination which
422 > enhances the tetrahedral geometry for hydrogen bonded structures),
423 > while $w^\prime$ is a purely empirical function.  A more detailed
424 > description of the functional parts and variables in this potential
425 > can be found in the original SSD
426 > articles.\cite{liu96:new_model,liu96:monte_carlo,chandra99:ssd_md,Ichiye03}
427 >
428 > Since SSD is a single-point {\it dipolar} model, the force
429 > calculations are simplified significantly relative to the standard
430 > {\it charged} multi-point models. In the original Monte Carlo
431 > simulations using this model, Ichiye {\it et al.} reported that using
432 > SSD decreased computer time by a factor of 6-7 compared to other
433 > models.\cite{liu96:new_model} What is most impressive is that these savings
434 > did not come at the expense of accurate depiction of the liquid state
435 > properties.  Indeed, SSD maintains reasonable agreement with the Soper
436 > diffraction data for the structural features of liquid
437 > water.\cite{Soper86,liu96:new_model} Additionally, the dynamical properties
438 > exhibited by SSD agree with experiment better than those of more
439 > computationally expensive models (like TIP3P and
440 > SPC/E).\cite{chandra99:ssd_md} The combination of speed and accurate depiction
441 > of solvent properties makes SSD a very attractive model for the
442 > simulation of large scale biochemical simulations.
443 >
444 > Recent constant pressure simulations revealed issues in the original
445 > SSD model that led to lower than expected densities at all target
446 > pressures.\cite{Ichiye03,Gezelter04} The default model in {\sc oopse}
447 > is therefore SSD/E, a density corrected derivative of SSD that
448 > exhibits improved liquid structure and transport behavior. If the use
449 > of a reaction field long-range interaction correction is desired, it
450 > is recommended that the parameters be modified to those of the SSD/RF
451 > model. Solvent parameters can be easily modified in an accompanying
452 > {\sc BASS} file as illustrated in the scheme below. A table of the
453 > parameter values and the drawbacks and benefits of the different
454 > density corrected SSD models can be found in reference
455 > \ref{Gezelter04}.
456 >
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
495 > capacity to simulate metallic systems, including some that have
496 > parallel computational abilities\cite{plimpton93}. Potentials that
497 > describe bonding transition metal
498 > systems\cite{Finnis84,Ercolessi88,Chen90,Qi99,Ercolessi02} have a
499 > attractive interaction which models  ``Embedding''
500 > a positively charged metal ion in the electron density due to the
501 > free valance ``sea'' of electrons created by the surrounding atoms in
502 > the system. A mostly repulsive pairwise part of the potential
503 > describes the interaction of the positively charged metal core ions
504 > with one another. A particular potential description called the
505 > Embedded Atom Method\cite{Daw84,FBD86,johnson89,Lu97}({\sc eam}) that has
506 > particularly wide adoption has been selected for inclusion in {\sc oopse}. A
507 > good review of {\sc eam} and other metallic potential formulations was done
508 > by Voter.\cite{voter}
509 >
510 > The {\sc eam} potential has the form:
511 > \begin{eqnarray}
512 > V & = & \sum_{i} F_{i}\left[\rho_{i}\right] + \sum_{i} \sum_{j \neq i}
513 > \phi_{ij}({\bf r}_{ij})  \\
514 > \rho_{i}  & = & \sum_{j \neq i} f_{j}({\bf r}_{ij})
515 > \end{eqnarray}S
516 >
517 > where $F_{i} $ is the embedding function that equates the energy required to embed a
518 > positively-charged core ion $i$ into a linear superposition of
519 > spherically averaged atomic electron densities given by
520 > $\rho_{i}$.  $\phi_{ij}$ is a primarily repulsive pairwise interaction
521 > between atoms $i$ and $j$. In the original formulation of
522 > {\sc eam} cite{Daw84}, $\phi_{ij}$ was an entirely repulsive term, however
523 > in later refinements to EAM have shown that non-uniqueness between $F$
524 > and $\phi$ allow for more general forms for $\phi$.\cite{Daw89}
525 > There is a cutoff distance, $r_{cut}$, which limits the
526 > summations in the {\sc eam} equation to the few dozen atoms
527 > surrounding atom $i$ for both the density $\rho$ and pairwise $\phi$
528 > 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}.
529 >
530 >
531 > \subsection{\label{Sec:pbc}Periodic Boundary Conditions}
532 >
533 > \newcommand{\roundme}{\operatorname{round}}
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 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 < w(\mathbf{r}_{ij},\boldsymbol{\Omega}_i,\boldsymbol{\Omega}_j)=\sin\theta_{ij}\sin2\theta_{ij}\cos2\phi_{ij},
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 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
559 + to its corresponding vector in box space first, \bigskip%
560   \begin{equation}
561 < 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,
561 > \mathbf{s}=\underline{\mathbf{H}}^{-1}\cdot\mathbf{r}%
562   \end{equation}
563 < where $w^0 = 0.07715$. The $w$ function is the tetrahedral attractive
564 < term that promotes hydrogen bonding orientations within the first
565 < solvation shell, and $w^\prime$ is a dipolar repulsion term that
566 < repels unrealistic dipolar arrangements within the first solvation
567 < 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}
563 > And then, each element of $\mathbf{s}$ is wrapped to lie between -0.5 to 0.5,
564 > \begin{equation}
565 > s_{i}^{\prime}=s_{i}-\roundme(s_{i})
566 > \end{equation}
567 > where
568  
569 < 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}
569 > %
570  
571 < Recent constant pressure simulations revealed issues in the original
572 < SSD model that led to lower than expected densities at all target
573 < pressures.\cite{Ichiye03,Gezelter04} Reparameterizations of the
574 < original SSD have resulted in improved density behavior, as well as
575 < alterations in the water structure and transport behavior. {\sc oopse} is
576 < easily modified to impliment these new potential parameter sets for
577 < the derivative water models: SSD1, SSD/E, and SSD/RF. All of the
578 < variable parameters are listed in the accompanying BASS file, and
275 < these parameters simply need to be changed to the updated values.
571 > \begin{equation}
572 > \roundme(x)=\left\{
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}
579  
580  
581 < \subsection{\label{sec:eam}Embedded Atom Model}
581 > For example, $\roundme(3.6)=4$,$\roundme(3.1)=3$, $\roundme(-3.6)=-4$, $\roundme(-3.1)=-3$.
582  
583 < here there be Monsters
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}
589 >

Diff Legend

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