| 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 signifigant, |
| 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 |
< |
{\sc give final breakdown of who wrote which section here.} |
| 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 |
> |
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 |
| 42 |
> |
talents in the development of the extended system integrators |
| 43 |
> |
(Sec.~\ref{oopseSec:noseHooverThermo}) as well as giving general |
| 44 |
> |
direction and oversight to the entire project. My responsibilities |
| 45 |
> |
covered the creation and specification of {\sc bass} |
| 46 |
> |
(Sec.~\ref{oopseSec:IOfiles}), the original development of the single |
| 47 |
> |
processor version of {\sc oopse}, contributions to the extended state |
| 48 |
> |
integrators (Sec.~\ref{oopseSec:noseHooverThermo}), the implementation |
| 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. {\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 |
|
|
| 61 |
|
might consider writing one's own programming code. However, as systems |
| 62 |
|
grow larger and more complex, building and maintaining code for the |
| 63 |
|
simulations becomes a time consuming task. In such cases it is usually |
| 64 |
< |
more convienent for a researcher to turn to pre-existing simulation |
| 64 |
> |
more convenient for a researcher to turn to pre-existing simulation |
| 65 |
|
packages. These packages, such as {\sc amber}\cite{pearlman:1995} and |
| 66 |
|
{\sc charmm}\cite{Brooks83}, provide powerful tools for researchers to |
| 67 |
|
conduct simulations of their systems without spending their time |
| 68 |
|
developing a code base to conduct their research. This then frees them |
| 69 |
< |
to perhaps explore experimental analouges to their models. |
| 69 |
> |
to perhaps explore experimental analogues to their models. |
| 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 emmbedded |
| 75 |
> |
interactions, rigid body dynamics, and metallic embedded |
| 76 |
|
potentials. When faced with these obstacles, a researcher must either |
| 77 |
|
develop their own code or license and extend one of the commercial |
| 78 |
|
packages. What we have elected to do, is develop a package of |
| 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 dcevelopment, and releaseing 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 emperical energy |
| 93 |
> |
We have structured this chapter to first discuss the empirical energy |
| 94 |
|
functions that {\sc oopse } implements in |
| 95 |
< |
Sec.~\ref{oopseSec:empericalEnergy}. Following that is a discusion of |
| 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 |
< |
mplements in the integration of the Newtonian equations of |
| 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, |
| 78 |
< |
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 |
|
|
| 112 |
|
atoms of an element, or be used for collections of atoms such as |
| 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 on atoms are not currently supported by {\sc oopse}. |
| 115 |
> |
dipoles). Charges, permanent dipoles, and Lennard-Jones parameters for |
| 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 the simple Ar molecule},label=sch:AtmMole] |
| 118 |
> |
\begin{lstlisting}[float,caption={[Specifier for molecules and atoms] A sample specification of an Ar molecule},label=sch:AtmMole] |
| 119 |
|
molecule{ |
| 120 |
|
name = "Ar"; |
| 121 |
|
nAtoms = 1; |
| 127 |
|
\end{lstlisting} |
| 128 |
|
|
| 129 |
|
|
| 130 |
< |
Atoms can be collected into secondary srtructures such as rigid bodies |
| 130 |
> |
Atoms can be collected into secondary structures such as rigid bodies |
| 131 |
|
or molecules. The molecule is a way for {\sc oopse} to keep track of |
| 132 |
|
the atoms in a simulation in logical manner. Molecular units store the |
| 133 |
< |
identities of all the atoms associated with themselves, and are |
| 134 |
< |
responsible for the evaluation of their own internal interactions |
| 135 |
< |
(\emph{i.e.}~bonds, bends, and torsions). Scheme \ref{sch:AtmMole} |
| 136 |
< |
shws how one creates a molecule in the \texttt{.mdl} files. The |
| 137 |
< |
position of the atoms given in the declaration are relative to the |
| 138 |
< |
origin of the molecule, and is used when creating a system containing |
| 139 |
< |
the molecule. |
| 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 ``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. |
| 140 |
|
|
| 141 |
|
As stated previously, one of the features that sets {\sc oopse} apart |
| 142 |
|
from most of the current molecular simulation packages is the ability |
| 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 |
| 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 |
| 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, |
| 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 |
| 189 |
< |
on this choice can be found in Sec.~\ref{sec:integrate}. An example |
| 190 |
< |
definition of a riged body can be seen in Scheme |
| 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 |
| 192 |
|
placements of the atoms relative to the origin of the rigid body, |
| 193 |
|
which itself has a position relative to the origin of the molecule. |
| 216 |
|
} |
| 217 |
|
\end{lstlisting} |
| 218 |
|
|
| 219 |
< |
\subsection{\label{sec:LJPot}The Lennard Jones Potential} |
| 219 |
> |
\subsection{\label{sec:LJPot}The Lennard Jones Force Field} |
| 220 |
|
|
| 221 |
|
The most basic force field implemented in {\sc oopse} is the |
| 222 |
< |
Lennard-Jones potential, which mimics the van der Waals interaction at |
| 222 |
> |
Lennard-Jones force field, which mimics the van der Waals interaction at |
| 223 |
|
long distances, and uses an empirical repulsion at short |
| 224 |
|
distances. The Lennard-Jones potential is given by: |
| 225 |
|
\begin{equation} |
| 233 |
|
Where $r_{ij}$ is the distance between particles $i$ and $j$, |
| 234 |
|
$\sigma_{ij}$ scales the length of the interaction, and |
| 235 |
|
$\epsilon_{ij}$ scales the well depth of the potential. Scheme |
| 236 |
< |
\ref{sch:LJFF} gives and example partial \texttt{.bass} file that |
| 237 |
< |
shows a system of 108 Ar particles simulated with the Lennard-Jones |
| 238 |
< |
force field. |
| 236 |
> |
\ref{sch:LJFF} gives and example \texttt{.bass} file that |
| 237 |
> |
sets up a system of 108 Ar particles to be simulated using the |
| 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}] |
| 241 |
|
|
| 214 |
– |
/* |
| 215 |
– |
* The Ar molecule is specified |
| 216 |
– |
* external to the.bass file |
| 217 |
– |
*/ |
| 218 |
– |
|
| 242 |
|
#include "argon.mdl" |
| 243 |
|
|
| 244 |
|
nComponents = 1; |
| 247 |
|
nMol = 108; |
| 248 |
|
} |
| 249 |
|
|
| 227 |
– |
/* |
| 228 |
– |
* The initial configuration is generated |
| 229 |
– |
* before the simulation is invoked. |
| 230 |
– |
*/ |
| 231 |
– |
|
| 250 |
|
initialConfig = "./argon.init"; |
| 251 |
|
|
| 252 |
|
forceField = "LJ"; |
| 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 |
| 280 |
|
\label{eq:epsilonMix} |
| 281 |
|
\end{equation} |
| 282 |
|
|
| 262 |
– |
|
| 263 |
– |
|
| 283 |
|
\subsection{\label{oopseSec:DUFF}Dipolar Unified-Atom Force Field} |
| 284 |
|
|
| 285 |
|
The dipolar unified-atom force field ({\sc duff}) was developed to |
| 286 |
|
simulate lipid bilayers. The simulations require a model capable of |
| 287 |
|
forming bilayers, while still being sufficiently computationally |
| 288 |
< |
efficient to allow large systems ($\approx$100's of phospholipids, |
| 289 |
< |
$\approx$1000's of waters) to be simulated for long times |
| 290 |
< |
($\approx$10's of nanoseconds). |
| 288 |
> |
efficient to allow large systems ($\sim$100's of phospholipids, |
| 289 |
> |
$\sim$1000's of waters) to be simulated for long times |
| 290 |
> |
($\sim$10's of nanoseconds). |
| 291 |
|
|
| 292 |
|
With this goal in mind, {\sc duff} has no point |
| 293 |
|
charges. Charge-neutral distributions were replaced with dipoles, |
| 294 |
|
while most atoms and groups of atoms were reduced to Lennard-Jones |
| 295 |
|
interaction sites. This simplification cuts the length scale of long |
| 296 |
< |
range interactions from $\frac{1}{r}$ to $\frac{1}{r^3}$, allowing us |
| 297 |
< |
to avoid the computationally expensive Ewald sum. Instead, we can use |
| 298 |
< |
neighbor-lists, reaction field, and cutoff radii for the dipolar |
| 299 |
< |
interactions. |
| 296 |
> |
range interactions from $\frac{1}{r}$ to $\frac{1}{r^3}$, and allows |
| 297 |
> |
us to avoid the computationally expensive Ewald sum. Instead, we can |
| 298 |
> |
use neighbor-lists and cutoff radii for the dipolar interactions, or |
| 299 |
> |
include a reaction field to mimic larger range interactions. |
| 300 |
|
|
| 301 |
|
As an example, lipid head-groups in {\sc duff} are represented as |
| 302 |
< |
point dipole interaction sites. By placing a dipole of 20.6~Debye at |
| 303 |
< |
the head group center of mass, our model mimics the head group of |
| 304 |
< |
phosphatidylcholine.\cite{Cevc87} Additionally, a large Lennard-Jones |
| 305 |
< |
site is located at the pseudoatom's center of mass. The model is |
| 306 |
< |
illustrated by the dark grey atom in Fig.~\ref{fig:lipidModel}. The |
| 307 |
< |
water model we use to complement the dipoles of the lipids is our |
| 308 |
< |
reparameterization of the soft sticky dipole (SSD) model of Ichiye |
| 302 |
> |
point dipole interaction sites. By placing a dipole at the head group |
| 303 |
> |
center of mass, our model mimics the charge separation found in common |
| 304 |
> |
phospholipids such as phosphatidylcholine.\cite{Cevc87} Additionally, |
| 305 |
> |
a large Lennard-Jones site is located at the pseudoatom's center of |
| 306 |
> |
mass. The model is illustrated by the red atom in |
| 307 |
> |
Fig.~\ref{oopseFig:lipidModel}. The water model we use to complement |
| 308 |
> |
the dipoles of the lipids is our reparameterization of the soft sticky |
| 309 |
> |
dipole (SSD) model of Ichiye |
| 310 |
|
\emph{et al.}\cite{liu96:new_model} |
| 311 |
|
|
| 312 |
|
\begin{figure} |
| 326 |
|
equilibria using Gibbs ensemble Monte Carlo simulation |
| 327 |
|
techniques.\cite{Siepmann1998} One of the advantages of TraPPE is that |
| 328 |
|
it generalizes the types of atoms in an alkyl chain to keep the number |
| 329 |
< |
of pseudoatoms to a minimum; the parameters for an atom such as |
| 329 |
> |
of pseudoatoms to a minimum; the parameters for a unified atom such as |
| 330 |
|
$\text{CH}_2$ do not change depending on what species are bonded to |
| 331 |
|
it. |
| 332 |
|
|
| 333 |
|
TraPPE also constrains all bonds to be of fixed length. Typically, |
| 334 |
|
bond vibrations are the fastest motions in a molecular dynamic |
| 335 |
|
simulation. Small time steps between force evaluations must be used to |
| 336 |
< |
ensure adequate sampling of the bond potential to ensure conservation |
| 337 |
< |
of energy. By constraining the bond lengths, larger time steps may be |
| 338 |
< |
used when integrating the equations of motion. A simulation using {\sc |
| 339 |
< |
duff} is illustrated in Scheme \ref{sch:DUFF}. |
| 336 |
> |
ensure adequate energy conservation in the bond degrees of freedom. By |
| 337 |
> |
constraining the bond lengths, larger time steps may be used when |
| 338 |
> |
integrating the equations of motion. A simulation using {\sc duff} is |
| 339 |
> |
illustrated in Scheme \ref{sch:DUFF}. |
| 340 |
|
|
| 341 |
|
\begin{lstlisting}[float,caption={[Invocation of {\sc duff}]Sample \texttt{.bass} file showing a simulation utilizing {\sc duff}},label={sch:DUFF}] |
| 342 |
|
|
| 365 |
|
The total potential energy function in {\sc duff} is |
| 366 |
|
\begin{equation} |
| 367 |
|
V = \sum^{N}_{I=1} V^{I}_{\text{Internal}} |
| 368 |
< |
+ \sum^{N}_{I=1} \sum_{J>I} V^{IJ}_{\text{Cross}} |
| 368 |
> |
+ \sum^{N-1}_{I=1} \sum_{J>I} V^{IJ}_{\text{Cross}} |
| 369 |
|
\label{eq:totalPotential} |
| 370 |
|
\end{equation} |
| 371 |
|
Where $V^{I}_{\text{Internal}}$ is the internal potential of molecule $I$: |
| 392 |
|
V_{\text{bend}}(\theta_{ijk}) = k_{\theta}( \theta_{ijk} - \theta_0 )^2 \label{eq:bendPot} |
| 393 |
|
\end{equation} |
| 394 |
|
Where $\theta_{ijk}$ is the angle defined by atoms $i$, $j$, and $k$ |
| 395 |
< |
(see Fig.~\ref{fig:lipidModel}), $\theta_0$ is the equilibrium |
| 395 |
> |
(see Fig.~\ref{oopseFig:lipidModel}), $\theta_0$ is the equilibrium |
| 396 |
|
bond angle, and $k_{\theta}$ is the force constant which determines the |
| 397 |
|
strength of the harmonic bend. The parameters for $k_{\theta}$ and |
| 398 |
|
$\theta_0$ are borrowed from those in TraPPE.\cite{Siepmann1998} |
| 405 |
|
+ c_3[1 + \cos(3\phi)] |
| 406 |
|
\label{eq:origTorsionPot} |
| 407 |
|
\end{equation} |
| 408 |
< |
Here $\phi$ is the angle defined by four bonded neighbors $i$, |
| 389 |
< |
$j$, $k$, and $l$ (again, see Fig.~\ref{fig:lipidModel}). For |
| 390 |
< |
computational efficiency, the torsion potential has been recast after |
| 391 |
< |
the method of CHARMM,\cite{charmm1983} in which the angle series is |
| 392 |
< |
converted to a power series of the form: |
| 408 |
> |
Where: |
| 409 |
|
\begin{equation} |
| 410 |
+ |
\cos\phi = (\hat{\mathbf{r}}_{ij} \times \hat{\mathbf{r}}_{jk}) \cdot |
| 411 |
+ |
(\hat{\mathbf{r}}_{jk} \times \hat{\mathbf{r}}_{kl}) |
| 412 |
+ |
\label{eq:torsPhi} |
| 413 |
+ |
\end{equation} |
| 414 |
+ |
Here, $\hat{\mathbf{r}}_{\alpha\beta}$ are the set of unit bond |
| 415 |
+ |
vectors between atoms $i$, $j$, $k$, and $l$. For computational |
| 416 |
+ |
efficiency, the torsion potential has been recast after the method of |
| 417 |
+ |
{\sc charmm},\cite{Brooks83} in which the angle series is converted to |
| 418 |
+ |
a power series of the form: |
| 419 |
+ |
\begin{equation} |
| 420 |
|
V_{\text{torsion}}(\phi) = |
| 421 |
|
k_3 \cos^3 \phi + k_2 \cos^2 \phi + k_1 \cos \phi + k_0 |
| 422 |
|
\label{eq:torsionPot} |
| 447 |
|
Where $V_{\text{LJ}}$ is the Lennard Jones potential, |
| 448 |
|
$V_{\text{dipole}}$ is the dipole dipole potential, and |
| 449 |
|
$V_{\text{sticky}}$ is the sticky potential defined by the SSD model |
| 450 |
< |
(Sec.~\ref{sec:SSD}). Note that not all atom types include all |
| 450 |
> |
(Sec.~\ref{oopseSec:SSD}). Note that not all atom types include all |
| 451 |
|
interactions. |
| 452 |
|
|
| 453 |
|
The dipole-dipole potential has the following form: |
| 456 |
|
\boldsymbol{\Omega}_{j}) = \frac{|\mu_i||\mu_j|}{4\pi\epsilon_{0}r_{ij}^{3}} \biggl[ |
| 457 |
|
\boldsymbol{\hat{u}}_{i} \cdot \boldsymbol{\hat{u}}_{j} |
| 458 |
|
- |
| 459 |
< |
\frac{3(\boldsymbol{\hat{u}}_i \cdot \mathbf{r}_{ij}) % |
| 460 |
< |
(\boldsymbol{\hat{u}}_j \cdot \mathbf{r}_{ij}) } |
| 435 |
< |
{r^{2}_{ij}} \biggr] |
| 459 |
> |
3(\boldsymbol{\hat{u}}_i \cdot \hat{\mathbf{r}}_{ij}) % |
| 460 |
> |
(\boldsymbol{\hat{u}}_j \cdot \hat{\mathbf{r}}_{ij}) \biggr] |
| 461 |
|
\label{eq:dipolePot} |
| 462 |
|
\end{equation} |
| 463 |
|
Here $\mathbf{r}_{ij}$ is the vector starting at atom $i$ pointing |
| 464 |
|
towards $j$, and $\boldsymbol{\Omega}_i$ and $\boldsymbol{\Omega}_j$ |
| 465 |
|
are the orientational degrees of freedom for atoms $i$ and $j$ |
| 466 |
|
respectively. $|\mu_i|$ is the magnitude of the dipole moment of atom |
| 467 |
< |
$i$, $\boldsymbol{\hat{u}}_i$ is the standard unit orientation |
| 468 |
< |
vector of $\boldsymbol{\Omega}_i$, and $\boldsymbol{\hat{r}}_{ij}$ is |
| 469 |
< |
the unit vector pointing along $\mathbf{r}_{ij}$. |
| 467 |
> |
$i$, $\boldsymbol{\hat{u}}_i$ is the standard unit orientation vector |
| 468 |
> |
of $\boldsymbol{\Omega}_i$, and $\boldsymbol{\hat{r}}_{ij}$ is the |
| 469 |
> |
unit vector pointing along $\mathbf{r}_{ij}$ |
| 470 |
> |
($\boldsymbol{\hat{r}}_{ij}=\mathbf{r}_{ij}/|\mathbf{r}_{ij}|$). |
| 471 |
|
|
| 472 |
+ |
To improve computational efficiency of the dipole-dipole interactions, |
| 473 |
+ |
{\sc oopse} employs an electrostatic cutoff radius. This parameter can |
| 474 |
+ |
be set in the \texttt{.bass} file, and controls the length scale over |
| 475 |
+ |
which dipole interactions are felt. To compensate for the |
| 476 |
+ |
discontinuity in the potential and the forces at the cutoff radius, we |
| 477 |
+ |
have implemented a switching function to smoothly scale the |
| 478 |
+ |
dipole-dipole interaction at the cutoff. |
| 479 |
+ |
\begin{equation} |
| 480 |
+ |
S(r_{ij}) = |
| 481 |
+ |
\begin{cases} |
| 482 |
+ |
1 & \text{if $r_{ij} \le r_t$},\\ |
| 483 |
+ |
\frac{(r_{\text{cut}} + 2r_{ij} - 3r_t)(r_{\text{cut}} - r_{ij})^2} |
| 484 |
+ |
{(r_{\text{cut}} - r_t)^2} |
| 485 |
+ |
& \text{if $r_t < r_{ij} \le r_{\text{cut}}$}, \\ |
| 486 |
+ |
0 & \text{if $r_{ij} > r_{\text{cut}}$.} |
| 487 |
+ |
\end{cases} |
| 488 |
+ |
\label{eq:dipoleSwitching} |
| 489 |
+ |
\end{equation} |
| 490 |
+ |
Here $S(r_{ij})$ scales the potential at a given $r_{ij}$, and $r_t$ |
| 491 |
+ |
is the taper radius some given thickness less than the electrostatic |
| 492 |
+ |
cutoff. The switching thickness can be set in the \texttt{.bass} file. |
| 493 |
|
|
| 494 |
< |
\subsubsection{\label{sec:SSD}The {\sc duff} Water Models: SSD/E and SSD/RF} |
| 494 |
> |
\subsection{\label{oopseSec:SSD}The {\sc duff} Water Models: SSD/E and SSD/RF} |
| 495 |
|
|
| 496 |
|
In the interest of computational efficiency, the default solvent used |
| 497 |
|
by {\sc oopse} is the extended Soft Sticky Dipole (SSD/E) water |
| 551 |
|
can be found in the original SSD |
| 552 |
|
articles.\cite{liu96:new_model,liu96:monte_carlo,chandra99:ssd_md,Ichiye03} |
| 553 |
|
|
| 554 |
< |
Since SSD is a single-point {\it dipolar} model, the force |
| 554 |
> |
Since SSD/E is a single-point {\it dipolar} model, the force |
| 555 |
|
calculations are simplified significantly relative to the standard |
| 556 |
|
{\it charged} multi-point models. In the original Monte Carlo |
| 557 |
|
simulations using this model, Ichiye {\it et al.} reported that using |
| 558 |
|
SSD decreased computer time by a factor of 6-7 compared to other |
| 559 |
|
models.\cite{liu96:new_model} What is most impressive is that these savings |
| 560 |
|
did not come at the expense of accurate depiction of the liquid state |
| 561 |
< |
properties. Indeed, SSD maintains reasonable agreement with the Soper |
| 561 |
> |
properties. Indeed, SSD/E maintains reasonable agreement with the Head-Gordon |
| 562 |
|
diffraction data for the structural features of liquid |
| 563 |
< |
water.\cite{Soper86,liu96:new_model} Additionally, the dynamical properties |
| 564 |
< |
exhibited by SSD agree with experiment better than those of more |
| 563 |
> |
water.\cite{hura00,liu96:new_model} Additionally, the dynamical properties |
| 564 |
> |
exhibited by SSD/E agree with experiment better than those of more |
| 565 |
|
computationally expensive models (like TIP3P and |
| 566 |
|
SPC/E).\cite{chandra99:ssd_md} The combination of speed and accurate depiction |
| 567 |
< |
of solvent properties makes SSD a very attractive model for the |
| 567 |
> |
of solvent properties makes SSD/E a very attractive model for the |
| 568 |
|
simulation of large scale biochemical simulations. |
| 569 |
|
|
| 570 |
|
Recent constant pressure simulations revealed issues in the original |
| 575 |
|
of a reaction field long-range interaction correction is desired, it |
| 576 |
|
is recommended that the parameters be modified to those of the SSD/RF |
| 577 |
|
model. Solvent parameters can be easily modified in an accompanying |
| 578 |
< |
{\sc BASS} file as illustrated in the scheme below. A table of the |
| 578 |
> |
\texttt{.bass} file as illustrated in the scheme below. A table of the |
| 579 |
|
parameter values and the drawbacks and benefits of the different |
| 580 |
< |
density corrected SSD models can be found in reference |
| 581 |
< |
\ref{Gezelter04}. |
| 580 |
> |
density corrected SSD models can be found in |
| 581 |
> |
reference~\cite{Gezelter04}. |
| 582 |
|
|
| 583 |
|
\begin{lstlisting}[float,caption={[A simulation of {\sc ssd} water]An example file showing a simulation including {\sc ssd} water.},label={sch:ssd}] |
| 584 |
|
|
| 595 |
|
forceField = "DUFF"; |
| 596 |
|
|
| 597 |
|
/* |
| 551 |
– |
* The reactionField flag toggles reaction |
| 552 |
– |
* field corrections. |
| 553 |
– |
*/ |
| 554 |
– |
|
| 555 |
– |
reactionField = false; // defaults to false |
| 556 |
– |
dielectric = 80.0; // dielectric for reaction field |
| 557 |
– |
|
| 558 |
– |
/* |
| 598 |
|
* The following two flags set the cutoff |
| 599 |
|
* radius for the electrostatic forces |
| 600 |
|
* as well as the skin thickness of the switching |
| 609 |
|
|
| 610 |
|
\subsection{\label{oopseSec:eam}Embedded Atom Method} |
| 611 |
|
|
| 612 |
< |
Several other molecular dynamics packages\cite{dynamo86} exist which have the |
| 612 |
> |
There are Molecular Dynamics packages which have the |
| 613 |
|
capacity to simulate metallic systems, including some that have |
| 614 |
|
parallel computational abilities\cite{plimpton93}. Potentials that |
| 615 |
|
describe bonding transition metal |
| 616 |
< |
systems\cite{Finnis84,Ercolessi88,Chen90,Qi99,Ercolessi02} have a |
| 616 |
> |
systems\cite{Finnis84,Ercolessi88,Chen90,Qi99,Ercolessi02} have an |
| 617 |
|
attractive interaction which models ``Embedding'' |
| 618 |
|
a positively charged metal ion in the electron density due to the |
| 619 |
|
free valance ``sea'' of electrons created by the surrounding atoms in |
| 620 |
< |
the system. A mostly repulsive pairwise part of the potential |
| 620 |
> |
the system. A mostly-repulsive pairwise part of the potential |
| 621 |
|
describes the interaction of the positively charged metal core ions |
| 622 |
|
with one another. A particular potential description called the |
| 623 |
|
Embedded Atom Method\cite{Daw84,FBD86,johnson89,Lu97}({\sc eam}) that has |
| 624 |
|
particularly wide adoption has been selected for inclusion in {\sc oopse}. A |
| 625 |
< |
good review of {\sc eam} and other metallic potential formulations was done |
| 625 |
> |
good review of {\sc eam} and other metallic potential formulations was written |
| 626 |
|
by Voter.\cite{voter} |
| 627 |
|
|
| 628 |
|
The {\sc eam} potential has the form: |
| 630 |
|
V & = & \sum_{i} F_{i}\left[\rho_{i}\right] + \sum_{i} \sum_{j \neq i} |
| 631 |
|
\phi_{ij}({\bf r}_{ij}) \\ |
| 632 |
|
\rho_{i} & = & \sum_{j \neq i} f_{j}({\bf r}_{ij}) |
| 633 |
< |
\end{eqnarray}S |
| 595 |
< |
|
| 633 |
> |
\end{eqnarray} |
| 634 |
|
where $F_{i} $ is the embedding function that equates the energy required to embed a |
| 635 |
|
positively-charged core ion $i$ into a linear superposition of |
| 636 |
|
spherically averaged atomic electron densities given by |
| 637 |
|
$\rho_{i}$. $\phi_{ij}$ is a primarily repulsive pairwise interaction |
| 638 |
|
between atoms $i$ and $j$. In the original formulation of |
| 639 |
< |
{\sc eam} cite{Daw84}, $\phi_{ij}$ was an entirely repulsive term, however |
| 639 |
> |
{\sc eam}\cite{Daw84}, $\phi_{ij}$ was an entirely repulsive term, however |
| 640 |
|
in later refinements to EAM have shown that non-uniqueness between $F$ |
| 641 |
|
and $\phi$ allow for more general forms for $\phi$.\cite{Daw89} |
| 642 |
|
There is a cutoff distance, $r_{cut}$, which limits the |
| 643 |
|
summations in the {\sc eam} equation to the few dozen atoms |
| 644 |
|
surrounding atom $i$ for both the density $\rho$ and pairwise $\phi$ |
| 645 |
< |
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}. |
| 645 |
> |
interactions. Foiles et al. fit EAM potentials for fcc metals Cu, Ag, Au, Ni, Pd, Pt and alloys of these metals\cite{FBD86}. These potential fits are in the DYNAMO 86 format and are included with {\sc oopse}. |
| 646 |
|
|
| 647 |
|
|
| 648 |
|
\subsection{\label{oopseSec:pbc}Periodic Boundary Conditions} |
| 649 |
|
|
| 650 |
|
\newcommand{\roundme}{\operatorname{round}} |
| 651 |
|
|
| 652 |
< |
\textit{Periodic boundary conditions} are widely used to simulate truly |
| 653 |
< |
macroscopic systems with a relatively small number of particles. The |
| 654 |
< |
simulation box is replicated throughout space to form an infinite lattice. |
| 655 |
< |
During the simulation, when a particle moves in the primary cell, its image in |
| 656 |
< |
other boxes move in exactly the same direction with exactly the same |
| 657 |
< |
orientation.Thus, as a particle leaves the primary cell, one of its images |
| 658 |
< |
will enter through the opposite face.If the simulation box is large enough to |
| 659 |
< |
avoid \textquotedblleft feeling\textquotedblright\ the symmetries of the |
| 660 |
< |
periodic lattice, surface effects can be ignored. Cubic, orthorhombic and |
| 661 |
< |
parallelepiped are the available periodic cells In OOPSE. We use a matrix to |
| 662 |
< |
describe the property of the simulation box. Therefore, both the size and |
| 625 |
< |
shape of the simulation box can be changed during the simulation. The |
| 626 |
< |
transformation from box space vector $\mathbf{s}$ to its corresponding real |
| 627 |
< |
space vector $\mathbf{r}$ is defined by |
| 652 |
> |
\textit{Periodic boundary conditions} are widely used to simulate bulk properties with a relatively small number of particles. The |
| 653 |
> |
simulation box is replicated throughout space to form an infinite |
| 654 |
> |
lattice. During the simulation, when a particle moves in the primary |
| 655 |
> |
cell, its image in other cells move in exactly the same direction with |
| 656 |
> |
exactly the same orientation. Thus, as a particle leaves the primary |
| 657 |
> |
cell, one of its images will enter through the opposite face. If the |
| 658 |
> |
simulation box is large enough to avoid ``feeling'' the symmetries of |
| 659 |
> |
the periodic lattice, surface effects can be ignored. The available |
| 660 |
> |
periodic cells in OOPSE are cubic, orthorhombic and parallelepiped. We |
| 661 |
> |
use a $3 \times 3$ matrix, $\mathbf{H}$, to describe the shape and |
| 662 |
> |
size of the simulation box. $\mathbf{H}$ is defined: |
| 663 |
|
\begin{equation} |
| 664 |
< |
\mathbf{r}=\underline{\mathbf{H}}\cdot\mathbf{s}% |
| 664 |
> |
\mathbf{H} = ( \mathbf{h}_x, \mathbf{h}_y, \mathbf{h}_z ) |
| 665 |
|
\end{equation} |
| 666 |
+ |
Where $\mathbf{h}_j$ is the column vector of the $j$th axis of the |
| 667 |
+ |
box. During the course of the simulation both the size and shape of |
| 668 |
+ |
the box can be changed to allow volume fluctations when constraining |
| 669 |
+ |
the pressure. |
| 670 |
|
|
| 671 |
< |
|
| 672 |
< |
where $H=(h_{x},h_{y},h_{z})$ is a transformation matrix made up of the three |
| 673 |
< |
box axis vectors. $h_{x},h_{y}$ and $h_{z}$ represent the three sides of the |
| 674 |
< |
simulation box respectively. |
| 675 |
< |
|
| 676 |
< |
To find the minimum image of a vector $\mathbf{r}$, we convert the real vector |
| 677 |
< |
to its corresponding vector in box space first, \bigskip% |
| 671 |
> |
A real space vector, $\mathbf{r}$ can be transformed in to a box space |
| 672 |
> |
vector, $\mathbf{s}$, and back through the following transformations: |
| 673 |
> |
\begin{align} |
| 674 |
> |
\mathbf{s} &= \mathbf{H}^{-1} \mathbf{r} \\ |
| 675 |
> |
\mathbf{r} &= \mathbf{H} \mathbf{s} |
| 676 |
> |
\end{align} |
| 677 |
> |
The vector $\mathbf{s}$ is now a vector expressed as the number of box |
| 678 |
> |
lengths in the $\mathbf{h}_x$, $\mathbf{h}_y$, and $\mathbf{h}_z$ |
| 679 |
> |
directions. To find the minimum image of a vector $\mathbf{r}$, we |
| 680 |
> |
first convert it to its corresponding vector in box space, and then, |
| 681 |
> |
cast each element to lie on the in the range $[-0.5,0.5]$: |
| 682 |
|
\begin{equation} |
| 640 |
– |
\mathbf{s}=\underline{\mathbf{H}}^{-1}\cdot\mathbf{r}% |
| 641 |
– |
\end{equation} |
| 642 |
– |
And then, each element of $\mathbf{s}$ is wrapped to lie between -0.5 to 0.5, |
| 643 |
– |
\begin{equation} |
| 683 |
|
s_{i}^{\prime}=s_{i}-\roundme(s_{i}) |
| 684 |
|
\end{equation} |
| 685 |
< |
where |
| 686 |
< |
|
| 648 |
< |
% |
| 649 |
< |
|
| 685 |
> |
Where $s_i$ is the $i$th element of $\mathbf{s}$, and |
| 686 |
> |
$\roundme(s_i)$is given by |
| 687 |
|
\begin{equation} |
| 688 |
< |
\roundme(x)=\left\{ |
| 689 |
< |
\begin{array}{cc}% |
| 690 |
< |
\lfloor{x+0.5}\rfloor & \text{if \ }x\geqslant 0 \\ |
| 691 |
< |
\lceil{x-0.5}\rceil & \text{otherwise}% |
| 692 |
< |
\end{array} |
| 656 |
< |
\right. |
| 688 |
> |
\roundme(x) = |
| 689 |
> |
\begin{cases} |
| 690 |
> |
\lfloor x+0.5 \rfloor & \text{if $x \ge 0$} \\ |
| 691 |
> |
\lceil x-0.5 \rceil & \text{if $x < 0$ } |
| 692 |
> |
\end{cases} |
| 693 |
|
\end{equation} |
| 694 |
+ |
Here $\lfloor x \rfloor$ is the floor operator, and gives the largest |
| 695 |
+ |
integer value that is not greater than $x$, and $\lceil x \rceil$ is |
| 696 |
+ |
the ceiling operator, and gives the smallest integer that is not less |
| 697 |
+ |
than $x$. For example, $\roundme(3.6)=4$, $\roundme(3.1)=3$, |
| 698 |
+ |
$\roundme(-3.6)=-4$, $\roundme(-3.1)=-3$. |
| 699 |
|
|
| 659 |
– |
|
| 660 |
– |
For example, $\roundme(3.6)=4$,$\roundme(3.1)=3$, $\roundme(-3.6)=-4$, $\roundme(-3.1)=-3$. |
| 661 |
– |
|
| 700 |
|
Finally, we obtain the minimum image coordinates $\mathbf{r}^{\prime}$ by |
| 701 |
< |
transforming back to real space,% |
| 664 |
< |
|
| 701 |
> |
transforming back to real space, |
| 702 |
|
\begin{equation} |
| 703 |
< |
\mathbf{r}^{\prime}=\underline{\mathbf{H}}^{-1}\cdot\mathbf{s}^{\prime}% |
| 703 |
> |
\mathbf{r}^{\prime}=\mathbf{H}^{-1}\mathbf{s}^{\prime}% |
| 704 |
|
\end{equation} |
| 705 |
+ |
In this way, particles are allowed to diffuse freely in $\mathbf{r}$, |
| 706 |
+ |
but their minimum images, $\mathbf{r}^{\prime}$ are used to compute |
| 707 |
+ |
the interatomic forces. |
| 708 |
|
|
| 709 |
|
|
| 710 |
|
\section{\label{oopseSec:IOfiles}Input and Output Files} |
| 711 |
|
|
| 712 |
|
\subsection{{\sc bass} and Model Files} |
| 713 |
|
|
| 714 |
< |
Every {\sc oopse} simuation begins with a {\sc bass} file. {\sc bass} |
| 715 |
< |
(\underline{B}izarre \underline{A}tom \underline{S}imulation |
| 716 |
< |
\underline{S}yntax) is a script syntax that is parsed by {\sc oopse} at |
| 717 |
< |
runtime. The {\sc bass} file allows for the user to completely describe the |
| 718 |
< |
system they are to simulate, as well as tailor {\sc oopse}'s behavior during |
| 719 |
< |
the simulation. {\sc bass} files are denoted with the extension |
| 714 |
> |
Every {\sc oopse} simulation begins with a Bizarre Atom Simulation |
| 715 |
> |
Syntax ({\sc bass}) file. {\sc bass} is a script syntax that is parsed |
| 716 |
> |
by {\sc oopse} at runtime. The {\sc bass} file allows for the user to |
| 717 |
> |
completely describe the system they wish to simulate, as well as tailor |
| 718 |
> |
{\sc oopse}'s behavior during the simulation. {\sc bass} files are |
| 719 |
> |
denoted with the extension |
| 720 |
|
\texttt{.bass}, an example file is shown in |
| 721 |
< |
Fig.~\ref{fig:bassExample}. |
| 721 |
> |
Scheme~\ref{sch:bassExample}. |
| 722 |
|
|
| 723 |
< |
\begin{figure} |
| 684 |
< |
\centering |
| 685 |
< |
\framebox[\linewidth]{\rule{0cm}{0.75\linewidth}I'm a {\sc bass} file!} |
| 686 |
< |
\caption{Here is an example \texttt{.bass} file} |
| 687 |
< |
\label{fig:bassExample} |
| 688 |
< |
\end{figure} |
| 723 |
> |
\begin{lstlisting}[float,caption={[An example of a complete {\sc bass} file] An example showing a complete {\sc bass} file.},label={sch:bassExample}] |
| 724 |
|
|
| 725 |
< |
Within the \texttt{.bass} file it is neccassary to provide a complete |
| 725 |
> |
molecule{ |
| 726 |
> |
name = "Ar"; |
| 727 |
> |
nAtoms = 1; |
| 728 |
> |
atom[0]{ |
| 729 |
> |
type="Ar"; |
| 730 |
> |
position( 0.0, 0.0, 0.0 ); |
| 731 |
> |
} |
| 732 |
> |
} |
| 733 |
> |
|
| 734 |
> |
nComponents = 1; |
| 735 |
> |
component{ |
| 736 |
> |
type = "Ar"; |
| 737 |
> |
nMol = 108; |
| 738 |
> |
} |
| 739 |
> |
|
| 740 |
> |
initialConfig = "./argon.init"; |
| 741 |
> |
|
| 742 |
> |
forceField = "LJ"; |
| 743 |
> |
ensemble = "NVE"; // specify the simulation enesemble |
| 744 |
> |
dt = 1.0; // the time step for integration |
| 745 |
> |
runTime = 1e3; // the total simulation run time |
| 746 |
> |
sampleTime = 100; // trajectory file frequency |
| 747 |
> |
statusTime = 50; // statistics file frequency |
| 748 |
> |
|
| 749 |
> |
\end{lstlisting} |
| 750 |
> |
|
| 751 |
> |
Within the \texttt{.bass} file it is necessary to provide a complete |
| 752 |
|
description of the molecule before it is actually placed in the |
| 753 |
< |
simulation. The {\sc bass} syntax was originally developed with this goal in |
| 754 |
< |
mind, and allows for the specification of all the atoms in a molecular |
| 755 |
< |
prototype, as well as any bonds, bends, or torsions. These |
| 753 |
> |
simulation. The {\sc bass} syntax was originally developed with this |
| 754 |
> |
goal in mind, and allows for the specification of all the atoms in a |
| 755 |
> |
molecular prototype, as well as any bonds, bends, or torsions. These |
| 756 |
|
descriptions can become lengthy for complex molecules, and it would be |
| 757 |
< |
inconvient to duplicate the simulation at the begining of each {\sc bass} |
| 758 |
< |
script. Addressing this issue {\sc bass} allows for the inclusion of model |
| 759 |
< |
files at the top of a \texttt{.bass} file. These model files, denoted |
| 760 |
< |
with the \texttt{.mdl} extension, allow the user to describe a |
| 761 |
< |
molecular prototype once, then simply include it into each simulation |
| 762 |
< |
containing that molecule. |
| 757 |
> |
inconvenient to duplicate the simulation at the beginning of each {\sc |
| 758 |
> |
bass} script. Addressing this issue {\sc bass} allows for the |
| 759 |
> |
inclusion of model files at the top of a \texttt{.bass} file. These |
| 760 |
> |
model files, denoted with the \texttt{.mdl} extension, allow the user |
| 761 |
> |
to describe a molecular prototype once, then simply include it into |
| 762 |
> |
each simulation containing that molecule. Returning to the example in |
| 763 |
> |
Scheme~\ref{sch:bassExample}, the \texttt{.mdl} file's contents would |
| 764 |
> |
be Scheme~\ref{sch:mdlExample}, and the new \texttt{.bass} file would |
| 765 |
> |
become Scheme~\ref{sch:bassExPrime}. |
| 766 |
|
|
| 767 |
+ |
\begin{lstlisting}[float,caption={An example \texttt{.mdl} file.},label={sch:mdlExample}] |
| 768 |
+ |
|
| 769 |
+ |
molecule{ |
| 770 |
+ |
name = "Ar"; |
| 771 |
+ |
nAtoms = 1; |
| 772 |
+ |
atom[0]{ |
| 773 |
+ |
type="Ar"; |
| 774 |
+ |
position( 0.0, 0.0, 0.0 ); |
| 775 |
+ |
} |
| 776 |
+ |
} |
| 777 |
+ |
|
| 778 |
+ |
\end{lstlisting} |
| 779 |
+ |
|
| 780 |
+ |
\begin{lstlisting}[float,caption={Revised {\sc bass} example.},label={sch:bassExPrime}] |
| 781 |
+ |
|
| 782 |
+ |
#include "argon.mdl" |
| 783 |
+ |
|
| 784 |
+ |
molecule{ |
| 785 |
+ |
name = "Ar"; |
| 786 |
+ |
nAtoms = 1; |
| 787 |
+ |
atom[0]{ |
| 788 |
+ |
type="Ar"; |
| 789 |
+ |
position( 0.0, 0.0, 0.0 ); |
| 790 |
+ |
} |
| 791 |
+ |
} |
| 792 |
+ |
|
| 793 |
+ |
nComponents = 1; |
| 794 |
+ |
component{ |
| 795 |
+ |
type = "Ar"; |
| 796 |
+ |
nMol = 108; |
| 797 |
+ |
} |
| 798 |
+ |
|
| 799 |
+ |
initialConfig = "./argon.init"; |
| 800 |
+ |
|
| 801 |
+ |
forceField = "LJ"; |
| 802 |
+ |
ensemble = "NVE"; |
| 803 |
+ |
dt = 1.0; |
| 804 |
+ |
runTime = 1e3; |
| 805 |
+ |
sampleTime = 100; |
| 806 |
+ |
statusTime = 50; |
| 807 |
+ |
|
| 808 |
+ |
\end{lstlisting} |
| 809 |
+ |
|
| 810 |
|
\subsection{\label{oopseSec:coordFiles}Coordinate Files} |
| 811 |
|
|
| 812 |
|
The standard format for storage of a systems coordinates is a modified |
| 813 |
|
xyz-file syntax, the exact details of which can be seen in |
| 814 |
< |
App.~\ref{appCoordFormat}. As all bonding and molecular information is |
| 815 |
< |
stored in the \texttt{.bass} and \texttt{.mdl} files, the coordinate |
| 816 |
< |
files are simply the complete set of coordinates for each atom at a |
| 817 |
< |
given simulation time. |
| 814 |
> |
Scheme~\ref{sch:dumpFormat}. As all bonding and molecular information |
| 815 |
> |
is stored in the \texttt{.bass} and \texttt{.mdl} files, the |
| 816 |
> |
coordinate files are simply the complete set of coordinates for each |
| 817 |
> |
atom at a given simulation time. One important note, although the |
| 818 |
> |
simulation propagates the complete rotation matrix, directional |
| 819 |
> |
entities are written out using quanternions, to save space in the |
| 820 |
> |
output files. |
| 821 |
|
|
| 822 |
< |
There are three major files used by {\sc oopse} written in the coordinate |
| 823 |
< |
format, they are as follows: the initialization file, the simulation |
| 824 |
< |
trajectory file, and the final coordinates of the simulation. The |
| 825 |
< |
initialization file is neccassary for {\sc oopse} to start the simulation |
| 826 |
< |
with the proper coordinates. It is typically denoted with the |
| 827 |
< |
extension \texttt{.init}. The trajectory file is created at the |
| 828 |
< |
beginning of the simulation, and is used to store snapshots of the |
| 829 |
< |
simulation at regular intervals. The first frame is a duplication of |
| 830 |
< |
the \texttt{.init} file, and each subsequent frame is appended to the |
| 831 |
< |
file at an interval specified in the \texttt{.bass} file. The |
| 832 |
< |
trajectory file is given the extension \texttt{.dump}. The final |
| 833 |
< |
coordinate file is the end of run or \texttt{.eor} file. The |
| 834 |
< |
\texttt{.eor} file stores the final configuration of teh system for a |
| 822 |
> |
\begin{lstlisting}[float,caption={[The format of the coordinate files]Shows the format of the coordinate files. The fist line is the number of atoms. The second line begins with the time stamp followed by the three $\mathbf{H}$ column vectors. The next lines are the atomic coordinates for all atoms in the system. First is the name followed by position, velocity, quanternions, and lastly angular momentum.},label=sch:dumpFormat] |
| 823 |
> |
|
| 824 |
> |
nAtoms |
| 825 |
> |
time; Hxx Hyx Hzx; Hxy Hyy Hzy; Hxz Hyz Hzz; |
| 826 |
> |
Name1 x y z vx vy vz q0 q1 q2 q3 jx jy jz |
| 827 |
> |
Name2 x y z vx vy vz q0 q1 q2 q3 jx jy jz |
| 828 |
> |
etc... |
| 829 |
> |
|
| 830 |
> |
\end{lstlisting} |
| 831 |
> |
|
| 832 |
> |
|
| 833 |
> |
There are three major files used by {\sc oopse} written in the |
| 834 |
> |
coordinate format, they are as follows: the initialization file |
| 835 |
> |
(\texttt{.init}), the simulation trajectory file (\texttt{.dump}), and |
| 836 |
> |
the final coordinates of the simulation. The initialization file is |
| 837 |
> |
necessary for {\sc oopse} to start the simulation with the proper |
| 838 |
> |
coordinates, and is generated before the simulation run. The |
| 839 |
> |
trajectory file is created at the beginning of the simulation, and is |
| 840 |
> |
used to store snapshots of the simulation at regular intervals. The |
| 841 |
> |
first frame is a duplication of the |
| 842 |
> |
\texttt{.init} file, and each subsequent frame is appended to the file |
| 843 |
> |
at an interval specified in the \texttt{.bass} file with the |
| 844 |
> |
\texttt{sampleTime} flag. The final coordinate file is the end of run file. The |
| 845 |
> |
\texttt{.eor} file stores the final configuration of the system for a |
| 846 |
|
given simulation. The file is updated at the same time as the |
| 847 |
< |
\texttt{.dump} file. However, it only contains the most recent |
| 847 |
> |
\texttt{.dump} file, however, it only contains the most recent |
| 848 |
|
frame. In this way, an \texttt{.eor} file may be used as the |
| 849 |
< |
initialization file to a second simulation in order to continue or |
| 850 |
< |
recover the previous simulation. |
| 849 |
> |
initialization file to a second simulation in order to continue a |
| 850 |
> |
simulation or recover one from a processor that has crashed during the |
| 851 |
> |
course of the run. |
| 852 |
|
|
| 853 |
< |
\subsection{Generation of Initial Coordinates} |
| 853 |
> |
\subsection{\label{oopseSec:initCoords}Generation of Initial Coordinates} |
| 854 |
|
|
| 855 |
< |
As was stated in Sec.~\ref{subSec:coordFiles}, an initialization file |
| 856 |
< |
is needed to provide the starting coordinates for a simulation. The |
| 857 |
< |
{\sc oopse} package provides a program called \texttt{sysBuilder} to aid in |
| 858 |
< |
the creation of the \texttt{.init} file. \texttt{sysBuilder} is {\sc bass} |
| 859 |
< |
aware, and will recognize arguments and parameters in the |
| 860 |
< |
\texttt{.bass} file that would otherwise be ignored by the |
| 861 |
< |
simulation. The program itself is under contiunual development, and is |
| 740 |
< |
offered here as a helper tool only. |
| 855 |
> |
As was stated in Sec.~\ref{oopseSec:coordFiles}, an initialization |
| 856 |
> |
file is needed to provide the starting coordinates for a |
| 857 |
> |
simulation. The {\sc oopse} package provides a program called |
| 858 |
> |
\texttt{sysBuilder} to aid in the creation of the \texttt{.init} |
| 859 |
> |
file. \texttt{sysBuilder} uses {\sc bass}, and will recognize |
| 860 |
> |
arguments and parameters in the \texttt{.bass} file that would |
| 861 |
> |
otherwise be ignored by the simulation. |
| 862 |
|
|
| 863 |
|
\subsection{The Statistics File} |
| 864 |
|
|
| 865 |
< |
The last output file generated by {\sc oopse} is the statistics file. This |
| 866 |
< |
file records such statistical quantities as the instantaneous |
| 867 |
< |
temperature, volume, pressure, etc. It is written out with the |
| 868 |
< |
frequency specified in the \texttt{.bass} file. The file allows the |
| 869 |
< |
user to observe the system variables as a function od simulation time |
| 870 |
< |
while the simulation is in progress. One useful function the |
| 871 |
< |
statistics file serves is to monitor the conserved quantity of a given |
| 872 |
< |
simulation ensemble, this allows the user to observe the stability of |
| 873 |
< |
the integrator. The statistics file is denoted with the \texttt{.stat} |
| 874 |
< |
file extension. |
| 865 |
> |
The last output file generated by {\sc oopse} is the statistics |
| 866 |
> |
file. This file records such statistical quantities as the |
| 867 |
> |
instantaneous temperature, volume, pressure, etc. It is written out |
| 868 |
> |
with the frequency specified in the \texttt{.bass} file with the |
| 869 |
> |
\texttt{statusTime} keyword. The file allows the user to observe the |
| 870 |
> |
system variables as a function of simulation time while the simulation |
| 871 |
> |
is in progress. One useful function the statistics file serves is to |
| 872 |
> |
monitor the conserved quantity of a given simulation ensemble, this |
| 873 |
> |
allows the user to observe the stability of the integrator. The |
| 874 |
> |
statistics file is denoted with the \texttt{.stat} file extension. |
| 875 |
|
|
| 876 |
|
\section{\label{oopseSec:mechanics}Mechanics} |
| 877 |
|
|
| 878 |
< |
\subsection{\label{integrate}Integrating the Equations of Motion: the Symplectic Step Integrator} |
| 878 |
> |
\subsection{\label{oopseSec:integrate}Integrating the Equations of Motion: the Symplectic Step Integrator} |
| 879 |
|
|
| 880 |
|
Integration of the equations of motion was carried out using the |
| 881 |
|
symplectic splitting method proposed by Dullweber \emph{et |
| 882 |
< |
al.}.\cite{Dullweber1997} The reason for this integrator selection |
| 883 |
< |
deals with poor energy conservation of rigid body systems using |
| 884 |
< |
quaternions. While quaternions work well for orientational motion in |
| 885 |
< |
alternate ensembles, the microcanonical ensemble has a constant energy |
| 886 |
< |
requirement that is quite sensitive to errors in the equations of |
| 887 |
< |
motion. The original implementation of this code utilized quaternions |
| 888 |
< |
for rotational motion propagation; however, a detailed investigation |
| 889 |
< |
showed that they resulted in a steady drift in the total energy, |
| 890 |
< |
something that has been observed by others.\cite{Laird97} |
| 882 |
> |
al.}.\cite{Dullweber1997} The reason for the selection of this |
| 883 |
> |
integrator, is the poor energy conservation of rigid body systems |
| 884 |
> |
using quaternion dynamics. While quaternions work well for |
| 885 |
> |
orientational motion in alternate ensembles, the microcanonical |
| 886 |
> |
ensemble has a constant energy requirement that is quite sensitive to |
| 887 |
> |
errors in the equations of motion. The original implementation of {\sc |
| 888 |
> |
oopse} utilized quaternions for rotational motion propagation; |
| 889 |
> |
however, a detailed investigation showed that they resulted in a |
| 890 |
> |
steady drift in the total energy, something that has been observed by |
| 891 |
> |
others.\cite{Laird97} |
| 892 |
|
|
| 893 |
|
The key difference in the integration method proposed by Dullweber |
| 894 |
< |
\emph{et al.} is that the entire rotation matrix is propagated from |
| 895 |
< |
one time step to the next. In the past, this would not have been as |
| 896 |
< |
feasible a option, being that the rotation matrix for a single body is |
| 897 |
< |
nine elements long as opposed to 3 or 4 elements for Euler angles and |
| 898 |
< |
quaternions respectively. System memory has become much less of an |
| 899 |
< |
issue in recent times, and this has resulted in substantial benefits |
| 900 |
< |
in energy conservation. There is still the issue of 5 or 6 additional |
| 779 |
< |
elements for describing the orientation of each particle, which will |
| 780 |
< |
increase dump files substantially. Simply translating the rotation |
| 781 |
< |
matrix into its component Euler angles or quaternions for storage |
| 782 |
< |
purposes relieves this burden. |
| 894 |
> |
\emph{et al}.~({\sc dlm}) is that the entire rotation matrix is propagated from |
| 895 |
> |
one time step to the next. In the past, this would not have been a |
| 896 |
> |
feasible option, since the rotation matrix for a single body is nine |
| 897 |
> |
elements long as opposed to three or four elements for Euler angles |
| 898 |
> |
and quaternions respectively. System memory has become much less of an |
| 899 |
> |
issue in recent times, and the {\sc dlm} method has used memory in |
| 900 |
> |
exchange for substantial benefits in energy conservation. |
| 901 |
|
|
| 902 |
< |
The symplectic splitting method allows for Verlet style integration of |
| 903 |
< |
both linear and angular motion of rigid bodies. In the integration |
| 904 |
< |
method, the orientational propagation involves a sequence of matrix |
| 902 |
> |
The {\sc dlm} method allows for Verlet style integration of both |
| 903 |
> |
linear and angular motion of rigid bodies. In the integration method, |
| 904 |
> |
the orientational propagation involves a sequence of matrix |
| 905 |
|
evaluations to update the rotation matrix.\cite{Dullweber1997} These |
| 906 |
< |
matrix rotations end up being more costly computationally than the |
| 907 |
< |
simpler arithmetic quaternion propagation. With the same time step, a |
| 908 |
< |
1000 SSD particle simulation shows an average 7\% increase in |
| 909 |
< |
computation time using the symplectic step method in place of |
| 910 |
< |
quaternions. This cost is more than justified when comparing the |
| 911 |
< |
energy conservation of the two methods as illustrated in figure |
| 794 |
< |
\ref{timestep}. |
| 906 |
> |
matrix rotations are more costly computationally than the simpler |
| 907 |
> |
arithmetic quaternion propagation. With the same time step, a 1000 SSD |
| 908 |
> |
particle simulation shows an average 7\% increase in computation time |
| 909 |
> |
using the {\sc dlm} method in place of quaternions. This cost is more |
| 910 |
> |
than justified when comparing the energy conservation of the two |
| 911 |
> |
methods as illustrated in Fig.~\ref{timestep}. |
| 912 |
|
|
| 913 |
|
\begin{figure} |
| 914 |
|
\centering |
| 915 |
|
\includegraphics[width=\linewidth]{timeStep.eps} |
| 916 |
< |
\caption{Energy conservation using quaternion based integration versus |
| 917 |
< |
the symplectic step method proposed by Dullweber \emph{et al.} with |
| 916 |
> |
\caption[Energy conservation for quaternion versus {\sc dlm} dynamics]{Energy conservation using quaternion based integration versus |
| 917 |
> |
the {\sc dlm} method with |
| 918 |
|
increasing time step. For each time step, the dotted line is total |
| 919 |
< |
energy using the symplectic step integrator, and the solid line comes |
| 919 |
> |
energy using the {\sc dlm} integrator, and the solid line comes |
| 920 |
|
from the quaternion integrator. The larger time step plots are shifted |
| 921 |
|
up from the true energy baseline for clarity.} |
| 922 |
|
\label{timestep} |
| 923 |
|
\end{figure} |
| 924 |
|
|
| 925 |
< |
In figure \ref{timestep}, the resulting energy drift at various time |
| 926 |
< |
steps for both the symplectic step and quaternion integration schemes |
| 925 |
> |
In Fig.~\ref{timestep}, the resulting energy drift at various time |
| 926 |
> |
steps for both the {\sc dlm} and quaternion integration schemes |
| 927 |
|
is compared. All of the 1000 SSD particle simulations started with the |
| 928 |
|
same configuration, and the only difference was the method for |
| 929 |
|
handling rotational motion. At time steps of 0.1 and 0.5 fs, both |
| 930 |
|
methods for propagating particle rotation conserve energy fairly well, |
| 931 |
|
with the quaternion method showing a slight energy drift over time in |
| 932 |
|
the 0.5 fs time step simulation. At time steps of 1 and 2 fs, the |
| 933 |
< |
energy conservation benefits of the symplectic step method are clearly |
| 933 |
> |
energy conservation benefits of the {\sc dlm} method are clearly |
| 934 |
|
demonstrated. Thus, while maintaining the same degree of energy |
| 935 |
|
conservation, one can take considerably longer time steps, leading to |
| 936 |
|
an overall reduction in computation time. |
| 939 |
|
time steps up to three femtoseconds. A slight energy drift on the |
| 940 |
|
order of 0.012 kcal/mol per nanosecond was observed at a time step of |
| 941 |
|
four femtoseconds, and as expected, this drift increases dramatically |
| 942 |
< |
with increasing time step. To insure accuracy in the constant energy |
| 826 |
< |
simulations, time steps were set at 2 fs and kept at this value for |
| 827 |
< |
constant pressure simulations as well. |
| 942 |
> |
with increasing time step. |
| 943 |
|
|
| 944 |
|
|
| 945 |
|
\subsection{\label{sec:extended}Extended Systems for other Ensembles} |
| 948 |
|
{\sc oopse} implements a |
| 949 |
|
|
| 950 |
|
|
| 951 |
< |
\subsubsection{\label{sec:noseHooverThermo}Nose-Hoover Thermostatting} |
| 951 |
> |
\subsection{\label{oopseSec:noseHooverThermo}Nose-Hoover Thermostatting} |
| 952 |
|
|
| 953 |
|
To mimic the effects of being in a constant temperature ({\sc nvt}) |
| 954 |
|
ensemble, {\sc oopse} uses the Nose-Hoover extended system |
| 975 |
|
to values of a few ps. Within a {\sc bass} file, $\tau_{T}$ could be |
| 976 |
|
set to 1 ps using the {\tt tauThermostat = 1000; } command. |
| 977 |
|
|
| 978 |
+ |
\subsection{\label{oopseSec:rattle}The {\sc rattle} Method for Bond |
| 979 |
+ |
Constraints} |
| 980 |
|
|
| 981 |
< |
\subsection{\label{Sec:zcons}Z-Constraint Method} |
| 981 |
> |
In order to satisfy the constraints of fixed bond lengths within {\sc |
| 982 |
> |
oopse}, we have implemented the {\sc rattle} algorithm of |
| 983 |
> |
Andersen.\cite{andersen83} The algorithm is a velocity verlet |
| 984 |
> |
formulation of the {\sc shake} method\cite{ryckaert77} of iteratively |
| 985 |
> |
solving the Lagrange multipliers of constraint. The system of lagrange |
| 986 |
> |
multipliers allows one to reformulate the equations of motion with |
| 987 |
> |
explicit constraint forces on the equations of |
| 988 |
> |
motion.\cite{fowles99:lagrange} |
| 989 |
|
|
| 990 |
< |
Based on fluctuatin-dissipation theorem,\bigskip\ force auto-correlation |
| 990 |
> |
Consider a system described by qoordinates $q_1$ and $q_2$ subject to an |
| 991 |
> |
equation of constraint: |
| 992 |
> |
\begin{equation} |
| 993 |
> |
\sigma(q_1, q_2,t) = 0 |
| 994 |
> |
\label{oopseEq:lm1} |
| 995 |
> |
\end{equation} |
| 996 |
> |
The Lagrange formulation of the equations of motion can be written: |
| 997 |
> |
\begin{equation} |
| 998 |
> |
\delta\int_{t_1}^{t_2}L\, dt = |
| 999 |
> |
\int_{t_1}^{t_2} \sum_i \biggl [ \frac{\partial L}{\partial q_i} |
| 1000 |
> |
- \frac{d}{dt}\biggl(\frac{\partial L}{\partial \dot{q}_i} |
| 1001 |
> |
\biggr ) \biggr] \delta q_i \, dt = 0 |
| 1002 |
> |
\label{oopseEq:lm2} |
| 1003 |
> |
\end{equation} |
| 1004 |
> |
Here, $\delta q_i$ is not independent for each $q$, as $q_1$ and $q_2$ |
| 1005 |
> |
are linked by $\sigma$. However, $\sigma$ is fixed at any given |
| 1006 |
> |
instant of time, giving: |
| 1007 |
> |
\begin{align} |
| 1008 |
> |
\delta\sigma &= \biggl( \frac{\partial\sigma}{\partial q_1} \delta q_1 % |
| 1009 |
> |
+ \frac{\partial\sigma}{\partial q_2} \delta q_2 \biggr) = 0 \\ |
| 1010 |
> |
% |
| 1011 |
> |
\frac{\partial\sigma}{\partial q_1} \delta q_1 &= % |
| 1012 |
> |
- \frac{\partial\sigma}{\partial q_2} \delta q_2 \\ |
| 1013 |
> |
% |
| 1014 |
> |
\delta q_2 &= - \biggl(\frac{\partial\sigma}{\partial q_1} \bigg / % |
| 1015 |
> |
\frac{\partial\sigma}{\partial q_2} \biggr) \delta q_1 |
| 1016 |
> |
\end{align} |
| 1017 |
> |
Substituted back into Eq.~\ref{oopseEq:lm2}, |
| 1018 |
> |
\begin{equation} |
| 1019 |
> |
\int_{t_1}^{t_2}\biggl [ \biggl(\frac{\partial L}{\partial q_1} |
| 1020 |
> |
- \frac{d}{dt}\,\frac{\partial L}{\partial \dot{q}_1} |
| 1021 |
> |
\biggr) |
| 1022 |
> |
- \biggl( \frac{\partial L}{\partial q_1} |
| 1023 |
> |
- \frac{d}{dt}\,\frac{\partial L}{\partial \dot{q}_1} |
| 1024 |
> |
\biggr) \biggl(\frac{\partial\sigma}{\partial q_1} \bigg / % |
| 1025 |
> |
\frac{\partial\sigma}{\partial q_2} \biggr)\biggr] \delta q_1 \, dt = 0 |
| 1026 |
> |
\label{oopseEq:lm3} |
| 1027 |
> |
\end{equation} |
| 1028 |
> |
Leading to, |
| 1029 |
> |
\begin{equation} |
| 1030 |
> |
\frac{\biggl(\frac{\partial L}{\partial q_1} |
| 1031 |
> |
- \frac{d}{dt}\,\frac{\partial L}{\partial \dot{q}_1} |
| 1032 |
> |
\biggr)}{\frac{\partial\sigma}{\partial q_1}} = |
| 1033 |
> |
\frac{\biggl(\frac{\partial L}{\partial q_2} |
| 1034 |
> |
- \frac{d}{dt}\,\frac{\partial L}{\partial \dot{q}_2} |
| 1035 |
> |
\biggr)}{\frac{\partial\sigma}{\partial q_2}} |
| 1036 |
> |
\label{oopseEq:lm4} |
| 1037 |
> |
\end{equation} |
| 1038 |
> |
This relation can only be statisfied, if both are equal to a single |
| 1039 |
> |
function $-\lambda(t)$, |
| 1040 |
> |
\begin{align} |
| 1041 |
> |
\frac{\biggl(\frac{\partial L}{\partial q_1} |
| 1042 |
> |
- \frac{d}{dt}\,\frac{\partial L}{\partial \dot{q}_1} |
| 1043 |
> |
\biggr)}{\frac{\partial\sigma}{\partial q_1}} &= -\lambda(t) \\ |
| 1044 |
> |
% |
| 1045 |
> |
\frac{\partial L}{\partial q_1} |
| 1046 |
> |
- \frac{d}{dt}\,\frac{\partial L}{\partial \dot{q}_1} &= |
| 1047 |
> |
-\lambda(t)\,\frac{\partial\sigma}{\partial q_1} \\ |
| 1048 |
> |
% |
| 1049 |
> |
\frac{\partial L}{\partial q_1} |
| 1050 |
> |
- \frac{d}{dt}\,\frac{\partial L}{\partial \dot{q}_1} |
| 1051 |
> |
+ \mathcal{G}_i &= 0 |
| 1052 |
> |
\end{align} |
| 1053 |
> |
Where $\mathcal{G}_i$, the force of constraint on $i$, is: |
| 1054 |
> |
\begin{equation} |
| 1055 |
> |
\mathcal{G}_i = \lambda(t)\,\frac{\partial\sigma}{\partial q_1} |
| 1056 |
> |
\label{oopseEq:lm5} |
| 1057 |
> |
\end{equation} |
| 1058 |
> |
|
| 1059 |
> |
In a simulation, this would involve the solution of a set of $(m + n)$ |
| 1060 |
> |
number of equations. Where $m$ is the number of constraints, and $n$ |
| 1061 |
> |
is the number of constrained coordinates. In practice, this is not |
| 1062 |
> |
done, as the matrix inversion neccassary to solve the system of |
| 1063 |
> |
equations would be very time consuming to solve. Additionally, the |
| 1064 |
> |
numerical error in the solution of the set of $\lambda$'s would be |
| 1065 |
> |
compounded by the error inherent in propagating by the Velocity Verlet |
| 1066 |
> |
algorithm ($\Delta t^4$). The verlet propagation error is negligible |
| 1067 |
> |
in an unconstrained system, as one is interested in the statisitics of |
| 1068 |
> |
the run, and not that the run be numerically exact to the ``true'' |
| 1069 |
> |
integration. This relates back to the ergodic hypothesis that a time |
| 1070 |
> |
integral of a valid trajectory will still give the correct enesemble |
| 1071 |
> |
average. However, in the case of constraints, if the equations of |
| 1072 |
> |
motion leave the ``true'' trajectory, they are departing from the |
| 1073 |
> |
constrained surface. The method that is used, is to iteratively solve |
| 1074 |
> |
for $\lambda(t)$ at each time step. |
| 1075 |
> |
|
| 1076 |
> |
In {\sc rattle} the equations of motion are modified subject to the |
| 1077 |
> |
following two constraints: |
| 1078 |
> |
\begin{align} |
| 1079 |
> |
\sigma_{ij}[\mathbf{r}(t)] \equiv |
| 1080 |
> |
[ \mathbf{r}_i(t) - \mathbf{r}_j(t)]^2 - d_{ij}^2 &= 0 % |
| 1081 |
> |
\label{oopseEq:c1} \\ |
| 1082 |
> |
% |
| 1083 |
> |
[\mathbf{\dot{r}}_i(t) - \mathbf{\dot{r}}_j(t)] \cdot |
| 1084 |
> |
[\mathbf{r}_i(t) - \mathbf{r}_j(t)] &= 0 \label{oopseEq:c2} |
| 1085 |
> |
\end{align} |
| 1086 |
> |
Eq.~\ref{oopseEq:c1} is the set of bond constraints, where $d_{ij}$ is |
| 1087 |
> |
the constrained distance between atom $i$ and |
| 1088 |
> |
$j$. Eq.~\ref{oopseEq:c2} constrains the velocities of $i$ and $j$ to |
| 1089 |
> |
be perpindicular to the bond vector, so that the bond can neither grow |
| 1090 |
> |
nor shrink. The constrained dynamics equations become: |
| 1091 |
> |
\begin{equation} |
| 1092 |
> |
m_i \mathbf{\ddot{r}}_i = \mathbf{F}_i + \mathbf{\mathcal{G}}_i |
| 1093 |
> |
\label{oopseEq:r1} |
| 1094 |
> |
\end{equation} |
| 1095 |
> |
Where, |
| 1096 |
> |
\begin{equation} |
| 1097 |
> |
\mathbf{\mathcal{G}}_i = - \sum_j \lambda_{ij}(t)\,\nabla \sigma_{ij} |
| 1098 |
> |
\label{oopseEq:r2} |
| 1099 |
> |
\end{equation} |
| 1100 |
> |
|
| 1101 |
> |
In Velocity Verlet, if $\Delta t = h$, the propagation can be written: |
| 1102 |
> |
\begin{align} |
| 1103 |
> |
\mathbf{r}_i(t+h) &= |
| 1104 |
> |
\mathbf{r}_i(t) + h\mathbf{\dot{r}}(t) + |
| 1105 |
> |
\frac{h^2}{2m_i}\,\Bigl[ \mathbf{F}_i(t) + |
| 1106 |
> |
\mathbf{\mathcal{G}}_{Ri}(t) \Bigr] \label{oopseEq:vv1} \\ |
| 1107 |
> |
% |
| 1108 |
> |
\mathbf{\dot{r}}_i(t+h) &= |
| 1109 |
> |
\mathbf{\dot{r}}_i(t) + \frac{h}{2m_i} |
| 1110 |
> |
\Bigl[ \mathbf{F}_i(t) + \mathbf{\mathcal{G}}_{Ri}(t) + |
| 1111 |
> |
\mathbf{F}_i(t+h) + \mathbf{\mathcal{G}}_{Vi}(t+h) \Bigr] % |
| 1112 |
> |
\label{oopseEq:vv2} |
| 1113 |
> |
\end{align} |
| 1114 |
> |
|
| 1115 |
> |
|
| 1116 |
> |
|
| 1117 |
> |
\subsection{\label{oopseSec:zcons}Z-Constraint Method} |
| 1118 |
> |
|
| 1119 |
> |
Based on fluctuation-dissipation theorem, a force auto-correlation |
| 1120 |
|
method was developed to investigate the dynamics of ions inside the ion |
| 1121 |
|
channels.\cite{Roux91} Time-dependent friction coefficient can be calculated |
| 1122 |
< |
from the deviation of the instaneous force from its mean force. |
| 1122 |
> |
from the deviation of the instantaneous force from its mean force. |
| 1123 |
|
|
| 1124 |
|
% |
| 1125 |
|
|
| 1155 |
|
resetting the coordinate, we reset the forces of z-constraint molecules as |
| 1156 |
|
well as subtract the total constraint forces from the rest of the system after |
| 1157 |
|
force calculation at each time step. |
| 1158 |
< |
\begin{verbatim} |
| 1159 |
< |
$F_{\alpha i}=0$ |
| 1160 |
< |
$V_{\alpha i}=V_{\alpha i}-\frac{\sum\limits_{i}M_{_{\alpha i}}V_{\alpha i}}{\sum\limits_{i}M_{_{\alpha i}}}$ |
| 1161 |
< |
$F_{\alpha i}=F_{\alpha i}-\frac{M_{_{\alpha i}}}{\sum\limits_{\alpha}\sum\limits_{i}M_{_{\alpha i}}}\sum\limits_{\beta}F_{\beta}$ |
| 1162 |
< |
$V_{\alpha i}=V_{\alpha i}-\frac{\sum\limits_{\alpha}\sum\limits_{i}M_{_{\alpha i}}V_{\alpha i}}{\sum\limits_{\alpha}\sum\limits_{i}M_{_{\alpha i}}}$ |
| 1163 |
< |
\end{verbatim} |
| 1158 |
> |
\begin{align} |
| 1159 |
> |
F_{\alpha i}&=0\\ |
| 1160 |
> |
V_{\alpha i}&=V_{\alpha i}-\frac{\sum\limits_{i}M_{_{\alpha i}}V_{\alpha i}}{\sum\limits_{i}M_{_{\alpha i}}}\\ |
| 1161 |
> |
F_{\alpha i}&=F_{\alpha i}-\frac{M_{_{\alpha i}}}{\sum\limits_{\alpha}\sum\limits_{i}M_{_{\alpha i}}}\sum\limits_{\beta}F_{\beta}\\ |
| 1162 |
> |
V_{\alpha i}&=V_{\alpha i}-\frac{\sum\limits_{\alpha}\sum\limits_{i}M_{_{\alpha i}}V_{\alpha i}}{\sum\limits_{\alpha}\sum\limits_{i}M_{_{\alpha i}}} |
| 1163 |
> |
\end{align} |
| 1164 |
|
|
| 1165 |
|
At the very beginning of the simulation, the molecules may not be at its |
| 1166 |
|
constraint position. To move the z-constraint molecule to the specified |
| 1310 |
|
between the frames contained within the two blocks are |
| 1311 |
|
calculated. Once completed, the memory blocks are incremented, and the |
| 1312 |
|
process is repeated. A diagram illustrating the process is shown in |
| 1313 |
< |
Fig.~\ref{fig:dynamicPropsMemory}. As was the case with \texttt{staticProps}, |
| 1314 |
< |
multiple properties may be calculated in a single run to avoid |
| 1315 |
< |
multiple reads on the same file. |
| 1313 |
> |
Fig.~\ref{oopseFig:dynamicPropsMemory}. As was the case with |
| 1314 |
> |
\texttt{staticProps}, multiple properties may be calculated in a |
| 1315 |
> |
single run to avoid multiple reads on the same file. |
| 1316 |
|
|
| 1064 |
– |
\begin{figure} |
| 1065 |
– |
\centering |
| 1066 |
– |
\includegraphics[width=\linewidth]{dynamicPropsMem.eps} |
| 1067 |
– |
\caption{This diagram illustrates the dynamic memory allocation used by \texttt{dynamicProps}, which follows the scheme: $\sum^{N_{\text{memory blocks}}}_{i=1}[ \operatorname{self}(i) + \sum^{N_{\text{memory blocks}}}_{j>i} \operatorname{cross}(i,j)]$. The shaded region represents the self correlation of the memory block, and the open blocks are read one at a time and the cross correlations between blocks are calculated.} |
| 1068 |
– |
\label{fig:dynamicPropsMemory} |
| 1069 |
– |
\end{figure} |
| 1317 |
|
|
| 1318 |
+ |
|
| 1319 |
|
\section{\label{oopseSec:design}Program Design} |
| 1320 |
|
|
| 1321 |
< |
\subsection{\label{sec:architecture} OOPSE Architecture} |
| 1321 |
> |
\subsection{\label{sec:architecture} {\sc oopse} Architecture} |
| 1322 |
|
|
| 1323 |
< |
The core of OOPSE is divided into two main object libraries: {\texttt |
| 1324 |
< |
libBASS} and {\texttt libmdtools}. {\texttt libBASS} is the library |
| 1325 |
< |
developed around the parseing engine and {\texttt libmdtools} is the |
| 1326 |
< |
software library developed around the simulation engine. |
| 1323 |
> |
The core of OOPSE is divided into two main object libraries: |
| 1324 |
> |
\texttt{libBASS} and \texttt{libmdtools}. \texttt{libBASS} is the |
| 1325 |
> |
library developed around the parsing engine and \texttt{libmdtools} |
| 1326 |
> |
is the software library developed around the simulation engine. These |
| 1327 |
> |
two libraries are designed to encompass all the basic functions and |
| 1328 |
> |
tools that {\sc oopse} provides. Utility programs, such as the |
| 1329 |
> |
property analyzers, need only link against the software libraries to |
| 1330 |
> |
gain access to parsing, force evaluation, and input / output |
| 1331 |
> |
routines. |
| 1332 |
|
|
| 1333 |
< |
\subsection{\label{sec:parallelization} Parallelization of OOPSE} |
| 1333 |
> |
Contained in \texttt{libBASS} are all the routines associated with |
| 1334 |
> |
reading and parsing the \texttt{.bass} input files. Given a |
| 1335 |
> |
\texttt{.bass} file, \texttt{libBASS} will open it and any associated |
| 1336 |
> |
\texttt{.mdl} files; then create structures in memory that are |
| 1337 |
> |
templates of all the molecules specified in the input files. In |
| 1338 |
> |
addition, any simulation parameters set in the \texttt{.bass} file |
| 1339 |
> |
will be placed in a structure for later query by the controlling |
| 1340 |
> |
program. |
| 1341 |
|
|
| 1342 |
< |
Although processor power is doubling roughly every 18 months according |
| 1343 |
< |
to the famous Moore's Law\cite{moore}, it is still unreasonable to |
| 1344 |
< |
simulate systems of more then a 1000 atoms on a single processor. To |
| 1345 |
< |
facilitate study of larger system sizes or smaller systems on long |
| 1346 |
< |
time scales in a reasonable period of time, parallel methods were |
| 1347 |
< |
developed allowing multiple CPU's to share the simulation |
| 1348 |
< |
workload. Three general categories of parallel decomposition method's |
| 1349 |
< |
have been developed including atomic, spatial and force decomposition |
| 1350 |
< |
methods. |
| 1342 |
> |
Located in \texttt{libmdtools} are all other routines necessary to a |
| 1343 |
> |
Molecular Dynamics simulation. The library uses the main data |
| 1344 |
> |
structures returned by \texttt{libBASS} to initialize the various |
| 1345 |
> |
parts of the simulation: the atom structures and positions, the force |
| 1346 |
> |
field, the integrator, \emph{et cetera}. After initialization, the |
| 1347 |
> |
library can be used to perform a variety of tasks: integrate a |
| 1348 |
> |
Molecular Dynamics trajectory, query phase space information from a |
| 1349 |
> |
specific frame of a completed trajectory, or even recalculate force or |
| 1350 |
> |
energetic information about specific frames from a completed |
| 1351 |
> |
trajectory. |
| 1352 |
|
|
| 1353 |
+ |
With these core libraries in place, several programs have been |
| 1354 |
+ |
developed to utilize the routines provided by \texttt{libBASS} and |
| 1355 |
+ |
\texttt{libmdtools}. The main program of the package is \texttt{oopse} |
| 1356 |
+ |
and the corresponding parallel version \texttt{oopse\_MPI}. These two |
| 1357 |
+ |
programs will take the \texttt{.bass} file, and create then integrate |
| 1358 |
+ |
the simulation specified in the script. The two analysis programs |
| 1359 |
+ |
\texttt{staticProps} and \texttt{dynamicProps} utilize the core |
| 1360 |
+ |
libraries to initialize and read in trajectories from previously |
| 1361 |
+ |
completed simulations, in addition to the ability to use functionality |
| 1362 |
+ |
from \texttt{libmdtools} to recalculate forces and energies at key |
| 1363 |
+ |
frames in the trajectories. Lastly, the family of system building |
| 1364 |
+ |
programs (Sec.~\ref{oopseSec:initCoords}) also use the libraries to |
| 1365 |
+ |
store and output the system configurations they create. |
| 1366 |
+ |
|
| 1367 |
+ |
\subsection{\label{oopseSec:parallelization} Parallelization of {\sc oopse}} |
| 1368 |
+ |
|
| 1369 |
+ |
Although processor power is continually growing month by month, it is |
| 1370 |
+ |
still unreasonable to simulate systems of more then a 1000 atoms on a |
| 1371 |
+ |
single processor. To facilitate study of larger system sizes or |
| 1372 |
+ |
smaller systems on long time scales in a reasonable period of time, |
| 1373 |
+ |
parallel methods were developed allowing multiple CPU's to share the |
| 1374 |
+ |
simulation workload. Three general categories of parallel |
| 1375 |
+ |
decomposition method's have been developed including atomic, spatial |
| 1376 |
+ |
and force decomposition methods. |
| 1377 |
+ |
|
| 1378 |
|
Algorithmically simplest of the three method's is atomic decomposition |
| 1379 |
|
where N particles in a simulation are split among P processors for the |
| 1380 |
|
duration of the simulation. Computational cost scales as an optimal |
| 1409 |
|
favorably then spatial decomposition up to 10,000 atoms and favorably |
| 1410 |
|
competes with spatial methods for up to 100,000 atoms. |
| 1411 |
|
|
| 1412 |
< |
\subsection{\label{openSource}Open Source and Distribution License} |
| 1412 |
> |
\subsection{\label{oopseSec:memAlloc}Memory Issues in Trajectory Analysis} |
| 1413 |
|
|
| 1414 |
+ |
For large simulations, the trajectory files can sometimes reach sizes |
| 1415 |
+ |
in excess of several gigabytes. In order to effectively analyze that |
| 1416 |
+ |
amount of data+, two memory management schemes have been devised for |
| 1417 |
+ |
\texttt{staticProps} and for \texttt{dynamicProps}. The first scheme, |
| 1418 |
+ |
developed for \texttt{staticProps}, is the simplest. As each frame's |
| 1419 |
+ |
statistics are calculated independent of each other, memory is |
| 1420 |
+ |
allocated for each frame, then freed once correlation calculations are |
| 1421 |
+ |
complete for the snapshot. To prevent multiple passes through a |
| 1422 |
+ |
potentially large file, \texttt{staticProps} is capable of calculating |
| 1423 |
+ |
all requested correlations per frame with only a single pair loop in |
| 1424 |
+ |
each frame and a single read through of the file. |
| 1425 |
|
|
| 1426 |
< |
\section{\label{oopseSec:conclusion}Conclusion} |
| 1426 |
> |
The second, more advanced memory scheme, is used by |
| 1427 |
> |
\texttt{dynamicProps}. Here, the program must have multiple frames in |
| 1428 |
> |
memory to calculate time dependent correlations. In order to prevent a |
| 1429 |
> |
situation where the program runs out of memory due to large |
| 1430 |
> |
trajectories, the user is able to specify that the trajectory be read |
| 1431 |
> |
in blocks. The number of frames in each block is specified by the |
| 1432 |
> |
user, and upon reading a block of the trajectory, |
| 1433 |
> |
\texttt{dynamicProps} will calculate all of the time correlation frame |
| 1434 |
> |
pairs within the block. After in block correlations are complete, a |
| 1435 |
> |
second block of the trajectory is read, and the cross correlations are |
| 1436 |
> |
calculated between the two blocks. this second block is then freed and |
| 1437 |
> |
then incremented and the process repeated until the end of the |
| 1438 |
> |
trajectory. Once the end is reached, the first block is freed then |
| 1439 |
> |
incremented, and the again the internal time correlations are |
| 1440 |
> |
calculated. The algorithm with the second block is then repeated with |
| 1441 |
> |
the new origin block, until all frame pairs have been correlated in |
| 1442 |
> |
time. This process is illustrated in |
| 1443 |
> |
Fig.~\ref{oopseFig:dynamicPropsMemory}. |
| 1444 |
|
|
| 1445 |
< |
\begin{itemize} |
| 1446 |
< |
|
| 1447 |
< |
\item Restate capabilities |
| 1445 |
> |
\begin{figure} |
| 1446 |
> |
\centering |
| 1447 |
> |
\includegraphics[width=\linewidth]{dynamicPropsMem.eps} |
| 1448 |
> |
\caption[A representation of the block correlations in \texttt{dynamicProps}]{This diagram illustrates the memory management used by \texttt{dynamicProps}, which follows the scheme: $\sum^{N_{\text{memory blocks}}}_{i=1}[ \operatorname{self}(i) + \sum^{N_{\text{memory blocks}}}_{j>i} \operatorname{cross}(i,j)]$. The shaded region represents the self correlation of the memory block, and the open blocks are read one at a time and the cross correlations between blocks are calculated.} |
| 1449 |
> |
\label{oopseFig:dynamicPropsMemory} |
| 1450 |
> |
\end{figure} |
| 1451 |
|
|
| 1452 |
< |
\item recap major structure / design choices |
| 1452 |
> |
\subsection{\label{openSource}Open Source and Distribution License} |
| 1453 |
|
|
| 1454 |
< |
\begin{itemize} |
| 1138 |
< |
|
| 1139 |
< |
\item parallel |
| 1140 |
< |
\item symplectic integration |
| 1141 |
< |
\item languages |
| 1454 |
> |
\section{\label{oopseSec:conclusion}Conclusion} |
| 1455 |
|
|
| 1456 |
< |
\end{itemize} |
| 1456 |
> |
We have presented the design and implementation of our open source |
| 1457 |
> |
simulation package {\sc oopse}. The package offers novel |
| 1458 |
> |
capabilities to the field of Molecular Dynamics simulation packages in |
| 1459 |
> |
the form of dipolar force fields, and symplectic integration of rigid |
| 1460 |
> |
body dynamics. It is capable of scaling across multiple processors |
| 1461 |
> |
through the use of MPI. It also implements several integration |
| 1462 |
> |
ensembles allowing the end user control over temperature and |
| 1463 |
> |
pressure. In addition, it is capable of integrating constrained |
| 1464 |
> |
dynamics through both the {\sc rattle} algorithm and the z-constraint |
| 1465 |
> |
method. |
| 1466 |
|
|
| 1467 |
< |
\item How well does it meet the primary goal |
| 1467 |
> |
These features are all brought together in a single open-source |
| 1468 |
> |
development package. This allows researchers to not only benefit from |
| 1469 |
> |
{\sc oopse}, but also contribute to {\sc oopse}'s development as |
| 1470 |
> |
well.Documentation and source code for {\sc oopse} can be downloaded |
| 1471 |
> |
from \texttt{http://www.openscience.org/oopse/}. |
| 1472 |
|
|
| 1147 |
– |
\end{itemize} |
