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

Comparing trunk/mattDisertation/oopse.tex (file contents):
Revision 1054 by mmeineke, Mon Feb 16 20:35:58 2004 UTC vs.
Revision 1068 by mmeineke, Wed Feb 25 21:48:44 2004 UTC

# Line 22 | Line 22 | Gezelter. Although my contributions to {\sc oopse} are
22   simulation package of this size and scope would not have been possible
23   without the collaborative efforts of my colleagues: Charles
24   F.~Vardeman II, Teng Lin, Christopher J.~Fennell and J.~Daniel
25 < Gezelter. Although my contributions to {\sc oopse} are significant,
26 < consideration of my work apart from the others, would not give a
25 > Gezelter. Although my contributions to {\sc oopse} are major,
26 > consideration of my work apart from the others would not give a
27   complete description to the package's capabilities. As such, all
28   contributions to {\sc oopse} to date are presented in this chapter.
29  
30 < Charles Vardeman is responsible for the parallelization of {\sc oopse}
31 < (Sec.~\ref{oopseSec:parallelization}) as well as the inclusion of the
32 < embedded-atom potential (Sec.~\ref{oopseSec:eam}). Teng Lin's
33 < contributions include refinement of the periodic boundary conditions
30 > Charles Vardeman is responsible for the parallelization of the long
31 > range forces in {\sc oopse} (Sec.~\ref{oopseSec:parallelization}) as
32 > well as the inclusion of the embedded-atom potential for transition
33 > metals (Sec.~\ref{oopseSec:eam}). Teng Lin's contributions include
34 > refinement of the periodic boundary conditions
35   (Sec.~\ref{oopseSec:pbc}), the z-constraint method
36   (Sec.~\ref{oopseSec:zcons}), refinement of the property analysis
37   programs (Sec.~\ref{oopseSec:props}), and development in the extended
38 < state integrators (Sec.~\ref{oopseSec:noseHooverThermo}). Christopher
38 > system integrators (Sec.~\ref{oopseSec:noseHooverThermo}). Christopher
39   Fennell worked on the symplectic integrator
40   (Sec.~\ref{oopseSec:integrate}) and the refinement of the {\sc ssd}
41   water model (Sec.~\ref{oopseSec:SSD}). Daniel Gezelter lent his
# Line 48 | Line 49 | initialization (Sec.~\ref{oopseSec:initCoords}) utilit
49   of the Lennard-Jones (Sec.~\ref{sec:LJPot}) and {\sc duff}
50   (Sec.~\ref{oopseSec:DUFF}) force fields, and initial implementation of
51   the property analysis (Sec.~\ref{oopseSec:props}) and system
52 < initialization (Sec.~\ref{oopseSec:initCoords}) utility programs.
52 > initialization (Sec.~\ref{oopseSec:initCoords}) utility programs. {\sc
53 > oopse}, like many other Molecular Dynamics programs, is a work in
54 > progress, and will continue to be so for many graduate student
55 > lifetimes.
56  
57   \section{\label{sec:intro}Introduction}
58  
# Line 66 | Line 70 | code was not originally designed to do. Examples of un
70  
71   Despite their utility, problems with these packages arise when
72   researchers try to develop techniques or energetic models that the
73 < code was not originally designed to do. Examples of uncommonly
73 > code was not originally designed to simulate. Examples of uncommonly
74   implemented techniques and energetics include; dipole-dipole
75   interactions, rigid body dynamics, and metallic embedded
76   potentials. When faced with these obstacles, a researcher must either
# Line 75 | Line 79 | Having written {\sc oopse} we are implementing the con
79   simulation code capable of implementing the types of models upon which
80   our research is based.
81  
82 < Having written {\sc oopse} we are implementing the concept of Open
83 < Source development, and releasing our source code into the public
84 < domain. It is our intent that by doing so, other researchers might
85 < benefit from our work, and add their own contributions to the
86 < package. The license under which {\sc oopse} is distributed allows any
87 < researcher to download and modify the source code for their own
88 < use. In this way further development of {\sc oopse} is not limited to
89 < only the models of interest to ourselves, but also those of the
90 < community of scientists who contribute back to the project.
82 > In developing {\sc oopse}, we have adhered to the precepts of Open
83 > Source development, and are releasing our source code with a
84 > permissive license. It is our intent that by doing so, other
85 > researchers might benefit from our work, and add their own
86 > contributions to the package. The license under which {\sc oopse} is
87 > distributed allows any researcher to download and modify the source
88 > code for their own use. In this way further development of {\sc oopse}
89 > is not limited to only the models of interest to ourselves, but also
90 > those of the community of scientists who contribute back to the
91 > project.
92  
93   We have structured this chapter to first discuss the empirical energy
94   functions that {\sc oopse } implements in
95   Sec.~\ref{oopseSec:empiricalEnergy}. Following that is a discussion of
96   the various input and output files associated with the package
97 < (Sec.~\ref{oopseSec:IOfiles}). In Sec.~\ref{oopseSec:mechanics}
97 > (Sec.~\ref{oopseSec:IOfiles}). Sec.~\ref{oopseSec:mechanics}
98   elucidates the various Molecular Dynamics algorithms {\sc oopse}
99   implements in the integration of the Newtonian equations of
100   motion. Basic analysis of the trajectories obtained from the
101   simulation is discussed in Sec.~\ref{oopseSec:props}. Program design
102 < considerations as well as the software distribution license is
103 < presented in Sec.~\ref{oopseSec:design}. And lastly,
99 < Sec.~\ref{oopseSec:conclusion} concludes the chapter.
102 > considerations are presented in Sec.~\ref{oopseSec:design}. And
103 > lastly, Sec.~\ref{oopseSec:conclusion} concludes the chapter.
104  
105   \section{\label{oopseSec:empiricalEnergy}The Empirical Energy Functions}
106  
# Line 109 | Line 113 | a given atom type are set in the force field parameter
113   methyl and carbonyl groups. The atoms are also capable of having
114   directional components associated with them (\emph{e.g.}~permanent
115   dipoles). Charges, permanent dipoles, and Lennard-Jones parameters for
116 < a given atom type are set in the force field parameter files
116 > a given atom type are set in the force field parameter files.
117  
118   \begin{lstlisting}[float,caption={[Specifier for molecules and atoms] A sample specification of an Ar molecule},label=sch:AtmMole]
119   molecule{
# Line 129 | Line 133 | internal interactions (\emph{i.e.}~bonds, bends, and t
133   identities of all the atoms and rigid bodies associated with
134   themselves, and are responsible for the evaluation of their own
135   internal interactions (\emph{i.e.}~bonds, bends, and torsions). Scheme
136 < \ref{sch:AtmMole} shows how one creates a molecule in a
136 > \ref{sch:AtmMole} shows how one creates a molecule in a ``model'' or
137   \texttt{.mdl} file. The position of the atoms given in the
138   declaration are relative to the origin of the molecule, and is used
139   when creating a system containing the molecule.
# Line 139 | Line 143 | included in most simulation packages because of the re
143   to handle rigid body dynamics. Rigid bodies are non-spherical
144   particles or collections of particles that have a constant internal
145   potential and move collectively.\cite{Goldstein01} They are not
146 < included in most simulation packages because of the requirement to
147 < propagate the orientational degrees of freedom. Until recently,
148 < integrators which propagate orientational motion have been lacking.
146 > included in most simulation packages because of the algorithmic
147 > complexity involved in propagating orientational degrees of
148 > freedom. Until recently, integrators which propagate orientational
149 > motion have been much worse than those available for translational
150 > motion.
151  
152   Moving a rigid body involves determination of both the force and
153   torque applied by the surroundings, which directly affect the
# Line 150 | Line 156 | than the force in that it is the torque applied on the
156   first be calculated for all the internal particles. The total force on
157   the rigid body is simply the sum of these external forces.
158   Accumulation of the total torque on the rigid body is more complex
159 < than the force in that it is the torque applied on the center of mass
160 < that dictates rotational motion. The torque on rigid body {\it i} is
159 > than the force because the torque is applied to the center of mass of
160 > the rigid body. The torque on rigid body $i$ is
161   \begin{equation}
162   \boldsymbol{\tau}_i=
163 <        \sum_{a}(\mathbf{r}_{ia}-\mathbf{r}_i)\times \mathbf{f}_{ia}
164 <        + \boldsymbol{\tau}_{ia},
163 >        \sum_{a}\biggl[(\mathbf{r}_{ia}-\mathbf{r}_i)\times \mathbf{f}_{ia}
164 >        + \boldsymbol{\tau}_{ia}\biggr]
165   \label{eq:torqueAccumulate}
166   \end{equation}
167   where $\boldsymbol{\tau}_i$ and $\mathbf{r}_i$ are the torque on and
# Line 164 | Line 170 | the rigid body. In order to move between the space fix
170   position of, and torque on the component particles of the rigid body.
171  
172   The summation of the total torque is done in the body fixed axis of
173 < the rigid body. In order to move between the space fixed and body
173 > each rigid body. In order to move between the space fixed and body
174   fixed coordinate axes, parameters describing the orientation must be
175   maintained for each rigid body. At a minimum, the rotation matrix
176   (\textbf{A}) can be described by the three Euler angles ($\phi,
# Line 179 | Line 185 | entire nine parameter rotation matrix internally. Furt
185   systems.\cite{Evans77}
186  
187   {\sc oopse} utilizes a relatively new scheme that propagates the
188 < entire nine parameter rotation matrix internally. Further discussion
188 > entire nine parameter rotation matrix. Further discussion
189   on this choice can be found in Sec.~\ref{oopseSec:integrate}. An example
190   definition of a rigid body can be seen in Scheme
191   \ref{sch:rigidBody}. The positions in the atom definitions are the
# Line 232 | Line 238 | Lennard-Jones force field.
238   Lennard-Jones force field.
239  
240   \begin{lstlisting}[float,caption={[Invocation of the Lennard-Jones force field] A sample system using the Lennard-Jones force field.},label={sch:LJFF}]
235 –
236 – /*
237 – * The Ar molecule is specified
238 – * external to the.bass file
239 – */
241  
242   #include "argon.mdl"
243  
# Line 246 | Line 247 | component{
247    nMol = 108;
248   }
249  
249 – /*
250 – * The initial configuration is generated
251 – * before the simulation is invoked.
252 – */
253 –
250   initialConfig = "./argon.init";
251  
252   forceField = "LJ";
# Line 259 | Line 255 | a cut-off radius.\cite{allen87:csl} The cutoff radius
255   Because this potential is calculated between all pairs, the force
256   evaluation can become computationally expensive for large systems. To
257   keep the pair evaluations to a manageable number, {\sc oopse} employs
258 < a cut-off radius.\cite{allen87:csl} The cutoff radius is set to be
258 > a cut-off radius.\cite{allen87:csl} The cutoff radius can either be
259 > specified in the \texttt{.bass} file, or left as its default value of
260   $2.5\sigma_{ii}$, where $\sigma_{ii}$ is the largest Lennard-Jones
261   length parameter present in the simulation. Truncating the calculation
262   at $r_{\text{cut}}$ introduces a discontinuity into the potential
263 < energy. To offset this discontinuity, the energy value at
264 < $r_{\text{cut}}$ is subtracted from the potential. This causes the
265 < potential to go to zero smoothly at the cut-off radius.
263 > energy and the force. To offset this discontinuity in the potential,
264 > the energy value at $r_{\text{cut}}$ is subtracted from the
265 > potential. This causes the potential to go to zero smoothly at the
266 > cut-off radius, and preserves conservation of energy in integrating
267 > the equations of motion.
268  
269   Interactions between dissimilar particles requires the generation of
270   cross term parameters for $\sigma$ and $\epsilon$. These are
# Line 296 | Line 295 | range interactions from $\frac{1}{r}$ to $\frac{1}{r^3
295   charges. Charge-neutral distributions were replaced with dipoles,
296   while most atoms and groups of atoms were reduced to Lennard-Jones
297   interaction sites. This simplification cuts the length scale of long
298 < range interactions from $\frac{1}{r}$ to $\frac{1}{r^3}$, allowing us
299 < to avoid the computationally expensive Ewald sum. Instead, we can use
300 < neighbor-lists, reaction field, and cutoff radii for the dipolar
301 < interactions.
298 > range interactions from $\frac{1}{r}$ to $\frac{1}{r^3}$, and allows
299 > us to avoid the computationally expensive Ewald sum. Instead, we can
300 > use neighbor-lists and cutoff radii for the dipolar interactions, or
301 > include a reaction field to mimic larger range interactions.
302  
303   As an example, lipid head-groups in {\sc duff} are represented as
304 < point dipole interaction sites. By placing a dipole of 20.6~Debye at
305 < the head group center of mass, our model mimics the head group of
306 < phosphatidylcholine.\cite{Cevc87} Additionally, a large Lennard-Jones
307 < site is located at the pseudoatom's center of mass. The model is
308 < illustrated by the dark grey atom in Fig.~\ref{oopseFig:lipidModel}. The
309 < water model we use to complement the dipoles of the lipids is our
310 < reparameterization of the soft sticky dipole (SSD) model of Ichiye
304 > point dipole interaction sites. By placing a dipole at the head group
305 > center of mass, our model mimics the charge separation found in common
306 > phospholipids such as phosphatidylcholine.\cite{Cevc87} Additionally,
307 > a large Lennard-Jones site is located at the pseudoatom's center of
308 > mass. The model is illustrated by the red atom in
309 > Fig.~\ref{oopseFig:lipidModel}. The water model we use to complement
310 > the dipoles of the lipids is our reparameterization of the soft sticky
311 > dipole (SSD) model of Ichiye
312   \emph{et al.}\cite{liu96:new_model}
313  
314   \begin{figure}
# Line 328 | Line 328 | of pseudoatoms to a minimum; the parameters for an ato
328   equilibria using Gibbs ensemble Monte Carlo simulation
329   techniques.\cite{Siepmann1998} One of the advantages of TraPPE is that
330   it generalizes the types of atoms in an alkyl chain to keep the number
331 < of pseudoatoms to a minimum; the parameters for an atom such as
331 > of pseudoatoms to a minimum; the parameters for a unified atom such as
332   $\text{CH}_2$ do not change depending on what species are bonded to
333   it.
334  
335   TraPPE also constrains all bonds to be of fixed length. Typically,
336   bond vibrations are the fastest motions in a molecular dynamic
337   simulation. Small time steps between force evaluations must be used to
338 < ensure adequate sampling of the bond potential to ensure conservation
339 < of energy. By constraining the bond lengths, larger time steps may be
340 < used when integrating the equations of motion. A simulation using {\sc
341 < duff} is illustrated in Scheme \ref{sch:DUFF}.
338 > ensure adequate energy conservation in the bond degrees of freedom. By
339 > constraining the bond lengths, larger time steps may be used when
340 > integrating the equations of motion. A simulation using {\sc duff} is
341 > illustrated in Scheme \ref{sch:DUFF}.
342  
343   \begin{lstlisting}[float,caption={[Invocation of {\sc duff}]Sample \texttt{.bass} file showing a simulation utilizing {\sc duff}},label={sch:DUFF}]
344  
# Line 407 | Line 407 | Here $\phi$ is the angle defined by four bonded neighb
407          + c_3[1 + \cos(3\phi)]
408   \label{eq:origTorsionPot}
409   \end{equation}
410 < Here $\phi$ is the angle defined by four bonded neighbors $i$,
411 < $j$, $k$, and $l$ (again, see Fig.~\ref{oopseFig:lipidModel}). For
412 < computational efficiency, the torsion potential has been recast after
413 < the method of {\sc charmm},\cite{Brooks83} in which the angle series is
414 < converted to a power series of the form:
410 > Where:
411   \begin{equation}
412 + \cos\phi = (\hat{\mathbf{r}}_{ij} \times \hat{\mathbf{r}}_{jk}) \cdot
413 +        (\hat{\mathbf{r}}_{jk} \times \hat{\mathbf{r}}_{kl})
414 + \label{eq:torsPhi}
415 + \end{equation}
416 + Here, $\hat{\mathbf{r}}_{\alpha\beta}$ are the set of unit bond
417 + vectors between atoms $i$, $j$, $k$, and $l$. For computational
418 + efficiency, the torsion potential has been recast after the method of
419 + {\sc charmm},\cite{Brooks83} in which the angle series is converted to
420 + a power series of the form:
421 + \begin{equation}
422   V_{\text{torsion}}(\phi) =  
423          k_3 \cos^3 \phi + k_2 \cos^2 \phi + k_1 \cos \phi + k_0
424   \label{eq:torsionPot}
# Line 452 | Line 458 | V_{\text{dipole}}(\mathbf{r}_{ij},\boldsymbol{\Omega}_
458          \boldsymbol{\Omega}_{j}) = \frac{|\mu_i||\mu_j|}{4\pi\epsilon_{0}r_{ij}^{3}} \biggl[
459          \boldsymbol{\hat{u}}_{i} \cdot \boldsymbol{\hat{u}}_{j}
460          -
461 <        \frac{3(\boldsymbol{\hat{u}}_i \cdot \mathbf{r}_{ij}) %
462 <                (\boldsymbol{\hat{u}}_j \cdot \mathbf{r}_{ij}) }
457 <                {r^{2}_{ij}} \biggr]
461 >        3(\boldsymbol{\hat{u}}_i \cdot \hat{\mathbf{r}}_{ij}) %
462 >                (\boldsymbol{\hat{u}}_j \cdot \hat{\mathbf{r}}_{ij}) \biggr]
463   \label{eq:dipolePot}
464   \end{equation}
465   Here $\mathbf{r}_{ij}$ is the vector starting at atom $i$ pointing
# Line 466 | Line 471 | unit vector pointing along $\mathbf{r}_{ij}$
471   unit vector pointing along $\mathbf{r}_{ij}$
472   ($\boldsymbol{\hat{r}}_{ij}=\mathbf{r}_{ij}/|\mathbf{r}_{ij}|$).
473  
474 + To improve computational efficiency of the dipole-dipole interactions,
475 + {\sc oopse} employs an electrostatic cutoff radius. This parameter can
476 + be set in the \texttt{.bass} file, and controls the length scale over
477 + which dipole interactions are felt. To compensate for the
478 + discontinuity in the potential and the forces at the cutoff radius, we
479 + have implemented a switching function to smoothly scale the
480 + dipole-dipole interaction at the cutoff.
481 + \begin{equation}
482 + S(r_{ij}) =
483 +        \begin{cases}
484 +        1 & \text{if $r_{ij} \le r_t$},\\
485 +        \frac{(r_{\text{cut}} + 2r_{ij} - 3r_t)(r_{\text{cut}} - r_{ij})^2}
486 +        {(r_{\text{cut}} - r_t)^2}
487 +        & \text{if $r_t < r_{ij} \le r_{\text{cut}}$}, \\
488 +        0 & \text{if $r_{ij} > r_{\text{cut}}$.}
489 +        \end{cases}
490 + \label{eq:dipoleSwitching}
491 + \end{equation}
492 + Here $S(r_{ij})$ scales the potential at a given $r_{ij}$, and $r_t$
493 + is the taper radius some given thickness less than the electrostatic
494 + cutoff. The switching thickness can be set in the \texttt{.bass} file.
495  
496 < \subsubsection{\label{oopseSec:SSD}The {\sc duff} Water Models: SSD/E and SSD/RF}
496 > \subsection{\label{oopseSec:SSD}The {\sc duff} Water Models: SSD/E and SSD/RF}
497  
498   In the interest of computational efficiency, the default solvent used
499   by {\sc oopse} is the extended Soft Sticky Dipole (SSD/E) water
# Line 527 | Line 553 | Since SSD is a single-point {\it dipolar} model, the f
553   can be found in the original SSD
554   articles.\cite{liu96:new_model,liu96:monte_carlo,chandra99:ssd_md,Ichiye03}
555  
556 < Since SSD is a single-point {\it dipolar} model, the force
556 > Since SSD/E is a single-point {\it dipolar} model, the force
557   calculations are simplified significantly relative to the standard
558   {\it charged} multi-point models. In the original Monte Carlo
559   simulations using this model, Ichiye {\it et al.} reported that using
560   SSD decreased computer time by a factor of 6-7 compared to other
561   models.\cite{liu96:new_model} What is most impressive is that these savings
562   did not come at the expense of accurate depiction of the liquid state
563 < properties.  Indeed, SSD maintains reasonable agreement with the Soper
563 > properties.  Indeed, SSD/E maintains reasonable agreement with the Head-Gordon
564   diffraction data for the structural features of liquid
565 < water.\cite{Soper86,liu96:new_model} Additionally, the dynamical properties
566 < exhibited by SSD agree with experiment better than those of more
565 > water.\cite{hura00,liu96:new_model} Additionally, the dynamical properties
566 > exhibited by SSD/E agree with experiment better than those of more
567   computationally expensive models (like TIP3P and
568   SPC/E).\cite{chandra99:ssd_md} The combination of speed and accurate depiction
569 < of solvent properties makes SSD a very attractive model for the
569 > of solvent properties makes SSD/E a very attractive model for the
570   simulation of large scale biochemical simulations.
571  
572   Recent constant pressure simulations revealed issues in the original
# Line 551 | Line 577 | model. Solvent parameters can be easily modified in an
577   of a reaction field long-range interaction correction is desired, it
578   is recommended that the parameters be modified to those of the SSD/RF
579   model. Solvent parameters can be easily modified in an accompanying
580 < {\sc BASS} file as illustrated in the scheme below. A table of the
580 > \texttt{.bass} file as illustrated in the scheme below. A table of the
581   parameter values and the drawbacks and benefits of the different
582   density corrected SSD models can be found in
583   reference~\cite{Gezelter04}.
# Line 571 | Line 597 | forceField = "DUFF";
597   forceField = "DUFF";
598  
599   /*
574 – * The reactionField flag toggles reaction
575 – * field corrections.
576 – */
577 –
578 – reactionField = false; // defaults to false
579 – dielectric = 80.0; // dielectric for reaction field
580 –
581 – /*
600   * The following two flags set the cutoff
601   * radius for the electrostatic forces
602   * as well as the skin thickness of the switching
# Line 597 | Line 615 | systems\cite{Finnis84,Ercolessi88,Chen90,Qi99,Ercoless
615   capacity to simulate metallic systems, including some that have
616   parallel computational abilities\cite{plimpton93}. Potentials that
617   describe bonding transition metal
618 < systems\cite{Finnis84,Ercolessi88,Chen90,Qi99,Ercolessi02} have a
618 > systems\cite{Finnis84,Ercolessi88,Chen90,Qi99,Ercolessi02} have an
619   attractive interaction which models  ``Embedding''
620   a positively charged metal ion in the electron density due to the
621   free valance ``sea'' of electrons created by the surrounding atoms in
622 < the system. A mostly repulsive pairwise part of the potential
622 > the system. A mostly-repulsive pairwise part of the potential
623   describes the interaction of the positively charged metal core ions
624   with one another. A particular potential description called the
625   Embedded Atom Method\cite{Daw84,FBD86,johnson89,Lu97}({\sc eam}) that has
626   particularly wide adoption has been selected for inclusion in {\sc oopse}. A
627 < good review of {\sc eam} and other metallic potential formulations was done
627 > good review of {\sc eam} and other metallic potential formulations was written
628   by Voter.\cite{voter}
629  
630   The {\sc eam} potential has the form:
# Line 614 | Line 632 | V & = & \sum_{i} F_{i}\left[\rho_{i}\right] + \sum_{i}
632   V & = & \sum_{i} F_{i}\left[\rho_{i}\right] + \sum_{i} \sum_{j \neq i}
633   \phi_{ij}({\bf r}_{ij})  \\
634   \rho_{i}  & = & \sum_{j \neq i} f_{j}({\bf r}_{ij})
635 < \end{eqnarray}S
618 <
635 > \end{eqnarray}
636   where $F_{i} $ is the embedding function that equates the energy required to embed a
637   positively-charged core ion $i$ into a linear superposition of
638   spherically averaged atomic electron densities given by
# Line 634 | Line 651 | interactions. Foiles et al. fit EAM potentials for fcc
651  
652   \newcommand{\roundme}{\operatorname{round}}
653  
654 < \textit{Periodic boundary conditions} are widely used to simulate truly
655 < macroscopic systems with a relatively small number of particles. The
656 < simulation box is replicated throughout space to form an infinite lattice.
657 < During the simulation, when a particle moves in the primary cell, its image in
658 < other boxes move in exactly the same direction with exactly the same
659 < orientation.Thus, as a particle leaves the primary cell, one of its images
660 < will enter through the opposite face.If the simulation box is large enough to
661 < avoid \textquotedblleft feeling\textquotedblright\ the symmetries of the
662 < periodic lattice, surface effects can be ignored. Cubic, orthorhombic and
663 < parallelepiped are the available periodic cells In OOPSE. We use a matrix to
664 < describe the property of the simulation box. Therefore, both the size and
648 < shape of the simulation box can be changed during the simulation. The
649 < transformation from box space vector $\mathbf{s}$ to its corresponding real
650 < space vector $\mathbf{r}$ is defined by
654 > \textit{Periodic boundary conditions} are widely used to simulate bulk properties with a relatively small number of particles. The
655 > simulation box is replicated throughout space to form an infinite
656 > lattice.  During the simulation, when a particle moves in the primary
657 > cell, its image in other cells move in exactly the same direction with
658 > exactly the same orientation. Thus, as a particle leaves the primary
659 > cell, one of its images will enter through the opposite face. If the
660 > simulation box is large enough to avoid ``feeling'' the symmetries of
661 > the periodic lattice, surface effects can be ignored. The available
662 > periodic cells in OOPSE are cubic, orthorhombic and parallelepiped. We
663 > use a $3 \times 3$ matrix, $\mathbf{H}$, to describe the shape and
664 > size of the simulation box. $\mathbf{H}$ is defined:
665   \begin{equation}
666 < \mathbf{r}=\underline{\mathbf{H}}\cdot\mathbf{s}%
666 > \mathbf{H} = ( \mathbf{h}_x, \mathbf{h}_y, \mathbf{h}_z )
667   \end{equation}
668 + Where $\mathbf{h}_j$ is the column vector of the $j$th axis of the
669 + box.  During the course of the simulation both the size and shape of
670 + the box can be changed to allow volume fluctations when constraining
671 + the pressure.
672  
673 <
674 < where $H=(h_{x},h_{y},h_{z})$ is a transformation matrix made up of the three
675 < box axis vectors. $h_{x},h_{y}$ and $h_{z}$ represent the three sides of the
676 < simulation box respectively.
677 <
678 < To find the minimum image of a vector $\mathbf{r}$, we convert the real vector
679 < to its corresponding vector in box space first, \bigskip%
673 > A real space vector, $\mathbf{r}$ can be transformed in to a box space
674 > vector, $\mathbf{s}$, and back through the following transformations:
675 > \begin{align}
676 > \mathbf{s} &= \mathbf{H}^{-1} \mathbf{r} \\
677 > \mathbf{r} &= \mathbf{H} \mathbf{s}
678 > \end{align}
679 > The vector $\mathbf{s}$ is now a vector expressed as the number of box
680 > lengths in the $\mathbf{h}_x$, $\mathbf{h}_y$, and $\mathbf{h}_z$
681 > directions. To find the minimum image of a vector $\mathbf{r}$, we
682 > first convert it to its corresponding vector in box space, and then,
683 > cast each element to lie on the in the range $[-0.5,0.5]$:
684   \begin{equation}
685 < \mathbf{s}=\underline{\mathbf{H}}^{-1}\cdot\mathbf{r}%
685 > s_{i}^{\prime}=s_{i}-\roundme(s_{i})
686   \end{equation}
687 < And then, each element of $\mathbf{s}$ is wrapped to lie between -0.5 to 0.5,
687 > Where $s_i$ is the $i$th element of $\mathbf{s}$, and
688 > $\roundme(s_i)$is given by
689   \begin{equation}
690 < s_{i}^{\prime}=s_{i}-\roundme(s_{i})
690 > \roundme(x) =
691 >        \begin{cases}
692 >        \lfloor x+0.5 \rfloor & \text{if $x \ge 0$} \\
693 >        \lceil x-0.5 \rceil & \text{if $x < 0$ }
694 >        \end{cases}
695   \end{equation}
696 < where
697 <
698 < %
699 <
700 < \begin{equation}
674 < \roundme(x)=\left\{
675 < \begin{array}{cc}%
676 < \lfloor{x+0.5}\rfloor & \text{if \ }x\geqslant 0 \\
677 < \lceil{x-0.5}\rceil & \text{otherwise}%
678 < \end{array}
679 < \right.
680 < \end{equation}
681 <
682 <
683 < For example, $\roundme(3.6)=4$,$\roundme(3.1)=3$, $\roundme(-3.6)=-4$, $\roundme(-3.1)=-3$.
696 > Here $\lfloor x \rfloor$ is the floor operator, and gives the largest
697 > integer value that is not greater than $x$, and $\lceil x \rceil$ is
698 > the ceiling operator, and gives the smallest integer that is not less
699 > than $x$.  For example, $\roundme(3.6)=4$, $\roundme(3.1)=3$,
700 > $\roundme(-3.6)=-4$, $\roundme(-3.1)=-3$.
701  
702   Finally, we obtain the minimum image coordinates $\mathbf{r}^{\prime}$ by
703 < transforming back to real space,%
687 <
703 > transforming back to real space,
704   \begin{equation}
705 < \mathbf{r}^{\prime}=\underline{\mathbf{H}}^{-1}\cdot\mathbf{s}^{\prime}%
705 > \mathbf{r}^{\prime}=\mathbf{H}^{-1}\mathbf{s}^{\prime}%
706   \end{equation}
707 + In this way, particles are allowed to diffuse freely in $\mathbf{r}$,
708 + but their minimum images, $\mathbf{r}^{\prime}$ are used to compute
709 + the interatomic forces.
710  
711  
712   \section{\label{oopseSec:IOfiles}Input and Output Files}

Diff Legend

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