| 1 |
|
\documentclass[11pt]{article} |
| 2 |
|
\usepackage{amsmath} |
| 3 |
< |
\usepackage{amssymb} |
| 3 |
> |
\usepackage{amssymb} |
| 4 |
|
\usepackage{endfloat} |
| 5 |
– |
%\usepackage{berkeley} |
| 5 |
|
\usepackage{listings} |
| 6 |
+ |
\usepackage{berkeley} |
| 7 |
|
\usepackage{graphicx} |
| 8 |
|
\usepackage[ref]{overcite} |
| 9 |
|
\usepackage{setspace} |
| 18 |
|
\renewcommand\citemid{\ } % no comma in optional reference note |
| 19 |
|
|
| 20 |
|
\begin{document} |
| 21 |
< |
\lstset{language=C,float,frame=tblr,frameround=tttt} |
| 21 |
> |
\lstset{language=C,frame=TB,basicstyle=\small,basicstyle=\ttfamily, % |
| 22 |
> |
xleftmargin=0.5in, xrightmargin=0.5in,captionpos=b, % |
| 23 |
> |
abovecaptionskip=0.5cm, belowcaptionskip=0.5cm} |
| 24 |
|
\renewcommand{\lstlistingname}{Scheme} |
| 25 |
< |
\title{{\sc oopse}: An Open Source Object-Oriented Parallel Simulation |
| 25 |
> |
\title{{\sc oopse}: An Object-Oriented Parallel Simulation |
| 26 |
|
Engine for Molecular Dynamics} |
| 27 |
|
|
| 28 |
< |
\author{Matthew A. Meineke, Charles F. Vardeman II, Teng Lin, Christopher J. Fennell and J. Daniel Gezelter\\ |
| 28 |
> |
\author{Matthew A. Meineke, Charles F. Vardeman II, Teng Lin,\\ |
| 29 |
> |
Christopher J. Fennell and J. Daniel Gezelter\\ |
| 30 |
|
Department of Chemistry and Biochemistry\\ |
| 31 |
|
University of Notre Dame\\ |
| 32 |
|
Notre Dame, Indiana 46556} |
| 35 |
|
\maketitle |
| 36 |
|
|
| 37 |
|
\begin{abstract} |
| 38 |
< |
We detail the capabilities of a new open-source parallel simulation |
| 39 |
< |
package ({\sc oopse}) that can perform molecular dynamics simulations |
| 40 |
< |
on atom types that are missing from other popular packages. In |
| 41 |
< |
particular, {\sc oopse} is capable of performing orientational |
| 42 |
< |
dynamics on dipolar systems, and it can handle simulations of metallic |
| 43 |
< |
systems using the embedded atom method ({\sc eam}). |
| 38 |
> |
{\sc oopse} is a new molecular dynamics simulation program which is |
| 39 |
> |
capable of efficiently integrating equations of motion for atom types |
| 40 |
> |
with orientational degrees of freedom (e.g. ``sticky'' atoms and point |
| 41 |
> |
dipoles). Transition metals can also be simulated using the embedded |
| 42 |
> |
atom method ({\sc eam}) potential included in the code. Parallel |
| 43 |
> |
simulations are carried out using the force-based decomposition |
| 44 |
> |
method. Simulations are specified using a very simple C-based |
| 45 |
> |
meta-data language. A number of advanced integrators are included, |
| 46 |
> |
and the basic integrator for orientational dynamics provides |
| 47 |
> |
substantial improvements over older quaternion-based schemes. |
| 48 |
|
\end{abstract} |
| 49 |
|
|
| 50 |
|
\section{\label{sec:intro}Introduction} |
| 51 |
|
|
| 52 |
< |
When choosing to simulate a chemical system with molecular dynamics, |
| 53 |
< |
there are a variety of options available. For simple systems, one |
| 54 |
< |
might consider writing one's own programming code. However, as systems |
| 55 |
< |
grow larger and more complex, building and maintaining code for the |
| 56 |
< |
simulations becomes a time consuming task. In such cases it is usually |
| 57 |
< |
more convenient for a researcher to turn to pre-existing simulation |
| 58 |
< |
packages. These packages, such as {\sc amber}\cite{pearlman:1995} and |
| 59 |
< |
{\sc charmm}\cite{Brooks83}, provide powerful tools for researchers to |
| 60 |
< |
conduct simulations of their systems without spending their time |
| 61 |
< |
developing a code base to conduct their research. This then frees them |
| 62 |
< |
to perhaps explore experimental analogues to their models. |
| 52 |
> |
There are a number of excellent molecular dynamics packages available |
| 53 |
> |
to the chemical physics |
| 54 |
> |
community.\cite{Brooks83,MacKerell98,pearlman:1995,Gromacs,Gromacs3,DL_POLY,Tinker,Paradyn,namd,macromodel} |
| 55 |
> |
All of these packages are stable, polished programs which solve many |
| 56 |
> |
problems of interest. Most are now capable of performing molecular |
| 57 |
> |
dynamics simulations on parallel computers. Some have source code |
| 58 |
> |
which is freely available to the entire scientific community. Few, |
| 59 |
> |
however, are capable of efficiently integrating the equations of |
| 60 |
> |
motion for atom types with orientational degrees of freedom |
| 61 |
> |
(e.g. point dipoles, and ``sticky'' atoms). And only one of the |
| 62 |
> |
programs referenced can handle transition metal force fields like the |
| 63 |
> |
Embedded Atom Method ({\sc eam}). The direction our research program |
| 64 |
> |
has taken us now involves the use of atoms with orientational degrees |
| 65 |
> |
of freedom as well as transition metals. Since these simulation |
| 66 |
> |
methods may be of some use to other researchers, we have decided to |
| 67 |
> |
release our program (and all related source code) to the scientific |
| 68 |
> |
community. |
| 69 |
|
|
| 70 |
< |
Despite their utility, problems with these packages arise when |
| 71 |
< |
researchers try to develop techniques or energetic models that the |
| 72 |
< |
code was not originally designed to simulate. Examples of techniques |
| 73 |
< |
and energetics not commonly implemented include; dipole-dipole |
| 74 |
< |
interactions, rigid body dynamics, and metallic potentials. When faced |
| 75 |
< |
with these obstacles, a researcher must either develop their own code |
| 76 |
< |
or license and extend one of the commercial packages. What we have |
| 77 |
< |
elected to do is develop a body of simulation code capable of |
| 78 |
< |
implementing the types of models upon which our research is based. |
| 70 |
> |
This paper communicates the algorithmic details of our program, which |
| 71 |
> |
we have been calling the Object-Oriented Parallel Simulation Engine |
| 72 |
> |
(i.e. {\sc oopse}). We have structured this paper to first discuss |
| 73 |
> |
the underlying concepts in this simulation package |
| 74 |
> |
(Sec. \ref{oopseSec:IOfiles}). The empirical energy functions |
| 75 |
> |
implemented are discussed in Sec.~\ref{oopseSec:empiricalEnergy}. |
| 76 |
> |
Sec.~\ref{oopseSec:mechanics} describes the various Molecular Dynamics |
| 77 |
> |
algorithms {\sc oopse} implements in the integration of Hamilton's |
| 78 |
> |
equations of motion. Program design considerations for parallel |
| 79 |
> |
computing are presented in |
| 80 |
> |
Sec.~\ref{oopseSec:parallelization}. Concluding remarks are presented |
| 81 |
> |
in Sec.~\ref{oopseSec:conclusion}. |
| 82 |
|
|
| 83 |
< |
In developing {\sc oopse}, we have adhered to the precepts of Open |
| 68 |
< |
Source development, and are releasing our source code with a |
| 69 |
< |
permissive license. It is our intent that by doing so, other |
| 70 |
< |
researchers might benefit from our work, and add their own |
| 71 |
< |
contributions to the package. The license under which {\sc oopse} is |
| 72 |
< |
distributed allows any researcher to download and modify the source |
| 73 |
< |
code for their own use. In this way further development of {\sc oopse} |
| 74 |
< |
is not limited to only the models of interest to ourselves, but also |
| 75 |
< |
those of the community of scientists who contribute back to the |
| 76 |
< |
project. |
| 83 |
> |
\section{\label{oopseSec:IOfiles}Concepts \& Files} |
| 84 |
|
|
| 85 |
< |
We have structured this chapter to first discuss the empirical energy |
| 86 |
< |
functions that {\sc oopse } implements in |
| 87 |
< |
Sec.~\ref{oopseSec:empiricalEnergy}. Following that is a discussion of |
| 88 |
< |
the various input and output files associated with the package |
| 89 |
< |
(Sec.~\ref{oopseSec:IOfiles}). Sec.~\ref{oopseSec:mechanics} |
| 90 |
< |
elucidates the various Molecular Dynamics algorithms {\sc oopse} |
| 84 |
< |
implements in the integration of the Newtonian equations of |
| 85 |
< |
motion. Basic analysis of the trajectories obtained from the |
| 86 |
< |
simulation is discussed in Sec.~\ref{oopseSec:props}. Program design |
| 87 |
< |
considerations are presented in Sec.~\ref{oopseSec:design}. And |
| 88 |
< |
lastly, Sec.~\ref{oopseSec:conclusion} concludes the chapter. |
| 85 |
> |
A simulation in {\sc oopse} is built using a few fundamental |
| 86 |
> |
conceptual building blocks most of which are chemically intuitive. |
| 87 |
> |
The basic unit of a simulation is an {\tt atom}. The parameters |
| 88 |
> |
describing an {\tt atom} have been generalized to make it as flexible |
| 89 |
> |
as possible; this means that in addition to translational degrees of |
| 90 |
> |
freedom, {\tt Atoms} may also have {\it orientational} degrees of freedom. |
| 91 |
|
|
| 92 |
< |
\section{\label{oopseSec:empiricalEnergy}The Empirical Energy Functions} |
| 92 |
> |
The fundamental (static) properties of {\tt atoms} are defined by the |
| 93 |
> |
{\tt forceField} chosen for the simulation. The atomic properties |
| 94 |
> |
specified by a {\tt forceField} might include (but are not limited to) |
| 95 |
> |
charge, $\sigma$ and $\epsilon$ values for Lennard-Jones interactions, |
| 96 |
> |
the strength of the dipole moment ($\mu$), the mass, and the moments |
| 97 |
> |
of inertia. Other more complicated properties of atoms might also be |
| 98 |
> |
specified by the {\tt forceField}. |
| 99 |
|
|
| 100 |
< |
\subsection{\label{oopseSec:atomsMolecules}Atoms, Molecules and Rigid Bodies} |
| 100 |
> |
{\tt Atoms} can be grouped together in many ways. A {\tt rigidBody} |
| 101 |
> |
contains atoms that exert no forces on one another and which move as a |
| 102 |
> |
single rigid unit. A {\tt cutoffGroup} may contain atoms which |
| 103 |
> |
function together as a (rigid {\it or} non-rigid) unit for potential |
| 104 |
> |
energy calculations, |
| 105 |
> |
\begin{equation} |
| 106 |
> |
V_{ab} = s(r_{ab}) \sum_{i \in a} \sum_{j \in b} V_{ij}(r_{ij}) |
| 107 |
> |
\end{equation} |
| 108 |
> |
Here, $a$ and $b$ are two {\tt cutoffGroups} containing multiple atoms |
| 109 |
> |
($a = \left\{i\right\}$ and $b = \left\{j\right\}$). $s(r_{ab})$ is a |
| 110 |
> |
generalized switching function which insures that the atoms in the two |
| 111 |
> |
{\tt cutoffGroups} are treated identically as the two groups enter or |
| 112 |
> |
leave an interaction region. |
| 113 |
|
|
| 114 |
< |
The basic unit of an {\sc oopse} simulation is the atom. The |
| 115 |
< |
parameters describing the atom are generalized to make the atom as |
| 116 |
< |
flexible a representation as possible. They may represent specific |
| 117 |
< |
atoms of an element, or be used for collections of atoms such as |
| 98 |
< |
methyl and carbonyl groups. The atoms are also capable of having |
| 99 |
< |
directional components associated with them (\emph{e.g.}~permanent |
| 100 |
< |
dipoles). Charges, permanent dipoles, and Lennard-Jones parameters for |
| 101 |
< |
a given atom type are set in the force field parameter files. |
| 114 |
> |
{\tt Atoms} may also be grouped in more traditional ways into {\tt |
| 115 |
> |
bonds}, {\tt bends}, and {\tt torsions}. These groupings allow the |
| 116 |
> |
correct choice of interaction parameters for short-range interactions |
| 117 |
> |
to be chosen from the definitions in the {\tt forceField}. |
| 118 |
|
|
| 119 |
< |
\begin{lstlisting}[float,caption={[Specifier for molecules and atoms] A sample specification of an Ar molecule},label=sch:AtmMole] |
| 119 |
> |
All of these groups of {\tt atoms} are brought together in the {\tt |
| 120 |
> |
molecule}, which is the fundamental structure for setting up and {\sc |
| 121 |
> |
oopse} simulation. {\tt Molecules} contain lists of {\tt atoms} |
| 122 |
> |
followed by listings of the other atomic groupings ({\tt bonds}, {\tt |
| 123 |
> |
bends}, {\tt torsions}, {\tt rigidBodies}, and {\tt cutoffGroups}) |
| 124 |
> |
which relate the atoms to one another. |
| 125 |
> |
|
| 126 |
> |
Simulations often involve heterogeneous collections of molecules. To |
| 127 |
> |
specify a mixture of {\tt molecule} types, {\sc oopse} uses {\tt |
| 128 |
> |
components}. Even simulations containing only one type of molecule |
| 129 |
> |
must specify a single {\tt component}. |
| 130 |
> |
|
| 131 |
> |
Starting a simulation requires two types of information: {\it |
| 132 |
> |
meta-data}, which describes the types of objects present in the |
| 133 |
> |
simulation, and {\it configuration} information, which describes the |
| 134 |
> |
initial state of these objects. The meta-data is given to {\sc oopse} |
| 135 |
> |
using a C-based syntax that is parsed at the beginning of the |
| 136 |
> |
simulation. Configuration information is specified using an extended |
| 137 |
> |
XYZ file format. Both the meta-data and configuration file formats |
| 138 |
> |
are described in the following sections. |
| 139 |
> |
|
| 140 |
> |
\subsection{Meta-data Files} |
| 141 |
> |
|
| 142 |
> |
{\sc oopse} uses a C-based script syntax to parse the meta-data files |
| 143 |
> |
at run time. These files allow the user to completely describe the |
| 144 |
> |
system they wish to simulate, as well as tailor {\sc oopse}'s behavior |
| 145 |
> |
during the simulation. Meta-data files are typically denoted with the |
| 146 |
> |
extension {\tt .md} (which can stand for Meta-Data or Molecular |
| 147 |
> |
Dynamics or Molecule Definition depending on the user's mood). An |
| 148 |
> |
example meta-data file is shown in Scheme~\ref{sch:mdExample}. |
| 149 |
> |
|
| 150 |
> |
\begin{lstlisting}[float,caption={[An example of a complete meta-data |
| 151 |
> |
file] An example showing a complete meta-data |
| 152 |
> |
file.},label={sch:mdExample}] |
| 153 |
> |
|
| 154 |
> |
molecule{ |
| 155 |
> |
name = "Ar"; |
| 156 |
> |
nAtoms = 1; |
| 157 |
> |
atom[0]{ |
| 158 |
> |
type="Ar"; |
| 159 |
> |
position( 0.0, 0.0, 0.0 ); |
| 160 |
> |
} |
| 161 |
> |
} |
| 162 |
> |
|
| 163 |
> |
nComponents = 1; |
| 164 |
> |
component{ |
| 165 |
> |
type = "Ar"; |
| 166 |
> |
nMol = 108; |
| 167 |
> |
} |
| 168 |
> |
|
| 169 |
> |
initialConfig = "./argon.in"; |
| 170 |
> |
|
| 171 |
> |
forceField = "LJ"; |
| 172 |
> |
ensemble = "NVE"; // specify the simulation ensemble |
| 173 |
> |
dt = 1.0; // the time step for integration |
| 174 |
> |
runTime = 1e3; // the total simulation run time |
| 175 |
> |
sampleTime = 100; // trajectory file frequency |
| 176 |
> |
statusTime = 50; // statistics file frequency |
| 177 |
> |
|
| 178 |
> |
\end{lstlisting} |
| 179 |
> |
|
| 180 |
> |
Within the meta-data file it is necessary to provide a complete |
| 181 |
> |
description of the molecule before it is actually placed in the |
| 182 |
> |
simulation. {\sc oopse}'s meta-data syntax was originally developed |
| 183 |
> |
with this goal in mind, and allows for the use of {\it include files} |
| 184 |
> |
to specify all atoms in a molecular prototype, as well as any bonds, |
| 185 |
> |
bends, or torsions. Include files allow the user to describe a |
| 186 |
> |
molecular prototype once, then simply include it into each simulation |
| 187 |
> |
containing that molecule. Returning to the example in |
| 188 |
> |
Scheme~\ref{sch:mdExample}, the include file's contents would be |
| 189 |
> |
Scheme~\ref{sch:mdIncludeExample}, and the new meta-data file would |
| 190 |
> |
become Scheme~\ref{sch:mdExPrime}. |
| 191 |
> |
|
| 192 |
> |
\begin{lstlisting}[float,caption={An example molecule definition in an |
| 193 |
> |
include file.},label={sch:mdIncludeExample}] |
| 194 |
> |
|
| 195 |
|
molecule{ |
| 196 |
|
name = "Ar"; |
| 197 |
|
nAtoms = 1; |
| 200 |
|
position( 0.0, 0.0, 0.0 ); |
| 201 |
|
} |
| 202 |
|
} |
| 203 |
+ |
|
| 204 |
|
\end{lstlisting} |
| 205 |
|
|
| 206 |
+ |
\begin{lstlisting}[float,caption={Revised meta-data example.},label={sch:mdExPrime}] |
| 207 |
|
|
| 208 |
< |
Atoms can be collected into secondary structures such as rigid bodies |
| 116 |
< |
or molecules. The molecule is a way for {\sc oopse} to keep track of |
| 117 |
< |
the atoms in a simulation in logical manner. Molecular units store the |
| 118 |
< |
identities of all the atoms and rigid bodies associated with |
| 119 |
< |
themselves, and are responsible for the evaluation of their own |
| 120 |
< |
internal interactions (\emph{i.e.}~bonds, bends, and torsions). Scheme |
| 121 |
< |
\ref{sch:AtmMole} shows how one creates a molecule in a ``model'' or |
| 122 |
< |
\texttt{.mdl} file. The position of the atoms given in the |
| 123 |
< |
declaration are relative to the origin of the molecule, and is used |
| 124 |
< |
when creating a system containing the molecule. |
| 208 |
> |
#include "argon.md" |
| 209 |
|
|
| 210 |
< |
As stated previously, one of the features that sets {\sc oopse} apart |
| 211 |
< |
from most of the current molecular simulation packages is the ability |
| 212 |
< |
to handle rigid body dynamics. Rigid bodies are non-spherical |
| 213 |
< |
particles or collections of particles that have a constant internal |
| 210 |
> |
nComponents = 1; |
| 211 |
> |
component{ |
| 212 |
> |
type = "Ar"; |
| 213 |
> |
nMol = 108; |
| 214 |
> |
} |
| 215 |
> |
|
| 216 |
> |
initialConfig = "./argon.in"; |
| 217 |
> |
|
| 218 |
> |
forceField = "LJ"; |
| 219 |
> |
ensemble = "NVE"; |
| 220 |
> |
dt = 1.0; |
| 221 |
> |
runTime = 1e3; |
| 222 |
> |
sampleTime = 100; |
| 223 |
> |
statusTime = 50; |
| 224 |
> |
|
| 225 |
> |
\end{lstlisting} |
| 226 |
> |
|
| 227 |
> |
\subsection{\label{oopseSec:atomsMolecules}Atoms, Molecules, and other |
| 228 |
> |
ways of grouping atoms} |
| 229 |
> |
|
| 230 |
> |
As mentioned above, the fundamental unit for an {\sc oopse} simulation |
| 231 |
> |
is the {\tt atom}. Atoms can be collected into secondary structures |
| 232 |
> |
such as {\tt rigidBodies}, {\tt cutoffGroups}, or {\tt molecules}. The |
| 233 |
> |
{\tt molecule} is a way for {\sc oopse} to keep track of the atoms in |
| 234 |
> |
a simulation in logical manner. Molecular units store the identities |
| 235 |
> |
of all the atoms and rigid bodies associated with themselves, and they |
| 236 |
> |
are responsible for the evaluation of their own internal interactions |
| 237 |
> |
(\emph{i.e.}~bonds, bends, and torsions). Scheme |
| 238 |
> |
\ref{sch:mdIncludeExample} shows how one creates a molecule in an |
| 239 |
> |
included meta-data file. The positions of the atoms given in the |
| 240 |
> |
declaration are relative to the origin of the molecule, and the origin |
| 241 |
> |
is used when creating a system containing the molecule. |
| 242 |
> |
|
| 243 |
> |
One of the features that sets {\sc oopse} apart from most of the |
| 244 |
> |
current molecular simulation packages is the ability to handle rigid |
| 245 |
> |
body dynamics. Rigid bodies are non-spherical particles or collections |
| 246 |
> |
of particles (e.g. $\mbox{C}_{60}$) that have a constant internal |
| 247 |
|
potential and move collectively.\cite{Goldstein01} They are not |
| 248 |
|
included in most simulation packages because of the algorithmic |
| 249 |
< |
complexity involved in propagating orientational degrees of |
| 250 |
< |
freedom. Until recently, integrators which propagate orientational |
| 251 |
< |
motion have been much worse than those available for translational |
| 252 |
< |
motion. |
| 249 |
> |
complexity involved in propagating orientational degrees of freedom. |
| 250 |
> |
Integrators which propagate orientational motion with an acceptable |
| 251 |
> |
level of energy conservation for molecular dynamics are relatively |
| 252 |
> |
new inventions. |
| 253 |
|
|
| 254 |
|
Moving a rigid body involves determination of both the force and |
| 255 |
|
torque applied by the surroundings, which directly affect the |
| 259 |
|
the rigid body is simply the sum of these external forces. |
| 260 |
|
Accumulation of the total torque on the rigid body is more complex |
| 261 |
|
than the force because the torque is applied to the center of mass of |
| 262 |
< |
the rigid body. The torque on rigid body $i$ is |
| 262 |
> |
the rigid body. The space-fixed torque on rigid body $i$ is |
| 263 |
|
\begin{equation} |
| 264 |
|
\boldsymbol{\tau}_i= |
| 265 |
|
\sum_{a}\biggl[(\mathbf{r}_{ia}-\mathbf{r}_i)\times \mathbf{f}_{ia} |
| 281 |
|
$\psi$.\cite{Goldstein01} In order to avoid numerical instabilities |
| 282 |
|
inherent in using the Euler angles, the four parameter ``quaternion'' |
| 283 |
|
scheme is often used. The elements of $\mathsf{A}$ can be expressed as |
| 284 |
< |
arithmetic operations involving the four quaternions ($q_0, q_1, q_2,$ |
| 285 |
< |
and $q_3$).\cite{allen87:csl} Use of quaternions also leads to |
| 284 |
> |
arithmetic operations involving the four quaternions ($q_w, q_x, q_y,$ |
| 285 |
> |
and $q_z$).\cite{allen87:csl} Use of quaternions also leads to |
| 286 |
|
performance enhancements, particularly for very small |
| 287 |
|
systems.\cite{Evans77} |
| 288 |
|
|
| 289 |
< |
{\sc oopse} utilizes a relatively new scheme that propagates the |
| 290 |
< |
entire nine parameter rotation matrix. Further discussion |
| 291 |
< |
on this choice can be found in Sec.~\ref{oopseSec:integrate}. An example |
| 292 |
< |
definition of a rigid body can be seen in Scheme |
| 293 |
< |
\ref{sch:rigidBody}. The positions in the atom definitions are the |
| 294 |
< |
placements of the atoms relative to the origin of the rigid body, |
| 178 |
< |
which itself has a position relative to the origin of the molecule. |
| 289 |
> |
Rather than use one of the previously stated methods, {\sc oopse} |
| 290 |
> |
utilizes a relatively new scheme that propagates the entire nine |
| 291 |
> |
parameter rotation matrix. Further discussion on this choice can be |
| 292 |
> |
found in Sec.~\ref{oopseSec:integrate}. An example definition of a |
| 293 |
> |
rigid body can be seen in Scheme |
| 294 |
> |
\ref{sch:rigidBody}. |
| 295 |
|
|
| 296 |
< |
\begin{lstlisting}[float,caption={[Defining rigid bodies]A sample definition of a rigid body},label={sch:rigidBody}] |
| 296 |
> |
\begin{lstlisting}[float,caption={[Defining rigid bodies]A sample |
| 297 |
> |
definition of a molecule containing a rigid body and a cutoff |
| 298 |
> |
group},label={sch:rigidBody}] |
| 299 |
|
molecule{ |
| 300 |
|
name = "TIP3P"; |
| 301 |
|
nAtoms = 3; |
| 317 |
|
nMembers = 3; |
| 318 |
|
members(0, 1, 2); |
| 319 |
|
} |
| 320 |
+ |
|
| 321 |
+ |
nCutoffGroups = 1; |
| 322 |
+ |
cutoffGroup[0]{ |
| 323 |
+ |
nMembers = 3; |
| 324 |
+ |
members(0, 1, 2); |
| 325 |
+ |
} |
| 326 |
|
} |
| 327 |
|
\end{lstlisting} |
| 328 |
|
|
| 329 |
+ |
\subsection{\label{sec:miscConcepts}Creating a Metadata File} |
| 330 |
+ |
|
| 331 |
+ |
The actual creation of a metadata file requires several key |
| 332 |
+ |
components. The first part of the file needs to be the declaration of |
| 333 |
+ |
all of the molecule prototypes used in the simulation. This is |
| 334 |
+ |
typically done through included meta-data files. Only the molecules |
| 335 |
+ |
actually present in the simulation need to be declared; however, {\sc |
| 336 |
+ |
oopse} allows for the declaration of more molecules than are |
| 337 |
+ |
needed. This gives the user the ability to build up a library of |
| 338 |
+ |
commonly used molecules into a single include file. |
| 339 |
+ |
|
| 340 |
+ |
Once all prototypes are declared, the ordering of the rest of the |
| 341 |
+ |
script is less stringent. The molecular composition of the simulation |
| 342 |
+ |
is specified with {\tt component} statements. Each different type of |
| 343 |
+ |
molecule present in the simulation is considered a separate |
| 344 |
+ |
component. The number of components must be declared before the first |
| 345 |
+ |
component block statement (an example is shown in |
| 346 |
+ |
Sch.~\ref{sch:mdExPrime}). The component blocks tell {\sc oopse} the |
| 347 |
+ |
number of molecules that will be in the simulation, and the order in |
| 348 |
+ |
which the components blocks are declared sets the ordering of the real |
| 349 |
+ |
atoms in the configuration file as well as in the output files. The |
| 350 |
+ |
remainder of the script then sets the various simulation parameters |
| 351 |
+ |
for the system of interest. |
| 352 |
+ |
|
| 353 |
+ |
The required set of parameters that must be present in all simulations |
| 354 |
+ |
is given in Table~\ref{table:reqParams}. Since the user can use {\sc |
| 355 |
+ |
oopse} to perform energy minimizations as well as molecular dynamics |
| 356 |
+ |
simulations, one of the {\tt minimizer} or {\tt ensemble} keywords |
| 357 |
+ |
must be present. The {\tt ensemble} keyword is responsible for |
| 358 |
+ |
selecting the integration method used for the calculation of the |
| 359 |
+ |
equations of motion. An in depth discussion of the various methods |
| 360 |
+ |
available in {\sc oopse} can be found in |
| 361 |
+ |
Sec.~\ref{oopseSec:mechanics}. The {\tt minimizer} keyword selects |
| 362 |
+ |
which minimization method to use, and more details on the choices of |
| 363 |
+ |
minimizer parameters can be found in |
| 364 |
+ |
Sec.~\ref{oopseSec:minimizer}. The {\tt forceField} statement is |
| 365 |
+ |
important for the selection of which forces will be used in the course |
| 366 |
+ |
of the simulation. {\sc oopse} supports several force fields, as |
| 367 |
+ |
outlined in Sec.~\ref{oopseSec:empiricalEnergy}. The force fields are |
| 368 |
+ |
interchangeable between simulations, with the only requirement being |
| 369 |
+ |
that all atoms needed by the simulation are defined within the |
| 370 |
+ |
selected force field. |
| 371 |
+ |
|
| 372 |
+ |
For molecular dynamics simulations, the time step between force |
| 373 |
+ |
evaluations is set with the {\tt dt} parameter, and {\tt runTime} will |
| 374 |
+ |
set the time length of the simulation. Note, that {\tt runTime} is an |
| 375 |
+ |
absolute time, meaning if the simulation is started at t = 10.0~ns |
| 376 |
+ |
with a {\tt runTime} of 25.0~ns, the simulation will only run for an |
| 377 |
+ |
additional 15.0~ns. |
| 378 |
+ |
|
| 379 |
+ |
For energy minimizations, it is not necessary to specify {\tt dt} or |
| 380 |
+ |
{\tt runTime}. |
| 381 |
+ |
|
| 382 |
+ |
The final required parameter is the {\tt initialConfig} |
| 383 |
+ |
statement. This will set the initial coordinates for the system, as |
| 384 |
+ |
well as the initial time if the {\tt useInitalTime} flag is set to |
| 385 |
+ |
{\tt true}. The format of the file specified in {\tt initialConfig}, |
| 386 |
+ |
is given in Sec.~\ref{oopseSec:coordFiles}. Additional parameters are |
| 387 |
+ |
summarized in Table~\ref{table:genParams}. |
| 388 |
+ |
|
| 389 |
+ |
It is important to note the fundamental units in all files which are |
| 390 |
+ |
read and written by {\sc oopse}. Energies are in $\mbox{kcal |
| 391 |
+ |
mol}^{-1}$, distances are in $\mbox{\AA}$, times are in $\mbox{fs}$, |
| 392 |
+ |
translational velocities are in $\mbox{\AA~fs}^{-1}$, and masses are |
| 393 |
+ |
in $\mbox{amu}$. Orientational degrees of freedom are described using |
| 394 |
+ |
quaternions (unitless, but $q_w^2 + q_x^2 + q_y^2 + q_z^2 = 1$), |
| 395 |
+ |
body-fixed angular momenta ($\mbox{amu \AA}^{2} \mbox{radians |
| 396 |
+ |
fs}^{-1}$), and body-fixed moments of inertia ($\mbox{amu \AA}^{2}$). |
| 397 |
+ |
|
| 398 |
+ |
\begin{table} |
| 399 |
+ |
\caption{Meta-data Keywords: Required Parameters} |
| 400 |
+ |
\label{table:reqParams} |
| 401 |
+ |
\begin{center} |
| 402 |
+ |
% Note when adding or removing columns, the \hsize numbers must add up to the total number |
| 403 |
+ |
% of columns. |
| 404 |
+ |
\begin{tabularx}{\linewidth}% |
| 405 |
+ |
{>{\setlength{\hsize}{1.00\hsize}}X% |
| 406 |
+ |
>{\setlength{\hsize}{0.4\hsize}}X% |
| 407 |
+ |
>{\setlength{\hsize}{1.2\hsize}}X% |
| 408 |
+ |
>{\setlength{\hsize}{1.4\hsize}}X} |
| 409 |
+ |
|
| 410 |
+ |
{\bf keyword} & {\bf units} & {\bf use} & {\bf remarks} \\ \hline |
| 411 |
+ |
|
| 412 |
+ |
{\tt forceField} & string & Sets the force field. & Possible force fields are DUFF, LJ, and EAM. \\ |
| 413 |
+ |
{\tt nComponents} & integer & Sets the number of components. & Needs to appear before the first {\tt Component} block. \\ |
| 414 |
+ |
{\tt initialConfig} & string & Sets the file containing the initial configuration. & Can point to any file containing the configuration in the correct order. \\ |
| 415 |
+ |
{\tt minimizer}& string & Chooses a minimizer & Possible minimizers |
| 416 |
+ |
are SD and CG. Either {\tt ensemble} or {\tt minimizer} must be specified. \\ |
| 417 |
+ |
{\tt ensemble} & string & Sets the ensemble. & Possible ensembles are |
| 418 |
+ |
NVE, NVT, NPTi, NPTf, and NPTxyz. Either {\tt ensemble} |
| 419 |
+ |
or {\tt minimizer} must be specified. \\ |
| 420 |
+ |
{\tt dt} & fs & Sets the time step. & Selection of {\tt dt} should be |
| 421 |
+ |
small enough to sample the fastest motion of the simulation. ({\tt |
| 422 |
+ |
dt} is required for molecular dynamics simulations)\\ |
| 423 |
+ |
{\tt runTime} & fs & Sets the time at which the simulation should |
| 424 |
+ |
end. & This is an absolute time, and will end the simulation when the |
| 425 |
+ |
current time meets or exceeds the {\tt runTime}. ({\tt runTime} is |
| 426 |
+ |
required for molecular dynamics simulations)\\ |
| 427 |
+ |
|
| 428 |
+ |
\end{tabularx} |
| 429 |
+ |
\end{center} |
| 430 |
+ |
\end{table} |
| 431 |
+ |
|
| 432 |
+ |
\begin{table} |
| 433 |
+ |
\caption{Meta-data Keywords: General Parameters} |
| 434 |
+ |
\label{table:genParams} |
| 435 |
+ |
\begin{center} |
| 436 |
+ |
% Note when adding or removing columns, the \hsize numbers must add up to the total number |
| 437 |
+ |
% of columns. |
| 438 |
+ |
\begin{tabularx}{\linewidth}% |
| 439 |
+ |
{>{\setlength{\hsize}{1.00\hsize}}X% |
| 440 |
+ |
>{\setlength{\hsize}{0.4\hsize}}X% |
| 441 |
+ |
>{\setlength{\hsize}{1.2\hsize}}X% |
| 442 |
+ |
>{\setlength{\hsize}{1.4\hsize}}X} |
| 443 |
+ |
|
| 444 |
+ |
{\bf keyword} & {\bf units} & {\bf use} & {\bf remarks} \\ \hline |
| 445 |
+ |
|
| 446 |
+ |
{\tt finalConfig} & string & Sets the name of the final |
| 447 |
+ |
output file. & Useful when stringing simulations together. Defaults to |
| 448 |
+ |
the root name of the initial meta-data file but with an {\tt .eor} |
| 449 |
+ |
extension. \\ |
| 450 |
+ |
{\tt useInitialTime} & logical & Sets whether the initial time is taken from the {\tt .in} file. & Useful when recovering a simulation from a crashed processor. Default is false. \\ |
| 451 |
+ |
{\tt sampleTime} & fs & Sets the frequency at which the {\tt .dump} |
| 452 |
+ |
file is written. & The default is equal to the {\tt runTime}. \\ |
| 453 |
+ |
{\tt statusTime} & fs & Sets the frequency at which the {\tt .stat} |
| 454 |
+ |
file is written. & The default is equal to the {\tt sampleTime}. \\ |
| 455 |
+ |
{\tt cutoffRadius} & $\mbox{\AA}$ & Manually sets the cutoff radius & |
| 456 |
+ |
Defaults to $15\mbox{\AA}$ for systems containing charges or dipoles or to $2.5 |
| 457 |
+ |
\sigma_{L}$, where $\sigma_{L}$ is the largest LJ $\sigma$ in the |
| 458 |
+ |
simulation. \\ |
| 459 |
+ |
{\tt switchingRadius} & $\mbox{\AA}$ & Manually sets the inner radius for the switching function. & Defaults to 95~\% of the {\tt cutoffRadius}. \\ |
| 460 |
+ |
{\tt useReactionField} & logical & Turns the reaction field correction on/off. & Default is false. \\ |
| 461 |
+ |
{\tt dielectric} & unitless & Sets the dielectric constant for reaction field. & If {\tt useReactionField} is true, then {\tt dielectric} must be set. \\ |
| 462 |
+ |
{\tt usePeriodicBoundaryConditions} & & & \\ |
| 463 |
+ |
& logical & Turns periodic boundary conditions on/off. & Default is true. \\ |
| 464 |
+ |
{\tt seed } & integer & Sets the seed value for the random number |
| 465 |
+ |
generator. & The seed needs to be at least 9 digits long. The default |
| 466 |
+ |
is to take the seed from the CPU clock. \\ |
| 467 |
+ |
{\tt forceFieldVariant} & string & Sets the name of the variant of the |
| 468 |
+ |
force field. & {\sc eam} has three variants: {\tt u3}, {\tt u6}, and |
| 469 |
+ |
{\tt VC}. |
| 470 |
+ |
|
| 471 |
+ |
\end{tabularx} |
| 472 |
+ |
\end{center} |
| 473 |
+ |
\end{table} |
| 474 |
+ |
|
| 475 |
+ |
|
| 476 |
+ |
\subsection{\label{oopseSec:coordFiles}Coordinate Files} |
| 477 |
+ |
|
| 478 |
+ |
The standard format for storage of a systems coordinates is a modified |
| 479 |
+ |
xyz-file syntax, the exact details of which can be seen in |
| 480 |
+ |
Scheme~\ref{sch:dumpFormat}. As all bonding and molecular information |
| 481 |
+ |
is stored in the meta-data files, the coordinate files contain only |
| 482 |
+ |
the coordinates of the objects which move independently during the |
| 483 |
+ |
simulation. It is important to note that {\it not all atoms} are |
| 484 |
+ |
capable of independent motion. Atoms which are part of rigid bodies |
| 485 |
+ |
are not ``integrable objects'' in the equations of motion; the rigid |
| 486 |
+ |
bodies themselves are the integrable objects. Therefore, the |
| 487 |
+ |
coordinate file contains coordinates of all the {\tt |
| 488 |
+ |
integrableObjects} in the system. For systems without rigid bodies, |
| 489 |
+ |
this is simply the coordinates of all the atoms. |
| 490 |
+ |
|
| 491 |
+ |
It is important to note that although the simulation propagates the |
| 492 |
+ |
complete rotation matrix, directional entities are written out using |
| 493 |
+ |
quaternions to save space in the output files. All objects (atoms, |
| 494 |
+ |
orientational atoms, and rigid bodies) are given quaternions and |
| 495 |
+ |
angular momenta in coordinate files which are output by OOPSE, but it |
| 496 |
+ |
is not necessary for the user to specify the quaternions or angular |
| 497 |
+ |
momenta for atoms without orientational degrees of freedom. |
| 498 |
+ |
|
| 499 |
+ |
\begin{lstlisting}[float,caption={[The format of the coordinate |
| 500 |
+ |
files] An example of the format of the coordinate files. The fist line |
| 501 |
+ |
is the number of {\tt integrableObjects} (freely-moving atoms and |
| 502 |
+ |
rigid bodies). The second line begins with the time stamp followed by |
| 503 |
+ |
the three $\mathsf{H}$ column vectors. It is important to note that |
| 504 |
+ |
for extended system ensembles, additional information pertinent to the |
| 505 |
+ |
integrators may be stored on this line as well. The next lines are the |
| 506 |
+ |
coordinates for all integrable objects in the system. The name of the |
| 507 |
+ |
integrable object is followed by position, velocity, quaternions, and |
| 508 |
+ |
lastly, body fixed angular momentum.},label=sch:dumpFormat] |
| 509 |
+ |
|
| 510 |
+ |
nIntegrable |
| 511 |
+ |
time; Hxx Hyx Hzx; Hxy Hyy Hzy; Hxz Hyz Hzz; |
| 512 |
+ |
Name1 x y z vx vy vz qw qx qy qz jx jy jz |
| 513 |
+ |
Name2 x y z vx vy vz qw qx qy qz jx jy jz |
| 514 |
+ |
etc... |
| 515 |
+ |
|
| 516 |
+ |
\end{lstlisting} |
| 517 |
+ |
|
| 518 |
+ |
The {\tt name} field for atoms is simply the atom type as specified in |
| 519 |
+ |
the meta-data file. The {\tt name} field for a rigid body is |
| 520 |
+ |
specified as {\tt MOLTYPE\_RB\_N}, to specify that this is {\tt |
| 521 |
+ |
rigidBody} N in a molecule of type MOLTYPE. In simulations with rigid |
| 522 |
+ |
body models of water, a sample coordinate line might be: |
| 523 |
+ |
|
| 524 |
+ |
\begin{tt} |
| 525 |
+ |
TIP3P\_RB\_0 x y z vx vy vz qw qx qy qz jx jy jz |
| 526 |
+ |
\end{tt} |
| 527 |
+ |
|
| 528 |
+ |
which tells the program that the rigid body representing a TIP3P |
| 529 |
+ |
molecule (rigid body \# 0) is listed on that line. |
| 530 |
+ |
|
| 531 |
+ |
There are three files used by {\sc oopse} which are written in the |
| 532 |
+ |
coordinate format. They are: the initial coordinate file |
| 533 |
+ |
(\texttt{.in}), the simulation trajectory file (\texttt{.dump}), and |
| 534 |
+ |
the final coordinates or ``end-of-run'' for the simulation |
| 535 |
+ |
(\texttt{.eor}). The initial coordinate file is necessary for {\sc |
| 536 |
+ |
oopse} to start the simulation with the proper coordinates, and this |
| 537 |
+ |
file must be generated by the user before the simulation run. The |
| 538 |
+ |
trajectory (or ``dump'') file is updated during simulation and is used |
| 539 |
+ |
to store snapshots of the coordinates at regular intervals. The first |
| 540 |
+ |
frame is a duplication of the |
| 541 |
+ |
\texttt{.in} file, and each subsequent frame is appended to the file |
| 542 |
+ |
at an interval specified in the meta-data file with the |
| 543 |
+ |
\texttt{sampleTime} flag. The final coordinate file is the |
| 544 |
+ |
``end-of-run'' file. The \texttt{.eor} file stores the final |
| 545 |
+ |
configuration of the system for a given simulation. The file is |
| 546 |
+ |
updated at the same time as the \texttt{.dump} file, but it only |
| 547 |
+ |
contains the most recent frame. In this way, an \texttt{.eor} file may |
| 548 |
+ |
be used to initialize a second simulation should it be necessary to |
| 549 |
+ |
recover from a crash or power outage. |
| 550 |
+ |
|
| 551 |
+ |
\subsection{\label{oopseSec:initCoords}Generation of Initial Coordinates} |
| 552 |
+ |
|
| 553 |
+ |
As was stated in Sec.~\ref{oopseSec:coordFiles}, an initial coordinate |
| 554 |
+ |
file is needed to provide the starting coordinates for a simulation. |
| 555 |
+ |
Since each simulation is different, system creation is left to the end |
| 556 |
+ |
user; however, we have included a few sample programs which make some |
| 557 |
+ |
specialized structures. The {\tt .in} file must list the integrable |
| 558 |
+ |
objects in the correct order. The ordering of the integrable objects |
| 559 |
+ |
relies on the ordering of molecules within the meta-data file. {\sc |
| 560 |
+ |
oopse} expects the order to comply with the following guidelines: |
| 561 |
+ |
\begin{enumerate} |
| 562 |
+ |
\item All of the molecules of the first declared component are given |
| 563 |
+ |
before proceeding to the molecules of the second component, and so on |
| 564 |
+ |
for all subsequently declared components. |
| 565 |
+ |
\item The ordering of the atoms for each molecule follows the order |
| 566 |
+ |
declared in the molecule's declaration within the model file. |
| 567 |
+ |
\item Only atoms which are not members of a {\tt rigidBody} are |
| 568 |
+ |
included. |
| 569 |
+ |
\item Rigid Body coordinates for a molecule are listed immediately |
| 570 |
+ |
after the the other atoms in a molecule. Some molecules may be |
| 571 |
+ |
entirely rigid, in which case, only the rigid body coordinates are |
| 572 |
+ |
given. |
| 573 |
+ |
\end{enumerate} |
| 574 |
+ |
An example is given in the meta-data file in Scheme~\ref{sch:initEx1} |
| 575 |
+ |
which results in the {\tt .in} file shown in Scheme~\ref{sch:initEx2}. |
| 576 |
+ |
|
| 577 |
+ |
\begin{lstlisting}[float,caption={Example declaration of the |
| 578 |
+ |
$\text{I}_2$ molecule and the HCl molecule. The two molecules are then |
| 579 |
+ |
included into a simulation.}, label=sch:initEx1] |
| 580 |
+ |
|
| 581 |
+ |
molecule{ |
| 582 |
+ |
name = "I2"; |
| 583 |
+ |
nAtoms = 2; |
| 584 |
+ |
atom[0]{ |
| 585 |
+ |
type = "I"; |
| 586 |
+ |
} |
| 587 |
+ |
atom[1]{ |
| 588 |
+ |
type = "I"; |
| 589 |
+ |
} |
| 590 |
+ |
nBonds = 1; |
| 591 |
+ |
bond[0]{ |
| 592 |
+ |
members( 0, 1); |
| 593 |
+ |
} |
| 594 |
+ |
} |
| 595 |
+ |
|
| 596 |
+ |
molecule{ |
| 597 |
+ |
name = "HCl" |
| 598 |
+ |
nAtoms = 2; |
| 599 |
+ |
atom[0]{ |
| 600 |
+ |
type = "H"; |
| 601 |
+ |
} |
| 602 |
+ |
atom[1]{ |
| 603 |
+ |
type = "Cl"; |
| 604 |
+ |
} |
| 605 |
+ |
nBonds = 1; |
| 606 |
+ |
bond[0]{ |
| 607 |
+ |
members( 0, 1); |
| 608 |
+ |
} |
| 609 |
+ |
} |
| 610 |
+ |
|
| 611 |
+ |
nComponents = 2; |
| 612 |
+ |
component{ |
| 613 |
+ |
type = "HCl"; |
| 614 |
+ |
nMol = 4; |
| 615 |
+ |
} |
| 616 |
+ |
component{ |
| 617 |
+ |
type = "I2"; |
| 618 |
+ |
nMol = 1; |
| 619 |
+ |
} |
| 620 |
+ |
|
| 621 |
+ |
initialConfig = "mixture.in"; |
| 622 |
+ |
|
| 623 |
+ |
\end{lstlisting} |
| 624 |
+ |
|
| 625 |
+ |
\begin{lstlisting}[float,caption={The contents of the {\tt |
| 626 |
+ |
mixture.in} file matching the declarations in |
| 627 |
+ |
Scheme~\ref{sch:initEx1}. Note that even though $\text{I}_2$ is |
| 628 |
+ |
declared before HCl, the {\tt .in} file follows the order {\it in |
| 629 |
+ |
which the components were included}.},label=sch:initEx2] |
| 630 |
+ |
|
| 631 |
+ |
10 |
| 632 |
+ |
0.0; 10.0 0.0 0.0; 0.0 10.0 0.0; 0.0 0.0 10.0; |
| 633 |
+ |
H ... |
| 634 |
+ |
Cl ... |
| 635 |
+ |
H ... |
| 636 |
+ |
Cl ... |
| 637 |
+ |
H ... |
| 638 |
+ |
Cl ... |
| 639 |
+ |
H ... |
| 640 |
+ |
Cl ... |
| 641 |
+ |
I ... |
| 642 |
+ |
I ... |
| 643 |
+ |
|
| 644 |
+ |
\end{lstlisting} |
| 645 |
+ |
|
| 646 |
+ |
|
| 647 |
+ |
\subsection{The Statistics File} |
| 648 |
+ |
|
| 649 |
+ |
The last output file generated by {\sc oopse} is the statistics |
| 650 |
+ |
file. This file records such statistical quantities as the |
| 651 |
+ |
instantaneous temperature (in $K$), volume (in $\mbox{\AA}^{3}$), |
| 652 |
+ |
pressure (in $\mbox{atm}$), etc. It is written out with the frequency |
| 653 |
+ |
specified in the meta-data file with the |
| 654 |
+ |
\texttt{statusTime} keyword. The file allows the user to observe the |
| 655 |
+ |
system variables as a function of simulation time while the simulation |
| 656 |
+ |
is in progress. One useful function the statistics file serves is to |
| 657 |
+ |
monitor the conserved quantity of a given simulation ensemble, |
| 658 |
+ |
allowing the user to gauge the stability of the integrator. The |
| 659 |
+ |
statistics file is denoted with the \texttt{.stat} file extension. |
| 660 |
+ |
|
| 661 |
+ |
\section{\label{oopseSec:empiricalEnergy}The Empirical Energy |
| 662 |
+ |
Functions} |
| 663 |
+ |
|
| 664 |
+ |
Like many simulation packages, {\sc oopse} splits the potential energy |
| 665 |
+ |
into the short-ranged (bonded) portion and a long-range (non-bonded) |
| 666 |
+ |
potential, |
| 667 |
+ |
\begin{equation} |
| 668 |
+ |
V = V_{\mathrm{short-range}} + V_{\mathrm{long-range}}. |
| 669 |
+ |
\end{equation} |
| 670 |
+ |
The short-ranged portion includes the explicit bonds, bends, and |
| 671 |
+ |
torsions which have been defined in the meta-data file for the |
| 672 |
+ |
molecules which are present in the simulation. The functional forms and |
| 673 |
+ |
parameters for these interactions are defined by the force field which |
| 674 |
+ |
is chosen. |
| 675 |
+ |
|
| 676 |
+ |
Calculating the long-range (non-bonded) potential involves a sum over |
| 677 |
+ |
all pairs of atoms (except for those atoms which are involved in a |
| 678 |
+ |
bond, bend, or torsion with each other). If done poorly, calculating |
| 679 |
+ |
the the long-range interactions for $N$ atoms would involve $N(N-1)/2$ |
| 680 |
+ |
evaluations of atomic distances. To reduce the number of distance |
| 681 |
+ |
evaluations between pairs of atoms, {\sc oopse} uses a switched cutoff |
| 682 |
+ |
with Verlet neighbor lists.\cite{allen87:csl} It is well known that |
| 683 |
+ |
neutral groups which contain charges will exhibit pathological forces |
| 684 |
+ |
unless the cutoff is applied to the neutral groups evenly instead of |
| 685 |
+ |
to the individual atoms.\cite{leach01:mm} {\sc oopse} allows users to |
| 686 |
+ |
specify cutoff groups which may contain an arbitrary number of atoms |
| 687 |
+ |
in the molecule. Atoms in a cutoff group are treated as a single unit |
| 688 |
+ |
for the evaluation of the switching function: |
| 689 |
+ |
\begin{equation} |
| 690 |
+ |
V_{\mathrm{long-range}} = \sum_{a} \sum_{b>a} s(r_{ab}) \sum_{i \in a} \sum_{j \in b} V_{ij}(r_{ij}), |
| 691 |
+ |
\end{equation} |
| 692 |
+ |
where $r_{ab}$ is the distance between the centers of mass of the two |
| 693 |
+ |
cutoff groups ($a$ and $b$). |
| 694 |
+ |
|
| 695 |
+ |
The sums over $a$ and $b$ are over the cutoff groups that are present |
| 696 |
+ |
in the simulation. Atoms which are not explicitly defined as members |
| 697 |
+ |
of a {\tt cutoffGroup} are treated as a group consisting of only one |
| 698 |
+ |
atom. The switching function, $s(r)$ is the standard cubic switching |
| 699 |
+ |
function, |
| 700 |
+ |
\begin{equation} |
| 701 |
+ |
S(r) = |
| 702 |
+ |
\begin{cases} |
| 703 |
+ |
1 & \text{if $r \le r_{\text{sw}}$},\\ |
| 704 |
+ |
\frac{(r_{\text{cut}} + 2r - 3r_{\text{sw}})(r_{\text{cut}} - r)^2} |
| 705 |
+ |
{(r_{\text{cut}} - r_{\text{sw}})^2} |
| 706 |
+ |
& \text{if $r_{\text{sw}} < r \le r_{\text{cut}}$}, \\ |
| 707 |
+ |
0 & \text{if $r > r_{\text{cut}}$.} |
| 708 |
+ |
\end{cases} |
| 709 |
+ |
\label{eq:dipoleSwitching} |
| 710 |
+ |
\end{equation} |
| 711 |
+ |
Here, $r_{\text{sw}}$ is the {\tt switchingRadius}, or the distance |
| 712 |
+ |
beyond which interactions are reduced, and $r_{\text{cut}}$ is the |
| 713 |
+ |
{\tt cutoffRadius}, or the distance at which interactions are |
| 714 |
+ |
truncated. |
| 715 |
+ |
|
| 716 |
+ |
Users of {\sc oopse} do not need to specify the {\tt cutoffRadius} or |
| 717 |
+ |
{\tt switchingRadius}. In simulations containing only Lennard-Jones |
| 718 |
+ |
atoms, the cutoff radius has a default value of $2.5\sigma_{ii}$, |
| 719 |
+ |
where $\sigma_{ii}$ is the largest Lennard-Jones length parameter |
| 720 |
+ |
present in the simulation. In simulations containing charged or |
| 721 |
+ |
dipolar atoms, the default cutoff radius is $15 \mbox{\AA}$. |
| 722 |
+ |
|
| 723 |
+ |
The {\tt switchingRadius} is set to a default value of 95\% of the |
| 724 |
+ |
{\tt cutoffRadius}. In the special case of a simulation containing |
| 725 |
+ |
{\it only} Lennard-Jones atoms, the default switching radius takes the |
| 726 |
+ |
same value as the cutoff radius, and {\sc oopse} will use a shifted |
| 727 |
+ |
potential to remove discontinuities in the potential at the cutoff. |
| 728 |
+ |
Both radii may be specified in the meta-data file. |
| 729 |
+ |
|
| 730 |
+ |
Force fields can be added to {\sc oopse}, although it comes with a few |
| 731 |
+ |
simple examples (Lennard-Jones, {\sc duff}, {\sc water}, and {\sc |
| 732 |
+ |
eam}) which are explained in the following sections. |
| 733 |
+ |
|
| 734 |
|
\subsection{\label{sec:LJPot}The Lennard Jones Force Field} |
| 735 |
|
|
| 736 |
|
The most basic force field implemented in {\sc oopse} is the |
| 737 |
< |
Lennard-Jones force field, which mimics the van der Waals interaction at |
| 738 |
< |
long distances, and uses an empirical repulsion at short |
| 737 |
> |
Lennard-Jones force field, which mimics the van der Waals interaction |
| 738 |
> |
at long distances and uses an empirical repulsion at short |
| 739 |
|
distances. The Lennard-Jones potential is given by: |
| 740 |
|
\begin{equation} |
| 741 |
|
V_{\text{LJ}}(r_{ij}) = |
| 748 |
|
where $r_{ij}$ is the distance between particles $i$ and $j$, |
| 749 |
|
$\sigma_{ij}$ scales the length of the interaction, and |
| 750 |
|
$\epsilon_{ij}$ scales the well depth of the potential. Scheme |
| 751 |
< |
\ref{sch:LJFF} gives an example \texttt{.bass} file that |
| 751 |
> |
\ref{sch:LJFF} gives an example meta-data file that |
| 752 |
|
sets up a system of 108 Ar particles to be simulated using the |
| 753 |
|
Lennard-Jones force field. |
| 754 |
|
|
| 755 |
< |
\begin{lstlisting}[float,caption={[Invocation of the Lennard-Jones force field] A sample system using the Lennard-Jones force field.},label={sch:LJFF}] |
| 755 |
> |
\begin{lstlisting}[float,caption={[Invocation of the Lennard-Jones |
| 756 |
> |
force field] A sample meta-data file for a small Lennard-Jones |
| 757 |
> |
simulation.},label={sch:LJFF}] |
| 758 |
|
|
| 759 |
< |
#include "argon.mdl" |
| 759 |
> |
#include "argon.md" |
| 760 |
|
|
| 761 |
|
nComponents = 1; |
| 762 |
|
component{ |
| 764 |
|
nMol = 108; |
| 765 |
|
} |
| 766 |
|
|
| 767 |
< |
initialConfig = "./argon.init"; |
| 767 |
> |
initialConfig = "./argon.in"; |
| 768 |
|
|
| 769 |
|
forceField = "LJ"; |
| 770 |
|
\end{lstlisting} |
| 771 |
|
|
| 241 |
– |
Because this potential is calculated between all pairs, the force |
| 242 |
– |
evaluation can become computationally expensive for large systems. To |
| 243 |
– |
keep the pair evaluations to a manageable number, {\sc oopse} employs |
| 244 |
– |
a cut-off radius.\cite{allen87:csl} The cutoff radius can either be |
| 245 |
– |
specified in the \texttt{.bass} file, or left as its default value of |
| 246 |
– |
$2.5\sigma_{ii}$, where $\sigma_{ii}$ is the largest Lennard-Jones |
| 247 |
– |
length parameter present in the simulation. Truncating the calculation |
| 248 |
– |
at $r_{\text{cut}}$ introduces a discontinuity into the potential |
| 249 |
– |
energy and the force. To offset this discontinuity in the potential, |
| 250 |
– |
the energy value at $r_{\text{cut}}$ is subtracted from the |
| 251 |
– |
potential. This causes the potential to go to zero smoothly at the |
| 252 |
– |
cut-off radius, and preserves conservation of energy in integrating |
| 253 |
– |
the equations of motion. There still remains a discontinuity in the derivative (the forces), however, this does not significantly affect the dynamics. |
| 254 |
– |
|
| 772 |
|
Interactions between dissimilar particles requires the generation of |
| 773 |
< |
cross term parameters for $\sigma$ and $\epsilon$. These are |
| 774 |
< |
calculated through the Lorentz-Berthelot mixing |
| 773 |
> |
cross term parameters for $\sigma$ and $\epsilon$. These parameters |
| 774 |
> |
are determined using the Lorentz-Berthelot mixing |
| 775 |
|
rules:\cite{allen87:csl} |
| 776 |
|
\begin{equation} |
| 777 |
|
\sigma_{ij} = \frac{1}{2}[\sigma_{ii} + \sigma_{jj}], |
| 786 |
|
\subsection{\label{oopseSec:DUFF}Dipolar Unified-Atom Force Field} |
| 787 |
|
|
| 788 |
|
The dipolar unified-atom force field ({\sc duff}) was developed to |
| 789 |
< |
simulate lipid bilayers. The simulations require a model capable of |
| 790 |
< |
forming bilayers, while still being sufficiently computationally |
| 791 |
< |
efficient to allow large systems ($\sim$100's of phospholipids, |
| 792 |
< |
$\sim$1000's of waters) to be simulated for long times |
| 793 |
< |
($\sim$10's of nanoseconds). |
| 789 |
> |
simulate lipid bilayers. These types of simulations require a model |
| 790 |
> |
capable of forming bilayers, while still being sufficiently |
| 791 |
> |
computationally efficient to allow large systems ($\sim$100's of |
| 792 |
> |
phospholipids, $\sim$1000's of waters) to be simulated for long times |
| 793 |
> |
($\sim$10's of nanoseconds). With this goal in mind, {\sc duff} has no |
| 794 |
> |
point charges. Charge-neutral distributions are replaced with dipoles, |
| 795 |
> |
while most atoms and groups of atoms are reduced to Lennard-Jones |
| 796 |
> |
interaction sites. This simplification reduces the length scale of |
| 797 |
> |
long range interactions from $\frac{1}{r}$ to $\frac{1}{r^3}$, |
| 798 |
> |
removing the need for the computationally expensive Ewald |
| 799 |
> |
sum. Instead, Verlet neighbor-lists and cutoff radii are used for the |
| 800 |
> |
dipolar interactions, and, if desired, a reaction field may be added |
| 801 |
> |
to mimic longer range interactions. |
| 802 |
|
|
| 278 |
– |
With this goal in mind, {\sc duff} has no point |
| 279 |
– |
charges. Charge-neutral distributions were replaced with dipoles, |
| 280 |
– |
while most atoms and groups of atoms were reduced to Lennard-Jones |
| 281 |
– |
interaction sites. This simplification cuts the length scale of long |
| 282 |
– |
range interactions from $\frac{1}{r}$ to $\frac{1}{r^3}$, and allows |
| 283 |
– |
us to avoid the computationally expensive Ewald sum. Instead, we can |
| 284 |
– |
use neighbor-lists and cutoff radii for the dipolar interactions, or |
| 285 |
– |
include a reaction field to mimic larger range interactions. |
| 286 |
– |
|
| 803 |
|
As an example, lipid head-groups in {\sc duff} are represented as |
| 804 |
< |
point dipole interaction sites. By placing a dipole at the head |
| 805 |
< |
group's center of mass, our model mimics the charge separation found |
| 806 |
< |
in common phospholipid head groups such as |
| 807 |
< |
phosphatidylcholine.\cite{Cevc87} Additionally, a large Lennard-Jones |
| 808 |
< |
site is located at the pseudoatom's center of mass. The model is |
| 809 |
< |
illustrated by the red atom in Fig.~\ref{oopseFig:lipidModel}. The |
| 810 |
< |
water model we use to complement the dipoles of the lipids is our |
| 811 |
< |
reparameterization of the soft sticky dipole (SSD) model of Ichiye |
| 804 |
> |
point dipole interaction sites. Placing a dipole at the head group's |
| 805 |
> |
center of mass mimics the charge separation found in common |
| 806 |
> |
phospholipid head groups such as phosphatidylcholine.\cite{Cevc87} |
| 807 |
> |
Additionally, a large Lennard-Jones site is located at the |
| 808 |
> |
pseudoatom's center of mass. The model is illustrated by the red atom |
| 809 |
> |
in Fig.~\ref{oopseFig:lipidModel}. The water model we use to |
| 810 |
> |
complement the dipoles of the lipids is a |
| 811 |
> |
reparameterization\cite{fennell04} of the soft sticky dipole (SSD) |
| 812 |
> |
model of Ichiye |
| 813 |
|
\emph{et al.}\cite{liu96:new_model} |
| 814 |
|
|
| 815 |
|
\begin{figure} |
| 816 |
|
\centering |
| 817 |
< |
\includegraphics[width=\linewidth]{twoChainFig.pdf} |
| 818 |
< |
\caption[A representation of a lipid model in {\sc duff}]{A representation of the lipid model. $\phi$ is the torsion angle, $\theta$ % |
| 819 |
< |
is the bend angle, and $\mu$ is the dipole moment of the head group.} |
| 817 |
> |
\includegraphics[width=\linewidth]{lipidModel.eps} |
| 818 |
> |
\caption[A representation of a lipid model in {\sc duff}]{A |
| 819 |
> |
representation of the lipid model. $\phi$ is the torsion angle, |
| 820 |
> |
$\theta$ is the bend angle, and $\mu$ is the dipole moment of the head |
| 821 |
> |
group.} |
| 822 |
|
\label{oopseFig:lipidModel} |
| 823 |
|
\end{figure} |
| 824 |
|
|
| 825 |
< |
We have used a set of scalable parameters to model the alkyl groups |
| 826 |
< |
with Lennard-Jones sites. For this, we have borrowed parameters from |
| 827 |
< |
the TraPPE force field of Siepmann |
| 828 |
< |
\emph{et al}.\cite{Siepmann1998} TraPPE is a unified-atom |
| 829 |
< |
representation of n-alkanes, which is parametrized against phase |
| 830 |
< |
equilibria using Gibbs ensemble Monte Carlo simulation |
| 831 |
< |
techniques.\cite{Siepmann1998} One of the advantages of TraPPE is that |
| 832 |
< |
it generalizes the types of atoms in an alkyl chain to keep the number |
| 833 |
< |
of pseudoatoms to a minimum; the parameters for a unified atom such as |
| 834 |
< |
$\text{CH}_2$ do not change depending on what species are bonded to |
| 316 |
< |
it. |
| 825 |
> |
A set of scalable parameters has been used to model the alkyl groups |
| 826 |
> |
with Lennard-Jones sites. For this, parameters from the TraPPE force |
| 827 |
> |
field of Siepmann \emph{et al.}\cite{Siepmann1998} have been |
| 828 |
> |
utilized. TraPPE is a unified-atom representation of n-alkanes which |
| 829 |
> |
is parametrized against phase equilibria using Gibbs ensemble Monte |
| 830 |
> |
Carlo simulation techniques.\cite{Siepmann1998} One of the advantages |
| 831 |
> |
of TraPPE is that it generalizes the types of atoms in an alkyl chain |
| 832 |
> |
to keep the number of pseudoatoms to a minimum; thus, the parameters |
| 833 |
> |
for a unified atom such as $\text{CH}_2$ do not change depending on |
| 834 |
> |
what species are bonded to it. |
| 835 |
|
|
| 836 |
< |
TraPPE also constrains all bonds to be of fixed length. Typically, |
| 837 |
< |
bond vibrations are the fastest motions in a molecular dynamic |
| 838 |
< |
simulation. Small time steps between force evaluations must be used to |
| 839 |
< |
ensure adequate energy conservation in the bond degrees of freedom. By |
| 840 |
< |
constraining the bond lengths, larger time steps may be used when |
| 841 |
< |
integrating the equations of motion. A simulation using {\sc duff} is |
| 842 |
< |
illustrated in Scheme \ref{sch:DUFF}. |
| 836 |
> |
As is required by TraPPE, {\sc duff} also constrains all bonds to be |
| 837 |
> |
of fixed length. Typically, bond vibrations are the fastest motions in |
| 838 |
> |
a molecular dynamic simulation. With these vibrations present, small |
| 839 |
> |
time steps between force evaluations must be used to ensure adequate |
| 840 |
> |
energy conservation in the bond degrees of freedom. By constraining |
| 841 |
> |
the bond lengths, larger time steps may be used when integrating the |
| 842 |
> |
equations of motion. A simulation using {\sc duff} is illustrated in |
| 843 |
> |
Scheme \ref{sch:DUFF}. |
| 844 |
|
|
| 845 |
< |
\begin{lstlisting}[float,caption={[Invocation of {\sc duff}]A portion of a \texttt{.bass} file showing a simulation utilizing {\sc duff}},label={sch:DUFF}] |
| 845 |
> |
\begin{lstlisting}[float,caption={[Invocation of {\sc duff}]A portion |
| 846 |
> |
of a meta-data file showing a simulation utilizing {\sc |
| 847 |
> |
duff}},label={sch:DUFF}] |
| 848 |
|
|
| 849 |
< |
#include "water.mdl" |
| 850 |
< |
#include "lipid.mdl" |
| 849 |
> |
#include "water.md" |
| 850 |
> |
#include "lipid.md" |
| 851 |
|
|
| 852 |
|
nComponents = 2; |
| 853 |
|
component{ |
| 860 |
|
nMol = 1936; |
| 861 |
|
} |
| 862 |
|
|
| 863 |
< |
initialConfig = "bilayer.init"; |
| 863 |
> |
initialConfig = "bilayer.in"; |
| 864 |
|
|
| 865 |
|
forceField = "DUFF"; |
| 866 |
|
|
| 867 |
|
\end{lstlisting} |
| 868 |
|
|
| 869 |
< |
\subsection{\label{oopseSec:energyFunctions}{\sc duff} Energy Functions} |
| 869 |
> |
\subsubsection{\label{oopseSec:energyFunctions}{\sc duff} Energy Functions} |
| 870 |
|
|
| 871 |
|
The total potential energy function in {\sc duff} is |
| 872 |
|
\begin{equation} |
| 886 |
|
\label{eq:internalPotential} |
| 887 |
|
\end{equation} |
| 888 |
|
Here $V_{\text{bend}}$ is the bend potential for all 1, 3 bonded pairs |
| 889 |
< |
within the molecule $I$, and $V_{\text{torsion}}$ is the torsion potential |
| 890 |
< |
for all 1, 4 bonded pairs. The pairwise portions of the internal |
| 891 |
< |
potential are excluded for atom pairs that are involved in the same bond, bend, or torsion. All other atom pairs within the molecule are subject to the LJ pair potential. |
| 889 |
> |
within the molecule $I$, and $V_{\text{torsion}}$ is the torsion |
| 890 |
> |
potential for all 1, 4 bonded pairs. The pairwise portions of the |
| 891 |
> |
non-bonded interactions are excluded for atom pairs that are involved |
| 892 |
> |
in the smae bond, bend, or torsion. All other atom pairs within a |
| 893 |
> |
molecule are subject to the LJ pair potential. |
| 894 |
|
|
| 372 |
– |
|
| 895 |
|
The bend potential of a molecule is represented by the following function: |
| 896 |
|
\begin{equation} |
| 897 |
< |
V_{\text{bend}}(\theta_{ijk}) = k_{\theta}( \theta_{ijk} - \theta_0 )^2, \label{eq:bendPot} |
| 897 |
> |
V_{\text{bend}}(\theta_{ijk}) = k_{\theta}( \theta_{ijk} - \theta_0 |
| 898 |
> |
)^2, \label{eq:bendPot} |
| 899 |
|
\end{equation} |
| 900 |
|
where $\theta_{ijk}$ is the angle defined by atoms $i$, $j$, and $k$ |
| 901 |
|
(see Fig.~\ref{oopseFig:lipidModel}), $\theta_0$ is the equilibrium |
| 935 |
|
k_3 &= 4c_3. |
| 936 |
|
\end{align*} |
| 937 |
|
By recasting the potential as a power series, repeated trigonometric |
| 938 |
< |
evaluations are avoided during the calculation of the potential energy. |
| 938 |
> |
evaluations are avoided during the calculation of the potential |
| 939 |
> |
energy. |
| 940 |
|
|
| 941 |
|
|
| 942 |
< |
The cross potential between molecules $I$ and $J$, $V^{IJ}_{\text{Cross}}$, is |
| 943 |
< |
as follows: |
| 942 |
> |
The cross potential between molecules $I$ and $J$, |
| 943 |
> |
$V^{IJ}_{\text{Cross}}$, is as follows: |
| 944 |
|
\begin{equation} |
| 945 |
|
V^{IJ}_{\text{Cross}} = |
| 946 |
|
\sum_{i \in I} \sum_{j \in J} |
| 970 |
|
Here $\mathbf{r}_{ij}$ is the vector starting at atom $i$ pointing |
| 971 |
|
towards $j$, and $\boldsymbol{\Omega}_i$ and $\boldsymbol{\Omega}_j$ |
| 972 |
|
are the orientational degrees of freedom for atoms $i$ and $j$ |
| 973 |
< |
respectively. $|\mu_i|$ is the magnitude of the dipole moment of atom |
| 974 |
< |
$i$, $\boldsymbol{\hat{u}}_i$ is the standard unit orientation vector |
| 975 |
< |
of $\boldsymbol{\Omega}_i$, and $\boldsymbol{\hat{r}}_{ij}$ is the |
| 976 |
< |
unit vector pointing along $\mathbf{r}_{ij}$ |
| 973 |
> |
respectively. The magnitude of the dipole moment of atom $i$ is |
| 974 |
> |
$|\mu_i|$, $\boldsymbol{\hat{u}}_i$ is the standard unit orientation |
| 975 |
> |
vector of $\boldsymbol{\Omega}_i$, and $\boldsymbol{\hat{r}}_{ij}$ is |
| 976 |
> |
the unit vector pointing along $\mathbf{r}_{ij}$ |
| 977 |
|
($\boldsymbol{\hat{r}}_{ij}=\mathbf{r}_{ij}/|\mathbf{r}_{ij}|$). |
| 978 |
|
|
| 979 |
< |
To improve computational efficiency of the dipole-dipole interactions, |
| 980 |
< |
{\sc oopse} employs an electrostatic cutoff radius. This parameter can |
| 457 |
< |
be set in the \texttt{.bass} file, and controls the length scale over |
| 458 |
< |
which dipole interactions are felt. To compensate for the |
| 459 |
< |
discontinuity in the potential and the forces at the cutoff radius, we |
| 460 |
< |
have implemented a switching function to smoothly scale the |
| 461 |
< |
dipole-dipole interaction at the cutoff. |
| 462 |
< |
\begin{equation} |
| 463 |
< |
S(r_{ij}) = |
| 464 |
< |
\begin{cases} |
| 465 |
< |
1 & \text{if $r_{ij} \le r_t$},\\ |
| 466 |
< |
\frac{(r_{\text{cut}} + 2r_{ij} - 3r_t)(r_{\text{cut}} - r_{ij})^2} |
| 467 |
< |
{(r_{\text{cut}} - r_t)^2} |
| 468 |
< |
& \text{if $r_t < r_{ij} \le r_{\text{cut}}$}, \\ |
| 469 |
< |
0 & \text{if $r_{ij} > r_{\text{cut}}$.} |
| 470 |
< |
\end{cases} |
| 471 |
< |
\label{eq:dipoleSwitching} |
| 472 |
< |
\end{equation} |
| 473 |
< |
Here $S(r_{ij})$ scales the potential at a given $r_{ij}$, and $r_t$ |
| 474 |
< |
is the taper radius some given thickness less than the electrostatic |
| 475 |
< |
cutoff. The switching thickness can be set in the \texttt{.bass} file. |
| 979 |
> |
\subsubsection{\label{oopseSec:SSD}The {\sc duff} Water Models: SSD/E |
| 980 |
> |
and SSD/RF} |
| 981 |
|
|
| 477 |
– |
\subsection{\label{oopseSec:SSD}The {\sc duff} Water Models: SSD/E and SSD/RF} |
| 478 |
– |
|
| 982 |
|
In the interest of computational efficiency, the default solvent used |
| 983 |
|
by {\sc oopse} is the extended Soft Sticky Dipole (SSD/E) water |
| 984 |
|
model.\cite{fennell04} The original SSD was developed by Ichiye |
| 1037 |
|
can be found in the original SSD |
| 1038 |
|
articles.\cite{liu96:new_model,liu96:monte_carlo,chandra99:ssd_md,Ichiye03} |
| 1039 |
|
|
| 1040 |
+ |
\begin{figure} |
| 1041 |
+ |
\centering |
| 1042 |
+ |
\includegraphics[width=\linewidth]{waterAngle.eps} |
| 1043 |
+ |
\caption[Coordinate definition for the SSD/E water model]{Coordinates |
| 1044 |
+ |
for the interaction between two SSD/E water molecules. $\theta_{ij}$ |
| 1045 |
+ |
is the angle that $r_{ij}$ makes with the $\hat{z}$ vector in the |
| 1046 |
+ |
body-fixed frame for molecule $i$. The $\hat{z}$ vector bisects the |
| 1047 |
+ |
HOH angle in each water molecule. } |
| 1048 |
+ |
\label{oopseFig:ssd} |
| 1049 |
+ |
\end{figure} |
| 1050 |
+ |
|
| 1051 |
+ |
|
| 1052 |
|
Since SSD/E is a single-point {\it dipolar} model, the force |
| 1053 |
|
calculations are simplified significantly relative to the standard |
| 1054 |
|
{\it charged} multi-point models. In the original Monte Carlo |
| 1055 |
|
simulations using this model, Ichiye {\it et al.} reported that using |
| 1056 |
|
SSD decreased computer time by a factor of 6-7 compared to other |
| 1057 |
< |
models.\cite{liu96:new_model} What is most impressive is that these savings |
| 1058 |
< |
did not come at the expense of accurate depiction of the liquid state |
| 1059 |
< |
properties. Indeed, SSD/E maintains reasonable agreement with the Head-Gordon |
| 1060 |
< |
diffraction data for the structural features of liquid |
| 1061 |
< |
water.\cite{hura00,liu96:new_model} Additionally, the dynamical properties |
| 1062 |
< |
exhibited by SSD/E agree with experiment better than those of more |
| 1063 |
< |
computationally expensive models (like TIP3P and |
| 1064 |
< |
SPC/E).\cite{chandra99:ssd_md} The combination of speed and accurate depiction |
| 1065 |
< |
of solvent properties makes SSD/E a very attractive model for the |
| 1066 |
< |
simulation of large scale biochemical simulations. |
| 1057 |
> |
models.\cite{liu96:new_model} What is most impressive is that these |
| 1058 |
> |
savings did not come at the expense of accurate depiction of the |
| 1059 |
> |
liquid state properties. Indeed, SSD/E maintains reasonable agreement |
| 1060 |
> |
with the Head-Gordon diffraction data for the structural features of |
| 1061 |
> |
liquid water.\cite{hura00,liu96:new_model} Additionally, the dynamical |
| 1062 |
> |
properties exhibited by SSD/E agree with experiment better than those |
| 1063 |
> |
of more computationally expensive models (like TIP3P and |
| 1064 |
> |
SPC/E).\cite{chandra99:ssd_md} The combination of speed and accurate |
| 1065 |
> |
depiction of solvent properties makes SSD/E a very attractive model |
| 1066 |
> |
for the simulation of large scale biochemical simulations. |
| 1067 |
|
|
| 1068 |
|
Recent constant pressure simulations revealed issues in the original |
| 1069 |
|
SSD model that led to lower than expected densities at all target |
| 1072 |
|
exhibits improved liquid structure and transport behavior. If the use |
| 1073 |
|
of a reaction field long-range interaction correction is desired, it |
| 1074 |
|
is recommended that the parameters be modified to those of the SSD/RF |
| 1075 |
< |
model (an SSD variant parameterized for reaction field). Solvent parameters can be easily modified in an accompanying |
| 1076 |
< |
\texttt{.bass} file as illustrated in the scheme below. A table of the |
| 1077 |
< |
parameter values and the drawbacks and benefits of the different |
| 1078 |
< |
density corrected SSD models can be found in |
| 1079 |
< |
reference~\cite{fennell04}. |
| 1075 |
> |
model (an SSD variant parameterized for reaction field). These solvent |
| 1076 |
> |
parameters are listed and can be easily modified in the {\sc duff} |
| 1077 |
> |
force field file ({\tt DUFF.frc}). A table of the parameter values |
| 1078 |
> |
and the drawbacks and benefits of the different density corrected SSD |
| 1079 |
> |
models can be found in reference~\citen{fennell04}. |
| 1080 |
|
|
| 1081 |
< |
\begin{lstlisting}[float,caption={[A simulation of {\sc ssd} water]A portion of a \texttt{.bass} file showing a simulation including {\sc ssd} water.},label={sch:ssd}] |
| 1081 |
> |
\subsection{\label{oopseSec:WATER}The {\sc water} Force Field} |
| 1082 |
|
|
| 1083 |
< |
#include "water.mdl" |
| 1083 |
> |
In addition to the {\sc duff} force field's solvent description, a |
| 1084 |
> |
separate {\sc water} force field has been included for simulating most |
| 1085 |
> |
of the common rigid-body water models. This force field includes the |
| 1086 |
> |
simple and point-dipolar models (SSD, SSD1, SSD/E, SSD/RF, and DPD |
| 1087 |
> |
water), as well as the common charge-based models (SPC, SPC/E, TIP3P, |
| 1088 |
> |
TIP4P, and |
| 1089 |
> |
TIP5P).\cite{liu96:new_model,Ichiye03,fennell04,Marrink01,Berendsen81,Berendsen87,Jorgensen83,Mahoney00} |
| 1090 |
> |
In order to handle these models, charge-charge interactions were |
| 1091 |
> |
included in the force-loop: |
| 1092 |
> |
\begin{equation} |
| 1093 |
> |
V_{\text{charge}}(r_{ij}) = \sum_{ij}\frac{q_iq_je^2}{r_{ij}}, |
| 1094 |
> |
\end{equation} |
| 1095 |
> |
where $q$ represents the charge on particle $i$ or $j$, and $e$ is the |
| 1096 |
> |
charge of an electron in Coulombs. As with the other pair |
| 1097 |
> |
interactions, charges can be simulated with a pure cutoff or a |
| 1098 |
> |
reaction field. The {\sc water} force field can be easily expanded |
| 1099 |
> |
through modification of the {\sc water} force field file ({\tt |
| 1100 |
> |
WATER.frc}). By adding atom types and inserting the appropriate |
| 1101 |
> |
parameters, it is possible to extend the force field to handle rigid |
| 1102 |
> |
molecules other than water. |
| 1103 |
|
|
| 570 |
– |
nComponents = 1; |
| 571 |
– |
component{ |
| 572 |
– |
type = "SSD_water"; |
| 573 |
– |
nMol = 864; |
| 574 |
– |
} |
| 575 |
– |
|
| 576 |
– |
initialConfig = "liquidWater.init"; |
| 577 |
– |
|
| 578 |
– |
forceField = "DUFF"; |
| 579 |
– |
|
| 580 |
– |
/* |
| 581 |
– |
* The following two flags set the cutoff |
| 582 |
– |
* radius for the electrostatic forces |
| 583 |
– |
* as well as the skin thickness of the switching |
| 584 |
– |
* function. |
| 585 |
– |
*/ |
| 586 |
– |
|
| 587 |
– |
electrostaticCutoffRadius = 9.2; |
| 588 |
– |
electrostaticSkinThickness = 1.38; |
| 589 |
– |
|
| 590 |
– |
\end{lstlisting} |
| 591 |
– |
|
| 592 |
– |
|
| 1104 |
|
\subsection{\label{oopseSec:eam}Embedded Atom Method} |
| 1105 |
|
|
| 1106 |
< |
There are Molecular Dynamics packages which have the |
| 1107 |
< |
capacity to simulate metallic systems, including some that have |
| 1108 |
< |
parallel computational abilities\cite{plimpton93}. Potentials that |
| 1109 |
< |
describe bonding transition metal |
| 1110 |
< |
systems\cite{Finnis84,Ercolessi88,Chen90,Qi99,Ercolessi02} have an |
| 600 |
< |
attractive interaction which models ``Embedding'' |
| 601 |
< |
a positively charged metal ion in the electron density due to the |
| 1106 |
> |
{\sc oopse} implements a potential that describes bonding in |
| 1107 |
> |
transition metal |
| 1108 |
> |
systems.~\cite{Finnis84,Ercolessi88,Chen90,Qi99,Ercolessi02} This |
| 1109 |
> |
potential has an attractive interaction which models ``Embedding'' a |
| 1110 |
> |
positively charged pseudo-atom core in the electron density due to the |
| 1111 |
|
free valance ``sea'' of electrons created by the surrounding atoms in |
| 1112 |
< |
the system. A mostly-repulsive pairwise part of the potential |
| 1113 |
< |
describes the interaction of the positively charged metal core ions |
| 1114 |
< |
with one another. A particular potential description called the |
| 1115 |
< |
Embedded Atom Method\cite{Daw84,FBD86,johnson89,Lu97}({\sc eam}) that has |
| 1116 |
< |
particularly wide adoption has been selected for inclusion in {\sc oopse}. A |
| 1117 |
< |
good review of {\sc eam} and other metallic potential formulations was written |
| 1118 |
< |
by Voter.\cite{voter} |
| 1112 |
> |
the system. A pairwise part of the potential (which is primarily |
| 1113 |
> |
repulsive) describes the interaction of the positively charged metal |
| 1114 |
> |
core ions with one another. The Embedded Atom Method ({\sc |
| 1115 |
> |
eam})~\cite{Daw84,FBD86,johnson89,Lu97} has been widely adopted in the |
| 1116 |
> |
materials science community and has been included in {\sc oopse}. A |
| 1117 |
> |
good review of {\sc eam} and other formulations of metallic potentials |
| 1118 |
> |
was given by Voter.\cite{Voter:95} |
| 1119 |
|
|
| 1120 |
|
The {\sc eam} potential has the form: |
| 1121 |
< |
\begin{eqnarray} |
| 1122 |
< |
V & = & \sum_{i} F_{i}\left[\rho_{i}\right] + \sum_{i} \sum_{j \neq i} |
| 1123 |
< |
\phi_{ij}({\bf r}_{ij}), \\ |
| 1124 |
< |
\rho_{i} & = & \sum_{j \neq i} f_{j}({\bf r}_{ij}), |
| 1125 |
< |
\end{eqnarray} |
| 617 |
< |
where $F_{i} $ is the embedding function that equates the energy |
| 1121 |
> |
\begin{equation} |
| 1122 |
> |
V = \sum_{i} F_{i}\left[\rho_{i}\right] + \sum_{i} \sum_{j \neq i} |
| 1123 |
> |
\phi_{ij}({\bf r}_{ij}) |
| 1124 |
> |
\end{equation} |
| 1125 |
> |
where $F_{i} $ is an embedding functional that approximates the energy |
| 1126 |
|
required to embed a positively-charged core ion $i$ into a linear |
| 1127 |
|
superposition of spherically averaged atomic electron densities given |
| 1128 |
< |
by $\rho_{i}$. $\phi_{ij}$ is a primarily repulsive pairwise |
| 1129 |
< |
interaction between atoms $i$ and $j$. In the original formulation of |
| 1130 |
< |
{\sc eam}\cite{Daw84}, $\phi_{ij}$ was an entirely repulsive term, |
| 1131 |
< |
however in later refinements to {\sc eam} have shown that non-uniqueness |
| 1132 |
< |
between $F$ and $\phi$ allow for more general forms for |
| 1133 |
< |
$\phi$.\cite{Daw89} There is a cutoff distance, $r_{cut}$, which |
| 1134 |
< |
limits the summations in the {\sc eam} equation to the few dozen atoms |
| 1128 |
> |
by $\rho_{i}$, |
| 1129 |
> |
\begin{equation} |
| 1130 |
> |
\rho_{i} = \sum_{j \neq i} f_{j}({\bf r}_{ij}), |
| 1131 |
> |
\end{equation} |
| 1132 |
> |
Since the density at site $i$ ($\rho_i$) must be computed before the |
| 1133 |
> |
embedding functional can be evaluated, {\sc eam} and the related |
| 1134 |
> |
transition metal potentials require two loops through the atom pairs |
| 1135 |
> |
to compute the inter-atomic forces. |
| 1136 |
> |
|
| 1137 |
> |
The pairwise portion of the potential, $\phi_{ij}$, is a primarily |
| 1138 |
> |
repulsive interaction between atoms $i$ and $j$. In the original |
| 1139 |
> |
formulation of {\sc eam}\cite{Daw84}, $\phi_{ij}$ was an entirely |
| 1140 |
> |
repulsive term; however later refinements to {\sc eam} allowed for |
| 1141 |
> |
more general forms for $\phi$.\cite{Daw89} The effective cutoff |
| 1142 |
> |
distance, $r_{{\text cut}}$ is the distance at which the values of |
| 1143 |
> |
$f(r)$ and $\phi(r)$ drop to zero for all atoms present in the |
| 1144 |
> |
simulation. In practice, this distance is fairly small, limiting the |
| 1145 |
> |
summations in the {\sc eam} equation to the few dozen atoms |
| 1146 |
|
surrounding atom $i$ for both the density $\rho$ and pairwise $\phi$ |
| 1147 |
< |
interactions. Foiles \emph{et al}.~fit {\sc eam} potentials for the fcc |
| 629 |
< |
metals Cu, Ag, Au, Ni, Pd, Pt and alloys of these metals.\cite{FBD86} |
| 630 |
< |
These fits are included in {\sc oopse}. |
| 1147 |
> |
interactions. |
| 1148 |
|
|
| 1149 |
+ |
In computing forces for alloys, mixing rules as outlined by |
| 1150 |
+ |
Johnson~\cite{johnson89} are used to compute the heterogenous pair |
| 1151 |
+ |
potential, |
| 1152 |
+ |
\begin{eqnarray} |
| 1153 |
+ |
\label{eq:johnson} |
| 1154 |
+ |
\phi_{ab}(r)=\frac{1}{2}\left( |
| 1155 |
+ |
\frac{f_{b}(r)}{f_{a}(r)}\phi_{aa}(r)+ |
| 1156 |
+ |
\frac{f_{a}(r)}{f_{b}(r)}\phi_{bb}(r) |
| 1157 |
+ |
\right). |
| 1158 |
+ |
\end{eqnarray} |
| 1159 |
+ |
No mixing rule is needed for the densities, since the density at site |
| 1160 |
+ |
$i$ is simply the linear sum of density contributions of all the other |
| 1161 |
+ |
atoms. |
| 1162 |
+ |
|
| 1163 |
+ |
The {\sc eam} force field illustrates an additional feature of {\sc |
| 1164 |
+ |
oopse}. Foiles, Baskes and Daw fit {\sc eam} potentials for Cu, Ag, |
| 1165 |
+ |
Au, Ni, Pd, Pt and alloys of these metals.\cite{FBD86} These fits are |
| 1166 |
+ |
included in {\sc oopse} as the {\tt u3} variant of the {\sc eam} force |
| 1167 |
+ |
field. Voter and Chen reparamaterized a set of {\sc eam} functions |
| 1168 |
+ |
which do a better job of predicting melting points.\cite{Voter:87} |
| 1169 |
+ |
These functions are included in {\sc oopse} as the {\tt VC} variant of |
| 1170 |
+ |
the {\sc eam} force field. An additional set of functions (the |
| 1171 |
+ |
``Universal 6'' functions) are included in {\sc oopse} as the {\tt u6} |
| 1172 |
+ |
variant of {\sc eam}. For example, to specify the Voter-Chen variant |
| 1173 |
+ |
of the {\sc eam} force field, the user would add the {\tt |
| 1174 |
+ |
forceFieldVariant = "VC";} line to the meta-data file. |
| 1175 |
+ |
|
| 1176 |
+ |
The potential files used by the {\sc eam} force field are in the |
| 1177 |
+ |
standard {\tt funcfl} format, which is the format utilized by a number |
| 1178 |
+ |
of other codes (e.g. ParaDyn~\cite{Paradyn}, {\sc dynamo 86}). It |
| 1179 |
+ |
should be noted that the energy units in these files are in eV, not |
| 1180 |
+ |
$\mbox{kcal mol}^{-1}$ as in the rest of the {\sc oopse} force field |
| 1181 |
+ |
files. |
| 1182 |
+ |
|
| 1183 |
|
\subsection{\label{oopseSec:pbc}Periodic Boundary Conditions} |
| 1184 |
|
|
| 1185 |
|
\newcommand{\roundme}{\operatorname{round}} |
| 1186 |
|
|
| 1187 |
< |
\textit{Periodic boundary conditions} are widely used to simulate bulk properties with a relatively small number of particles. The |
| 1188 |
< |
simulation box is replicated throughout space to form an infinite |
| 1187 |
> |
\textit{Periodic boundary conditions} are widely used to simulate bulk |
| 1188 |
> |
properties with a relatively small number of particles. In this method |
| 1189 |
> |
the simulation box is replicated throughout space to form an infinite |
| 1190 |
|
lattice. During the simulation, when a particle moves in the primary |
| 1191 |
|
cell, its image in other cells move in exactly the same direction with |
| 1192 |
|
exactly the same orientation. Thus, as a particle leaves the primary |
| 1193 |
|
cell, one of its images will enter through the opposite face. If the |
| 1194 |
|
simulation box is large enough to avoid ``feeling'' the symmetries of |
| 1195 |
|
the periodic lattice, surface effects can be ignored. The available |
| 1196 |
< |
periodic cells in OOPSE are cubic, orthorhombic and parallelepiped. We |
| 1197 |
< |
use a $3 \times 3$ matrix, $\mathsf{H}$, to describe the shape and |
| 1198 |
< |
size of the simulation box. $\mathsf{H}$ is defined: |
| 1196 |
> |
periodic cells in {\sc oopse} are cubic, orthorhombic and |
| 1197 |
> |
parallelepiped. {\sc oopse} use a $3 \times 3$ matrix, $\mathsf{H}$, |
| 1198 |
> |
to describe the shape and size of the simulation box. $\mathsf{H}$ is |
| 1199 |
> |
defined: |
| 1200 |
|
\begin{equation} |
| 1201 |
|
\mathsf{H} = ( \mathbf{h}_x, \mathbf{h}_y, \mathbf{h}_z ), |
| 1202 |
|
\end{equation} |
| 1213 |
|
\end{align} |
| 1214 |
|
The vector $\mathbf{s}$ is now a vector expressed as the number of box |
| 1215 |
|
lengths in the $\mathbf{h}_x$, $\mathbf{h}_y$, and $\mathbf{h}_z$ |
| 1216 |
< |
directions. To find the minimum image of a vector $\mathbf{r}$, we |
| 1217 |
< |
first convert it to its corresponding vector in box space, and then, |
| 1218 |
< |
cast each element to lie in the range $[-0.5,0.5]$: |
| 1216 |
> |
directions. To find the minimum image of a vector $\mathbf{r}$, {\sc |
| 1217 |
> |
oopse} first converts it to its corresponding vector in box space, and |
| 1218 |
> |
then casts each element to lie in the range $[-0.5,0.5]$: |
| 1219 |
|
\begin{equation} |
| 1220 |
|
s_{i}^{\prime}=s_{i}-\roundme(s_{i}), |
| 1221 |
|
\end{equation} |
| 1231 |
|
Here $\lfloor x \rfloor$ is the floor operator, and gives the largest |
| 1232 |
|
integer value that is not greater than $x$, and $\lceil x \rceil$ is |
| 1233 |
|
the ceiling operator, and gives the smallest integer that is not less |
| 1234 |
< |
than $x$. For example, $\roundme(3.6)=4$, $\roundme(3.1)=3$, |
| 682 |
< |
$\roundme(-3.6)=-4$, $\roundme(-3.1)=-3$. |
| 1234 |
> |
than $x$. |
| 1235 |
|
|
| 1236 |
< |
Finally, we obtain the minimum image coordinates $\mathbf{r}^{\prime}$ by |
| 1237 |
< |
transforming back to real space, |
| 1236 |
> |
Finally, the minimum image coordinates $\mathbf{r}^{\prime}$ are |
| 1237 |
> |
obtained by transforming back to real space, |
| 1238 |
|
\begin{equation} |
| 1239 |
|
\mathbf{r}^{\prime}=\mathsf{H}^{-1}\mathbf{s}^{\prime}.% |
| 1240 |
|
\end{equation} |
| 1241 |
|
In this way, particles are allowed to diffuse freely in $\mathbf{r}$, |
| 1242 |
< |
but their minimum images, $\mathbf{r}^{\prime}$ are used to compute |
| 1242 |
> |
but their minimum images, or $\mathbf{r}^{\prime}$, are used to compute |
| 1243 |
|
the inter-atomic forces. |
| 1244 |
|
|
| 1245 |
|
|
| 694 |
– |
\section{\label{oopseSec:IOfiles}Input and Output Files} |
| 1246 |
|
|
| 696 |
– |
\subsection{{\sc bass} and Model Files} |
| 697 |
– |
|
| 698 |
– |
Every {\sc oopse} simulation begins with a Bizarre Atom Simulation |
| 699 |
– |
Syntax ({\sc bass}) file. {\sc bass} is a script syntax that is parsed |
| 700 |
– |
by {\sc oopse} at runtime. The {\sc bass} file allows for the user to |
| 701 |
– |
completely describe the system they wish to simulate, as well as tailor |
| 702 |
– |
{\sc oopse}'s behavior during the simulation. {\sc bass} files are |
| 703 |
– |
denoted with the extension |
| 704 |
– |
\texttt{.bass}, an example file is shown in |
| 705 |
– |
Scheme~\ref{sch:bassExample}. |
| 706 |
– |
|
| 707 |
– |
\begin{lstlisting}[float,caption={[An example of a complete {\sc bass} file] An example showing a complete {\sc bass} file.},label={sch:bassExample}] |
| 708 |
– |
|
| 709 |
– |
molecule{ |
| 710 |
– |
name = "Ar"; |
| 711 |
– |
nAtoms = 1; |
| 712 |
– |
atom[0]{ |
| 713 |
– |
type="Ar"; |
| 714 |
– |
position( 0.0, 0.0, 0.0 ); |
| 715 |
– |
} |
| 716 |
– |
} |
| 717 |
– |
|
| 718 |
– |
nComponents = 1; |
| 719 |
– |
component{ |
| 720 |
– |
type = "Ar"; |
| 721 |
– |
nMol = 108; |
| 722 |
– |
} |
| 723 |
– |
|
| 724 |
– |
initialConfig = "./argon.init"; |
| 725 |
– |
|
| 726 |
– |
forceField = "LJ"; |
| 727 |
– |
ensemble = "NVE"; // specify the simulation ensemble |
| 728 |
– |
dt = 1.0; // the time step for integration |
| 729 |
– |
runTime = 1e3; // the total simulation run time |
| 730 |
– |
sampleTime = 100; // trajectory file frequency |
| 731 |
– |
statusTime = 50; // statistics file frequency |
| 732 |
– |
|
| 733 |
– |
\end{lstlisting} |
| 734 |
– |
|
| 735 |
– |
Within the \texttt{.bass} file it is necessary to provide a complete |
| 736 |
– |
description of the molecule before it is actually placed in the |
| 737 |
– |
simulation. The {\sc bass} syntax was originally developed with this |
| 738 |
– |
goal in mind, and allows for the specification of all the atoms in a |
| 739 |
– |
molecular prototype, as well as any bonds, bends, or torsions. These |
| 740 |
– |
descriptions can become lengthy for complex molecules, and it would be |
| 741 |
– |
inconvenient to duplicate the simulation at the beginning of each {\sc |
| 742 |
– |
bass} script. Addressing this issue {\sc bass} allows for the |
| 743 |
– |
inclusion of model files at the top of a \texttt{.bass} file. These |
| 744 |
– |
model files, denoted with the \texttt{.mdl} extension, allow the user |
| 745 |
– |
to describe a molecular prototype once, then simply include it into |
| 746 |
– |
each simulation containing that molecule. Returning to the example in |
| 747 |
– |
Scheme~\ref{sch:bassExample}, the \texttt{.mdl} file's contents would |
| 748 |
– |
be Scheme~\ref{sch:mdlExample}, and the new \texttt{.bass} file would |
| 749 |
– |
become Scheme~\ref{sch:bassExPrime}. |
| 750 |
– |
|
| 751 |
– |
\begin{lstlisting}[float,caption={An example \texttt{.mdl} file.},label={sch:mdlExample}] |
| 752 |
– |
|
| 753 |
– |
molecule{ |
| 754 |
– |
name = "Ar"; |
| 755 |
– |
nAtoms = 1; |
| 756 |
– |
atom[0]{ |
| 757 |
– |
type="Ar"; |
| 758 |
– |
position( 0.0, 0.0, 0.0 ); |
| 759 |
– |
} |
| 760 |
– |
} |
| 761 |
– |
|
| 762 |
– |
\end{lstlisting} |
| 763 |
– |
|
| 764 |
– |
\begin{lstlisting}[float,caption={Revised {\sc bass} example.},label={sch:bassExPrime}] |
| 765 |
– |
|
| 766 |
– |
#include "argon.mdl" |
| 767 |
– |
|
| 768 |
– |
nComponents = 1; |
| 769 |
– |
component{ |
| 770 |
– |
type = "Ar"; |
| 771 |
– |
nMol = 108; |
| 772 |
– |
} |
| 773 |
– |
|
| 774 |
– |
initialConfig = "./argon.init"; |
| 775 |
– |
|
| 776 |
– |
forceField = "LJ"; |
| 777 |
– |
ensemble = "NVE"; |
| 778 |
– |
dt = 1.0; |
| 779 |
– |
runTime = 1e3; |
| 780 |
– |
sampleTime = 100; |
| 781 |
– |
statusTime = 50; |
| 782 |
– |
|
| 783 |
– |
\end{lstlisting} |
| 784 |
– |
|
| 785 |
– |
\subsection{\label{oopseSec:coordFiles}Coordinate Files} |
| 786 |
– |
|
| 787 |
– |
The standard format for storage of a systems coordinates is a modified |
| 788 |
– |
xyz-file syntax, the exact details of which can be seen in |
| 789 |
– |
Scheme~\ref{sch:dumpFormat}. As all bonding and molecular information |
| 790 |
– |
is stored in the \texttt{.bass} and \texttt{.mdl} files, the |
| 791 |
– |
coordinate files are simply the complete set of coordinates for each |
| 792 |
– |
atom at a given simulation time. One important note, although the |
| 793 |
– |
simulation propagates the complete rotation matrix, directional |
| 794 |
– |
entities are written out using quanternions, to save space in the |
| 795 |
– |
output files. |
| 796 |
– |
|
| 797 |
– |
\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 $\mathsf{H}$ column vectors. It is important to note, that for extended system ensembles, additional information pertinent to the integrators may be stored on this line as well. 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 velocities.},label=sch:dumpFormat] |
| 798 |
– |
|
| 799 |
– |
nAtoms |
| 800 |
– |
time; Hxx Hyx Hzx; Hxy Hyy Hzy; Hxz Hyz Hzz; |
| 801 |
– |
Name1 x y z vx vy vz q0 q1 q2 q3 jx jy jz |
| 802 |
– |
Name2 x y z vx vy vz q0 q1 q2 q3 jx jy jz |
| 803 |
– |
etc... |
| 804 |
– |
|
| 805 |
– |
\end{lstlisting} |
| 806 |
– |
|
| 807 |
– |
|
| 808 |
– |
There are three major files used by {\sc oopse} written in the |
| 809 |
– |
coordinate format, they are as follows: the initialization file |
| 810 |
– |
(\texttt{.init}), the simulation trajectory file (\texttt{.dump}), and |
| 811 |
– |
the final coordinates of the simulation. The initialization file is |
| 812 |
– |
necessary for {\sc oopse} to start the simulation with the proper |
| 813 |
– |
coordinates, and is generated before the simulation run. The |
| 814 |
– |
trajectory file is created at the beginning of the simulation, and is |
| 815 |
– |
used to store snapshots of the simulation at regular intervals. The |
| 816 |
– |
first frame is a duplication of the |
| 817 |
– |
\texttt{.init} file, and each subsequent frame is appended to the file |
| 818 |
– |
at an interval specified in the \texttt{.bass} file with the |
| 819 |
– |
\texttt{sampleTime} flag. The final coordinate file is the end of run file. The |
| 820 |
– |
\texttt{.eor} file stores the final configuration of the system for a |
| 821 |
– |
given simulation. The file is updated at the same time as the |
| 822 |
– |
\texttt{.dump} file, however, it only contains the most recent |
| 823 |
– |
frame. In this way, an \texttt{.eor} file may be used as the |
| 824 |
– |
initialization file to a second simulation in order to continue a |
| 825 |
– |
simulation or recover one from a processor that has crashed during the |
| 826 |
– |
course of the run. |
| 827 |
– |
|
| 828 |
– |
\subsection{\label{oopseSec:initCoords}Generation of Initial Coordinates} |
| 829 |
– |
|
| 830 |
– |
As was stated in Sec.~\ref{oopseSec:coordFiles}, an initialization |
| 831 |
– |
file is needed to provide the starting coordinates for a |
| 832 |
– |
simulation. The {\sc oopse} package provides several system building |
| 833 |
– |
programs to aid in the creation of the \texttt{.init} |
| 834 |
– |
file. The programs use {\sc bass}, and will recognize |
| 835 |
– |
arguments and parameters in the \texttt{.bass} file that would |
| 836 |
– |
otherwise be ignored by the simulation. |
| 837 |
– |
|
| 838 |
– |
\subsection{The Statistics File} |
| 839 |
– |
|
| 840 |
– |
The last output file generated by {\sc oopse} is the statistics |
| 841 |
– |
file. This file records such statistical quantities as the |
| 842 |
– |
instantaneous temperature, volume, pressure, etc. It is written out |
| 843 |
– |
with the frequency specified in the \texttt{.bass} file with the |
| 844 |
– |
\texttt{statusTime} keyword. The file allows the user to observe the |
| 845 |
– |
system variables as a function of simulation time while the simulation |
| 846 |
– |
is in progress. One useful function the statistics file serves is to |
| 847 |
– |
monitor the conserved quantity of a given simulation ensemble, this |
| 848 |
– |
allows the user to observe the stability of the integrator. The |
| 849 |
– |
statistics file is denoted with the \texttt{.stat} file extension. |
| 850 |
– |
|
| 1247 |
|
\section{\label{oopseSec:mechanics}Mechanics} |
| 1248 |
|
|
| 1249 |
|
\subsection{\label{oopseSec:integrate}Integrating the Equations of Motion: the |
| 1250 |
< |
DLM method} |
| 1250 |
> |
{\sc dlm} method} |
| 1251 |
|
|
| 1252 |
|
The default method for integrating the equations of motion in {\sc |
| 1253 |
|
oopse} is a velocity-Verlet version of the symplectic splitting method |
| 1254 |
|
proposed by Dullweber, Leimkuhler and McLachlan |
| 1255 |
< |
(DLM).\cite{Dullweber1997} When there are no directional atoms or |
| 1255 |
> |
({\sc dlm}).\cite{Dullweber1997} When there are no directional atoms or |
| 1256 |
|
rigid bodies present in the simulation, this integrator becomes the |
| 1257 |
|
standard velocity-Verlet integrator which is known to sample the |
| 1258 |
|
microcanonical (NVE) ensemble.\cite{Frenkel1996} |
| 1259 |
|
|
| 1260 |
|
Previous integration methods for orientational motion have problems |
| 1261 |
< |
that are avoided in the DLM method. Direct propagation of the Euler |
| 1261 |
> |
that are avoided in the {\sc dlm} method. Direct propagation of the Euler |
| 1262 |
|
angles has a known $1/\sin\theta$ divergence in the equations of |
| 1263 |
< |
motion for $\phi$ and $\psi$,\cite{allen87:csl} leading to |
| 1264 |
< |
numerical instabilities any time one of the directional atoms or rigid |
| 1265 |
< |
bodies has an orientation near $\theta=0$ or $\theta=\pi$. More |
| 1266 |
< |
modern quaternion-based integration methods have relatively poor |
| 1267 |
< |
energy conservation. While quaternions work well for orientational |
| 1268 |
< |
motion in other ensembles, the microcanonical ensemble has a |
| 1269 |
< |
constant energy requirement that is quite sensitive to errors in the |
| 1270 |
< |
equations of motion. An earlier implementation of {\sc oopse} |
| 1271 |
< |
utilized quaternions for propagation of rotational motion; however, a |
| 1272 |
< |
detailed investigation showed that they resulted in a steady drift in |
| 877 |
< |
the total energy, something that has been observed by |
| 878 |
< |
Laird {\it et al.}\cite{Laird97} |
| 1263 |
> |
motion for $\phi$ and $\psi$,\cite{allen87:csl} leading to numerical |
| 1264 |
> |
instabilities any time one of the directional atoms or rigid bodies |
| 1265 |
> |
has an orientation near $\theta=0$ or $\theta=\pi$. Quaternion-based |
| 1266 |
> |
integration methods work well for propagating orientational motion; |
| 1267 |
> |
however, energy conservation concerns arise when using the |
| 1268 |
> |
microcanonical (NVE) ensemble. An earlier implementation of {\sc |
| 1269 |
> |
oopse} utilized quaternions for propagation of rotational motion; |
| 1270 |
> |
however, a detailed investigation showed that they resulted in a |
| 1271 |
> |
steady drift in the total energy, something that has been observed by |
| 1272 |
> |
Laird {\it et al.}\cite{Laird97} |
| 1273 |
|
|
| 1274 |
|
The key difference in the integration method proposed by Dullweber |
| 1275 |
|
\emph{et al.} is that the entire $3 \times 3$ rotation matrix is |
| 1343 |
|
represented by ${\bf j}$. This equation of motion for angular momenta |
| 1344 |
|
is equivalent to the more familiar body-fixed forms, |
| 1345 |
|
\begin{eqnarray} |
| 1346 |
< |
\dot{j_{x}} & = & \tau^b_x(t) + |
| 1347 |
< |
\left(\overleftrightarrow{\mathsf{I}}_{yy} - \overleftrightarrow{\mathsf{I}}_{zz} \right) j_y j_z, \\ |
| 1348 |
< |
\dot{j_{y}} & = & \tau^b_y(t) + |
| 1349 |
< |
\left(\overleftrightarrow{\mathsf{I}}_{zz} - \overleftrightarrow{\mathsf{I}}_{xx} \right) j_z j_x,\\ |
| 1350 |
< |
\dot{j_{z}} & = & \tau^b_z(t) + |
| 1351 |
< |
\left(\overleftrightarrow{\mathsf{I}}_{xx} - \overleftrightarrow{\mathsf{I}}_{yy} \right) j_x j_y, |
| 1346 |
> |
\dot{j_{x}} & = & \tau^b_x(t) - |
| 1347 |
> |
\left(\overleftrightarrow{\mathsf{I}}_{yy}^{-1} - \overleftrightarrow{\mathsf{I}}_{zz}^{-1} \right) j_y j_z, \\ |
| 1348 |
> |
\dot{j_{y}} & = & \tau^b_y(t) - |
| 1349 |
> |
\left(\overleftrightarrow{\mathsf{I}}_{zz}^{-1} - \overleftrightarrow{\mathsf{I}}_{xx}^{-1} \right) j_z j_x,\\ |
| 1350 |
> |
\dot{j_{z}} & = & \tau^b_z(t) - |
| 1351 |
> |
\left(\overleftrightarrow{\mathsf{I}}_{xx}^{-1} - \overleftrightarrow{\mathsf{I}}_{yy}^{-1} \right) j_x j_y, |
| 1352 |
|
\end{eqnarray} |
| 1353 |
|
which utilize the body-fixed torques, ${\bf \tau}^b$. Torques are |
| 1354 |
|
most easily derived in the space-fixed frame, |
| 1366 |
|
Here $\hat{\bf u}$ is a unit vector pointing along the principal axis |
| 1367 |
|
of the particle in the space-fixed frame. |
| 1368 |
|
|
| 1369 |
< |
The DLM method uses a Trotter factorization of the orientational |
| 1369 |
> |
The {\sc dlm} method uses a Trotter factorization of the orientational |
| 1370 |
|
propagator. This has three effects: |
| 1371 |
|
\begin{enumerate} |
| 1372 |
|
\item the integrator is area-preserving in phase space (i.e. it is |
| 1453 |
|
|
| 1454 |
|
{\tt moveB:} |
| 1455 |
|
\begin{align*} |
| 1456 |
< |
{\bf v}\left(t + h \right) &\leftarrow {\bf v}\left(t + h / 2 \right) |
| 1457 |
< |
+ \frac{h}{2} \left( {\bf f}(t + h) / m \right), \\ |
| 1458 |
< |
% |
| 1459 |
< |
{\bf j}\left(t + h \right) &\leftarrow {\bf j}\left(t + h / 2 \right) |
| 1460 |
< |
+ \frac{h}{2} {\bf \tau}^b(t + h) . |
| 1461 |
< |
\end{align*} |
| 1456 |
> |
{\bf v}\left(t + h \right) &\leftarrow {\bf v}\left(t + h / 2 \right) |
| 1457 |
> |
+ \frac{h}{2} \left( {\bf f}(t + h) / m \right), \\ |
| 1458 |
> |
% |
| 1459 |
> |
{\bf j}\left(t + h \right) &\leftarrow {\bf j}\left(t + h / 2 \right) |
| 1460 |
> |
+ \frac{h}{2} {\bf \tau}^b(t + h) . |
| 1461 |
> |
\end{align*} |
| 1462 |
> |
|
| 1463 |
> |
The matrix rotations used in the {\sc dlm} method end up being more |
| 1464 |
> |
costly computationally than the simpler arithmetic quaternion |
| 1465 |
> |
propagation. With the same time step, a 1024-molecule water simulation |
| 1466 |
> |
incurs an average 12\% increase in computation time using the {\sc |
| 1467 |
> |
dlm} method in place of quaternions. This cost is more than justified |
| 1468 |
> |
when comparing the energy conservation achieved by the two |
| 1469 |
> |
methods. Figure ~\ref{quatdlm} provides a comparative analysis of the |
| 1470 |
> |
{\sc dlm} method versus the traditional quaternion scheme. |
| 1471 |
> |
|
| 1472 |
> |
\begin{figure} |
| 1473 |
> |
\centering |
| 1474 |
> |
\includegraphics[width=\linewidth]{quatvsdlm.eps} |
| 1475 |
> |
\caption[Energy conservation analysis of the {\sc dlm} and quaternion |
| 1476 |
> |
integration methods]{Analysis of the energy conservation of the {\sc |
| 1477 |
> |
dlm} and quaternion integration methods. $\delta \mathrm{E}_1$ is the |
| 1478 |
> |
linear drift in energy over time and $\delta \mathrm{E}_0$ is the |
| 1479 |
> |
standard deviation of energy fluctuations around this drift. All |
| 1480 |
> |
simulations were of a 1024-molecule simulation of SSD water at 298 K |
| 1481 |
> |
starting from the same initial configuration. Note that the {\sc dlm} |
| 1482 |
> |
method provides more than an order of magnitude improvement in both |
| 1483 |
> |
the energy drift and the size of the energy fluctuations when compared |
| 1484 |
> |
with the quaternion method at any given time step. At time steps |
| 1485 |
> |
larger than 4 fs, the quaternion scheme resulted in rapidly rising |
| 1486 |
> |
energies which eventually lead to simulation failure. Using the {\sc |
| 1487 |
> |
dlm} method, time steps up to 8 fs can be taken before this behavior |
| 1488 |
> |
is evident.} |
| 1489 |
> |
\label{quatdlm} |
| 1490 |
> |
\end{figure} |
| 1491 |
> |
|
| 1492 |
> |
In Fig.~\ref{quatdlm}, $\delta \mbox{E}_1$ is a measure of the linear |
| 1493 |
> |
energy drift in units of $\mbox{kcal mol}^{-1}$ per particle over a |
| 1494 |
> |
nanosecond of simulation time, and $\delta \mbox{E}_0$ is the standard |
| 1495 |
> |
deviation of the energy fluctuations in units of $\mbox{kcal |
| 1496 |
> |
mol}^{-1}$ per particle. In the top plot, it is apparent that the |
| 1497 |
> |
energy drift is reduced by a significant amount (2 to 3 orders of |
| 1498 |
> |
magnitude improvement at all tested time steps) by chosing the {\sc |
| 1499 |
> |
dlm} method over the simple non-symplectic quaternion integration |
| 1500 |
> |
method. In addition to this improvement in energy drift, the |
| 1501 |
> |
fluctuations in the total energy are also dampened by 1 to 2 orders of |
| 1502 |
> |
magnitude by utilizing the {\sc dlm} method. |
| 1503 |
> |
|
| 1504 |
> |
Although the {\sc dlm} method is more computationally expensive than |
| 1505 |
> |
the traditional quaternion scheme for propagating a single time step, |
| 1506 |
> |
consideration of the computational cost for a long simulation with a |
| 1507 |
> |
particular level of energy conservation is in order. A plot of energy |
| 1508 |
> |
drift versus computational cost was generated |
| 1509 |
> |
(Fig.~\ref{cpuCost}). This figure provides an estimate of the CPU time |
| 1510 |
> |
required under the two integration schemes for 1 nanosecond of |
| 1511 |
> |
simulation time for the model 1024-molecule system. By chosing a |
| 1512 |
> |
desired energy drift value it is possible to determine the CPU time |
| 1513 |
> |
required for both methods. If a $\delta \mbox{E}_1$ of $1 \times |
| 1514 |
> |
10^{-3} \mbox{kcal mol}^{-1}$ per particle is desired, a nanosecond of |
| 1515 |
> |
simulation time will require ~19 hours of CPU time with the {\sc dlm} |
| 1516 |
> |
integrator, while the quaternion scheme will require ~154 hours of CPU |
| 1517 |
> |
time. This demonstrates the computational advantage of the integration |
| 1518 |
> |
scheme utilized in {\sc oopse}. |
| 1519 |
|
|
| 1069 |
– |
The matrix rotations used in the DLM method end up being more costly |
| 1070 |
– |
computationally than the simpler arithmetic quaternion |
| 1071 |
– |
propagation. With the same time step, a 1000-molecule water simulation |
| 1072 |
– |
shows an average 7\% increase in computation time using the DLM method |
| 1073 |
– |
in place of quaternions. This cost is more than justified when |
| 1074 |
– |
comparing the energy conservation of the two methods as illustrated in |
| 1075 |
– |
Fig.~\ref{timestep}. |
| 1076 |
– |
|
| 1520 |
|
\begin{figure} |
| 1521 |
|
\centering |
| 1522 |
< |
\includegraphics[width=\linewidth]{timeStep.pdf} |
| 1523 |
< |
\caption[Energy conservation for quaternion versus DLM dynamics]{Energy conservation using quaternion based integration versus |
| 1524 |
< |
the method proposed by Dullweber \emph{et al.} with increasing time |
| 1525 |
< |
step. For each time step, the dotted line is total energy using the |
| 1526 |
< |
DLM integrator, and the solid line comes from the quaternion |
| 1527 |
< |
integrator. The larger time step plots are shifted up from the true |
| 1528 |
< |
energy baseline for clarity.} |
| 1529 |
< |
\label{timestep} |
| 1522 |
> |
\includegraphics[width=\linewidth]{compCost.eps} |
| 1523 |
> |
\caption[Energy drift as a function of required simulation run |
| 1524 |
> |
time]{Energy drift as a function of required simulation run time. |
| 1525 |
> |
$\delta \mathrm{E}_1$ is the linear drift in energy over time. |
| 1526 |
> |
Simulations were performed on a single 2.5 GHz Pentium 4 |
| 1527 |
> |
processor. Simulation time comparisons can be made by tracing |
| 1528 |
> |
horizontally from one curve to the other. For example, a simulation |
| 1529 |
> |
that takes ~24 hours using the {\sc dlm} method will take roughly 210 |
| 1530 |
> |
hours using the simple quaternion method if the same degree of energy |
| 1531 |
> |
conservation is desired.} |
| 1532 |
> |
\label{cpuCost} |
| 1533 |
|
\end{figure} |
| 1534 |
|
|
| 1089 |
– |
In Fig.~\ref{timestep}, the resulting energy drift at various time |
| 1090 |
– |
steps for both the DLM and quaternion integration schemes is |
| 1091 |
– |
compared. All of the 1000 molecule water simulations started with the |
| 1092 |
– |
same configuration, and the only difference was the method for |
| 1093 |
– |
handling rotational motion. At time steps of 0.1 and 0.5 fs, both |
| 1094 |
– |
methods for propagating molecule rotation conserve energy fairly well, |
| 1095 |
– |
with the quaternion method showing a slight energy drift over time in |
| 1096 |
– |
the 0.5 fs time step simulation. At time steps of 1 and 2 fs, the |
| 1097 |
– |
energy conservation benefits of the DLM method are clearly |
| 1098 |
– |
demonstrated. Thus, while maintaining the same degree of energy |
| 1099 |
– |
conservation, one can take considerably longer time steps, leading to |
| 1100 |
– |
an overall reduction in computation time. |
| 1101 |
– |
|
| 1535 |
|
There is only one specific keyword relevant to the default integrator, |
| 1536 |
|
and that is the time step for integrating the equations of motion. |
| 1537 |
|
|
| 1538 |
|
\begin{center} |
| 1539 |
|
\begin{tabular}{llll} |
| 1540 |
< |
{\bf variable} & {\bf {\tt .bass} keyword} & {\bf units} & {\bf |
| 1540 |
> |
{\bf variable} & {\bf Meta-data keyword} & {\bf units} & {\bf |
| 1541 |
|
default value} \\ |
| 1542 |
|
$h$ & {\tt dt = 2.0;} & fs & none |
| 1543 |
|
\end{tabular} |
| 1547 |
|
|
| 1548 |
|
{\sc oopse} implements a number of extended system integrators for |
| 1549 |
|
sampling from other ensembles relevant to chemical physics. The |
| 1550 |
< |
integrator can selected with the {\tt ensemble} keyword in the |
| 1551 |
< |
{\tt .bass} file: |
| 1550 |
> |
integrator can be selected with the {\tt ensemble} keyword in the |
| 1551 |
> |
meta-data file: |
| 1552 |
|
|
| 1553 |
|
\begin{center} |
| 1554 |
|
\begin{tabular}{lll} |
| 1555 |
< |
{\bf Integrator} & {\bf Ensemble} & {\bf {\tt .bass} line} \\ |
| 1555 |
> |
{\bf Integrator} & {\bf Ensemble} & {\bf Meta-data instruction} \\ |
| 1556 |
|
NVE & microcanonical & {\tt ensemble = NVE; } \\ |
| 1557 |
|
NVT & canonical & {\tt ensemble = NVT; } \\ |
| 1558 |
|
NPTi & isobaric-isothermal & {\tt ensemble = NPTi;} \\ |
| 1567 |
|
The relatively well-known Nos\'e-Hoover thermostat\cite{Hoover85} is |
| 1568 |
|
implemented in {\sc oopse}'s NVT integrator. This method couples an |
| 1569 |
|
extra degree of freedom (the thermostat) to the kinetic energy of the |
| 1570 |
< |
system, and has been shown to sample the canonical distribution in the |
| 1571 |
< |
system degrees of freedom while conserving a quantity that is, to |
| 1570 |
> |
system and it has been shown to sample the canonical distribution in |
| 1571 |
> |
the system degrees of freedom while conserving a quantity that is, to |
| 1572 |
|
within a constant, the Helmholtz free energy.\cite{melchionna93} |
| 1573 |
|
|
| 1574 |
|
NPT algorithms attempt to maintain constant pressure in the system by |
| 1592 |
|
|
| 1593 |
|
\begin{center} |
| 1594 |
|
\begin{tabular}{llll} |
| 1595 |
< |
{\bf variable} & {\bf {\tt .bass} keyword} & {\bf units} & {\bf |
| 1595 |
> |
{\bf variable} & {\bf Meta-data instruction} & {\bf units} & {\bf |
| 1596 |
|
default value} \\ |
| 1597 |
|
$T_{\mathrm{target}}$ & {\tt targetTemperature = 300;} & K & none \\ |
| 1598 |
|
$P_{\mathrm{target}}$ & {\tt targetPressure = 1;} & atm & none \\ |
| 1638 |
|
\end{equation} |
| 1639 |
|
Here, $f$ is the total number of degrees of freedom in the system, |
| 1640 |
|
\begin{equation} |
| 1641 |
< |
f = 3 N + 3 N_{\mathrm{orient}} - N_{\mathrm{constraints}}, |
| 1641 |
> |
f = 3 N + 2 N_{\mathrm{linear}} + 3 N_{\mathrm{non-linear}} - N_{\mathrm{constraints}}, |
| 1642 |
|
\end{equation} |
| 1643 |
|
and $K$ is the total kinetic energy, |
| 1644 |
|
\begin{equation} |
| 1645 |
|
K = \sum_{i=1}^{N} \frac{1}{2} m_i {\bf v}_i^T \cdot {\bf v}_i + |
| 1646 |
< |
\sum_{i=1}^{N_{\mathrm{orient}}} \frac{1}{2} {\bf j}_i^T \cdot |
| 1646 |
> |
\sum_{i=1}^{N_{\mathrm{linear}}+N_{\mathrm{non-linear}}} \frac{1}{2} {\bf j}_i^T \cdot |
| 1647 |
|
\overleftrightarrow{\mathsf{I}}_i^{-1} \cdot {\bf j}_i. |
| 1648 |
|
\end{equation} |
| 1649 |
+ |
$N_{\mathrm{linear}}$ is the number of linear rotors (i.e. with two |
| 1650 |
+ |
non-zero moments of inertia), and $N_{\mathrm{non-linear}}$ is the |
| 1651 |
+ |
number of non-linear rotors (i.e. with three non-zero moments of |
| 1652 |
+ |
inertia). |
| 1653 |
|
|
| 1654 |
|
In eq.(\ref{eq:nosehooverext}), $\tau_T$ is the time constant for |
| 1655 |
|
relaxation of the temperature to the target value. To set values for |
| 1656 |
|
$\tau_T$ or $T_{\mathrm{target}}$ in a simulation, one would use the |
| 1657 |
< |
{\tt tauThermostat} and {\tt targetTemperature} keywords in the {\tt |
| 1658 |
< |
.bass} file. The units for {\tt tauThermostat} are fs, and the units |
| 1659 |
< |
for the {\tt targetTemperature} are degrees K. The integration of |
| 1660 |
< |
the equations of motion is carried out in a velocity-Verlet style 2 |
| 1657 |
> |
{\tt tauThermostat} and {\tt targetTemperature} keywords in the |
| 1658 |
> |
meta-data file. The units for {\tt tauThermostat} are fs, and the |
| 1659 |
> |
units for the {\tt targetTemperature} are degrees K. The integration |
| 1660 |
> |
of the equations of motion is carried out in a velocity-Verlet style 2 |
| 1661 |
|
part algorithm: |
| 1662 |
|
|
| 1663 |
|
{\tt moveA:} |
| 1687 |
|
Here $\mathrm{rotate}(h * {\bf j} |
| 1688 |
|
\overleftrightarrow{\mathsf{I}}^{-1})$ is the same symplectic Trotter |
| 1689 |
|
factorization of the three rotation operations that was discussed in |
| 1690 |
< |
the section on the DLM integrator. Note that this operation modifies |
| 1690 |
> |
the section on the {\sc dlm} integrator. Note that this operation modifies |
| 1691 |
|
both the rotation matrix $\mathsf{A}$ and the angular momentum ${\bf |
| 1692 |
|
j}$. {\tt moveA} propagates velocities by a half time step, and |
| 1693 |
|
positional degrees of freedom by a full time step. The new positions |
| 1694 |
|
(and orientations) are then used to calculate a new set of forces and |
| 1695 |
|
torques in exactly the same way they are calculated in the {\tt |
| 1696 |
< |
doForces} portion of the DLM integrator. |
| 1696 |
> |
doForces} portion of the {\sc dlm} integrator. |
| 1697 |
|
|
| 1698 |
|
Once the forces and torques have been obtained at the new time step, |
| 1699 |
|
the temperature, velocities, and the extended system variable can be |
| 1719 |
|
\chi(t + h) \right) . |
| 1720 |
|
\end{align*} |
| 1721 |
|
|
| 1722 |
< |
Since ${\bf v}(t + h)$ and ${\bf j}(t + h)$ are required to caclculate |
| 1722 |
> |
Since ${\bf v}(t + h)$ and ${\bf j}(t + h)$ are required to calculate |
| 1723 |
|
$T(t + h)$ as well as $\chi(t + h)$, they indirectly depend on their |
| 1724 |
|
own values at time $t + h$. {\tt moveB} is therefore done in an |
| 1725 |
|
iterative fashion until $\chi(t + h)$ becomes self-consistent. The |
| 1747 |
|
\subsection{\label{sec:NPTi}Constant-pressure integration with |
| 1748 |
|
isotropic box deformations (NPTi)} |
| 1749 |
|
|
| 1750 |
< |
To carry out isobaric-isothermal ensemble calculations {\sc oopse} |
| 1750 |
> |
To carry out isobaric-isothermal ensemble calculations, {\sc oopse} |
| 1751 |
|
implements the Melchionna modifications to the Nos\'e-Hoover-Andersen |
| 1752 |
< |
equations of motion,\cite{melchionna93} |
| 1752 |
> |
equations of motion.\cite{melchionna93} The equations of motion are |
| 1753 |
> |
the same as NVT with the following exceptions: |
| 1754 |
|
|
| 1755 |
|
\begin{eqnarray} |
| 1756 |
|
\dot{{\bf r}} & = & {\bf v} + \eta \left( {\bf r} - {\bf R}_0 \right), \\ |
| 1757 |
|
\dot{{\bf v}} & = & \frac{{\bf f}}{m} - (\eta + \chi) {\bf v}, \\ |
| 1320 |
– |
\dot{\mathsf{A}} & = & \mathsf{A} \cdot |
| 1321 |
– |
\mbox{ skew}\left(\overleftrightarrow{I}^{-1} \cdot {\bf j}\right),\\ |
| 1322 |
– |
\dot{{\bf j}} & = & {\bf j} \times \left( \overleftrightarrow{I}^{-1} |
| 1323 |
– |
\cdot {\bf j} \right) - \mbox{ rot}\left(\mathsf{A}^{T} \cdot \frac{\partial |
| 1324 |
– |
V}{\partial \mathsf{A}} \right) - \chi {\bf j}, \\ |
| 1325 |
– |
\dot{\chi} & = & \frac{1}{\tau_{T}^2} \left( |
| 1326 |
– |
\frac{T}{T_{\mathrm{target}}} - 1 \right) ,\\ |
| 1758 |
|
\dot{\eta} & = & \frac{1}{\tau_{B}^2 f k_B T_{\mathrm{target}}} V \left( P - |
| 1759 |
|
P_{\mathrm{target}} \right), \\ |
| 1760 |
|
\dot{\mathcal{V}} & = & 3 \mathcal{V} \eta . |
| 1781 |
|
\overleftrightarrow{\mathsf{W}}(t). |
| 1782 |
|
\end{equation} |
| 1783 |
|
The kinetic contribution to the pressure tensor utilizes the {\it |
| 1784 |
< |
outer} product of the velocities denoted by the $\otimes$ symbol. The |
| 1784 |
> |
outer} product of the velocities, denoted by the $\otimes$ symbol. The |
| 1785 |
|
stress tensor is calculated from another outer product of the |
| 1786 |
|
inter-atomic separation vectors (${\bf r}_{ij} = {\bf r}_j - {\bf |
| 1787 |
|
r}_i$) with the forces between the same two atoms, |
| 1789 |
|
\overleftrightarrow{\mathsf{W}}(t) = \sum_{i} \sum_{j>i} {\bf r}_{ij}(t) |
| 1790 |
|
\otimes {\bf f}_{ij}(t). |
| 1791 |
|
\end{equation} |
| 1792 |
+ |
In systems containing cutoff groups, the stress tensor is computed |
| 1793 |
+ |
between the centers-of-mass of the cutoff groups: |
| 1794 |
+ |
\begin{equation} |
| 1795 |
+ |
\overleftrightarrow{\mathsf{W}}(t) = \sum_{a} \sum_{b} {\bf r}_{ab}(t) |
| 1796 |
+ |
\otimes {\bf f}_{ab}(t). |
| 1797 |
+ |
\end{equation} |
| 1798 |
+ |
where ${\bf r}_{ab}$ is the distance between the centers of mass, and |
| 1799 |
+ |
\begin{equation} |
| 1800 |
+ |
{\bf f}_{ab} = s(r_{ab}) \sum_{i \in a} \sum_{j \in b} {\bf f}_{ij} + |
| 1801 |
+ |
s^{\prime}(r_{ab}) \frac{{\bf r}_{ab}}{|r_{ab}|} \sum_{i \in a} \sum_{j |
| 1802 |
+ |
\in b} V_{ij}({\bf r}_{ij}). |
| 1803 |
+ |
\end{equation} |
| 1804 |
+ |
|
| 1805 |
|
The instantaneous pressure is then simply obtained from the trace of |
| 1806 |
< |
the Pressure tensor, |
| 1806 |
> |
the pressure tensor, |
| 1807 |
|
\begin{equation} |
| 1808 |
< |
P(t) = \frac{1}{3} \mathrm{Tr} \left( \overleftrightarrow{\mathsf{P}}(t). |
| 1809 |
< |
\right) |
| 1808 |
> |
P(t) = \frac{1}{3} \mathrm{Tr} \left( \overleftrightarrow{\mathsf{P}}(t) |
| 1809 |
> |
\right). |
| 1810 |
|
\end{equation} |
| 1811 |
|
|
| 1812 |
|
In eq.(\ref{eq:melchionna1}), $\tau_B$ is the time constant for |
| 1813 |
|
relaxation of the pressure to the target value. To set values for |
| 1814 |
|
$\tau_B$ or $P_{\mathrm{target}}$ in a simulation, one would use the |
| 1815 |
< |
{\tt tauBarostat} and {\tt targetPressure} keywords in the {\tt .bass} |
| 1815 |
> |
{\tt tauBarostat} and {\tt targetPressure} keywords in the meta-data |
| 1816 |
|
file. The units for {\tt tauBarostat} are fs, and the units for the |
| 1817 |
|
{\tt targetPressure} are atmospheres. Like in the NVT integrator, the |
| 1818 |
|
integration of the equations of motion is carried out in a |
| 1819 |
< |
velocity-Verlet style 2 part algorithm: |
| 1819 |
> |
velocity-Verlet style two part algorithm with only the following |
| 1820 |
> |
differences: |
| 1821 |
|
|
| 1822 |
|
{\tt moveA:} |
| 1823 |
|
\begin{align*} |
| 1379 |
– |
T(t) &\leftarrow \left\{{\bf v}(t)\right\}, \left\{{\bf j}(t)\right\} ,\\ |
| 1380 |
– |
% |
| 1824 |
|
P(t) &\leftarrow \left\{{\bf r}(t)\right\}, \left\{{\bf v}(t)\right\} ,\\ |
| 1825 |
|
% |
| 1826 |
|
{\bf v}\left(t + h / 2\right) &\leftarrow {\bf v}(t) |
| 1827 |
|
+ \frac{h}{2} \left( \frac{{\bf f}(t)}{m} - {\bf v}(t) |
| 1828 |
|
\left(\chi(t) + \eta(t) \right) \right), \\ |
| 1386 |
– |
% |
| 1387 |
– |
{\bf j}\left(t + h / 2 \right) &\leftarrow {\bf j}(t) |
| 1388 |
– |
+ \frac{h}{2} \left( {\bf \tau}^b(t) - {\bf j}(t) |
| 1389 |
– |
\chi(t) \right), \\ |
| 1390 |
– |
% |
| 1391 |
– |
\mathsf{A}(t + h) &\leftarrow \mathrm{rotate}\left(h * |
| 1392 |
– |
{\bf j}(t + h / 2) \overleftrightarrow{\mathsf{I}}^{-1} |
| 1393 |
– |
\right) ,\\ |
| 1829 |
|
% |
| 1395 |
– |
\chi\left(t + h / 2 \right) &\leftarrow \chi(t) + |
| 1396 |
– |
\frac{h}{2 \tau_T^2} \left( \frac{T(t)}{T_{\mathrm{target}}} - 1 |
| 1397 |
– |
\right) ,\\ |
| 1398 |
– |
% |
| 1830 |
|
\eta(t + h / 2) &\leftarrow \eta(t) + \frac{h |
| 1831 |
|
\mathcal{V}(t)}{2 N k_B T(t) \tau_B^2} \left( P(t) |
| 1832 |
|
- P_{\mathrm{target}} \right), \\ |
| 1840 |
|
\mathsf{H}(t). |
| 1841 |
|
\end{align*} |
| 1842 |
|
|
| 1843 |
< |
Most of these equations are identical to their counterparts in the NVT |
| 1413 |
< |
integrator, but the propagation of positions to time $t + h$ |
| 1843 |
> |
The propagation of positions to time $t + h$ |
| 1844 |
|
depends on the positions at the same time. {\sc oopse} carries out |
| 1845 |
|
this step iteratively (with a limit of 5 passes through the iterative |
| 1846 |
|
loop). Also, the simulation box $\mathsf{H}$ is scaled uniformly for |
| 1849 |
|
h / 2$. Reshaping the box uniformly also scales the volume of |
| 1850 |
|
the box by |
| 1851 |
|
\begin{equation} |
| 1852 |
< |
\mathcal{V}(t + h) \leftarrow e^{ - 3 h \eta(t + h /2)}. |
| 1853 |
< |
\mathcal{V}(t) |
| 1852 |
> |
\mathcal{V}(t + h) \leftarrow e^{ - 3 h \eta(t + h /2)} \times |
| 1853 |
> |
\mathcal{V}(t). |
| 1854 |
|
\end{equation} |
| 1855 |
|
|
| 1856 |
|
The {\tt doForces} step for the NPTi integrator is exactly the same as |
| 1857 |
< |
in both the DLM and NVT integrators. Once the forces and torques have |
| 1857 |
> |
in both the {\sc dlm} and NVT integrators. Once the forces and torques have |
| 1858 |
|
been obtained at the new time step, the velocities can be advanced to |
| 1859 |
|
the same time value. |
| 1860 |
|
|
| 1861 |
|
{\tt moveB:} |
| 1862 |
|
\begin{align*} |
| 1433 |
– |
T(t + h) &\leftarrow \left\{{\bf v}(t + h)\right\}, |
| 1434 |
– |
\left\{{\bf j}(t + h)\right\} ,\\ |
| 1435 |
– |
% |
| 1863 |
|
P(t + h) &\leftarrow \left\{{\bf r}(t + h)\right\}, |
| 1864 |
|
\left\{{\bf v}(t + h)\right\}, \\ |
| 1865 |
|
% |
| 1439 |
– |
\chi\left(t + h \right) &\leftarrow \chi\left(t + h / |
| 1440 |
– |
2 \right) + \frac{h}{2 \tau_T^2} \left( \frac{T(t+h)} |
| 1441 |
– |
{T_{\mathrm{target}}} - 1 \right), \\ |
| 1442 |
– |
% |
| 1866 |
|
\eta(t + h) &\leftarrow \eta(t + h / 2) + |
| 1867 |
|
\frac{h \mathcal{V}(t + h)}{2 N k_B T(t + h) |
| 1868 |
|
\tau_B^2} \left( P(t + h) - P_{\mathrm{target}} \right), \\ |
| 1879 |
|
\end{align*} |
| 1880 |
|
|
| 1881 |
|
Once again, since ${\bf v}(t + h)$ and ${\bf j}(t + h)$ are required |
| 1882 |
< |
to caclculate $T(t + h)$, $P(t + h)$, $\chi(t + h)$, and $\eta(t + |
| 1882 |
> |
to calculate $T(t + h)$, $P(t + h)$, $\chi(t + h)$, and $\eta(t + |
| 1883 |
|
h)$, they indirectly depend on their own values at time $t + h$. {\tt |
| 1884 |
|
moveB} is therefore done in an iterative fashion until $\chi(t + h)$ |
| 1885 |
|
and $\eta(t + h)$ become self-consistent. The relative tolerance for |
| 1919 |
|
{\it shape} as well as in the volume of the box. This method utilizes |
| 1920 |
|
the full $3 \times 3$ pressure tensor and introduces a tensor of |
| 1921 |
|
extended variables ($\overleftrightarrow{\eta}$) to control changes to |
| 1922 |
< |
the box shape. The equations of motion for this method are |
| 1922 |
> |
the box shape. The equations of motion for this method differ from |
| 1923 |
> |
those of NPTi as follows: |
| 1924 |
|
\begin{eqnarray} |
| 1925 |
|
\dot{{\bf r}} & = & {\bf v} + \overleftrightarrow{\eta} \cdot \left( {\bf r} - {\bf R}_0 \right), \\ |
| 1926 |
|
\dot{{\bf v}} & = & \frac{{\bf f}}{m} - (\overleftrightarrow{\eta} + |
| 1927 |
|
\chi \cdot \mathsf{1}) {\bf v}, \\ |
| 1504 |
– |
\dot{\mathsf{A}} & = & \mathsf{A} \cdot |
| 1505 |
– |
\mbox{ skew}\left(\overleftrightarrow{I}^{-1} \cdot {\bf j}\right) ,\\ |
| 1506 |
– |
\dot{{\bf j}} & = & {\bf j} \times \left( \overleftrightarrow{I}^{-1} |
| 1507 |
– |
\cdot {\bf j} \right) - \mbox{ rot}\left(\mathsf{A}^{T} \cdot \frac{\partial |
| 1508 |
– |
V}{\partial \mathsf{A}} \right) - \chi {\bf j} ,\\ |
| 1509 |
– |
\dot{\chi} & = & \frac{1}{\tau_{T}^2} \left( |
| 1510 |
– |
\frac{T}{T_{\mathrm{target}}} - 1 \right) ,\\ |
| 1928 |
|
\dot{\overleftrightarrow{\eta}} & = & \frac{1}{\tau_{B}^2 f k_B |
| 1929 |
|
T_{\mathrm{target}}} V \left( \overleftrightarrow{\mathsf{P}} - P_{\mathrm{target}}\mathsf{1} \right) ,\\ |
| 1930 |
|
\dot{\mathsf{H}} & = & \overleftrightarrow{\eta} \cdot \mathsf{H} . |
| 1940 |
|
|
| 1941 |
|
{\tt moveA:} |
| 1942 |
|
\begin{align*} |
| 1526 |
– |
T(t) &\leftarrow \left\{{\bf v}(t)\right\}, \left\{{\bf j}(t)\right\} ,\\ |
| 1527 |
– |
% |
| 1943 |
|
\overleftrightarrow{\mathsf{P}}(t) &\leftarrow \left\{{\bf r}(t)\right\}, |
| 1944 |
|
\left\{{\bf v}(t)\right\} ,\\ |
| 1945 |
|
% |
| 1948 |
|
\left(\chi(t)\mathsf{1} + \overleftrightarrow{\eta}(t) \right) \cdot |
| 1949 |
|
{\bf v}(t) \right), \\ |
| 1950 |
|
% |
| 1536 |
– |
{\bf j}\left(t + h / 2 \right) &\leftarrow {\bf j}(t) |
| 1537 |
– |
+ \frac{h}{2} \left( {\bf \tau}^b(t) - {\bf j}(t) |
| 1538 |
– |
\chi(t) \right), \\ |
| 1539 |
– |
% |
| 1540 |
– |
\mathsf{A}(t + h) &\leftarrow \mathrm{rotate}\left(h * |
| 1541 |
– |
{\bf j}(t + h / 2) \overleftrightarrow{\mathsf{I}}^{-1} |
| 1542 |
– |
\right), \\ |
| 1543 |
– |
% |
| 1544 |
– |
\chi\left(t + h / 2 \right) &\leftarrow \chi(t) + |
| 1545 |
– |
\frac{h}{2 \tau_T^2} \left( \frac{T(t)}{T_{\mathrm{target}}} |
| 1546 |
– |
- 1 \right), \\ |
| 1547 |
– |
% |
| 1951 |
|
\overleftrightarrow{\eta}(t + h / 2) &\leftarrow |
| 1952 |
|
\overleftrightarrow{\eta}(t) + \frac{h \mathcal{V}(t)}{2 N k_B |
| 1953 |
|
T(t) \tau_B^2} \left( \overleftrightarrow{\mathsf{P}}(t) |
| 1969 |
|
|
| 1970 |
|
{\tt moveB:} |
| 1971 |
|
\begin{align*} |
| 1569 |
– |
T(t + h) &\leftarrow \left\{{\bf v}(t + h)\right\}, |
| 1570 |
– |
\left\{{\bf j}(t + h)\right\}, \\ |
| 1571 |
– |
% |
| 1972 |
|
\overleftrightarrow{\mathsf{P}}(t + h) &\leftarrow \left\{{\bf r} |
| 1973 |
|
(t + h)\right\}, \left\{{\bf v}(t |
| 1974 |
|
+ h)\right\}, \left\{{\bf f}(t + h)\right\} ,\\ |
| 1975 |
|
% |
| 1576 |
– |
\chi\left(t + h \right) &\leftarrow \chi\left(t + h / |
| 1577 |
– |
2 \right) + \frac{h}{2 \tau_T^2} \left( \frac{T(t+ |
| 1578 |
– |
h)}{T_{\mathrm{target}}} - 1 \right), \\ |
| 1579 |
– |
% |
| 1976 |
|
\overleftrightarrow{\eta}(t + h) &\leftarrow |
| 1977 |
|
\overleftrightarrow{\eta}(t + h / 2) + |
| 1978 |
|
\frac{h \mathcal{V}(t + h)}{2 N k_B T(t + h) |
| 1984 |
|
\frac{{\bf f}(t + h)}{m} - |
| 1985 |
|
(\chi(t + h)\mathsf{1} + \overleftrightarrow{\eta}(t |
| 1986 |
|
+ h)) \right) \cdot {\bf v}(t + h), \\ |
| 1591 |
– |
% |
| 1592 |
– |
{\bf j}\left(t + h \right) &\leftarrow {\bf j}\left(t |
| 1593 |
– |
+ h / 2 \right) + \frac{h}{2} \left( {\bf \tau}^b(t |
| 1594 |
– |
+ h) - {\bf j}(t + h) \chi(t + h) \right) . |
| 1987 |
|
\end{align*} |
| 1988 |
|
|
| 1989 |
|
The iterative schemes for both {\tt moveA} and {\tt moveB} are |
| 2001 |
|
This integrator must be used with care, particularly in liquid |
| 2002 |
|
simulations. Liquids have very small restoring forces in the |
| 2003 |
|
off-diagonal directions, and the simulation box can very quickly form |
| 2004 |
< |
elongated and sheared geometries which become smaller than the |
| 2005 |
< |
electrostatic or Lennard-Jones cutoff radii. The NPTf integrator |
| 2006 |
< |
finds most use in simulating crystals or liquid crystals which assume |
| 1615 |
< |
non-orthorhombic geometries. |
| 2004 |
> |
elongated and sheared geometries which become smaller than the cutoff |
| 2005 |
> |
radius. The NPTf integrator finds most use in simulating crystals or |
| 2006 |
> |
liquid crystals which assume non-orthorhombic geometries. |
| 2007 |
|
|
| 2008 |
|
\subsection{\label{nptxyz}Constant pressure in 3 axes (NPTxyz)} |
| 2009 |
|
|
| 2026 |
|
orientational anisotropy in the system (i.e. in lipid bilayer |
| 2027 |
|
simulations). |
| 2028 |
|
|
| 2029 |
< |
\subsection{\label{oopseSec:rattle}The {\sc rattle} Method for Bond |
| 2029 |
> |
\subsection{\label{sec:constraints}Constraint Methods} |
| 2030 |
> |
|
| 2031 |
> |
\subsubsection{\label{oopseSec:rattle}The {\sc rattle} Method for Bond |
| 2032 |
|
Constraints} |
| 2033 |
|
|
| 2034 |
|
In order to satisfy the constraints of fixed bond lengths within {\sc |
| 2035 |
|
oopse}, we have implemented the {\sc rattle} algorithm of |
| 2036 |
< |
Andersen.\cite{andersen83} The algorithm is a velocity verlet |
| 2037 |
< |
formulation of the {\sc shake} method\cite{ryckaert77} of iteratively |
| 2038 |
< |
solving the Lagrange multipliers of constraint. The system of Lagrange |
| 2039 |
< |
multipliers allows one to reformulate the equations of motion with |
| 2040 |
< |
explicit constraint forces.\cite{fowles99:lagrange} |
| 2036 |
> |
Andersen.\cite{andersen83} {\sc rattle} is a velocity-Verlet |
| 2037 |
> |
formulation of the {\sc shake} method\cite{ryckaert77} for iteratively |
| 2038 |
> |
solving the Lagrange multipliers which maintain the holonomic |
| 2039 |
> |
constraints. Both methods are covered in depth in the |
| 2040 |
> |
literature,\cite{leach01:mm,allen87:csl} and a detailed description of |
| 2041 |
> |
this method would be redundant. |
| 2042 |
|
|
| 2043 |
< |
Consider a system described by coordinates $q_1$ and $q_2$ subject to an |
| 1650 |
< |
equation of constraint: |
| 1651 |
< |
\begin{equation} |
| 1652 |
< |
\sigma(q_1, q_2,t) = 0 |
| 1653 |
< |
\label{oopseEq:lm1} |
| 1654 |
< |
\end{equation} |
| 1655 |
< |
The Lagrange formulation of the equations of motion can be written: |
| 1656 |
< |
\begin{equation} |
| 1657 |
< |
\delta\int_{t_1}^{t_2}L\, dt = |
| 1658 |
< |
\int_{t_1}^{t_2} \sum_i \biggl [ \frac{\partial L}{\partial q_i} |
| 1659 |
< |
- \frac{d}{dt}\biggl(\frac{\partial L}{\partial \dot{q}_i} |
| 1660 |
< |
\biggr ) \biggr] \delta q_i \, dt = 0. |
| 1661 |
< |
\label{oopseEq:lm2} |
| 1662 |
< |
\end{equation} |
| 1663 |
< |
Here, $\delta q_i$ is not independent for each $q$, as $q_1$ and $q_2$ |
| 1664 |
< |
are linked by $\sigma$. However, $\sigma$ is fixed at any given |
| 1665 |
< |
instant of time, giving: |
| 1666 |
< |
\begin{align} |
| 1667 |
< |
\delta\sigma &= \biggl( \frac{\partial\sigma}{\partial q_1} \delta q_1 % |
| 1668 |
< |
+ \frac{\partial\sigma}{\partial q_2} \delta q_2 \biggr) = 0 ,\\ |
| 1669 |
< |
% |
| 1670 |
< |
\frac{\partial\sigma}{\partial q_1} \delta q_1 &= % |
| 1671 |
< |
- \frac{\partial\sigma}{\partial q_2} \delta q_2, \\ |
| 1672 |
< |
% |
| 1673 |
< |
\delta q_2 &= - \biggl(\frac{\partial\sigma}{\partial q_1} \bigg / % |
| 1674 |
< |
\frac{\partial\sigma}{\partial q_2} \biggr) \delta q_1. |
| 1675 |
< |
\end{align} |
| 1676 |
< |
Substituted back into Eq.~\ref{oopseEq:lm2}, |
| 1677 |
< |
\begin{equation} |
| 1678 |
< |
\int_{t_1}^{t_2}\biggl [ \biggl(\frac{\partial L}{\partial q_1} |
| 1679 |
< |
- \frac{d}{dt}\,\frac{\partial L}{\partial \dot{q}_1} |
| 1680 |
< |
\biggr) |
| 1681 |
< |
- \biggl( \frac{\partial L}{\partial q_1} |
| 1682 |
< |
- \frac{d}{dt}\,\frac{\partial L}{\partial \dot{q}_1} |
| 1683 |
< |
\biggr) \biggl(\frac{\partial\sigma}{\partial q_1} \bigg / % |
| 1684 |
< |
\frac{\partial\sigma}{\partial q_2} \biggr)\biggr] \delta q_1 \, dt = 0. |
| 1685 |
< |
\label{oopseEq:lm3} |
| 1686 |
< |
\end{equation} |
| 1687 |
< |
Leading to, |
| 1688 |
< |
\begin{equation} |
| 1689 |
< |
\frac{\biggl(\frac{\partial L}{\partial q_1} |
| 1690 |
< |
- \frac{d}{dt}\,\frac{\partial L}{\partial \dot{q}_1} |
| 1691 |
< |
\biggr)}{\frac{\partial\sigma}{\partial q_1}} = |
| 1692 |
< |
\frac{\biggl(\frac{\partial L}{\partial q_2} |
| 1693 |
< |
- \frac{d}{dt}\,\frac{\partial L}{\partial \dot{q}_2} |
| 1694 |
< |
\biggr)}{\frac{\partial\sigma}{\partial q_2}}. |
| 1695 |
< |
\label{oopseEq:lm4} |
| 1696 |
< |
\end{equation} |
| 1697 |
< |
This relation can only be statisfied, if both are equal to a single |
| 1698 |
< |
function $-\lambda(t)$, |
| 1699 |
< |
\begin{align} |
| 1700 |
< |
\frac{\biggl(\frac{\partial L}{\partial q_1} |
| 1701 |
< |
- \frac{d}{dt}\,\frac{\partial L}{\partial \dot{q}_1} |
| 1702 |
< |
\biggr)}{\frac{\partial\sigma}{\partial q_1}} &= -\lambda(t), \\ |
| 1703 |
< |
% |
| 1704 |
< |
\frac{\partial L}{\partial q_1} |
| 1705 |
< |
- \frac{d}{dt}\,\frac{\partial L}{\partial \dot{q}_1} &= |
| 1706 |
< |
-\lambda(t)\,\frac{\partial\sigma}{\partial q_1} ,\\ |
| 1707 |
< |
% |
| 1708 |
< |
\frac{\partial L}{\partial q_1} |
| 1709 |
< |
- \frac{d}{dt}\,\frac{\partial L}{\partial \dot{q}_1} |
| 1710 |
< |
+ \mathcal{G}_i &= 0, |
| 1711 |
< |
\end{align} |
| 1712 |
< |
where $\mathcal{G}_i$, the force of constraint on $i$, is: |
| 1713 |
< |
\begin{equation} |
| 1714 |
< |
\mathcal{G}_i = \lambda(t)\,\frac{\partial\sigma}{\partial q_1}. |
| 1715 |
< |
\label{oopseEq:lm5} |
| 1716 |
< |
\end{equation} |
| 2043 |
> |
\subsubsection{\label{oopseSec:zcons}The Z-Constraint Method} |
| 2044 |
|
|
| 2045 |
< |
In a simulation, this would involve the solution of a set of $(m + n)$ |
| 2046 |
< |
number of equations. Where $m$ is the number of constraints, and $n$ |
| 1720 |
< |
is the number of constrained coordinates. In practice, this is not |
| 1721 |
< |
done, as the matrix inversion necessary to solve the system of |
| 1722 |
< |
equations would be very time consuming to solve. Additionally, the |
| 1723 |
< |
numerical error in the solution of the set of $\lambda$'s would be |
| 1724 |
< |
compounded by the error inherent in propagating by the Velocity Verlet |
| 1725 |
< |
algorithm ($\Delta t^4$). The Verlet propagation error is negligible |
| 1726 |
< |
in an unconstrained system, as one is interested in the statistics of |
| 1727 |
< |
the run, and not that the run be numerically exact to the ``true'' |
| 1728 |
< |
integration. This relates back to the ergodic hypothesis that a time |
| 1729 |
< |
integral of a valid trajectory will still give the correct ensemble |
| 1730 |
< |
average. However, in the case of constraints, if the equations of |
| 1731 |
< |
motion leave the ``true'' trajectory, they are departing from the |
| 1732 |
< |
constrained surface. The method that is used, is to iteratively solve |
| 1733 |
< |
for $\lambda(t)$ at each time step. |
| 1734 |
< |
|
| 1735 |
< |
In {\sc rattle} the equations of motion are modified subject to the |
| 1736 |
< |
following two constraints: |
| 1737 |
< |
\begin{align} |
| 1738 |
< |
\sigma_{ij}[\mathbf{r}(t)] \equiv |
| 1739 |
< |
[ \mathbf{r}_i(t) - \mathbf{r}_j(t)]^2 - d_{ij}^2 &= 0 % |
| 1740 |
< |
\label{oopseEq:c1}, \\ |
| 1741 |
< |
% |
| 1742 |
< |
[\mathbf{\dot{r}}_i(t) - \mathbf{\dot{r}}_j(t)] \cdot |
| 1743 |
< |
[\mathbf{r}_i(t) - \mathbf{r}_j(t)] &= 0 .\label{oopseEq:c2} |
| 1744 |
< |
\end{align} |
| 1745 |
< |
Eq.~\ref{oopseEq:c1} is the set of bond constraints, where $d_{ij}$ is |
| 1746 |
< |
the constrained distance between atom $i$ and |
| 1747 |
< |
$j$. Eq.~\ref{oopseEq:c2} constrains the velocities of $i$ and $j$ to |
| 1748 |
< |
be perpendicular to the bond vector, so that the bond can neither grow |
| 1749 |
< |
nor shrink. The constrained dynamics equations become: |
| 1750 |
< |
\begin{equation} |
| 1751 |
< |
m_i \mathbf{\ddot{r}}_i = \mathbf{F}_i + \mathbf{\mathcal{G}}_i, |
| 1752 |
< |
\label{oopseEq:r1} |
| 1753 |
< |
\end{equation} |
| 1754 |
< |
where,$\mathbf{\mathcal{G}}_i$ are the forces of constraint on $i$, |
| 1755 |
< |
and are defined: |
| 1756 |
< |
\begin{equation} |
| 1757 |
< |
\mathbf{\mathcal{G}}_i = - \sum_j \lambda_{ij}(t)\,\nabla \sigma_{ij}. |
| 1758 |
< |
\label{oopseEq:r2} |
| 1759 |
< |
\end{equation} |
| 1760 |
< |
|
| 1761 |
< |
In Velocity Verlet, if $\Delta t = h$, the propagation can be written: |
| 1762 |
< |
\begin{align} |
| 1763 |
< |
\mathbf{r}_i(t+h) &= |
| 1764 |
< |
\mathbf{r}_i(t) + h\mathbf{\dot{r}}(t) + |
| 1765 |
< |
\frac{h^2}{2m_i}\,\Bigl[ \mathbf{F}_i(t) + |
| 1766 |
< |
\mathbf{\mathcal{G}}_{Ri}(t) \Bigr] \label{oopseEq:vv1}, \\ |
| 1767 |
< |
% |
| 1768 |
< |
\mathbf{\dot{r}}_i(t+h) &= |
| 1769 |
< |
\mathbf{\dot{r}}_i(t) + \frac{h}{2m_i} |
| 1770 |
< |
\Bigl[ \mathbf{F}_i(t) + \mathbf{\mathcal{G}}_{Ri}(t) + |
| 1771 |
< |
\mathbf{F}_i(t+h) + \mathbf{\mathcal{G}}_{Vi}(t+h) \Bigr] ,% |
| 1772 |
< |
\label{oopseEq:vv2} |
| 1773 |
< |
\end{align} |
| 1774 |
< |
where: |
| 1775 |
< |
\begin{align} |
| 1776 |
< |
\mathbf{\mathcal{G}}_{Ri}(t) &= |
| 1777 |
< |
-2 \sum_j \lambda_{Rij}(t) \mathbf{r}_{ij}(t) ,\\ |
| 1778 |
< |
% |
| 1779 |
< |
\mathbf{\mathcal{G}}_{Vi}(t+h) &= |
| 1780 |
< |
-2 \sum_j \lambda_{Vij}(t+h) \mathbf{r}(t+h). |
| 1781 |
< |
\end{align} |
| 1782 |
< |
Next, define: |
| 1783 |
< |
\begin{align} |
| 1784 |
< |
g_{ij} &= h \lambda_{Rij}(t) ,\\ |
| 1785 |
< |
k_{ij} &= h \lambda_{Vij}(t+h), \\ |
| 1786 |
< |
\mathbf{q}_i &= \mathbf{\dot{r}}_i(t) + \frac{h}{2m_i} \mathbf{F}_i(t) |
| 1787 |
< |
- \frac{1}{m_i}\sum_j g_{ij}\mathbf{r}_{ij}(t). |
| 1788 |
< |
\end{align} |
| 1789 |
< |
Using these definitions, Eq.~\ref{oopseEq:vv1} and \ref{oopseEq:vv2} |
| 1790 |
< |
can be rewritten as, |
| 1791 |
< |
\begin{align} |
| 1792 |
< |
\mathbf{r}_i(t+h) &= \mathbf{r}_i(t) + h \mathbf{q}_i ,\\ |
| 1793 |
< |
% |
| 1794 |
< |
\mathbf{\dot{r}}(t+h) &= \mathbf{q}_i + \frac{h}{2m_i}\mathbf{F}_i(t+h) |
| 1795 |
< |
-\frac{1}{m_i}\sum_j k_{ij} \mathbf{r}_{ij}(t+h). |
| 1796 |
< |
\end{align} |
| 1797 |
< |
|
| 1798 |
< |
To integrate the equations of motion, the {\sc rattle} algorithm first |
| 1799 |
< |
solves for $\mathbf{r}(t+h)$. Let, |
| 1800 |
< |
\begin{equation} |
| 1801 |
< |
\mathbf{q}_i = \mathbf{\dot{r}}(t) + \frac{h}{2m_i}\mathbf{F}_i(t). |
| 1802 |
< |
\end{equation} |
| 1803 |
< |
Here $\mathbf{q}_i$ corresponds to an initial unconstrained move. Next |
| 1804 |
< |
pick a constraint $j$, and let, |
| 1805 |
< |
\begin{equation} |
| 1806 |
< |
\mathbf{s} = \mathbf{r}_i(t) + h\mathbf{q}_i(t) |
| 1807 |
< |
- \mathbf{r}_j(t) + h\mathbf{q}_j(t). |
| 1808 |
< |
\label{oopseEq:ra1} |
| 1809 |
< |
\end{equation} |
| 1810 |
< |
If |
| 1811 |
< |
\begin{equation} |
| 1812 |
< |
\Big| |\mathbf{s}|^2 - d_{ij}^2 \Big| > \text{tolerance}, |
| 1813 |
< |
\end{equation} |
| 1814 |
< |
then the constraint is unsatisfied, and corrections are made to the |
| 1815 |
< |
positions. First we define a test corrected configuration as, |
| 1816 |
< |
\begin{align} |
| 1817 |
< |
\mathbf{r}_i^T(t+h) = \mathbf{r}_i(t) + h\biggl[\mathbf{q}_i - |
| 1818 |
< |
g_{ij}\,\frac{\mathbf{r}_{ij}(t)}{m_i} \biggr] ,\\ |
| 1819 |
< |
% |
| 1820 |
< |
\mathbf{r}_j^T(t+h) = \mathbf{r}_j(t) + h\biggl[\mathbf{q}_j + |
| 1821 |
< |
g_{ij}\,\frac{\mathbf{r}_{ij}(t)}{m_j} \biggr]. |
| 1822 |
< |
\end{align} |
| 1823 |
< |
And we chose $g_{ij}$ such that, $|\mathbf{r}_i^T - \mathbf{r}_j^T|^2 |
| 1824 |
< |
= d_{ij}^2$. Solving the quadratic for $g_{ij}$ we obtain the |
| 1825 |
< |
approximation, |
| 1826 |
< |
\begin{equation} |
| 1827 |
< |
g_{ij} = \frac{(s^2 - d^2)}{2h[\mathbf{s}\cdot\mathbf{r}_{ij}(t)] |
| 1828 |
< |
(\frac{1}{m_i} + \frac{1}{m_j})}. |
| 1829 |
< |
\end{equation} |
| 1830 |
< |
Although not an exact solution for $g_{ij}$, as this is an iterative |
| 1831 |
< |
scheme overall, the eventual solution will converge. With a trial |
| 1832 |
< |
$g_{ij}$, the new $\mathbf{q}$'s become, |
| 1833 |
< |
\begin{align} |
| 1834 |
< |
\mathbf{q}_i &= \mathbf{q}^{\text{old}}_i - g_{ij}\, |
| 1835 |
< |
\frac{\mathbf{r}_{ij}(t)}{m_i} ,\\ |
| 1836 |
< |
% |
| 1837 |
< |
\mathbf{q}_j &= \mathbf{q}^{\text{old}}_j + g_{ij}\, |
| 1838 |
< |
\frac{\mathbf{r}_{ij}(t)}{m_j} . |
| 1839 |
< |
\end{align} |
| 1840 |
< |
The whole algorithm is then repeated from Eq.~\ref{oopseEq:ra1} until |
| 1841 |
< |
all constraints are satisfied. |
| 1842 |
< |
|
| 1843 |
< |
The second step of {\sc rattle}, is to then update the velocities. The |
| 1844 |
< |
step starts with, |
| 1845 |
< |
\begin{equation} |
| 1846 |
< |
\mathbf{\dot{r}}_i(t+h) = \mathbf{q}_i + \frac{h}{2m_i}\mathbf{F}_i(t+h). |
| 1847 |
< |
\end{equation} |
| 1848 |
< |
Next we pick a constraint $j$, and calculate the dot product $\ell$. |
| 1849 |
< |
\begin{equation} |
| 1850 |
< |
\ell = \mathbf{r}_{ij}(t+h) \cdot \mathbf{\dot{r}}_{ij}(t+h). |
| 1851 |
< |
\label{oopseEq:rv1} |
| 1852 |
< |
\end{equation} |
| 1853 |
< |
Here if constraint Eq.~\ref{oopseEq:c2} holds, $\ell$ should be |
| 1854 |
< |
zero. Therefore if $\ell$ is greater than some tolerance, then |
| 1855 |
< |
corrections are made to the $i$ and $j$ velocities. |
| 1856 |
< |
\begin{align} |
| 1857 |
< |
\mathbf{\dot{r}}_i^T &= \mathbf{\dot{r}}_i(t+h) - k_{ij} |
| 1858 |
< |
\frac{\mathbf{\dot{r}}_{ij}(t+h)}{m_i}, \\ |
| 1859 |
< |
% |
| 1860 |
< |
\mathbf{\dot{r}}_j^T &= \mathbf{\dot{r}}_j(t+h) + k_{ij} |
| 1861 |
< |
\frac{\mathbf{\dot{r}}_{ij}(t+h)}{m_j}. |
| 1862 |
< |
\end{align} |
| 1863 |
< |
Like in the previous step, we select a value for $k_{ij}$ such that |
| 1864 |
< |
$\ell$ is zero. |
| 1865 |
< |
\begin{equation} |
| 1866 |
< |
k_{ij} = \frac{\ell}{d^2_{ij}(\frac{1}{m_i} + \frac{1}{m_j})}. |
| 1867 |
< |
\end{equation} |
| 1868 |
< |
The test velocities, $\mathbf{\dot{r}}^T_i$ and |
| 1869 |
< |
$\mathbf{\dot{r}}^T_j$, then replace their respective velocities, and |
| 1870 |
< |
the algorithm is iterated from Eq.~\ref{oopseEq:rv1} until all |
| 1871 |
< |
constraints are satisfied. |
| 1872 |
< |
|
| 1873 |
< |
|
| 1874 |
< |
\subsection{\label{oopseSec:zcons}Z-Constraint Method} |
| 1875 |
< |
|
| 1876 |
< |
Based on the fluctuation-dissipation theorem, a force auto-correlation |
| 1877 |
< |
method was developed by Roux and Karplus to investigate the dynamics |
| 2045 |
> |
A force auto-correlation method based on the fluctuation-dissipation |
| 2046 |
> |
theorem was developed by Roux and Karplus to investigate the dynamics |
| 2047 |
|
of ions inside ion channels.\cite{Roux91} The time-dependent friction |
| 2048 |
|
coefficient can be calculated from the deviation of the instantaneous |
| 2049 |
< |
force from its mean force. |
| 2049 |
> |
force from its mean value: |
| 2050 |
|
\begin{equation} |
| 2051 |
|
\xi(z,t)=\langle\delta F(z,t)\delta F(z,0)\rangle/k_{B}T, |
| 2052 |
|
\end{equation} |
| 2055 |
|
\delta F(z,t)=F(z,t)-\langle F(z,t)\rangle. |
| 2056 |
|
\end{equation} |
| 2057 |
|
|
| 1889 |
– |
|
| 2058 |
|
If the time-dependent friction decays rapidly, the static friction |
| 2059 |
|
coefficient can be approximated by |
| 2060 |
|
\begin{equation} |
| 2061 |
|
\xi_{\text{static}}(z)=\int_{0}^{\infty}\langle\delta F(z,t)\delta F(z,0)\rangle dt. |
| 2062 |
|
\end{equation} |
| 2063 |
< |
Allowing diffusion constant to then be calculated through the |
| 2063 |
> |
|
| 2064 |
> |
This allows the diffusion constant to then be calculated through the |
| 2065 |
|
Einstein relation:\cite{Marrink94} |
| 2066 |
|
\begin{equation} |
| 2067 |
|
D(z)=\frac{k_{B}T}{\xi_{\text{static}}(z)}=\frac{(k_{B}T)^{2}}{\int_{0}^{\infty |
| 2068 |
|
}\langle\delta F(z,t)\delta F(z,0)\rangle dt}.% |
| 2069 |
|
\end{equation} |
| 2070 |
|
|
| 2071 |
< |
The Z-Constraint method, which fixes the z coordinates of the |
| 2072 |
< |
molecules with respect to the center of the mass of the system, has |
| 2073 |
< |
been a method suggested to obtain the forces required for the force |
| 2074 |
< |
auto-correlation calculation.\cite{Marrink94} However, simply resetting the |
| 2075 |
< |
coordinate will move the center of the mass of the whole system. To |
| 2076 |
< |
avoid this problem, a new method was used in {\sc oopse}. Instead of |
| 2077 |
< |
resetting the coordinate, we reset the forces of z-constrained |
| 2078 |
< |
molecules as well as subtract the total constraint forces from the |
| 2079 |
< |
rest of the system after the force calculation at each time step. |
| 2071 |
> |
The Z-Constraint method, which fixes the $z$ coordinates of a few |
| 2072 |
> |
``tagged'' molecules with respect to the center of the mass of the |
| 2073 |
> |
system is a technique that was proposed to obtain the forces required |
| 2074 |
> |
for the force auto-correlation calculation.\cite{Marrink94} However, |
| 2075 |
> |
simply resetting the coordinate will move the center of the mass of |
| 2076 |
> |
the whole system. To avoid this problem, we have developed a new |
| 2077 |
> |
method that is utilized in {\sc oopse}. Instead of resetting the |
| 2078 |
> |
coordinates, we reset the forces of $z$-constrained molecules and |
| 2079 |
> |
subtract the total constraint forces from the rest of the system after |
| 2080 |
> |
the force calculation at each time step. |
| 2081 |
|
|
| 2082 |
< |
After the force calculation, define $G_\alpha$ as |
| 2082 |
> |
After the force calculation, the total force on molecule $\alpha$ is: |
| 2083 |
|
\begin{equation} |
| 2084 |
|
G_{\alpha} = \sum_i F_{\alpha i}, |
| 2085 |
|
\label{oopseEq:zc1} |
| 2086 |
|
\end{equation} |
| 2087 |
< |
where $F_{\alpha i}$ is the force in the z direction of atom $i$ in |
| 2088 |
< |
z-constrained molecule $\alpha$. The forces of the z constrained |
| 2089 |
< |
molecule are then set to: |
| 2087 |
> |
where $F_{\alpha i}$ is the force in the $z$ direction on atom $i$ in |
| 2088 |
> |
$z$-constrained molecule $\alpha$. The forces on the atoms in the |
| 2089 |
> |
$z$-constrained molecule are then adjusted to remove the total force |
| 2090 |
> |
on molecule $\alpha$: |
| 2091 |
|
\begin{equation} |
| 2092 |
|
F_{\alpha i} = F_{\alpha i} - |
| 2093 |
|
\frac{m_{\alpha i} G_{\alpha}}{\sum_i m_{\alpha i}}. |
| 2094 |
|
\end{equation} |
| 2095 |
< |
Here, $m_{\alpha i}$ is the mass of atom $i$ in the z-constrained |
| 2096 |
< |
molecule. Having rescaled the forces, the velocities must also be |
| 2097 |
< |
rescaled to subtract out any center of mass velocity in the z |
| 2098 |
< |
direction. |
| 2095 |
> |
Here, $m_{\alpha i}$ is the mass of atom $i$ in the $z$-constrained |
| 2096 |
> |
molecule. After the forces have been adjusted, the velocities must |
| 2097 |
> |
also be modified to subtract out molecule $\alpha$'s center-of-mass |
| 2098 |
> |
velocity in the $z$ direction. |
| 2099 |
|
\begin{equation} |
| 2100 |
|
v_{\alpha i} = v_{\alpha i} - |
| 2101 |
|
\frac{\sum_i m_{\alpha i} v_{\alpha i}}{\sum_i m_{\alpha i}}, |
| 2102 |
|
\end{equation} |
| 2103 |
< |
where $v_{\alpha i}$ is the velocity of atom $i$ in the z direction. |
| 2104 |
< |
Lastly, all of the accumulated z constrained forces must be subtracted |
| 2105 |
< |
from the system to keep the system center of mass from drifting. |
| 2103 |
> |
where $v_{\alpha i}$ is the velocity of atom $i$ in the $z$ direction. |
| 2104 |
> |
Lastly, all of the accumulated constraint forces must be subtracted |
| 2105 |
> |
from the rest of the unconstrained system to keep the system center of |
| 2106 |
> |
mass of the entire system from drifting. |
| 2107 |
|
\begin{equation} |
| 2108 |
|
F_{\beta i} = F_{\beta i} - \frac{m_{\beta i} \sum_{\alpha} G_{\alpha}} |
| 2109 |
|
{\sum_{\beta}\sum_i m_{\beta i}}, |
| 2110 |
|
\end{equation} |
| 2111 |
< |
where $\beta$ are all of the unconstrained molecules in the |
| 2111 |
> |
where $\beta$ denotes all {\it unconstrained} molecules in the |
| 2112 |
|
system. Similarly, the velocities of the unconstrained molecules must |
| 2113 |
< |
also be scaled. |
| 2113 |
> |
also be scaled: |
| 2114 |
|
\begin{equation} |
| 2115 |
< |
v_{\beta i} = v_{\beta i} + \sum_{\alpha} |
| 2116 |
< |
\frac{\sum_i m_{\alpha i} v_{\alpha i}}{\sum_i m_{\alpha i}}. |
| 2115 |
> |
v_{\beta i} = v_{\beta i} + \sum_{\alpha} \frac{\sum_i m_{\alpha i} |
| 2116 |
> |
v_{\alpha i}}{\sum_i m_{\alpha i}}. |
| 2117 |
|
\end{equation} |
| 2118 |
|
|
| 2119 |
< |
At the very beginning of the simulation, the molecules may not be at their |
| 2120 |
< |
constrained positions. To move a z-constrained molecule to its specified |
| 2121 |
< |
position, a simple harmonic potential is used |
| 2119 |
> |
This method will pin down the centers-of-mass of all of the |
| 2120 |
> |
$z$-constrained molecules, and will also keep the entire system fixed |
| 2121 |
> |
at the original system center-of-mass location. |
| 2122 |
> |
|
| 2123 |
> |
At the very beginning of the simulation, the molecules may not be at |
| 2124 |
> |
their desired positions. To steer a $z$-constrained molecule to its |
| 2125 |
> |
specified position, a simple harmonic potential is used: |
| 2126 |
|
\begin{equation} |
| 2127 |
|
U(t)=\frac{1}{2}k_{\text{Harmonic}}(z(t)-z_{\text{cons}})^{2},% |
| 2128 |
|
\end{equation} |
| 2129 |
< |
where $k_{\text{Harmonic}}$ is the harmonic force constant, $z(t)$ is the |
| 2130 |
< |
current $z$ coordinate of the center of mass of the constrained molecule, and |
| 2131 |
< |
$z_{\text{cons}}$ is the constrained position. The harmonic force operating |
| 2132 |
< |
on the z-constrained molecule at time $t$ can be calculated by |
| 2129 |
> |
where $k_{\text{Harmonic}}$ is an harmonic force constant, $z(t)$ is |
| 2130 |
> |
the current $z$ coordinate of the center of mass of the constrained |
| 2131 |
> |
molecule, and $z_{\text{cons}}$ is the desired constraint |
| 2132 |
> |
position. The harmonic force operating on the $z$-constrained molecule |
| 2133 |
> |
at time $t$ can be calculated by |
| 2134 |
|
\begin{equation} |
| 2135 |
|
F_{z_{\text{Harmonic}}}(t)=-\frac{\partial U(t)}{\partial z(t)}= |
| 2136 |
|
-k_{\text{Harmonic}}(z(t)-z_{\text{cons}}). |
| 2137 |
|
\end{equation} |
| 2138 |
|
|
| 2139 |
< |
\section{\label{oopseSec:props}Trajectory Analysis} |
| 2139 |
> |
The user may also specify the use of a constant velocity method |
| 2140 |
> |
(steered molecular dynamics) to move the molecules to their desired |
| 2141 |
> |
initial positions. Based on concepts from atomic force microscopy, |
| 2142 |
> |
{\sc smd} has been used to study many processes which occur via rare |
| 2143 |
> |
events on the time scale of a few hundreds of picoseconds. For |
| 2144 |
> |
example,{\sc smd} has been used to observe the dissociation of |
| 2145 |
> |
Streptavidin-biotin Complex.\cite{smd} |
| 2146 |
|
|
| 2147 |
< |
\subsection{\label{oopseSec:staticProps}Static Property Analysis} |
| 2147 |
> |
To use of the $z$-constraint method in an {\sc oopse} simulation, the |
| 2148 |
> |
molecules must be specified using the {\tt nZconstraints} keyword in |
| 2149 |
> |
the meta-data file. The other parameters for modifying the behavior |
| 2150 |
> |
of the $z$-constraint method are listed in table~\ref{table:zconParams}. |
| 2151 |
|
|
| 1966 |
– |
The static properties of the trajectories are analyzed with the |
| 1967 |
– |
program \texttt{staticProps}. The code is capable of calculating a |
| 1968 |
– |
number of pair correlations between species A and B. Some of which |
| 1969 |
– |
only apply to directional entities. The summary of pair correlations |
| 1970 |
– |
can be found in Table~\ref{oopseTb:gofrs} |
| 1971 |
– |
|
| 2152 |
|
\begin{table} |
| 2153 |
< |
\caption{THE DIFFERENT PAIR CORRELATIONS IN \texttt{staticProps}} |
| 2154 |
< |
\label{oopseTb:gofrs} |
| 2153 |
> |
\caption{Meta-data Keywords: Z-Constraint Parameters} |
| 2154 |
> |
\label{table:zconParams} |
| 2155 |
|
\begin{center} |
| 2156 |
< |
\begin{tabular}{|l|c|c|} |
| 2157 |
< |
\hline |
| 2158 |
< |
Name & Equation & Directional Atom \\ \hline |
| 2159 |
< |
$g_{\text{AB}}(r)$ & Eq.~\ref{eq:gofr} & neither \\ \hline |
| 2160 |
< |
$g_{\text{AB}}(r, \cos \theta)$ & Eq.~\ref{eq:gofrCosTheta} & A \\ \hline |
| 2161 |
< |
$g_{\text{AB}}(r, \cos \omega)$ & Eq.~\ref{eq:gofrCosOmega} & both \\ \hline |
| 2162 |
< |
$g_{\text{AB}}(x, y, z)$ & Eq.~\ref{eq:gofrXYZ} & neither \\ \hline |
| 2163 |
< |
$\langle \cos \omega \rangle_{\text{AB}}(r)$ & Eq.~\ref{eq:cosOmegaOfR} &% |
| 2164 |
< |
both \\ \hline |
| 2165 |
< |
\end{tabular} |
| 2166 |
< |
\begin{minipage}{\linewidth} |
| 2167 |
< |
\centering |
| 2168 |
< |
\vspace{2mm} |
| 2169 |
< |
The third column specifies which atom, if any, need be a directional entity. |
| 2170 |
< |
\end{minipage} |
| 2156 |
> |
% Note when adding or removing columns, the \hsize numbers must add up to the total number |
| 2157 |
> |
% of columns. |
| 2158 |
> |
\begin{tabularx}{\linewidth}% |
| 2159 |
> |
{>{\setlength{\hsize}{1.00\hsize}}X% |
| 2160 |
> |
>{\setlength{\hsize}{0.4\hsize}}X% |
| 2161 |
> |
>{\setlength{\hsize}{1.2\hsize}}X% |
| 2162 |
> |
>{\setlength{\hsize}{1.4\hsize}}X} |
| 2163 |
> |
|
| 2164 |
> |
{\bf keyword} & {\bf units} & {\bf use} & {\bf remarks} \\ \hline |
| 2165 |
> |
|
| 2166 |
> |
{\tt nZconstraints} & integer & The number of $z$-constrained |
| 2167 |
> |
molecules & If using the $z$-constraint method, {\tt nZconstraints} |
| 2168 |
> |
must be set \\ |
| 2169 |
> |
{\tt zconsTime} & fs & Sets the frequency at which the {\tt .fz} file |
| 2170 |
> |
is written & \\ |
| 2171 |
> |
{\tt zconsForcePolicy} & string & The strategy for subtracting |
| 2172 |
> |
the $z$-constraint force from the {\it unconstrained} molecules & Possible |
| 2173 |
> |
strategies are {\tt BYMASS} and {\tt BYNUMBER}. The default |
| 2174 |
> |
strategy is {\tt BYMASS}\\ |
| 2175 |
> |
{\tt zconsGap} & $\mbox{\AA}$ & Sets the distance between two adjacent |
| 2176 |
> |
constraint positions&Used mainly to move molecules through a |
| 2177 |
> |
simulation to estimate potentials of mean force. \\ |
| 2178 |
> |
{\tt zconsFixtime} & fs & Sets the length of time the $z$-constraint |
| 2179 |
> |
molecule is held fixed & {\tt zconsFixtime} must be set if {\tt |
| 2180 |
> |
zconsGap} is set\\ |
| 2181 |
> |
{\tt zconsUsingSMD} & logical & Flag for using Steered Molecular |
| 2182 |
> |
Dynamics to move the molecules to the correct constrained positions & |
| 2183 |
> |
Harmonic Forces are used by default\\ |
| 2184 |
> |
|
| 2185 |
> |
\end{tabularx} |
| 2186 |
|
\end{center} |
| 2187 |
|
\end{table} |
| 2188 |
|
|
| 2189 |
< |
The first pair correlation, $g_{\text{AB}}(r)$, is defined as follows: |
| 2189 |
> |
|
| 2190 |
> |
\section{\label{oopseSec:minimizer}Energy Minimization} |
| 2191 |
> |
|
| 2192 |
> |
As one of the basic procedures of molecular modeling, energy |
| 2193 |
> |
minimization is used to identify local configurations that are stable |
| 2194 |
> |
points on the potential energy surface. There is a vast literature on |
| 2195 |
> |
energy minimization algorithms have been developed to search for the |
| 2196 |
> |
global energy minimum as well as to find local structures which are |
| 2197 |
> |
stable fixed points on the surface. We have included two simple |
| 2198 |
> |
minimization algorithms: steepest descent, ({\sc sd}) and conjugate |
| 2199 |
> |
gradient ({\sc cg}) to help users find reasonable local minima from |
| 2200 |
> |
their initial configurations. Since {\sc oopse} handles atoms and |
| 2201 |
> |
rigid bodies which have orientational coordinates as well as |
| 2202 |
> |
translational coordinates, there is some subtlety to the choice of |
| 2203 |
> |
parameters for minimization algorithms. |
| 2204 |
> |
|
| 2205 |
> |
Given a coordinate set $x_{k}$ and a search direction $d_{k}$, a line |
| 2206 |
> |
search algorithm is performed along $d_{k}$ to produce |
| 2207 |
> |
$x_{k+1}=x_{k}+$ $\lambda _{k}d_{k}$. In the steepest descent ({\sc |
| 2208 |
> |
sd}) algorithm,% |
| 2209 |
|
\begin{equation} |
| 2210 |
< |
g_{\text{AB}}(r) = \frac{V}{N_{\text{A}}N_{\text{B}}}\langle %% |
| 1997 |
< |
\sum_{i \in \text{A}} \sum_{j \in \text{B}} %% |
| 1998 |
< |
\delta( r - |\mathbf{r}_{ij}|) \rangle, \label{eq:gofr} |
| 2210 |
> |
d_{k}=-\nabla V(x_{k}). |
| 2211 |
|
\end{equation} |
| 2212 |
< |
where $\mathbf{r}_{ij}$ is the vector |
| 2213 |
< |
\begin{equation*} |
| 2214 |
< |
\mathbf{r}_{ij} = \mathbf{r}_j - \mathbf{r}_i, \notag |
| 2215 |
< |
\end{equation*} |
| 2216 |
< |
and $\frac{V}{N_{\text{A}}N_{\text{B}}}$ normalizes the average over |
| 2005 |
< |
the expected pair density at a given $r$. |
| 2006 |
< |
|
| 2007 |
< |
The next two pair correlations, $g_{\text{AB}}(r, \cos \theta)$ and |
| 2008 |
< |
$g_{\text{AB}}(r, \cos \omega)$, are similar in that they are both two |
| 2009 |
< |
dimensional histograms. Both use $r$ for the primary axis then a |
| 2010 |
< |
$\cos$ for the secondary axis ($\cos \theta$ for |
| 2011 |
< |
Eq.~\ref{eq:gofrCosTheta} and $\cos \omega$ for |
| 2012 |
< |
Eq.~\ref{eq:gofrCosOmega}). This allows for the investigator to |
| 2013 |
< |
correlate alignment on directional entities. $g_{\text{AB}}(r, \cos |
| 2014 |
< |
\theta)$ is defined as follows: |
| 2212 |
> |
The gradient and the direction of next step are always orthogonal. |
| 2213 |
> |
This may cause oscillatory behavior in narrow valleys. To overcome |
| 2214 |
> |
this problem, the Fletcher-Reeves variant~\cite{FletcherReeves} of the |
| 2215 |
> |
conjugate gradient ({\sc cg}) algorithm is used to generate $d_{k+1}$ |
| 2216 |
> |
via simple recursion: |
| 2217 |
|
\begin{equation} |
| 2218 |
< |
g_{\text{AB}}(r, \cos \theta) = \frac{V}{N_{\text{A}}N_{\text{B}}}\langle |
| 2017 |
< |
\sum_{i \in \text{A}} \sum_{j \in \text{B}} |
| 2018 |
< |
\delta( \cos \theta - \cos \theta_{ij}) |
| 2019 |
< |
\delta( r - |\mathbf{r}_{ij}|) \rangle. |
| 2020 |
< |
\label{eq:gofrCosTheta} |
| 2218 |
> |
d_{k+1} =-\nabla V(x_{k+1})+\gamma_{k}d_{k} |
| 2219 |
|
\end{equation} |
| 2220 |
< |
Here |
| 2023 |
< |
\begin{equation*} |
| 2024 |
< |
\cos \theta_{ij} = \mathbf{\hat{i}} \cdot \mathbf{\hat{r}}_{ij}, |
| 2025 |
< |
\end{equation*} |
| 2026 |
< |
where $\mathbf{\hat{i}}$ is the unit directional vector of species $i$ |
| 2027 |
< |
and $\mathbf{\hat{r}}_{ij}$ is the unit vector associated with vector |
| 2028 |
< |
$\mathbf{r}_{ij}$. |
| 2029 |
< |
|
| 2030 |
< |
The second two dimensional histogram is of the form: |
| 2220 |
> |
where |
| 2221 |
|
\begin{equation} |
| 2222 |
< |
g_{\text{AB}}(r, \cos \omega) = |
| 2223 |
< |
\frac{V}{N_{\text{A}}N_{\text{B}}}\langle |
| 2034 |
< |
\sum_{i \in \text{A}} \sum_{j \in \text{B}} |
| 2035 |
< |
\delta( \cos \omega - \cos \omega_{ij}) |
| 2036 |
< |
\delta( r - |\mathbf{r}_{ij}|) \rangle. \label{eq:gofrCosOmega} |
| 2222 |
> |
\gamma_{k} =\frac{\nabla V(x_{k+1})^{T}\nabla V(x_{k+1})}{\nabla |
| 2223 |
> |
V(x_{k})^{T}\nabla V(x_{k})}. |
| 2224 |
|
\end{equation} |
| 2038 |
– |
Here |
| 2039 |
– |
\begin{equation*} |
| 2040 |
– |
\cos \omega_{ij} = \mathbf{\hat{i}} \cdot \mathbf{\hat{j}}. |
| 2041 |
– |
\end{equation*} |
| 2042 |
– |
Again, $\mathbf{\hat{i}}$ and $\mathbf{\hat{j}}$ are the unit |
| 2043 |
– |
directional vectors of species $i$ and $j$. |
| 2225 |
|
|
| 2226 |
< |
The static analysis code is also cable of calculating a three |
| 2227 |
< |
dimensional pair correlation of the form: |
| 2047 |
< |
\begin{equation}\label{eq:gofrXYZ} |
| 2048 |
< |
g_{\text{AB}}(x, y, z) = |
| 2049 |
< |
\frac{V}{N_{\text{A}}N_{\text{B}}}\langle |
| 2050 |
< |
\sum_{i \in \text{A}} \sum_{j \in \text{B}} |
| 2051 |
< |
\delta( x - x_{ij}) |
| 2052 |
< |
\delta( y - y_{ij}) |
| 2053 |
< |
\delta( z - z_{ij}) \rangle, |
| 2054 |
< |
\end{equation} |
| 2055 |
< |
where $x_{ij}$, $y_{ij}$, and $z_{ij}$ are the $x$, $y$, and $z$ |
| 2056 |
< |
components respectively of vector $\mathbf{r}_{ij}$. |
| 2057 |
< |
|
| 2058 |
< |
The final pair correlation is similar to |
| 2059 |
< |
Eq.~\ref{eq:gofrCosOmega}. $\langle \cos \omega |
| 2060 |
< |
\rangle_{\text{AB}}(r)$ is calculated in the following way: |
| 2061 |
< |
\begin{equation}\label{eq:cosOmegaOfR} |
| 2062 |
< |
\langle \cos \omega \rangle_{\text{AB}}(r) = |
| 2063 |
< |
\langle \sum_{i \in \text{A}} \sum_{j \in \text{B}} |
| 2064 |
< |
(\cos \omega_{ij}) \delta( r - |\mathbf{r}_{ij}|) \rangle. |
| 2065 |
< |
\end{equation} |
| 2066 |
< |
Here $\cos \omega_{ij}$ is defined in the same way as in |
| 2067 |
< |
Eq.~\ref{eq:gofrCosOmega}. This equation is a single dimensional pair |
| 2068 |
< |
correlation that gives the average correlation of two directional |
| 2069 |
< |
entities as a function of their distance from each other. |
| 2070 |
< |
|
| 2071 |
< |
\subsection{\label{dynamicProps}Dynamic Property Analysis} |
| 2072 |
< |
|
| 2073 |
< |
The dynamic properties of a trajectory are calculated with the program |
| 2074 |
< |
\texttt{dynamicProps}. The program calculates the following properties: |
| 2075 |
< |
\begin{gather} |
| 2076 |
< |
\langle | \mathbf{r}(t) - \mathbf{r}(0) |^2 \rangle, \label{eq:rms}\\ |
| 2077 |
< |
\langle \mathbf{v}(t) \cdot \mathbf{v}(0) \rangle, \label{eq:velCorr} \\ |
| 2078 |
< |
\langle \mathbf{j}(t) \cdot \mathbf{j}(0) \rangle. \label{eq:angularVelCorr} |
| 2079 |
< |
\end{gather} |
| 2080 |
< |
|
| 2081 |
< |
Eq.~\ref{eq:rms} is the root mean square displacement function. Which |
| 2082 |
< |
allows one to observe the average displacement of an atom as a |
| 2083 |
< |
function of time. The quantity is useful when calculating diffusion |
| 2084 |
< |
coefficients because of the Einstein Relation, which is valid at long |
| 2085 |
< |
times.\cite{allen87:csl} |
| 2226 |
> |
The Polak-Ribiere variant~\cite{PolakRibiere} of the conjugate |
| 2227 |
> |
gradient ($\gamma_{k}$) is defined as% |
| 2228 |
|
\begin{equation} |
| 2229 |
< |
2tD = \langle | \mathbf{r}(t) - \mathbf{r}(0) |^2 \rangle. |
| 2230 |
< |
\label{oopseEq:einstein} |
| 2229 |
> |
\gamma_{k}=\frac{[\nabla V(x_{k+1})-\nabla V(x)]^{T}\nabla V(x_{k+1})}{\nabla |
| 2230 |
> |
V(x_{k})^{T}\nabla V(x_{k})}% |
| 2231 |
|
\end{equation} |
| 2232 |
+ |
It is widely agreed that the Polak-Ribiere variant gives better |
| 2233 |
+ |
convergence than the Fletcher-Reeves variant, so the conjugate |
| 2234 |
+ |
gradient approach implemented in {\sc oopse} is the Polak-Ribiere |
| 2235 |
+ |
variant. |
| 2236 |
|
|
| 2237 |
< |
Eq.~\ref{eq:velCorr} and \ref{eq:angularVelCorr} are the translational |
| 2238 |
< |
velocity and angular velocity correlation functions respectively. The |
| 2239 |
< |
latter is only applicable to directional species in the |
| 2240 |
< |
simulation. The velocity autocorrelation functions are useful when |
| 2241 |
< |
determining vibrational information about the system of interest. |
| 2237 |
> |
The conjugate gradient method assumes that the conformation is close |
| 2238 |
> |
enough to a local minimum that the potential energy surface is very |
| 2239 |
> |
nearly quadratic. When the initial structure is far from the minimum, |
| 2240 |
> |
the steepest descent method can be superior to the conjugate gradient |
| 2241 |
> |
method. Hence, the steepest descent method is often used for the first |
| 2242 |
> |
10-100 steps of minimization. Another useful feature of minimization |
| 2243 |
> |
methods in {\sc oopse} is that a modified {\sc shake} algorithm can be |
| 2244 |
> |
applied during the minimization to constraint the bond lengths if this |
| 2245 |
> |
is required by the force field. Meta-data parameters concerning the |
| 2246 |
> |
minimizer are given in Table~\ref{table:minimizeParams} |
| 2247 |
|
|
| 2248 |
< |
\section{\label{oopseSec:design}Program Design} |
| 2248 |
> |
\begin{table} |
| 2249 |
> |
\caption{Meta-data Keywords: Energy Minimizer Parameters} |
| 2250 |
> |
\label{table:minimizeParams} |
| 2251 |
> |
\begin{center} |
| 2252 |
> |
% Note when adding or removing columns, the \hsize numbers must add up to the total number |
| 2253 |
> |
% of columns. |
| 2254 |
> |
\begin{tabularx}{\linewidth}% |
| 2255 |
> |
{>{\setlength{\hsize}{1.2\hsize}}X% |
| 2256 |
> |
>{\setlength{\hsize}{0.6\hsize}}X% |
| 2257 |
> |
>{\setlength{\hsize}{1.1\hsize}}X% |
| 2258 |
> |
>{\setlength{\hsize}{1.1\hsize}}X} |
| 2259 |
|
|
| 2260 |
< |
\subsection{\label{sec:architecture} {\sc oopse} Architecture} |
| 2260 |
> |
{\bf keyword} & {\bf units} & {\bf use} & {\bf remarks} \\ \hline |
| 2261 |
|
|
| 2262 |
< |
The core of OOPSE is divided into two main object libraries: |
| 2263 |
< |
\texttt{libBASS} and \texttt{libmdtools}. \texttt{libBASS} is the |
| 2264 |
< |
library developed around the parsing engine and \texttt{libmdtools} |
| 2265 |
< |
is the software library developed around the simulation engine. These |
| 2266 |
< |
two libraries are designed to encompass all the basic functions and |
| 2267 |
< |
tools that {\sc oopse} provides. Utility programs, such as the |
| 2268 |
< |
property analyzers, need only link against the software libraries to |
| 2269 |
< |
gain access to parsing, force evaluation, and input / output |
| 2270 |
< |
routines. |
| 2262 |
> |
{\tt minimizer} & string & selects the minimization method to be used |
| 2263 |
> |
& either {\tt CG} (conjugate gradient) or {\tt SD} (steepest |
| 2264 |
> |
descent) \\ |
| 2265 |
> |
{\tt minimizerMaxIter} & steps & Sets the maximum number of iterations |
| 2266 |
> |
for the energy minimization & The default value is 200\\ |
| 2267 |
> |
{\tt minimizerWriteFreq} & steps & Sets the frequency with which the {\tt .dump} and {\tt .stat} files are writtern during energy minimization & \\ |
| 2268 |
> |
{\tt minimizerStepSize} & $\mbox{\AA}$ & Sets the step size for the |
| 2269 |
> |
line search & The default value is 0.01\\ |
| 2270 |
> |
{\tt minimizerFTol} & $\mbox{kcal mol}^{-1}$ & Sets the energy tolerance |
| 2271 |
> |
for stopping the minimziation. & The default value is $10^{-8}$\\ |
| 2272 |
> |
{\tt minimizerGTol} & $\mbox{kcal mol}^{-1}\mbox{\AA}^{-1}$ & Sets the |
| 2273 |
> |
gradient tolerance for stopping the minimization. & The default value |
| 2274 |
> |
is $10^{-8}$\\ |
| 2275 |
> |
{\tt minimizerLSTol} & $\mbox{kcal mol}^{-1}$ & Sets line search |
| 2276 |
> |
tolerance for terminating each step of the minimization. & The default |
| 2277 |
> |
value is $10^{-8}$\\ |
| 2278 |
> |
{\tt minimizerLSMaxIter} & steps & Sets the maximum number of |
| 2279 |
> |
iterations for each line search & The default value is 50\\ |
| 2280 |
|
|
| 2281 |
< |
Contained in \texttt{libBASS} are all the routines associated with |
| 2282 |
< |
reading and parsing the \texttt{.bass} input files. Given a |
| 2283 |
< |
\texttt{.bass} file, \texttt{libBASS} will open it and any associated |
| 2114 |
< |
\texttt{.mdl} files; then create structures in memory that are |
| 2115 |
< |
templates of all the molecules specified in the input files. In |
| 2116 |
< |
addition, any simulation parameters set in the \texttt{.bass} file |
| 2117 |
< |
will be placed in a structure for later query by the controlling |
| 2118 |
< |
program. |
| 2281 |
> |
\end{tabularx} |
| 2282 |
> |
\end{center} |
| 2283 |
> |
\end{table} |
| 2284 |
|
|
| 2285 |
< |
Located in \texttt{libmdtools} are all other routines necessary to a |
| 2121 |
< |
Molecular Dynamics simulation. The library uses the main data |
| 2122 |
< |
structures returned by \texttt{libBASS} to initialize the various |
| 2123 |
< |
parts of the simulation: the atom structures and positions, the force |
| 2124 |
< |
field, the integrator, \emph{et cetera}. After initialization, the |
| 2125 |
< |
library can be used to perform a variety of tasks: integrate a |
| 2126 |
< |
Molecular Dynamics trajectory, query phase space information from a |
| 2127 |
< |
specific frame of a completed trajectory, or even recalculate force or |
| 2128 |
< |
energetic information about specific frames from a completed |
| 2129 |
< |
trajectory. |
| 2285 |
> |
\section{\label{oopseSec:parallelization} Parallel Simulation Implementation} |
| 2286 |
|
|
| 2287 |
< |
With these core libraries in place, several programs have been |
| 2288 |
< |
developed to utilize the routines provided by \texttt{libBASS} and |
| 2289 |
< |
\texttt{libmdtools}. The main program of the package is \texttt{oopse} |
| 2290 |
< |
and the corresponding parallel version \texttt{oopse\_MPI}. These two |
| 2291 |
< |
programs will take the \texttt{.bass} file, and create (and integrate) |
| 2292 |
< |
the simulation specified in the script. The two analysis programs |
| 2293 |
< |
\texttt{staticProps} and \texttt{dynamicProps} utilize the core |
| 2294 |
< |
libraries to initialize and read in trajectories from previously |
| 2139 |
< |
completed simulations, in addition to the ability to use functionality |
| 2140 |
< |
from \texttt{libmdtools} to recalculate forces and energies at key |
| 2141 |
< |
frames in the trajectories. Lastly, the family of system building |
| 2142 |
< |
programs (Sec.~\ref{oopseSec:initCoords}) also use the libraries to |
| 2143 |
< |
store and output the system configurations they create. |
| 2287 |
> |
Although processor power is continually improving, it is still |
| 2288 |
> |
unreasonable to simulate systems of more than 10,000 atoms on a single |
| 2289 |
> |
processor. To facilitate study of larger system sizes or smaller |
| 2290 |
> |
systems for longer time scales, parallel methods were developed to |
| 2291 |
> |
allow multiple CPU's to share the simulation workload. Three general |
| 2292 |
> |
categories of parallel decomposition methods have been developed: |
| 2293 |
> |
these are the atomic,\cite{Fox88} spatial~\cite{plimpton95} and |
| 2294 |
> |
force~\cite{Paradyn} decomposition methods. |
| 2295 |
|
|
| 2296 |
< |
\subsection{\label{oopseSec:parallelization} Parallelization of {\sc oopse}} |
| 2297 |
< |
|
| 2298 |
< |
Although processor power is continually growing roughly following |
| 2299 |
< |
Moore's Law, it is still unreasonable to simulate systems of more then |
| 2149 |
< |
a 1000 atoms on a single processor. To facilitate study of larger |
| 2150 |
< |
system sizes or smaller systems on long time scales in a reasonable |
| 2151 |
< |
period of time, parallel methods were developed allowing multiple |
| 2152 |
< |
CPU's to share the simulation workload. Three general categories of |
| 2153 |
< |
parallel decomposition methods have been developed including atomic, |
| 2154 |
< |
spatial and force decomposition methods. |
| 2155 |
< |
|
| 2156 |
< |
Algorithmically simplest of the three methods is atomic decomposition |
| 2157 |
< |
where N particles in a simulation are split among P processors for the |
| 2158 |
< |
duration of the simulation. Computational cost scales as an optimal |
| 2159 |
< |
$\mathcal{O}(N/P)$ for atomic decomposition. Unfortunately all |
| 2296 |
> |
Algorithmically simplest of the three methods is atomic decomposition, |
| 2297 |
> |
where $N$ particles in a simulation are split among $P$ processors for |
| 2298 |
> |
the duration of the simulation. Computational cost scales as an |
| 2299 |
> |
optimal $\mathcal{O}(N/P)$ for atomic decomposition. Unfortunately, all |
| 2300 |
|
processors must communicate positions and forces with all other |
| 2301 |
< |
processors at every force evaluation, leading communication costs to |
| 2302 |
< |
scale as an unfavorable $\mathcal{O}(N)$, \emph{independent of the |
| 2301 |
> |
processors at every force evaluation, leading the communication costs |
| 2302 |
> |
to scale as an unfavorable $\mathcal{O}(N)$, \emph{independent of the |
| 2303 |
|
number of processors}. This communication bottleneck led to the |
| 2304 |
< |
development of spatial and force decomposition methods in which |
| 2304 |
> |
development of spatial and force decomposition methods, in which |
| 2305 |
|
communication among processors scales much more favorably. Spatial or |
| 2306 |
|
domain decomposition divides the physical spatial domain into 3D boxes |
| 2307 |
|
in which each processor is responsible for calculation of forces and |
| 2308 |
|
positions of particles located in its box. Particles are reassigned to |
| 2309 |
|
different processors as they move through simulation space. To |
| 2310 |
< |
calculate forces on a given particle, a processor must know the |
| 2310 |
> |
calculate forces on a given particle, a processor must simply know the |
| 2311 |
|
positions of particles within some cutoff radius located on nearby |
| 2312 |
< |
processors instead of the positions of particles on all |
| 2312 |
> |
processors rather than the positions of particles on all |
| 2313 |
|
processors. Both communication between processors and computation |
| 2314 |
|
scale as $\mathcal{O}(N/P)$ in the spatial method. However, spatial |
| 2315 |
|
decomposition adds algorithmic complexity to the simulation code and |
| 2316 |
< |
is not very efficient for small N since the overall communication |
| 2316 |
> |
is not very efficient for small $N$, since the overall communication |
| 2317 |
|
scales as the surface to volume ratio $\mathcal{O}(N/P)^{2/3}$ in |
| 2318 |
|
three dimensions. |
| 2319 |
|
|
| 2320 |
|
The parallelization method used in {\sc oopse} is the force |
| 2321 |
< |
decomposition method. Force decomposition assigns particles to |
| 2322 |
< |
processors based on a block decomposition of the force |
| 2321 |
> |
decomposition method.\cite{hendrickson:95} Force decomposition assigns |
| 2322 |
> |
particles to processors based on a block decomposition of the force |
| 2323 |
|
matrix. Processors are split into an optimally square grid forming row |
| 2324 |
|
and column processor groups. Forces are calculated on particles in a |
| 2325 |
< |
given row by particles located in that processors column |
| 2326 |
< |
assignment. Force decomposition is less complex to implement than the |
| 2327 |
< |
spatial method but still scales computationally as $\mathcal{O}(N/P)$ |
| 2328 |
< |
and scales as $\mathcal{O}(N/\sqrt{P})$ in communication |
| 2329 |
< |
cost. Plimpton has also found that force decompositions scale more |
| 2330 |
< |
favorably than spatial decompositions for systems up to 10,000 atoms |
| 2331 |
< |
and favorably compete with spatial methods up to 100,000 |
| 2332 |
< |
atoms.\cite{plimpton95} |
| 2333 |
< |
|
| 2334 |
< |
\subsection{\label{oopseSec:memAlloc}Memory Issues in Trajectory Analysis} |
| 2335 |
< |
|
| 2336 |
< |
For large simulations, the trajectory files can sometimes reach sizes |
| 2197 |
< |
in excess of several gigabytes. In order to effectively analyze that |
| 2198 |
< |
amount of data, two memory management schemes have been devised for |
| 2199 |
< |
\texttt{staticProps} and for \texttt{dynamicProps}. The first scheme, |
| 2200 |
< |
developed for \texttt{staticProps}, is the simplest. As each frame's |
| 2201 |
< |
statistics are calculated independent of each other, memory is |
| 2202 |
< |
allocated for each frame, then freed once correlation calculations are |
| 2203 |
< |
complete for the snapshot. To prevent multiple passes through a |
| 2204 |
< |
potentially large file, \texttt{staticProps} is capable of calculating |
| 2205 |
< |
all requested correlations per frame with only a single pair loop in |
| 2206 |
< |
each frame and a single read of the file. |
| 2207 |
< |
|
| 2208 |
< |
The second, more advanced memory scheme, is used by |
| 2209 |
< |
\texttt{dynamicProps}. Here, the program must have multiple frames in |
| 2210 |
< |
memory to calculate time dependent correlations. In order to prevent a |
| 2211 |
< |
situation where the program runs out of memory due to large |
| 2212 |
< |
trajectories, the user is able to specify that the trajectory be read |
| 2213 |
< |
in blocks. The number of frames in each block is specified by the |
| 2214 |
< |
user, and upon reading a block of the trajectory, |
| 2215 |
< |
\texttt{dynamicProps} will calculate all of the time correlation frame |
| 2216 |
< |
pairs within the block. After in-block correlations are complete, a |
| 2217 |
< |
second block of the trajectory is read, and the cross correlations are |
| 2218 |
< |
calculated between the two blocks. This second block is then freed and |
| 2219 |
< |
then incremented and the process repeated until the end of the |
| 2220 |
< |
trajectory. Once the end is reached, the first block is freed then |
| 2221 |
< |
incremented, and the again the internal time correlations are |
| 2222 |
< |
calculated. The algorithm with the second block is then repeated with |
| 2223 |
< |
the new origin block, until all frame pairs have been correlated in |
| 2224 |
< |
time. This process is illustrated in |
| 2225 |
< |
Fig.~\ref{oopseFig:dynamicPropsMemory}. |
| 2226 |
< |
|
| 2227 |
< |
%\begin{figure} |
| 2228 |
< |
%\centering |
| 2229 |
< |
%\includegraphics[width=\linewidth]{dynamicPropsMem.eps} |
| 2230 |
< |
%\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.} |
| 2231 |
< |
%\label{oopseFig:dynamicPropsMemory} |
| 2232 |
< |
%\end{figure} |
| 2233 |
< |
|
| 2325 |
> |
given row by particles located in that processor's column |
| 2326 |
> |
assignment. One deviation from the algorithm described by Hendrickson |
| 2327 |
> |
{\it et al.} is the use of column ordering based on the row indexes |
| 2328 |
> |
preventing the need for a transpose operation necessitating a second |
| 2329 |
> |
communication step when gathering the final force components. Force |
| 2330 |
> |
decomposition is less complex to implement than the spatial method but |
| 2331 |
> |
still scales computationally as $\mathcal{O}(N/P)$ and scales as |
| 2332 |
> |
$\mathcal{O}(N/\sqrt{P})$ in communication cost. Plimpton has also |
| 2333 |
> |
found that force decompositions scale more favorably than spatial |
| 2334 |
> |
decompositions for systems up to 10,000 atoms and favorably compete |
| 2335 |
> |
with spatial methods up to 100,000 atoms.\cite{plimpton95} |
| 2336 |
> |
|
| 2337 |
|
\section{\label{oopseSec:conclusion}Conclusion} |
| 2338 |
|
|
| 2339 |
< |
We have presented the design and implementation of our open source |
| 2340 |
< |
simulation package {\sc oopse}. The package offers novel capabilities |
| 2341 |
< |
to the field of Molecular Dynamics simulation packages in the form of |
| 2342 |
< |
dipolar force fields, and symplectic integration of rigid body |
| 2343 |
< |
dynamics. It is capable of scaling across multiple processors through |
| 2344 |
< |
the use of force based decomposition using MPI. It also implements |
| 2345 |
< |
several advanced integrators allowing the end user control over |
| 2346 |
< |
temperature and pressure. In addition, it is capable of integrating |
| 2347 |
< |
constrained dynamics through both the {\sc rattle} algorithm and the |
| 2348 |
< |
z-constraint method. |
| 2339 |
> |
We have presented a new parallel simulation program called {\sc |
| 2340 |
> |
oopse}. This program offers some novel capabilities, but mostly makes |
| 2341 |
> |
available a library of modern object-oriented code for the scientific |
| 2342 |
> |
community to use freely. Notably, {\sc oopse} can handle symplectic |
| 2343 |
> |
integration of objects (atoms and rigid bodies) which have |
| 2344 |
> |
orientational degrees of freedom. It can also work with transition |
| 2345 |
> |
metal force fields and point-dipoles. It is capable of scaling across |
| 2346 |
> |
multiple processors through the use of force based decomposition. It |
| 2347 |
> |
also implements several advanced integrators allowing the end user |
| 2348 |
> |
control over temperature and pressure. In addition, it is capable of |
| 2349 |
> |
integrating constrained dynamics through both the {\sc rattle} |
| 2350 |
> |
algorithm and the $z$-constraint method. |
| 2351 |
|
|
| 2352 |
< |
These features are all brought together in a single open-source |
| 2353 |
< |
program. This allows researchers to not only benefit from |
| 2354 |
< |
{\sc oopse}, but also contribute to {\sc oopse}'s development as |
| 2355 |
< |
well. |
| 2352 |
> |
We encourage other researchers to download and apply this program to |
| 2353 |
> |
their own research problems. By making the code available, we hope to |
| 2354 |
> |
encourage other researchers to contribute their own code and make it a |
| 2355 |
> |
more powerful package for everyone in the molecular dynamics community |
| 2356 |
> |
to use. All source code for {\sc oopse} is available for download at |
| 2357 |
> |
{\tt http://oopse.org}. |
| 2358 |
|
|
| 2252 |
– |
|
| 2359 |
|
\newpage |
| 2360 |
|
\section{Acknowledgments} |
| 2255 |
– |
The authors would like to thank espresso for fueling this work, and |
| 2256 |
– |
would also like to send a special acknowledgement to single malt |
| 2257 |
– |
scotch for its wonderful calming effects and its ability to make the |
| 2258 |
– |
troubles of the world float away. |
| 2361 |
|
|
| 2362 |
< |
\bibliographystyle{achemso} |
| 2362 |
> |
Development of {\sc oopse} was funded by a New Faculty Award from the |
| 2363 |
> |
Camille and Henry Dreyfus Foundation and by the National Science |
| 2364 |
> |
Foundation under grant CHE-0134881. Computation time was provided by |
| 2365 |
> |
the Notre Dame Bunch-of-Boxes (B.o.B) computer cluster under NSF grant |
| 2366 |
> |
DMR-0079647. |
| 2367 |
> |
|
| 2368 |
> |
\bibliographystyle{jcc} |
| 2369 |
|
\bibliography{oopsePaper} |
| 2370 |
|
|
| 2371 |
|
\end{document} |
