1 |
mmeineke |
1121 |
\documentclass[11pt]{article} |
2 |
|
|
\usepackage{amsmath} |
3 |
gezelter |
1439 |
\usepackage{amssymb} |
4 |
mmeineke |
1121 |
\usepackage{endfloat} |
5 |
|
|
\usepackage{listings} |
6 |
gezelter |
1428 |
\usepackage{berkeley} |
7 |
mmeineke |
1121 |
\usepackage{graphicx} |
8 |
|
|
\usepackage[ref]{overcite} |
9 |
|
|
\usepackage{setspace} |
10 |
|
|
\usepackage{tabularx} |
11 |
|
|
\pagestyle{plain} |
12 |
|
|
\pagenumbering{arabic} |
13 |
|
|
\oddsidemargin 0.0cm \evensidemargin 0.0cm |
14 |
|
|
\topmargin -21pt \headsep 10pt |
15 |
|
|
\textheight 9.0in \textwidth 6.5in |
16 |
|
|
\brokenpenalty=10000 |
17 |
|
|
\renewcommand{\baselinestretch}{1.2} |
18 |
|
|
\renewcommand\citemid{\ } % no comma in optional reference note |
19 |
|
|
|
20 |
|
|
\begin{document} |
21 |
mmeineke |
1123 |
\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 |
mmeineke |
1121 |
\renewcommand{\lstlistingname}{Scheme} |
25 |
gezelter |
1439 |
\title{{\sc oopse}: An Object-Oriented Parallel Simulation |
26 |
mmeineke |
1121 |
Engine for Molecular Dynamics} |
27 |
|
|
|
28 |
mmeineke |
1155 |
\author{Matthew A. Meineke, Charles F. Vardeman II, Teng Lin,\\ |
29 |
|
|
Christopher J. Fennell and J. Daniel Gezelter\\ |
30 |
mmeineke |
1121 |
Department of Chemistry and Biochemistry\\ |
31 |
|
|
University of Notre Dame\\ |
32 |
|
|
Notre Dame, Indiana 46556} |
33 |
|
|
|
34 |
|
|
\date{\today} |
35 |
|
|
\maketitle |
36 |
|
|
|
37 |
|
|
\begin{abstract} |
38 |
gezelter |
1428 |
{\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 |
gezelter |
1439 |
and the basic integrator for orientational dynamics provides |
47 |
|
|
substantial improvements over older quaternion-based schemes. |
48 |
mmeineke |
1121 |
\end{abstract} |
49 |
|
|
|
50 |
|
|
\section{\label{sec:intro}Introduction} |
51 |
|
|
|
52 |
gezelter |
1428 |
There are a number of excellent molecular dynamics packages available |
53 |
|
|
to the chemical physics |
54 |
gezelter |
1439 |
community.\cite{Brooks83,MacKerell98,pearlman:1995,Gromacs,Gromacs3,DL_POLY,Tinker,Paradyn,namd,macromodel} |
55 |
gezelter |
1428 |
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 |
gezelter |
1439 |
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 |
mmeineke |
1121 |
|
70 |
gezelter |
1428 |
This paper communicates the algorithmic details of our program, which |
71 |
gezelter |
1439 |
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 |
gezelter |
1428 |
(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 |
gezelter |
1425 |
Sec.~\ref{oopseSec:parallelization}. Concluding remarks are presented |
81 |
|
|
in Sec.~\ref{oopseSec:conclusion}. |
82 |
mmeineke |
1121 |
|
83 |
mmeineke |
1155 |
\section{\label{oopseSec:IOfiles}Concepts \& Files} |
84 |
mmeineke |
1121 |
|
85 |
gezelter |
1425 |
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 |
mmeineke |
1155 |
|
92 |
gezelter |
1425 |
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 |
mmeineke |
1155 |
|
100 |
gezelter |
1425 |
{\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 |
mmeineke |
1155 |
|
114 |
gezelter |
1425 |
{\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 |
|
|
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 |
mmeineke |
1155 |
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 |
gezelter |
1425 |
initialConfig = "./argon.in"; |
170 |
mmeineke |
1155 |
|
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 |
gezelter |
1425 |
Within the meta-data file it is necessary to provide a complete |
181 |
mmeineke |
1155 |
description of the molecule before it is actually placed in the |
182 |
gezelter |
1425 |
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 |
mmeineke |
1155 |
|
192 |
gezelter |
1425 |
\begin{lstlisting}[float,caption={An example molecule definition in an |
193 |
|
|
include file.},label={sch:mdIncludeExample}] |
194 |
mmeineke |
1155 |
|
195 |
|
|
molecule{ |
196 |
|
|
name = "Ar"; |
197 |
|
|
nAtoms = 1; |
198 |
|
|
atom[0]{ |
199 |
|
|
type="Ar"; |
200 |
|
|
position( 0.0, 0.0, 0.0 ); |
201 |
|
|
} |
202 |
|
|
} |
203 |
|
|
|
204 |
|
|
\end{lstlisting} |
205 |
|
|
|
206 |
gezelter |
1425 |
\begin{lstlisting}[float,caption={Revised meta-data example.},label={sch:mdExPrime}] |
207 |
mmeineke |
1155 |
|
208 |
gezelter |
1425 |
#include "argon.md" |
209 |
mmeineke |
1155 |
|
210 |
|
|
nComponents = 1; |
211 |
|
|
component{ |
212 |
|
|
type = "Ar"; |
213 |
|
|
nMol = 108; |
214 |
|
|
} |
215 |
|
|
|
216 |
gezelter |
1425 |
initialConfig = "./argon.in"; |
217 |
mmeineke |
1155 |
|
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 |
gezelter |
1425 |
\subsection{\label{oopseSec:atomsMolecules}Atoms, Molecules, and other |
228 |
|
|
ways of grouping atoms} |
229 |
mmeineke |
1121 |
|
230 |
gezelter |
1425 |
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 |
mmeineke |
1121 |
|
243 |
gezelter |
1425 |
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 |
mmeineke |
1121 |
potential and move collectively.\cite{Goldstein01} They are not |
248 |
|
|
included in most simulation packages because of the algorithmic |
249 |
gezelter |
1425 |
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 |
mmeineke |
1121 |
|
254 |
|
|
Moving a rigid body involves determination of both the force and |
255 |
|
|
torque applied by the surroundings, which directly affect the |
256 |
|
|
translational and rotational motion in turn. In order to accumulate |
257 |
|
|
the total force on a rigid body, the external forces and torques must |
258 |
|
|
first be calculated for all the internal particles. The total force on |
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 |
gezelter |
1425 |
the rigid body. The space-fixed torque on rigid body $i$ is |
263 |
mmeineke |
1121 |
\begin{equation} |
264 |
|
|
\boldsymbol{\tau}_i= |
265 |
|
|
\sum_{a}\biggl[(\mathbf{r}_{ia}-\mathbf{r}_i)\times \mathbf{f}_{ia} |
266 |
|
|
+ \boldsymbol{\tau}_{ia}\biggr], |
267 |
|
|
\label{eq:torqueAccumulate} |
268 |
|
|
\end{equation} |
269 |
|
|
where $\boldsymbol{\tau}_i$ and $\mathbf{r}_i$ are the torque on and |
270 |
|
|
position of the center of mass respectively, while $\mathbf{f}_{ia}$, |
271 |
|
|
$\mathbf{r}_{ia}$, and $\boldsymbol{\tau}_{ia}$ are the force on, |
272 |
|
|
position of, and torque on the component particles of the rigid body. |
273 |
|
|
|
274 |
|
|
The summation of the total torque is done in the body fixed axis of |
275 |
|
|
each rigid body. In order to move between the space fixed and body |
276 |
|
|
fixed coordinate axes, parameters describing the orientation must be |
277 |
|
|
maintained for each rigid body. At a minimum, the rotation matrix |
278 |
|
|
($\mathsf{A}$) can be described by the three Euler angles ($\phi, |
279 |
|
|
\theta,$ and $\psi$), where the elements of $\mathsf{A}$ are composed of |
280 |
|
|
trigonometric operations involving $\phi, \theta,$ and |
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 |
gezelter |
1439 |
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 |
mmeineke |
1121 |
performance enhancements, particularly for very small |
287 |
|
|
systems.\cite{Evans77} |
288 |
|
|
|
289 |
gezelter |
1425 |
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 |
mmeineke |
1168 |
\ref{sch:rigidBody}. |
295 |
mmeineke |
1121 |
|
296 |
gezelter |
1425 |
\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 |
mmeineke |
1121 |
molecule{ |
300 |
|
|
name = "TIP3P"; |
301 |
|
|
nAtoms = 3; |
302 |
|
|
atom[0]{ |
303 |
|
|
type = "O_TIP3P"; |
304 |
|
|
position( 0.0, 0.0, -0.06556 ); |
305 |
|
|
} |
306 |
|
|
atom[1]{ |
307 |
|
|
type = "H_TIP3P"; |
308 |
|
|
position( 0.0, 0.75695, 0.52032 ); |
309 |
|
|
} |
310 |
|
|
atom[2]{ |
311 |
|
|
type = "H_TIP3P"; |
312 |
|
|
position( 0.0, -0.75695, 0.52032 ); |
313 |
|
|
} |
314 |
|
|
|
315 |
|
|
nRigidBodies = 1; |
316 |
|
|
rigidBody[0]{ |
317 |
|
|
nMembers = 3; |
318 |
|
|
members(0, 1, 2); |
319 |
|
|
} |
320 |
gezelter |
1425 |
|
321 |
|
|
nCutoffGroups = 1; |
322 |
|
|
cutoffGroup[0]{ |
323 |
|
|
nMembers = 3; |
324 |
|
|
members(0, 1, 2); |
325 |
|
|
} |
326 |
mmeineke |
1121 |
} |
327 |
|
|
\end{lstlisting} |
328 |
|
|
|
329 |
gezelter |
1425 |
\subsection{\label{sec:miscConcepts}Creating a Metadata File} |
330 |
mmeineke |
1155 |
|
331 |
gezelter |
1425 |
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 |
mmeineke |
1155 |
|
340 |
gezelter |
1425 |
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 |
mmeineke |
1155 |
|
353 |
gezelter |
1425 |
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 |
gezelter |
1428 |
outlined in Sec.~\ref{oopseSec:empiricalEnergy}. The force fields are |
368 |
gezelter |
1425 |
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 |
mmeineke |
1168 |
|
372 |
gezelter |
1425 |
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 |
gezelter |
1439 |
translational velocities are in $\mbox{\AA~fs}^{-1}$, and masses are |
393 |
gezelter |
1425 |
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 |
mmeineke |
1168 |
\begin{table} |
399 |
gezelter |
1425 |
\caption{Meta-data Keywords: Required Parameters} |
400 |
mmeineke |
1168 |
\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 |
gezelter |
1439 |
{\tt forceField} & string & Sets the force field. & Possible force fields are DUFF, LJ, and EAM. \\ |
413 |
mmeineke |
1168 |
{\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 |
gezelter |
1425 |
{\tt minimizer}& string & Chooses a minimizer & Possible minimizers |
416 |
gezelter |
1439 |
are SD and CG. Either {\tt ensemble} or {\tt minimizer} must be specified. \\ |
417 |
gezelter |
1425 |
{\tt ensemble} & string & Sets the ensemble. & Possible ensembles are |
418 |
gezelter |
1439 |
NVE, NVT, NPTi, NPTf, and NPTxyz. Either {\tt ensemble} |
419 |
gezelter |
1425 |
or {\tt minimizer} must be specified. \\ |
420 |
|
|
{\tt dt} & fs & Sets the time step. & Selection of {\tt dt} should be |
421 |
gezelter |
1439 |
small enough to sample the fastest motion of the simulation. ({\tt |
422 |
|
|
dt} is required for molecular dynamics simulations)\\ |
423 |
gezelter |
1425 |
{\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 |
gezelter |
1439 |
current time meets or exceeds the {\tt runTime}. ({\tt runTime} is |
426 |
|
|
required for molecular dynamics simulations)\\ |
427 |
mmeineke |
1168 |
|
428 |
|
|
\end{tabularx} |
429 |
|
|
\end{center} |
430 |
|
|
\end{table} |
431 |
|
|
|
432 |
|
|
\begin{table} |
433 |
gezelter |
1425 |
\caption{Meta-data Keywords: General Parameters} |
434 |
mmeineke |
1168 |
\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 |
gezelter |
1425 |
{\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 |
gezelter |
1439 |
{\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 |
gezelter |
1425 |
simulation. \\ |
459 |
|
|
{\tt switchingRadius} & $\mbox{\AA}$ & Manually sets the inner radius for the switching function. & Defaults to 95~\% of the {\tt cutoffRadius}. \\ |
460 |
gezelter |
1439 |
{\tt useReactionField} & logical & Turns the reaction field correction on/off. & Default is false. \\ |
461 |
mmeineke |
1168 |
{\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 |
gezelter |
1439 |
& logical & Turns periodic boundary conditions on/off. & Default is true. \\ |
464 |
gezelter |
1425 |
{\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 |
gezelter |
1439 |
force field. & {\sc eam} has three variants: {\tt u3}, {\tt u6}, and |
469 |
gezelter |
1425 |
{\tt VC}. |
470 |
mmeineke |
1168 |
|
471 |
|
|
\end{tabularx} |
472 |
|
|
\end{center} |
473 |
|
|
\end{table} |
474 |
|
|
|
475 |
|
|
|
476 |
mmeineke |
1155 |
\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 |
gezelter |
1425 |
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 |
mmeineke |
1155 |
|
491 |
gezelter |
1425 |
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 |
mmeineke |
1155 |
|
499 |
gezelter |
1425 |
\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 |
mmeineke |
1155 |
time; Hxx Hyx Hzx; Hxy Hyy Hzy; Hxz Hyz Hzz; |
512 |
gezelter |
1425 |
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 |
mmeineke |
1155 |
etc... |
515 |
|
|
|
516 |
|
|
\end{lstlisting} |
517 |
|
|
|
518 |
gezelter |
1425 |
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 |
mmeineke |
1155 |
|
524 |
gezelter |
1425 |
\begin{tt} |
525 |
|
|
TIP3P\_RB\_0 x y z vx vy vz qw qx qy qz jx jy jz |
526 |
|
|
\end{tt} |
527 |
mmeineke |
1155 |
|
528 |
gezelter |
1425 |
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 |
mmeineke |
1155 |
\subsection{\label{oopseSec:initCoords}Generation of Initial Coordinates} |
552 |
|
|
|
553 |
gezelter |
1425 |
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 |
mmeineke |
1155 |
\begin{enumerate} |
562 |
gezelter |
1425 |
\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 |
gezelter |
1439 |
included. |
569 |
gezelter |
1425 |
\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 |
mmeineke |
1155 |
\end{enumerate} |
574 |
gezelter |
1425 |
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 |
mmeineke |
1155 |
|
577 |
gezelter |
1425 |
\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 |
mmeineke |
1155 |
|
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 |
gezelter |
1425 |
initialConfig = "mixture.in"; |
622 |
mmeineke |
1155 |
|
623 |
|
|
\end{lstlisting} |
624 |
|
|
|
625 |
gezelter |
1425 |
\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 |
mmeineke |
1155 |
|
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 |
gezelter |
1425 |
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 |
mmeineke |
1155 |
\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 |
gezelter |
1425 |
monitor the conserved quantity of a given simulation ensemble, |
658 |
|
|
allowing the user to gauge the stability of the integrator. The |
659 |
mmeineke |
1155 |
statistics file is denoted with the \texttt{.stat} file extension. |
660 |
|
|
|
661 |
gezelter |
1425 |
\section{\label{oopseSec:empiricalEnergy}The Empirical Energy |
662 |
|
|
Functions} |
663 |
mmeineke |
1155 |
|
664 |
gezelter |
1425 |
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 |
gezelter |
1439 |
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 |
mmeineke |
1155 |
|
676 |
gezelter |
1439 |
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 |
gezelter |
1425 |
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 |
gezelter |
1439 |
The sums over $a$ and $b$ are over the cutoff groups that are present |
696 |
gezelter |
1425 |
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 |
gezelter |
1439 |
dipolar atoms, the default cutoff radius is $15 \mbox{\AA}$. |
722 |
gezelter |
1425 |
|
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 |
gezelter |
1439 |
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 |
gezelter |
1425 |
|
734 |
mmeineke |
1121 |
\subsection{\label{sec:LJPot}The Lennard Jones Force Field} |
735 |
|
|
|
736 |
|
|
The most basic force field implemented in {\sc oopse} is the |
737 |
gezelter |
1425 |
Lennard-Jones force field, which mimics the van der Waals interaction |
738 |
|
|
at long distances and uses an empirical repulsion at short |
739 |
mmeineke |
1121 |
distances. The Lennard-Jones potential is given by: |
740 |
|
|
\begin{equation} |
741 |
|
|
V_{\text{LJ}}(r_{ij}) = |
742 |
|
|
4\epsilon_{ij} \biggl[ |
743 |
|
|
\biggl(\frac{\sigma_{ij}}{r_{ij}}\biggr)^{12} |
744 |
|
|
- \biggl(\frac{\sigma_{ij}}{r_{ij}}\biggr)^{6} |
745 |
|
|
\biggr], |
746 |
|
|
\label{eq:lennardJonesPot} |
747 |
|
|
\end{equation} |
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 |
gezelter |
1425 |
\ref{sch:LJFF} gives an example meta-data file that |
752 |
mmeineke |
1121 |
sets up a system of 108 Ar particles to be simulated using the |
753 |
|
|
Lennard-Jones force field. |
754 |
|
|
|
755 |
gezelter |
1425 |
\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 |
mmeineke |
1121 |
|
759 |
gezelter |
1425 |
#include "argon.md" |
760 |
mmeineke |
1121 |
|
761 |
|
|
nComponents = 1; |
762 |
|
|
component{ |
763 |
|
|
type = "Ar"; |
764 |
|
|
nMol = 108; |
765 |
|
|
} |
766 |
|
|
|
767 |
gezelter |
1425 |
initialConfig = "./argon.in"; |
768 |
mmeineke |
1121 |
|
769 |
|
|
forceField = "LJ"; |
770 |
|
|
\end{lstlisting} |
771 |
|
|
|
772 |
|
|
Interactions between dissimilar particles requires the generation of |
773 |
gezelter |
1425 |
cross term parameters for $\sigma$ and $\epsilon$. These parameters |
774 |
|
|
are determined using the Lorentz-Berthelot mixing |
775 |
mmeineke |
1121 |
rules:\cite{allen87:csl} |
776 |
|
|
\begin{equation} |
777 |
|
|
\sigma_{ij} = \frac{1}{2}[\sigma_{ii} + \sigma_{jj}], |
778 |
|
|
\label{eq:sigmaMix} |
779 |
|
|
\end{equation} |
780 |
|
|
and |
781 |
|
|
\begin{equation} |
782 |
|
|
\epsilon_{ij} = \sqrt{\epsilon_{ii} \epsilon_{jj}}. |
783 |
|
|
\label{eq:epsilonMix} |
784 |
|
|
\end{equation} |
785 |
|
|
|
786 |
|
|
\subsection{\label{oopseSec:DUFF}Dipolar Unified-Atom Force Field} |
787 |
|
|
|
788 |
|
|
The dipolar unified-atom force field ({\sc duff}) was developed to |
789 |
gezelter |
1425 |
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 |
mmeineke |
1121 |
|
803 |
|
|
As an example, lipid head-groups in {\sc duff} are represented as |
804 |
gezelter |
1425 |
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 |
mmeineke |
1121 |
\emph{et al.}\cite{liu96:new_model} |
814 |
|
|
|
815 |
|
|
\begin{figure} |
816 |
|
|
\centering |
817 |
gezelter |
1425 |
\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 |
mmeineke |
1121 |
\label{oopseFig:lipidModel} |
823 |
|
|
\end{figure} |
824 |
|
|
|
825 |
gezelter |
1425 |
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 |
mmeineke |
1121 |
|
836 |
gezelter |
1425 |
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 |
mmeineke |
1121 |
|
845 |
gezelter |
1425 |
\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 |
mmeineke |
1121 |
|
849 |
gezelter |
1425 |
#include "water.md" |
850 |
|
|
#include "lipid.md" |
851 |
mmeineke |
1121 |
|
852 |
|
|
nComponents = 2; |
853 |
|
|
component{ |
854 |
|
|
type = "simpleLipid_16"; |
855 |
|
|
nMol = 60; |
856 |
|
|
} |
857 |
|
|
|
858 |
|
|
component{ |
859 |
|
|
type = "SSD_water"; |
860 |
|
|
nMol = 1936; |
861 |
|
|
} |
862 |
|
|
|
863 |
gezelter |
1425 |
initialConfig = "bilayer.in"; |
864 |
mmeineke |
1121 |
|
865 |
|
|
forceField = "DUFF"; |
866 |
|
|
|
867 |
|
|
\end{lstlisting} |
868 |
|
|
|
869 |
mmeineke |
1168 |
\subsubsection{\label{oopseSec:energyFunctions}{\sc duff} Energy Functions} |
870 |
mmeineke |
1121 |
|
871 |
|
|
The total potential energy function in {\sc duff} is |
872 |
|
|
\begin{equation} |
873 |
|
|
V = \sum^{N}_{I=1} V^{I}_{\text{Internal}} |
874 |
|
|
+ \sum^{N-1}_{I=1} \sum_{J>I} V^{IJ}_{\text{Cross}}, |
875 |
|
|
\label{eq:totalPotential} |
876 |
|
|
\end{equation} |
877 |
|
|
where $V^{I}_{\text{Internal}}$ is the internal potential of molecule $I$: |
878 |
|
|
\begin{equation} |
879 |
|
|
V^{I}_{\text{Internal}} = |
880 |
|
|
\sum_{\theta_{ijk} \in I} V_{\text{bend}}(\theta_{ijk}) |
881 |
|
|
+ \sum_{\phi_{ijkl} \in I} V_{\text{torsion}}(\phi_{ijkl}) |
882 |
|
|
+ \sum_{i \in I} \sum_{(j>i+4) \in I} |
883 |
|
|
\biggl[ V_{\text{LJ}}(r_{ij}) + V_{\text{dipole}} |
884 |
|
|
(\mathbf{r}_{ij},\boldsymbol{\Omega}_{i},\boldsymbol{\Omega}_{j}) |
885 |
|
|
\biggr]. |
886 |
|
|
\label{eq:internalPotential} |
887 |
|
|
\end{equation} |
888 |
|
|
Here $V_{\text{bend}}$ is the bend potential for all 1, 3 bonded pairs |
889 |
gezelter |
1425 |
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 |
mmeineke |
1121 |
|
895 |
|
|
The bend potential of a molecule is represented by the following function: |
896 |
|
|
\begin{equation} |
897 |
gezelter |
1425 |
V_{\text{bend}}(\theta_{ijk}) = k_{\theta}( \theta_{ijk} - \theta_0 |
898 |
|
|
)^2, \label{eq:bendPot} |
899 |
mmeineke |
1121 |
\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 |
902 |
|
|
bond angle, and $k_{\theta}$ is the force constant which determines the |
903 |
|
|
strength of the harmonic bend. The parameters for $k_{\theta}$ and |
904 |
|
|
$\theta_0$ are borrowed from those in TraPPE.\cite{Siepmann1998} |
905 |
|
|
|
906 |
|
|
The torsion potential and parameters are also borrowed from TraPPE. It is |
907 |
|
|
of the form: |
908 |
|
|
\begin{equation} |
909 |
|
|
V_{\text{torsion}}(\phi) = c_1[1 + \cos \phi] |
910 |
|
|
+ c_2[1 + \cos(2\phi)] |
911 |
|
|
+ c_3[1 + \cos(3\phi)], |
912 |
|
|
\label{eq:origTorsionPot} |
913 |
|
|
\end{equation} |
914 |
|
|
where: |
915 |
|
|
\begin{equation} |
916 |
|
|
\cos\phi = (\hat{\mathbf{r}}_{ij} \times \hat{\mathbf{r}}_{jk}) \cdot |
917 |
|
|
(\hat{\mathbf{r}}_{jk} \times \hat{\mathbf{r}}_{kl}). |
918 |
|
|
\label{eq:torsPhi} |
919 |
|
|
\end{equation} |
920 |
|
|
Here, $\hat{\mathbf{r}}_{\alpha\beta}$ are the set of unit bond |
921 |
|
|
vectors between atoms $i$, $j$, $k$, and $l$. For computational |
922 |
|
|
efficiency, the torsion potential has been recast after the method of |
923 |
|
|
{\sc charmm},\cite{Brooks83} in which the angle series is converted to |
924 |
|
|
a power series of the form: |
925 |
|
|
\begin{equation} |
926 |
|
|
V_{\text{torsion}}(\phi) = |
927 |
|
|
k_3 \cos^3 \phi + k_2 \cos^2 \phi + k_1 \cos \phi + k_0, |
928 |
|
|
\label{eq:torsionPot} |
929 |
|
|
\end{equation} |
930 |
|
|
where: |
931 |
|
|
\begin{align*} |
932 |
|
|
k_0 &= c_1 + c_3, \\ |
933 |
|
|
k_1 &= c_1 - 3c_3, \\ |
934 |
|
|
k_2 &= 2 c_2, \\ |
935 |
|
|
k_3 &= 4c_3. |
936 |
|
|
\end{align*} |
937 |
|
|
By recasting the potential as a power series, repeated trigonometric |
938 |
gezelter |
1425 |
evaluations are avoided during the calculation of the potential |
939 |
|
|
energy. |
940 |
mmeineke |
1121 |
|
941 |
|
|
|
942 |
gezelter |
1425 |
The cross potential between molecules $I$ and $J$, |
943 |
|
|
$V^{IJ}_{\text{Cross}}$, is as follows: |
944 |
mmeineke |
1121 |
\begin{equation} |
945 |
|
|
V^{IJ}_{\text{Cross}} = |
946 |
|
|
\sum_{i \in I} \sum_{j \in J} |
947 |
|
|
\biggl[ V_{\text{LJ}}(r_{ij}) + V_{\text{dipole}} |
948 |
|
|
(\mathbf{r}_{ij},\boldsymbol{\Omega}_{i},\boldsymbol{\Omega}_{j}) |
949 |
|
|
+ V_{\text{sticky}} |
950 |
|
|
(\mathbf{r}_{ij},\boldsymbol{\Omega}_{i},\boldsymbol{\Omega}_{j}) |
951 |
|
|
\biggr], |
952 |
|
|
\label{eq:crossPotentail} |
953 |
|
|
\end{equation} |
954 |
|
|
where $V_{\text{LJ}}$ is the Lennard Jones potential, |
955 |
|
|
$V_{\text{dipole}}$ is the dipole dipole potential, and |
956 |
|
|
$V_{\text{sticky}}$ is the sticky potential defined by the SSD model |
957 |
|
|
(Sec.~\ref{oopseSec:SSD}). Note that not all atom types include all |
958 |
|
|
interactions. |
959 |
|
|
|
960 |
|
|
The dipole-dipole potential has the following form: |
961 |
|
|
\begin{equation} |
962 |
|
|
V_{\text{dipole}}(\mathbf{r}_{ij},\boldsymbol{\Omega}_{i}, |
963 |
|
|
\boldsymbol{\Omega}_{j}) = \frac{|\mu_i||\mu_j|}{4\pi\epsilon_{0}r_{ij}^{3}} \biggl[ |
964 |
|
|
\boldsymbol{\hat{u}}_{i} \cdot \boldsymbol{\hat{u}}_{j} |
965 |
|
|
- |
966 |
|
|
3(\boldsymbol{\hat{u}}_i \cdot \hat{\mathbf{r}}_{ij}) % |
967 |
|
|
(\boldsymbol{\hat{u}}_j \cdot \hat{\mathbf{r}}_{ij}) \biggr]. |
968 |
|
|
\label{eq:dipolePot} |
969 |
|
|
\end{equation} |
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 |
gezelter |
1425 |
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 |
mmeineke |
1121 |
($\boldsymbol{\hat{r}}_{ij}=\mathbf{r}_{ij}/|\mathbf{r}_{ij}|$). |
978 |
|
|
|
979 |
gezelter |
1425 |
\subsubsection{\label{oopseSec:SSD}The {\sc duff} Water Models: SSD/E |
980 |
|
|
and SSD/RF} |
981 |
mmeineke |
1121 |
|
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 |
985 |
|
|
\emph{et al.}\cite{liu96:new_model} as a modified form of the hard-sphere |
986 |
|
|
water model proposed by Bratko, Blum, and |
987 |
|
|
Luzar.\cite{Bratko85,Bratko95} It consists of a single point dipole |
988 |
|
|
with a Lennard-Jones core and a sticky potential that directs the |
989 |
|
|
particles to assume the proper hydrogen bond orientation in the first |
990 |
|
|
solvation shell. Thus, the interaction between two SSD water molecules |
991 |
|
|
\emph{i} and \emph{j} is given by the potential |
992 |
|
|
\begin{equation} |
993 |
|
|
V_{ij} = |
994 |
|
|
V_{ij}^{LJ} (r_{ij})\ + V_{ij}^{dp} |
995 |
|
|
(\mathbf{r}_{ij},\boldsymbol{\Omega}_i,\boldsymbol{\Omega}_j)\ + |
996 |
|
|
V_{ij}^{sp} |
997 |
|
|
(\mathbf{r}_{ij},\boldsymbol{\Omega}_i,\boldsymbol{\Omega}_j), |
998 |
|
|
\label{eq:ssdPot} |
999 |
|
|
\end{equation} |
1000 |
|
|
where the $\mathbf{r}_{ij}$ is the position vector between molecules |
1001 |
|
|
\emph{i} and \emph{j} with magnitude equal to the distance $r_{ij}$, and |
1002 |
|
|
$\boldsymbol{\Omega}_i$ and $\boldsymbol{\Omega}_j$ represent the |
1003 |
|
|
orientations of the respective molecules. The Lennard-Jones and dipole |
1004 |
|
|
parts of the potential are given by equations \ref{eq:lennardJonesPot} |
1005 |
|
|
and \ref{eq:dipolePot} respectively. The sticky part is described by |
1006 |
|
|
the following, |
1007 |
|
|
\begin{equation} |
1008 |
|
|
u_{ij}^{sp}(\mathbf{r}_{ij},\boldsymbol{\Omega}_i,\boldsymbol{\Omega}_j)= |
1009 |
|
|
\frac{\nu_0}{2}[s(r_{ij})w(\mathbf{r}_{ij}, |
1010 |
|
|
\boldsymbol{\Omega}_i,\boldsymbol{\Omega}_j) + |
1011 |
|
|
s^\prime(r_{ij})w^\prime(\mathbf{r}_{ij}, |
1012 |
|
|
\boldsymbol{\Omega}_i,\boldsymbol{\Omega}_j)]\ , |
1013 |
|
|
\label{eq:stickyPot} |
1014 |
|
|
\end{equation} |
1015 |
|
|
where $\nu_0$ is a strength parameter for the sticky potential, and |
1016 |
|
|
$s$ and $s^\prime$ are cubic switching functions which turn off the |
1017 |
|
|
sticky interaction beyond the first solvation shell. The $w$ function |
1018 |
|
|
can be thought of as an attractive potential with tetrahedral |
1019 |
|
|
geometry: |
1020 |
|
|
\begin{equation} |
1021 |
|
|
w({\bf r}_{ij},{\bf \Omega}_i,{\bf \Omega}_j)= |
1022 |
|
|
\sin\theta_{ij}\sin2\theta_{ij}\cos2\phi_{ij}, |
1023 |
|
|
\label{eq:stickyW} |
1024 |
|
|
\end{equation} |
1025 |
|
|
while the $w^\prime$ function counters the normal aligned and |
1026 |
|
|
anti-aligned structures favored by point dipoles: |
1027 |
|
|
\begin{equation} |
1028 |
|
|
w^\prime({\bf r}_{ij},{\bf \Omega}_i,{\bf \Omega}_j)= |
1029 |
|
|
(\cos\theta_{ij}-0.6)^2(\cos\theta_{ij}+0.8)^2-w^0, |
1030 |
|
|
\label{eq:stickyWprime} |
1031 |
|
|
\end{equation} |
1032 |
|
|
It should be noted that $w$ is proportional to the sum of the $Y_3^2$ |
1033 |
|
|
and $Y_3^{-2}$ spherical harmonics (a linear combination which |
1034 |
|
|
enhances the tetrahedral geometry for hydrogen bonded structures), |
1035 |
|
|
while $w^\prime$ is a purely empirical function. A more detailed |
1036 |
|
|
description of the functional parts and variables in this potential |
1037 |
|
|
can be found in the original SSD |
1038 |
|
|
articles.\cite{liu96:new_model,liu96:monte_carlo,chandra99:ssd_md,Ichiye03} |
1039 |
|
|
|
1040 |
gezelter |
1425 |
\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 |
mmeineke |
1121 |
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 |
gezelter |
1425 |
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 |
mmeineke |
1121 |
|
1068 |
|
|
Recent constant pressure simulations revealed issues in the original |
1069 |
|
|
SSD model that led to lower than expected densities at all target |
1070 |
|
|
pressures.\cite{Ichiye03,fennell04} The default model in {\sc oopse} |
1071 |
|
|
is therefore SSD/E, a density corrected derivative of SSD that |
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 |
gezelter |
1425 |
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 |
mmeineke |
1121 |
|
1081 |
chrisfen |
1434 |
\subsection{\label{oopseSec:WATER}The {\sc water} Force Field} |
1082 |
|
|
|
1083 |
|
|
In addition to the {\sc duff} force field's solvent description, a |
1084 |
gezelter |
1439 |
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 |
chrisfen |
1434 |
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 |
gezelter |
1439 |
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 |
chrisfen |
1434 |
|
1104 |
mmeineke |
1121 |
\subsection{\label{oopseSec:eam}Embedded Atom Method} |
1105 |
|
|
|
1106 |
gezelter |
1425 |
{\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 |
mmeineke |
1121 |
free valance ``sea'' of electrons created by the surrounding atoms in |
1112 |
gezelter |
1425 |
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 |
mmeineke |
1121 |
|
1120 |
|
|
The {\sc eam} potential has the form: |
1121 |
gezelter |
1425 |
\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 |
mmeineke |
1121 |
required to embed a positively-charged core ion $i$ into a linear |
1127 |
|
|
superposition of spherically averaged atomic electron densities given |
1128 |
gezelter |
1425 |
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 |
mmeineke |
1121 |
surrounding atom $i$ for both the density $\rho$ and pairwise $\phi$ |
1147 |
gezelter |
1425 |
interactions. |
1148 |
mmeineke |
1121 |
|
1149 |
gezelter |
1425 |
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 |
mmeineke |
1121 |
\subsection{\label{oopseSec:pbc}Periodic Boundary Conditions} |
1184 |
|
|
|
1185 |
|
|
\newcommand{\roundme}{\operatorname{round}} |
1186 |
|
|
|
1187 |
gezelter |
1425 |
\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 |
mmeineke |
1121 |
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 |
gezelter |
1425 |
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 |
mmeineke |
1121 |
\begin{equation} |
1201 |
|
|
\mathsf{H} = ( \mathbf{h}_x, \mathbf{h}_y, \mathbf{h}_z ), |
1202 |
|
|
\end{equation} |
1203 |
|
|
where $\mathbf{h}_{\alpha}$ is the column vector of the $\alpha$ axis of the |
1204 |
|
|
box. During the course of the simulation both the size and shape of |
1205 |
|
|
the box can be changed to allow volume fluctuations when constraining |
1206 |
|
|
the pressure. |
1207 |
|
|
|
1208 |
|
|
A real space vector, $\mathbf{r}$ can be transformed in to a box space |
1209 |
|
|
vector, $\mathbf{s}$, and back through the following transformations: |
1210 |
|
|
\begin{align} |
1211 |
|
|
\mathbf{s} &= \mathsf{H}^{-1} \mathbf{r}, \\ |
1212 |
|
|
\mathbf{r} &= \mathsf{H} \mathbf{s}. |
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 |
gezelter |
1425 |
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 |
mmeineke |
1121 |
\begin{equation} |
1220 |
|
|
s_{i}^{\prime}=s_{i}-\roundme(s_{i}), |
1221 |
|
|
\end{equation} |
1222 |
|
|
where $s_i$ is the $i$th element of $\mathbf{s}$, and |
1223 |
|
|
$\roundme(s_i)$ is given by |
1224 |
|
|
\begin{equation} |
1225 |
|
|
\roundme(x) = |
1226 |
|
|
\begin{cases} |
1227 |
|
|
\lfloor x+0.5 \rfloor & \text{if $x \ge 0$,} \\ |
1228 |
|
|
\lceil x-0.5 \rceil & \text{if $x < 0$.} |
1229 |
|
|
\end{cases} |
1230 |
|
|
\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 |
gezelter |
1425 |
than $x$. |
1235 |
mmeineke |
1121 |
|
1236 |
gezelter |
1425 |
Finally, the minimum image coordinates $\mathbf{r}^{\prime}$ are |
1237 |
|
|
obtained by transforming back to real space, |
1238 |
mmeineke |
1121 |
\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 |
gezelter |
1425 |
but their minimum images, or $\mathbf{r}^{\prime}$, are used to compute |
1243 |
mmeineke |
1121 |
the inter-atomic forces. |
1244 |
|
|
|
1245 |
|
|
|
1246 |
|
|
|
1247 |
|
|
\section{\label{oopseSec:mechanics}Mechanics} |
1248 |
|
|
|
1249 |
|
|
\subsection{\label{oopseSec:integrate}Integrating the Equations of Motion: the |
1250 |
gezelter |
1439 |
{\sc dlm} method} |
1251 |
mmeineke |
1121 |
|
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 |
gezelter |
1439 |
({\sc dlm}).\cite{Dullweber1997} When there are no directional atoms or |
1256 |
mmeineke |
1121 |
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 |
gezelter |
1439 |
that are avoided in the {\sc dlm} method. Direct propagation of the Euler |
1262 |
mmeineke |
1121 |
angles has a known $1/\sin\theta$ divergence in the equations of |
1263 |
gezelter |
1425 |
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 |
mmeineke |
1121 |
|
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 |
1276 |
|
|
propagated from one time step to the next. In the past, this would not |
1277 |
|
|
have been feasible, since the rotation matrix for a single body has |
1278 |
|
|
nine elements compared with the more memory-efficient methods (using |
1279 |
|
|
three Euler angles or 4 quaternions). Computer memory has become much |
1280 |
|
|
less costly in recent years, and this can be translated into |
1281 |
|
|
substantial benefits in energy conservation. |
1282 |
|
|
|
1283 |
|
|
The basic equations of motion being integrated are derived from the |
1284 |
|
|
Hamiltonian for conservative systems containing rigid bodies, |
1285 |
|
|
\begin{equation} |
1286 |
|
|
H = \sum_{i} \left( \frac{1}{2} m_i {\bf v}_i^T \cdot {\bf v}_i + |
1287 |
|
|
\frac{1}{2} {\bf j}_i^T \cdot \overleftrightarrow{\mathsf{I}}_i^{-1} \cdot |
1288 |
|
|
{\bf j}_i \right) + |
1289 |
|
|
V\left(\left\{{\bf r}\right\}, \left\{\mathsf{A}\right\}\right), |
1290 |
|
|
\end{equation} |
1291 |
|
|
where ${\bf r}_i$ and ${\bf v}_i$ are the cartesian position vector |
1292 |
|
|
and velocity of the center of mass of particle $i$, and ${\bf j}_i$, |
1293 |
|
|
$\overleftrightarrow{\mathsf{I}}_i$ are the body-fixed angular |
1294 |
|
|
momentum and moment of inertia tensor respectively, and the |
1295 |
|
|
superscript $T$ denotes the transpose of the vector. $\mathsf{A}_i$ |
1296 |
|
|
is the $3 \times 3$ rotation matrix describing the instantaneous |
1297 |
|
|
orientation of the particle. $V$ is the potential energy function |
1298 |
|
|
which may depend on both the positions $\left\{{\bf r}\right\}$ and |
1299 |
|
|
orientations $\left\{\mathsf{A}\right\}$ of all particles. The |
1300 |
|
|
equations of motion for the particle centers of mass are derived from |
1301 |
|
|
Hamilton's equations and are quite simple, |
1302 |
|
|
\begin{eqnarray} |
1303 |
|
|
\dot{{\bf r}} & = & {\bf v}, \\ |
1304 |
|
|
\dot{{\bf v}} & = & \frac{{\bf f}}{m}, |
1305 |
|
|
\end{eqnarray} |
1306 |
|
|
where ${\bf f}$ is the instantaneous force on the center of mass |
1307 |
|
|
of the particle, |
1308 |
|
|
\begin{equation} |
1309 |
|
|
{\bf f} = - \frac{\partial}{\partial |
1310 |
|
|
{\bf r}} V(\left\{{\bf r}(t)\right\}, \left\{\mathsf{A}(t)\right\}). |
1311 |
|
|
\end{equation} |
1312 |
|
|
|
1313 |
|
|
The equations of motion for the orientational degrees of freedom are |
1314 |
|
|
\begin{eqnarray} |
1315 |
|
|
\dot{\mathsf{A}} & = & \mathsf{A} \cdot |
1316 |
|
|
\mbox{ skew}\left(\overleftrightarrow{\mathsf{I}}^{-1} \cdot {\bf j}\right),\\ |
1317 |
|
|
\dot{{\bf j}} & = & {\bf j} \times \left( \overleftrightarrow{\mathsf{I}}^{-1} |
1318 |
|
|
\cdot {\bf j} \right) - \mbox{ rot}\left(\mathsf{A}^{T} \cdot \frac{\partial |
1319 |
|
|
V}{\partial \mathsf{A}} \right). |
1320 |
|
|
\end{eqnarray} |
1321 |
|
|
In these equations of motion, the $\mbox{skew}$ matrix of a vector |
1322 |
|
|
${\bf v} = \left( v_1, v_2, v_3 \right)$ is defined: |
1323 |
|
|
\begin{equation} |
1324 |
|
|
\mbox{skew}\left( {\bf v} \right) := \left( |
1325 |
|
|
\begin{array}{ccc} |
1326 |
|
|
0 & v_3 & - v_2 \\ |
1327 |
|
|
-v_3 & 0 & v_1 \\ |
1328 |
|
|
v_2 & -v_1 & 0 |
1329 |
|
|
\end{array} |
1330 |
|
|
\right). |
1331 |
|
|
\end{equation} |
1332 |
|
|
The $\mbox{rot}$ notation refers to the mapping of the $3 \times 3$ |
1333 |
|
|
rotation matrix to a vector of orientations by first computing the |
1334 |
|
|
skew-symmetric part $\left(\mathsf{A} - \mathsf{A}^{T}\right)$ and |
1335 |
|
|
then associating this with a length 3 vector by inverting the |
1336 |
|
|
$\mbox{skew}$ function above: |
1337 |
|
|
\begin{equation} |
1338 |
|
|
\mbox{rot}\left(\mathsf{A}\right) := \mbox{ skew}^{-1}\left(\mathsf{A} |
1339 |
|
|
- \mathsf{A}^{T} \right). |
1340 |
|
|
\end{equation} |
1341 |
|
|
Written this way, the $\mbox{rot}$ operation creates a set of |
1342 |
|
|
conjugate angle coordinates to the body-fixed angular momenta |
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 |
gezelter |
1425 |
\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 |
mmeineke |
1121 |
\end{eqnarray} |
1353 |
|
|
which utilize the body-fixed torques, ${\bf \tau}^b$. Torques are |
1354 |
|
|
most easily derived in the space-fixed frame, |
1355 |
|
|
\begin{equation} |
1356 |
|
|
{\bf \tau}^b(t) = \mathsf{A}(t) \cdot {\bf \tau}^s(t), |
1357 |
|
|
\end{equation} |
1358 |
|
|
where the torques are either derived from the forces on the |
1359 |
|
|
constituent atoms of the rigid body, or for directional atoms, |
1360 |
|
|
directly from derivatives of the potential energy, |
1361 |
|
|
\begin{equation} |
1362 |
|
|
{\bf \tau}^s(t) = - \hat{\bf u}(t) \times \left( \frac{\partial} |
1363 |
|
|
{\partial \hat{\bf u}} V\left(\left\{ {\bf r}(t) \right\}, \left\{ |
1364 |
|
|
\mathsf{A}(t) \right\}\right) \right). |
1365 |
|
|
\end{equation} |
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 |
gezelter |
1439 |
The {\sc dlm} method uses a Trotter factorization of the orientational |
1370 |
mmeineke |
1121 |
propagator. This has three effects: |
1371 |
|
|
\begin{enumerate} |
1372 |
|
|
\item the integrator is area-preserving in phase space (i.e. it is |
1373 |
|
|
{\it symplectic}), |
1374 |
|
|
\item the integrator is time-{\it reversible}, making it suitable for Hybrid |
1375 |
|
|
Monte Carlo applications, and |
1376 |
|
|
\item the error for a single time step is of order $\mathcal{O}\left(h^4\right)$ |
1377 |
|
|
for timesteps of length $h$. |
1378 |
|
|
\end{enumerate} |
1379 |
|
|
|
1380 |
|
|
The integration of the equations of motion is carried out in a |
1381 |
|
|
velocity-Verlet style 2-part algorithm, where $h= \delta t$: |
1382 |
|
|
|
1383 |
|
|
{\tt moveA:} |
1384 |
|
|
\begin{align*} |
1385 |
|
|
{\bf v}\left(t + h / 2\right) &\leftarrow {\bf v}(t) |
1386 |
|
|
+ \frac{h}{2} \left( {\bf f}(t) / m \right), \\ |
1387 |
|
|
% |
1388 |
|
|
{\bf r}(t + h) &\leftarrow {\bf r}(t) |
1389 |
|
|
+ h {\bf v}\left(t + h / 2 \right), \\ |
1390 |
|
|
% |
1391 |
|
|
{\bf j}\left(t + h / 2 \right) &\leftarrow {\bf j}(t) |
1392 |
|
|
+ \frac{h}{2} {\bf \tau}^b(t), \\ |
1393 |
|
|
% |
1394 |
|
|
\mathsf{A}(t + h) &\leftarrow \mathrm{rotate}\left( h {\bf j} |
1395 |
|
|
(t + h / 2) \cdot \overleftrightarrow{\mathsf{I}}^{-1} \right). |
1396 |
|
|
\end{align*} |
1397 |
|
|
|
1398 |
|
|
In this context, the $\mathrm{rotate}$ function is the reversible product |
1399 |
|
|
of the three body-fixed rotations, |
1400 |
|
|
\begin{equation} |
1401 |
|
|
\mathrm{rotate}({\bf a}) = \mathsf{G}_x(a_x / 2) \cdot |
1402 |
|
|
\mathsf{G}_y(a_y / 2) \cdot \mathsf{G}_z(a_z) \cdot \mathsf{G}_y(a_y / |
1403 |
|
|
2) \cdot \mathsf{G}_x(a_x /2), |
1404 |
|
|
\end{equation} |
1405 |
|
|
where each rotational propagator, $\mathsf{G}_\alpha(\theta)$, rotates |
1406 |
|
|
both the rotation matrix ($\mathsf{A}$) and the body-fixed angular |
1407 |
|
|
momentum (${\bf j}$) by an angle $\theta$ around body-fixed axis |
1408 |
|
|
$\alpha$, |
1409 |
|
|
\begin{equation} |
1410 |
|
|
\mathsf{G}_\alpha( \theta ) = \left\{ |
1411 |
|
|
\begin{array}{lcl} |
1412 |
|
|
\mathsf{A}(t) & \leftarrow & \mathsf{A}(0) \cdot \mathsf{R}_\alpha(\theta)^T, \\ |
1413 |
|
|
{\bf j}(t) & \leftarrow & \mathsf{R}_\alpha(\theta) \cdot {\bf j}(0). |
1414 |
|
|
\end{array} |
1415 |
|
|
\right. |
1416 |
|
|
\end{equation} |
1417 |
|
|
$\mathsf{R}_\alpha$ is a quadratic approximation to |
1418 |
|
|
the single-axis rotation matrix. For example, in the small-angle |
1419 |
|
|
limit, the rotation matrix around the body-fixed x-axis can be |
1420 |
|
|
approximated as |
1421 |
|
|
\begin{equation} |
1422 |
|
|
\mathsf{R}_x(\theta) \approx \left( |
1423 |
|
|
\begin{array}{ccc} |
1424 |
|
|
1 & 0 & 0 \\ |
1425 |
|
|
0 & \frac{1-\theta^2 / 4}{1 + \theta^2 / 4} & -\frac{\theta}{1+ |
1426 |
|
|
\theta^2 / 4} \\ |
1427 |
|
|
0 & \frac{\theta}{1+ |
1428 |
|
|
\theta^2 / 4} & \frac{1-\theta^2 / 4}{1 + \theta^2 / 4} |
1429 |
|
|
\end{array} |
1430 |
|
|
\right). |
1431 |
|
|
\end{equation} |
1432 |
|
|
All other rotations follow in a straightforward manner. |
1433 |
|
|
|
1434 |
|
|
After the first part of the propagation, the forces and body-fixed |
1435 |
|
|
torques are calculated at the new positions and orientations |
1436 |
|
|
|
1437 |
|
|
{\tt doForces:} |
1438 |
|
|
\begin{align*} |
1439 |
|
|
{\bf f}(t + h) &\leftarrow |
1440 |
|
|
- \left(\frac{\partial V}{\partial {\bf r}}\right)_{{\bf r}(t + h)}, \\ |
1441 |
|
|
% |
1442 |
|
|
{\bf \tau}^{s}(t + h) &\leftarrow {\bf u}(t + h) |
1443 |
|
|
\times \frac{\partial V}{\partial {\bf u}}, \\ |
1444 |
|
|
% |
1445 |
|
|
{\bf \tau}^{b}(t + h) &\leftarrow \mathsf{A}(t + h) |
1446 |
|
|
\cdot {\bf \tau}^s(t + h). |
1447 |
|
|
\end{align*} |
1448 |
|
|
|
1449 |
|
|
{\sc oopse} automatically updates ${\bf u}$ when the rotation matrix |
1450 |
|
|
$\mathsf{A}$ is calculated in {\tt moveA}. Once the forces and |
1451 |
|
|
torques have been obtained at the new time step, the velocities can be |
1452 |
|
|
advanced to the same time value. |
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*} |
1462 |
|
|
|
1463 |
gezelter |
1439 |
The matrix rotations used in the {\sc dlm} method end up being more |
1464 |
|
|
costly computationally than the simpler arithmetic quaternion |
1465 |
chrisfen |
1434 |
propagation. With the same time step, a 1024-molecule water simulation |
1466 |
gezelter |
1439 |
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 |
mmeineke |
1121 |
|
1472 |
|
|
\begin{figure} |
1473 |
|
|
\centering |
1474 |
chrisfen |
1434 |
\includegraphics[width=\linewidth]{quatvsdlm.eps} |
1475 |
|
|
\caption[Energy conservation analysis of the {\sc dlm} and quaternion |
1476 |
gezelter |
1439 |
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 |
chrisfen |
1434 |
\label{quatdlm} |
1490 |
mmeineke |
1121 |
\end{figure} |
1491 |
|
|
|
1492 |
gezelter |
1439 |
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 |
mmeineke |
1121 |
|
1504 |
gezelter |
1439 |
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 |
chrisfen |
1434 |
|
1520 |
|
|
\begin{figure} |
1521 |
|
|
\centering |
1522 |
|
|
\includegraphics[width=\linewidth]{compCost.eps} |
1523 |
|
|
\caption[Energy drift as a function of required simulation run |
1524 |
gezelter |
1439 |
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 |
chrisfen |
1434 |
\label{cpuCost} |
1533 |
|
|
\end{figure} |
1534 |
|
|
|
1535 |
mmeineke |
1121 |
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 |
gezelter |
1425 |
{\bf variable} & {\bf Meta-data keyword} & {\bf units} & {\bf |
1541 |
mmeineke |
1121 |
default value} \\ |
1542 |
|
|
$h$ & {\tt dt = 2.0;} & fs & none |
1543 |
|
|
\end{tabular} |
1544 |
|
|
\end{center} |
1545 |
|
|
|
1546 |
|
|
\subsection{\label{sec:extended}Extended Systems for other Ensembles} |
1547 |
|
|
|
1548 |
|
|
{\sc oopse} implements a number of extended system integrators for |
1549 |
|
|
sampling from other ensembles relevant to chemical physics. The |
1550 |
gezelter |
1425 |
integrator can be selected with the {\tt ensemble} keyword in the |
1551 |
|
|
meta-data file: |
1552 |
mmeineke |
1121 |
|
1553 |
|
|
\begin{center} |
1554 |
|
|
\begin{tabular}{lll} |
1555 |
gezelter |
1425 |
{\bf Integrator} & {\bf Ensemble} & {\bf Meta-data instruction} \\ |
1556 |
mmeineke |
1121 |
NVE & microcanonical & {\tt ensemble = NVE; } \\ |
1557 |
|
|
NVT & canonical & {\tt ensemble = NVT; } \\ |
1558 |
|
|
NPTi & isobaric-isothermal & {\tt ensemble = NPTi;} \\ |
1559 |
|
|
& (with isotropic volume changes) & \\ |
1560 |
|
|
NPTf & isobaric-isothermal & {\tt ensemble = NPTf;} \\ |
1561 |
|
|
& (with changes to box shape) & \\ |
1562 |
|
|
NPTxyz & approximate isobaric-isothermal & {\tt ensemble = NPTxyz;} \\ |
1563 |
|
|
& (with separate barostats on each box dimension) & \\ |
1564 |
|
|
\end{tabular} |
1565 |
|
|
\end{center} |
1566 |
|
|
|
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 |
gezelter |
1425 |
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 |
mmeineke |
1121 |
within a constant, the Helmholtz free energy.\cite{melchionna93} |
1573 |
|
|
|
1574 |
|
|
NPT algorithms attempt to maintain constant pressure in the system by |
1575 |
|
|
coupling the volume of the system to a barostat. {\sc oopse} contains |
1576 |
|
|
three different constant pressure algorithms. The first two, NPTi and |
1577 |
|
|
NPTf have been shown to conserve a quantity that is, to within a |
1578 |
|
|
constant, the Gibbs free energy.\cite{melchionna93} The Melchionna |
1579 |
|
|
modification to the Hoover barostat is implemented in both NPTi and |
1580 |
|
|
NPTf. NPTi allows only isotropic changes in the simulation box, while |
1581 |
|
|
box {\it shape} variations are allowed in NPTf. The NPTxyz integrator |
1582 |
|
|
has {\it not} been shown to sample from the isobaric-isothermal |
1583 |
|
|
ensemble. It is useful, however, in that it maintains orthogonality |
1584 |
|
|
for the axes of the simulation box while attempting to equalize |
1585 |
|
|
pressure along the three perpendicular directions in the box. |
1586 |
|
|
|
1587 |
|
|
Each of the extended system integrators requires additional keywords |
1588 |
|
|
to set target values for the thermodynamic state variables that are |
1589 |
|
|
being held constant. Keywords are also required to set the |
1590 |
|
|
characteristic decay times for the dynamics of the extended |
1591 |
|
|
variables. |
1592 |
|
|
|
1593 |
|
|
\begin{center} |
1594 |
|
|
\begin{tabular}{llll} |
1595 |
gezelter |
1425 |
{\bf variable} & {\bf Meta-data instruction} & {\bf units} & {\bf |
1596 |
mmeineke |
1121 |
default value} \\ |
1597 |
|
|
$T_{\mathrm{target}}$ & {\tt targetTemperature = 300;} & K & none \\ |
1598 |
|
|
$P_{\mathrm{target}}$ & {\tt targetPressure = 1;} & atm & none \\ |
1599 |
|
|
$\tau_T$ & {\tt tauThermostat = 1e3;} & fs & none \\ |
1600 |
|
|
$\tau_B$ & {\tt tauBarostat = 5e3;} & fs & none \\ |
1601 |
|
|
& {\tt resetTime = 200;} & fs & none \\ |
1602 |
|
|
& {\tt useInitialExtendedSystemState = true;} & logical & |
1603 |
|
|
true |
1604 |
|
|
\end{tabular} |
1605 |
|
|
\end{center} |
1606 |
|
|
|
1607 |
|
|
Two additional keywords can be used to either clear the extended |
1608 |
|
|
system variables periodically ({\tt resetTime}), or to maintain the |
1609 |
|
|
state of the extended system variables between simulations ({\tt |
1610 |
|
|
useInitialExtendedSystemState}). More details on these variables |
1611 |
|
|
and their use in the integrators follows below. |
1612 |
|
|
|
1613 |
|
|
\subsection{\label{oopseSec:noseHooverThermo}Nos\'{e}-Hoover Thermostatting} |
1614 |
|
|
|
1615 |
|
|
The Nos\'e-Hoover equations of motion are given by\cite{Hoover85} |
1616 |
|
|
\begin{eqnarray} |
1617 |
|
|
\dot{{\bf r}} & = & {\bf v}, \\ |
1618 |
|
|
\dot{{\bf v}} & = & \frac{{\bf f}}{m} - \chi {\bf v} ,\\ |
1619 |
|
|
\dot{\mathsf{A}} & = & \mathsf{A} \cdot |
1620 |
|
|
\mbox{ skew}\left(\overleftrightarrow{\mathsf{I}}^{-1} \cdot {\bf j}\right), \\ |
1621 |
|
|
\dot{{\bf j}} & = & {\bf j} \times \left( \overleftrightarrow{\mathsf{I}}^{-1} |
1622 |
|
|
\cdot {\bf j} \right) - \mbox{ rot}\left(\mathsf{A}^{T} \cdot \frac{\partial |
1623 |
|
|
V}{\partial \mathsf{A}} \right) - \chi {\bf j}. |
1624 |
|
|
\label{eq:nosehoovereom} |
1625 |
|
|
\end{eqnarray} |
1626 |
|
|
|
1627 |
|
|
$\chi$ is an ``extra'' variable included in the extended system, and |
1628 |
|
|
it is propagated using the first order equation of motion |
1629 |
|
|
\begin{equation} |
1630 |
|
|
\dot{\chi} = \frac{1}{\tau_{T}^2} \left( \frac{T}{T_{\mathrm{target}}} - 1 \right). |
1631 |
|
|
\label{eq:nosehooverext} |
1632 |
|
|
\end{equation} |
1633 |
|
|
|
1634 |
|
|
The instantaneous temperature $T$ is proportional to the total kinetic |
1635 |
|
|
energy (both translational and orientational) and is given by |
1636 |
|
|
\begin{equation} |
1637 |
|
|
T = \frac{2 K}{f k_B} |
1638 |
|
|
\end{equation} |
1639 |
|
|
Here, $f$ is the total number of degrees of freedom in the system, |
1640 |
|
|
\begin{equation} |
1641 |
gezelter |
1439 |
f = 3 N + 2 N_{\mathrm{linear}} + 3 N_{\mathrm{non-linear}} - N_{\mathrm{constraints}}, |
1642 |
mmeineke |
1121 |
\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 |
gezelter |
1439 |
\sum_{i=1}^{N_{\mathrm{linear}}+N_{\mathrm{non-linear}}} \frac{1}{2} {\bf j}_i^T \cdot |
1647 |
mmeineke |
1121 |
\overleftrightarrow{\mathsf{I}}_i^{-1} \cdot {\bf j}_i. |
1648 |
|
|
\end{equation} |
1649 |
gezelter |
1439 |
$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 |
mmeineke |
1121 |
|
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 |
gezelter |
1425 |
{\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 |
mmeineke |
1121 |
part algorithm: |
1662 |
|
|
|
1663 |
|
|
{\tt moveA:} |
1664 |
|
|
\begin{align*} |
1665 |
|
|
T(t) &\leftarrow \left\{{\bf v}(t)\right\}, \left\{{\bf j}(t)\right\} ,\\ |
1666 |
|
|
% |
1667 |
|
|
{\bf v}\left(t + h / 2\right) &\leftarrow {\bf v}(t) |
1668 |
|
|
+ \frac{h}{2} \left( \frac{{\bf f}(t)}{m} - {\bf v}(t) |
1669 |
|
|
\chi(t)\right), \\ |
1670 |
|
|
% |
1671 |
|
|
{\bf r}(t + h) &\leftarrow {\bf r}(t) |
1672 |
|
|
+ h {\bf v}\left(t + h / 2 \right) ,\\ |
1673 |
|
|
% |
1674 |
|
|
{\bf j}\left(t + h / 2 \right) &\leftarrow {\bf j}(t) |
1675 |
|
|
+ \frac{h}{2} \left( {\bf \tau}^b(t) - {\bf j}(t) |
1676 |
|
|
\chi(t) \right) ,\\ |
1677 |
|
|
% |
1678 |
|
|
\mathsf{A}(t + h) &\leftarrow \mathrm{rotate} |
1679 |
|
|
\left(h * {\bf j}(t + h / 2) |
1680 |
|
|
\overleftrightarrow{\mathsf{I}}^{-1} \right) ,\\ |
1681 |
|
|
% |
1682 |
|
|
\chi\left(t + h / 2 \right) &\leftarrow \chi(t) |
1683 |
|
|
+ \frac{h}{2 \tau_T^2} \left( \frac{T(t)} |
1684 |
|
|
{T_{\mathrm{target}}} - 1 \right) . |
1685 |
|
|
\end{align*} |
1686 |
|
|
|
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 |
gezelter |
1439 |
the section on the {\sc dlm} integrator. Note that this operation modifies |
1691 |
mmeineke |
1121 |
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 |
gezelter |
1439 |
doForces} portion of the {\sc dlm} integrator. |
1697 |
mmeineke |
1121 |
|
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 |
1700 |
|
|
advanced to the same time value. |
1701 |
|
|
|
1702 |
|
|
{\tt moveB:} |
1703 |
|
|
\begin{align*} |
1704 |
|
|
T(t + h) &\leftarrow \left\{{\bf v}(t + h)\right\}, |
1705 |
|
|
\left\{{\bf j}(t + h)\right\}, \\ |
1706 |
|
|
% |
1707 |
|
|
\chi\left(t + h \right) &\leftarrow \chi\left(t + h / |
1708 |
|
|
2 \right) + \frac{h}{2 \tau_T^2} \left( \frac{T(t+h)} |
1709 |
|
|
{T_{\mathrm{target}}} - 1 \right), \\ |
1710 |
|
|
% |
1711 |
|
|
{\bf v}\left(t + h \right) &\leftarrow {\bf v}\left(t |
1712 |
|
|
+ h / 2 \right) + \frac{h}{2} \left( |
1713 |
|
|
\frac{{\bf f}(t + h)}{m} - {\bf v}(t + h) |
1714 |
|
|
\chi(t h)\right) ,\\ |
1715 |
|
|
% |
1716 |
|
|
{\bf j}\left(t + h \right) &\leftarrow {\bf j}\left(t |
1717 |
|
|
+ h / 2 \right) + \frac{h}{2} |
1718 |
|
|
\left( {\bf \tau}^b(t + h) - {\bf j}(t + h) |
1719 |
|
|
\chi(t + h) \right) . |
1720 |
|
|
\end{align*} |
1721 |
|
|
|
1722 |
gezelter |
1425 |
Since ${\bf v}(t + h)$ and ${\bf j}(t + h)$ are required to calculate |
1723 |
mmeineke |
1121 |
$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 |
1726 |
|
|
relative tolerance for the self-consistency check defaults to a value |
1727 |
|
|
of $\mbox{10}^{-6}$, but {\sc oopse} will terminate the iteration |
1728 |
|
|
after 4 loops even if the consistency check has not been satisfied. |
1729 |
|
|
|
1730 |
|
|
The Nos\'e-Hoover algorithm is known to conserve a Hamiltonian for the |
1731 |
|
|
extended system that is, to within a constant, identical to the |
1732 |
|
|
Helmholtz free energy,\cite{melchionna93} |
1733 |
|
|
\begin{equation} |
1734 |
|
|
H_{\mathrm{NVT}} = V + K + f k_B T_{\mathrm{target}} \left( |
1735 |
|
|
\frac{\tau_{T}^2 \chi^2(t)}{2} + \int_{0}^{t} \chi(t^\prime) dt^\prime |
1736 |
|
|
\right). |
1737 |
|
|
\end{equation} |
1738 |
|
|
Poor choices of $h$ or $\tau_T$ can result in non-conservation |
1739 |
|
|
of $H_{\mathrm{NVT}}$, so the conserved quantity is maintained in the |
1740 |
|
|
last column of the {\tt .stat} file to allow checks on the quality of |
1741 |
|
|
the integration. |
1742 |
|
|
|
1743 |
|
|
Bond constraints are applied at the end of both the {\tt moveA} and |
1744 |
|
|
{\tt moveB} portions of the algorithm. Details on the constraint |
1745 |
|
|
algorithms are given in section \ref{oopseSec:rattle}. |
1746 |
|
|
|
1747 |
|
|
\subsection{\label{sec:NPTi}Constant-pressure integration with |
1748 |
|
|
isotropic box deformations (NPTi)} |
1749 |
|
|
|
1750 |
gezelter |
1425 |
To carry out isobaric-isothermal ensemble calculations, {\sc oopse} |
1751 |
mmeineke |
1121 |
implements the Melchionna modifications to the Nos\'e-Hoover-Andersen |
1752 |
gezelter |
1425 |
equations of motion.\cite{melchionna93} The equations of motion are |
1753 |
|
|
the same as NVT with the following exceptions: |
1754 |
mmeineke |
1121 |
|
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}, \\ |
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 . |
1761 |
|
|
\label{eq:melchionna1} |
1762 |
|
|
\end{eqnarray} |
1763 |
|
|
|
1764 |
|
|
$\chi$ and $\eta$ are the ``extra'' degrees of freedom in the extended |
1765 |
|
|
system. $\chi$ is a thermostat, and it has the same function as it |
1766 |
|
|
does in the Nos\'e-Hoover NVT integrator. $\eta$ is a barostat which |
1767 |
|
|
controls changes to the volume of the simulation box. ${\bf R}_0$ is |
1768 |
|
|
the location of the center of mass for the entire system, and |
1769 |
|
|
$\mathcal{V}$ is the volume of the simulation box. At any time, the |
1770 |
|
|
volume can be calculated from the determinant of the matrix which |
1771 |
|
|
describes the box shape: |
1772 |
|
|
\begin{equation} |
1773 |
|
|
\mathcal{V} = \det(\mathsf{H}). |
1774 |
|
|
\end{equation} |
1775 |
|
|
|
1776 |
|
|
The NPTi integrator requires an instantaneous pressure. This quantity |
1777 |
|
|
is calculated via the pressure tensor, |
1778 |
|
|
\begin{equation} |
1779 |
|
|
\overleftrightarrow{\mathsf{P}}(t) = \frac{1}{\mathcal{V}(t)} \left( |
1780 |
|
|
\sum_{i=1}^{N} m_i {\bf v}_i(t) \otimes {\bf v}_i(t) \right) + |
1781 |
|
|
\overleftrightarrow{\mathsf{W}}(t). |
1782 |
|
|
\end{equation} |
1783 |
|
|
The kinetic contribution to the pressure tensor utilizes the {\it |
1784 |
gezelter |
1425 |
outer} product of the velocities, denoted by the $\otimes$ symbol. The |
1785 |
mmeineke |
1121 |
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, |
1788 |
|
|
\begin{equation} |
1789 |
|
|
\overleftrightarrow{\mathsf{W}}(t) = \sum_{i} \sum_{j>i} {\bf r}_{ij}(t) |
1790 |
|
|
\otimes {\bf f}_{ij}(t). |
1791 |
|
|
\end{equation} |
1792 |
gezelter |
1425 |
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 |
gezelter |
1439 |
s^{\prime}(r_{ab}) \frac{{\bf r}_{ab}}{|r_{ab}|} \sum_{i \in a} \sum_{j |
1802 |
gezelter |
1425 |
\in b} V_{ij}({\bf r}_{ij}). |
1803 |
|
|
\end{equation} |
1804 |
|
|
|
1805 |
mmeineke |
1121 |
The instantaneous pressure is then simply obtained from the trace of |
1806 |
gezelter |
1425 |
the pressure tensor, |
1807 |
mmeineke |
1121 |
\begin{equation} |
1808 |
gezelter |
1439 |
P(t) = \frac{1}{3} \mathrm{Tr} \left( \overleftrightarrow{\mathsf{P}}(t) |
1809 |
|
|
\right). |
1810 |
mmeineke |
1121 |
\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 |
gezelter |
1425 |
{\tt tauBarostat} and {\tt targetPressure} keywords in the meta-data |
1816 |
mmeineke |
1121 |
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 |
gezelter |
1439 |
velocity-Verlet style two part algorithm with only the following |
1820 |
|
|
differences: |
1821 |
mmeineke |
1121 |
|
1822 |
|
|
{\tt moveA:} |
1823 |
|
|
\begin{align*} |
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), \\ |
1829 |
|
|
% |
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), \\ |
1833 |
|
|
% |
1834 |
|
|
{\bf r}(t + h) &\leftarrow {\bf r}(t) + h |
1835 |
|
|
\left\{ {\bf v}\left(t + h / 2 \right) |
1836 |
|
|
+ \eta(t + h / 2)\left[ {\bf r}(t + h) |
1837 |
|
|
- {\bf R}_0 \right] \right\} ,\\ |
1838 |
|
|
% |
1839 |
|
|
\mathsf{H}(t + h) &\leftarrow e^{-h \eta(t + h / 2)} |
1840 |
|
|
\mathsf{H}(t). |
1841 |
|
|
\end{align*} |
1842 |
|
|
|
1843 |
mmeineke |
1179 |
The propagation of positions to time $t + h$ |
1844 |
mmeineke |
1121 |
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 |
1847 |
|
|
one full time step by an exponential factor that depends on the value |
1848 |
|
|
of $\eta$ at time $t + |
1849 |
|
|
h / 2$. Reshaping the box uniformly also scales the volume of |
1850 |
|
|
the box by |
1851 |
|
|
\begin{equation} |
1852 |
gezelter |
1425 |
\mathcal{V}(t + h) \leftarrow e^{ - 3 h \eta(t + h /2)} \times |
1853 |
gezelter |
1439 |
\mathcal{V}(t). |
1854 |
mmeineke |
1121 |
\end{equation} |
1855 |
|
|
|
1856 |
|
|
The {\tt doForces} step for the NPTi integrator is exactly the same as |
1857 |
gezelter |
1439 |
in both the {\sc dlm} and NVT integrators. Once the forces and torques have |
1858 |
mmeineke |
1121 |
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*} |
1863 |
|
|
P(t + h) &\leftarrow \left\{{\bf r}(t + h)\right\}, |
1864 |
|
|
\left\{{\bf v}(t + h)\right\}, \\ |
1865 |
|
|
% |
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), \\ |
1869 |
|
|
% |
1870 |
|
|
{\bf v}\left(t + h \right) &\leftarrow {\bf v}\left(t |
1871 |
|
|
+ h / 2 \right) + \frac{h}{2} \left( |
1872 |
|
|
\frac{{\bf f}(t + h)}{m} - {\bf v}(t + h) |
1873 |
|
|
(\chi(t + h) + \eta(t + h)) \right) ,\\ |
1874 |
|
|
% |
1875 |
|
|
{\bf j}\left(t + h \right) &\leftarrow {\bf j}\left(t |
1876 |
|
|
+ h / 2 \right) + \frac{h}{2} \left( {\bf |
1877 |
|
|
\tau}^b(t + h) - {\bf j}(t + h) |
1878 |
|
|
\chi(t + h) \right) . |
1879 |
|
|
\end{align*} |
1880 |
|
|
|
1881 |
|
|
Once again, since ${\bf v}(t + h)$ and ${\bf j}(t + h)$ are required |
1882 |
gezelter |
1425 |
to calculate $T(t + h)$, $P(t + h)$, $\chi(t + h)$, and $\eta(t + |
1883 |
mmeineke |
1121 |
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 |
1886 |
|
|
the self-consistency check defaults to a value of $\mbox{10}^{-6}$, |
1887 |
|
|
but {\sc oopse} will terminate the iteration after 4 loops even if the |
1888 |
|
|
consistency check has not been satisfied. |
1889 |
|
|
|
1890 |
|
|
The Melchionna modification of the Nos\'e-Hoover-Andersen algorithm is |
1891 |
|
|
known to conserve a Hamiltonian for the extended system that is, to |
1892 |
|
|
within a constant, identical to the Gibbs free energy, |
1893 |
|
|
\begin{equation} |
1894 |
|
|
H_{\mathrm{NPTi}} = V + K + f k_B T_{\mathrm{target}} \left( |
1895 |
|
|
\frac{\tau_{T}^2 \chi^2(t)}{2} + \int_{0}^{t} \chi(t^\prime) dt^\prime |
1896 |
|
|
\right) + P_{\mathrm{target}} \mathcal{V}(t). |
1897 |
|
|
\end{equation} |
1898 |
|
|
Poor choices of $\delta t$, $\tau_T$, or $\tau_B$ can result in |
1899 |
|
|
non-conservation of $H_{\mathrm{NPTi}}$, so the conserved quantity is |
1900 |
|
|
maintained in the last column of the {\tt .stat} file to allow checks |
1901 |
|
|
on the quality of the integration. It is also known that this |
1902 |
|
|
algorithm samples the equilibrium distribution for the enthalpy |
1903 |
|
|
(including contributions for the thermostat and barostat), |
1904 |
|
|
\begin{equation} |
1905 |
|
|
H_{\mathrm{NPTi}} = V + K + \frac{f k_B T_{\mathrm{target}}}{2} \left( |
1906 |
|
|
\chi^2 \tau_T^2 + \eta^2 \tau_B^2 \right) + P_{\mathrm{target}} |
1907 |
|
|
\mathcal{V}(t). |
1908 |
|
|
\end{equation} |
1909 |
|
|
|
1910 |
|
|
Bond constraints are applied at the end of both the {\tt moveA} and |
1911 |
|
|
{\tt moveB} portions of the algorithm. Details on the constraint |
1912 |
|
|
algorithms are given in section \ref{oopseSec:rattle}. |
1913 |
|
|
|
1914 |
|
|
\subsection{\label{sec:NPTf}Constant-pressure integration with a |
1915 |
|
|
flexible box (NPTf)} |
1916 |
|
|
|
1917 |
|
|
There is a relatively simple generalization of the |
1918 |
|
|
Nos\'e-Hoover-Andersen method to include changes in the simulation box |
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 |
gezelter |
1425 |
the box shape. The equations of motion for this method differ from |
1923 |
|
|
those of NPTi as follows: |
1924 |
mmeineke |
1121 |
\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}, \\ |
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} . |
1931 |
|
|
\label{eq:melchionna2} |
1932 |
|
|
\end{eqnarray} |
1933 |
|
|
|
1934 |
|
|
Here, $\mathsf{1}$ is the unit matrix and $\overleftrightarrow{\mathsf{P}}$ |
1935 |
|
|
is the pressure tensor. Again, the volume, $\mathcal{V} = \det |
1936 |
|
|
\mathsf{H}$. |
1937 |
|
|
|
1938 |
|
|
The propagation of the equations of motion is nearly identical to the |
1939 |
|
|
NPTi integration: |
1940 |
|
|
|
1941 |
|
|
{\tt moveA:} |
1942 |
|
|
\begin{align*} |
1943 |
|
|
\overleftrightarrow{\mathsf{P}}(t) &\leftarrow \left\{{\bf r}(t)\right\}, |
1944 |
|
|
\left\{{\bf v}(t)\right\} ,\\ |
1945 |
|
|
% |
1946 |
|
|
{\bf v}\left(t + h / 2\right) &\leftarrow {\bf v}(t) |
1947 |
|
|
+ \frac{h}{2} \left( \frac{{\bf f}(t)}{m} - |
1948 |
|
|
\left(\chi(t)\mathsf{1} + \overleftrightarrow{\eta}(t) \right) \cdot |
1949 |
|
|
{\bf v}(t) \right), \\ |
1950 |
|
|
% |
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) |
1954 |
|
|
- P_{\mathrm{target}}\mathsf{1} \right), \\ |
1955 |
|
|
% |
1956 |
|
|
{\bf r}(t + h) &\leftarrow {\bf r}(t) + h \left\{ {\bf v} |
1957 |
|
|
\left(t + h / 2 \right) + \overleftrightarrow{\eta}(t + |
1958 |
|
|
h / 2) \cdot \left[ {\bf r}(t + h) |
1959 |
|
|
- {\bf R}_0 \right] \right\}, \\ |
1960 |
|
|
% |
1961 |
|
|
\mathsf{H}(t + h) &\leftarrow \mathsf{H}(t) \cdot e^{-h |
1962 |
|
|
\overleftrightarrow{\eta}(t + h / 2)} . |
1963 |
|
|
\end{align*} |
1964 |
|
|
{\sc oopse} uses a power series expansion truncated at second order |
1965 |
|
|
for the exponential operation which scales the simulation box. |
1966 |
|
|
|
1967 |
|
|
The {\tt moveB} portion of the algorithm is largely unchanged from the |
1968 |
|
|
NPTi integrator: |
1969 |
|
|
|
1970 |
|
|
{\tt moveB:} |
1971 |
|
|
\begin{align*} |
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 |
|
|
% |
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) |
1979 |
|
|
\tau_B^2} \left( \overleftrightarrow{P}(t + h) |
1980 |
|
|
- P_{\mathrm{target}}\mathsf{1} \right) ,\\ |
1981 |
|
|
% |
1982 |
|
|
{\bf v}\left(t + h \right) &\leftarrow {\bf v}\left(t |
1983 |
|
|
+ h / 2 \right) + \frac{h}{2} \left( |
1984 |
|
|
\frac{{\bf f}(t + h)}{m} - |
1985 |
|
|
(\chi(t + h)\mathsf{1} + \overleftrightarrow{\eta}(t |
1986 |
|
|
+ h)) \right) \cdot {\bf v}(t + h), \\ |
1987 |
|
|
\end{align*} |
1988 |
|
|
|
1989 |
|
|
The iterative schemes for both {\tt moveA} and {\tt moveB} are |
1990 |
|
|
identical to those described for the NPTi integrator. |
1991 |
|
|
|
1992 |
|
|
The NPTf integrator is known to conserve the following Hamiltonian: |
1993 |
|
|
\begin{equation} |
1994 |
|
|
H_{\mathrm{NPTf}} = V + K + f k_B T_{\mathrm{target}} \left( |
1995 |
|
|
\frac{\tau_{T}^2 \chi^2(t)}{2} + \int_{0}^{t} \chi(t^\prime) dt^\prime |
1996 |
|
|
\right) + P_{\mathrm{target}} \mathcal{V}(t) + \frac{f k_B |
1997 |
|
|
T_{\mathrm{target}}}{2} |
1998 |
|
|
\mathrm{Tr}\left[\overleftrightarrow{\eta}(t)\right]^2 \tau_B^2. |
1999 |
|
|
\end{equation} |
2000 |
|
|
|
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 |
gezelter |
1425 |
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 |
mmeineke |
1121 |
|
2008 |
|
|
\subsection{\label{nptxyz}Constant pressure in 3 axes (NPTxyz)} |
2009 |
|
|
|
2010 |
|
|
There is one additional extended system integrator which is somewhat |
2011 |
|
|
simpler than the NPTf method described above. In this case, the three |
2012 |
|
|
axes have independent barostats which each attempt to preserve the |
2013 |
|
|
target pressure along the box walls perpendicular to that particular |
2014 |
|
|
axis. The lengths of the box axes are allowed to fluctuate |
2015 |
|
|
independently, but the angle between the box axes does not change. |
2016 |
|
|
The equations of motion are identical to those described above, but |
2017 |
|
|
only the {\it diagonal} elements of $\overleftrightarrow{\eta}$ are |
2018 |
|
|
computed. The off-diagonal elements are set to zero (even when the |
2019 |
|
|
pressure tensor has non-zero off-diagonal elements). |
2020 |
|
|
|
2021 |
|
|
It should be noted that the NPTxyz integrator is {\it not} known to |
2022 |
|
|
preserve any Hamiltonian of interest to the chemical physics |
2023 |
|
|
community. The integrator is extremely useful, however, in generating |
2024 |
|
|
initial conditions for other integration methods. It {\it is} suitable |
2025 |
|
|
for use with liquid simulations, or in cases where there is |
2026 |
|
|
orientational anisotropy in the system (i.e. in lipid bilayer |
2027 |
|
|
simulations). |
2028 |
|
|
|
2029 |
mmeineke |
1134 |
\subsection{\label{sec:constraints}Constraint Methods} |
2030 |
|
|
|
2031 |
|
|
\subsubsection{\label{oopseSec:rattle}The {\sc rattle} Method for Bond |
2032 |
mmeineke |
1121 |
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 |
gezelter |
1425 |
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 |
mmeineke |
1121 |
|
2043 |
gezelter |
1425 |
\subsubsection{\label{oopseSec:zcons}The Z-Constraint Method} |
2044 |
mmeineke |
1121 |
|
2045 |
gezelter |
1425 |
A force auto-correlation method based on the fluctuation-dissipation |
2046 |
|
|
theorem was developed by Roux and Karplus to investigate the dynamics |
2047 |
mmeineke |
1121 |
of ions inside ion channels.\cite{Roux91} The time-dependent friction |
2048 |
|
|
coefficient can be calculated from the deviation of the instantaneous |
2049 |
gezelter |
1425 |
force from its mean value: |
2050 |
mmeineke |
1121 |
\begin{equation} |
2051 |
|
|
\xi(z,t)=\langle\delta F(z,t)\delta F(z,0)\rangle/k_{B}T, |
2052 |
|
|
\end{equation} |
2053 |
|
|
where% |
2054 |
|
|
\begin{equation} |
2055 |
|
|
\delta F(z,t)=F(z,t)-\langle F(z,t)\rangle. |
2056 |
|
|
\end{equation} |
2057 |
|
|
|
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 |
gezelter |
1425 |
|
2064 |
|
|
This allows the diffusion constant to then be calculated through the |
2065 |
mmeineke |
1121 |
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 |
gezelter |
1425 |
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 |
mmeineke |
1121 |
|
2082 |
gezelter |
1439 |
After the force calculation, the total force on molecule $\alpha$ is: |
2083 |
mmeineke |
1121 |
\begin{equation} |
2084 |
|
|
G_{\alpha} = \sum_i F_{\alpha i}, |
2085 |
|
|
\label{oopseEq:zc1} |
2086 |
|
|
\end{equation} |
2087 |
gezelter |
1425 |
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 |
mmeineke |
1121 |
\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 |
gezelter |
1425 |
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 |
mmeineke |
1121 |
\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 |
gezelter |
1439 |
where $v_{\alpha i}$ is the velocity of atom $i$ in the $z$ direction. |
2104 |
gezelter |
1425 |
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 |
mmeineke |
1121 |
\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 |
gezelter |
1425 |
where $\beta$ denotes all {\it unconstrained} molecules in the |
2112 |
mmeineke |
1121 |
system. Similarly, the velocities of the unconstrained molecules must |
2113 |
gezelter |
1425 |
also be scaled: |
2114 |
mmeineke |
1121 |
\begin{equation} |
2115 |
gezelter |
1425 |
v_{\beta i} = v_{\beta i} + \sum_{\alpha} \frac{\sum_i m_{\alpha i} |
2116 |
|
|
v_{\alpha i}}{\sum_i m_{\alpha i}}. |
2117 |
mmeineke |
1121 |
\end{equation} |
2118 |
|
|
|
2119 |
gezelter |
1425 |
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 |
mmeineke |
1121 |
\begin{equation} |
2127 |
|
|
U(t)=\frac{1}{2}k_{\text{Harmonic}}(z(t)-z_{\text{cons}})^{2},% |
2128 |
|
|
\end{equation} |
2129 |
gezelter |
1425 |
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 |
mmeineke |
1121 |
\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 |
gezelter |
1425 |
The user may also specify the use of a constant velocity method |
2140 |
|
|
(steered molecular dynamics) to move the molecules to their desired |
2141 |
gezelter |
1439 |
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 |
gezelter |
1425 |
|
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 |
|
|
|
2152 |
mmeineke |
1179 |
\begin{table} |
2153 |
gezelter |
1439 |
\caption{Meta-data Keywords: Z-Constraint Parameters} |
2154 |
mmeineke |
1179 |
\label{table:zconParams} |
2155 |
|
|
\begin{center} |
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 |
gezelter |
1439 |
{\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 |
mmeineke |
1179 |
|
2185 |
|
|
\end{tabularx} |
2186 |
|
|
\end{center} |
2187 |
|
|
\end{table} |
2188 |
|
|
|
2189 |
|
|
|
2190 |
gezelter |
1428 |
\section{\label{oopseSec:minimizer}Energy Minimization} |
2191 |
mmeineke |
1179 |
|
2192 |
gezelter |
1425 |
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 |
gezelter |
1439 |
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 |
mmeineke |
1179 |
|
2205 |
gezelter |
1425 |
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 |
gezelter |
1439 |
$x_{k+1}=x_{k}+$ $\lambda _{k}d_{k}$. In the steepest descent ({\sc |
2208 |
|
|
sd}) algorithm,% |
2209 |
mmeineke |
1179 |
\begin{equation} |
2210 |
gezelter |
1439 |
d_{k}=-\nabla V(x_{k}). |
2211 |
mmeineke |
1179 |
\end{equation} |
2212 |
gezelter |
1425 |
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 |
gezelter |
1439 |
\begin{equation} |
2218 |
|
|
d_{k+1} =-\nabla V(x_{k+1})+\gamma_{k}d_{k} |
2219 |
|
|
\end{equation} |
2220 |
|
|
where |
2221 |
|
|
\begin{equation} |
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} |
2225 |
mmeineke |
1179 |
|
2226 |
gezelter |
1425 |
The Polak-Ribiere variant~\cite{PolakRibiere} of the conjugate |
2227 |
|
|
gradient ($\gamma_{k}$) is defined as% |
2228 |
mmeineke |
1179 |
\begin{equation} |
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 |
gezelter |
1439 |
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 |
mmeineke |
1179 |
|
2237 |
gezelter |
1425 |
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 |
mmeineke |
1179 |
|
2248 |
|
|
\begin{table} |
2249 |
gezelter |
1439 |
\caption{Meta-data Keywords: Energy Minimizer Parameters} |
2250 |
mmeineke |
1179 |
\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 |
gezelter |
1425 |
{>{\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 |
mmeineke |
1179 |
|
2260 |
|
|
{\bf keyword} & {\bf units} & {\bf use} & {\bf remarks} \\ \hline |
2261 |
|
|
|
2262 |
gezelter |
1425 |
{\tt minimizer} & string & selects the minimization method to be used |
2263 |
|
|
& either {\tt CG} (conjugate gradient) or {\tt SD} (steepest |
2264 |
|
|
descent) \\ |
2265 |
gezelter |
1439 |
{\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 |
mmeineke |
1179 |
|
2281 |
|
|
\end{tabularx} |
2282 |
|
|
\end{center} |
2283 |
|
|
\end{table} |
2284 |
|
|
|
2285 |
gezelter |
1425 |
\section{\label{oopseSec:parallelization} Parallel Simulation Implementation} |
2286 |
mmeineke |
1179 |
|
2287 |
gezelter |
1425 |
Although processor power is continually improving, it is still |
2288 |
gezelter |
1439 |
unreasonable to simulate systems of more than 10,000 atoms on a single |
2289 |
gezelter |
1425 |
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 |
mmeineke |
1121 |
|
2296 |
gezelter |
1425 |
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 |
gezelter |
1439 |
optimal $\mathcal{O}(N/P)$ for atomic decomposition. Unfortunately, all |
2300 |
mmeineke |
1121 |
processors must communicate positions and forces with all other |
2301 |
gezelter |
1425 |
processors at every force evaluation, leading the communication costs |
2302 |
|
|
to scale as an unfavorable $\mathcal{O}(N)$, \emph{independent of the |
2303 |
mmeineke |
1121 |
number of processors}. This communication bottleneck led to the |
2304 |
gezelter |
1425 |
development of spatial and force decomposition methods, in which |
2305 |
mmeineke |
1121 |
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 |
gezelter |
1425 |
calculate forces on a given particle, a processor must simply know the |
2311 |
mmeineke |
1121 |
positions of particles within some cutoff radius located on nearby |
2312 |
gezelter |
1425 |
processors rather than the positions of particles on all |
2313 |
mmeineke |
1121 |
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 |
gezelter |
1425 |
is not very efficient for small $N$, since the overall communication |
2317 |
mmeineke |
1121 |
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 |
gezelter |
1439 |
decomposition method.\cite{hendrickson:95} Force decomposition assigns |
2322 |
|
|
particles to processors based on a block decomposition of the force |
2323 |
mmeineke |
1121 |
matrix. Processors are split into an optimally square grid forming row |
2324 |
|
|
and column processor groups. Forces are calculated on particles in a |
2325 |
gezelter |
1425 |
given row by particles located in that processor's column |
2326 |
gezelter |
1439 |
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 |
mmeineke |
1121 |
\section{\label{oopseSec:conclusion}Conclusion} |
2338 |
|
|
|
2339 |
gezelter |
1439 |
We have presented a new parallel simulation program called {\sc |
2340 |
gezelter |
1425 |
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 |
mmeineke |
1121 |
|
2352 |
gezelter |
1425 |
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 |
mmeineke |
1121 |
|
2359 |
|
|
\newpage |
2360 |
|
|
\section{Acknowledgments} |
2361 |
|
|
|
2362 |
gezelter |
1425 |
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 |
gezelter |
1439 |
\bibliographystyle{jcc} |
2369 |
mmeineke |
1121 |
\bibliography{oopsePaper} |
2370 |
|
|
|
2371 |
|
|
\end{document} |