1 |
mmeineke |
1045 |
\chapter{\label{chapt:oopse}OOPSE: AN OPEN SOURCE OBJECT-ORIENTED PARALLEL SIMULATION ENGINE FOR MOLECULAR DYNAMICS} |
2 |
mmeineke |
1044 |
|
3 |
|
|
|
4 |
|
|
|
5 |
mmeineke |
1045 |
%% \begin{abstract} |
6 |
|
|
%% We detail the capabilities of a new open-source parallel simulation |
7 |
|
|
%% package ({\sc oopse}) that can perform molecular dynamics simulations |
8 |
|
|
%% on atom types that are missing from other popular packages. In |
9 |
|
|
%% particular, {\sc oopse} is capable of performing orientational |
10 |
|
|
%% dynamics on dipolar systems, and it can handle simulations of metallic |
11 |
|
|
%% systems using the embedded atom method ({\sc eam}). |
12 |
|
|
%% \end{abstract} |
13 |
mmeineke |
1044 |
|
14 |
mmeineke |
1045 |
\lstset{language=C,frame=TB,basicstyle=\small,basicstyle=\ttfamily, % |
15 |
|
|
xleftmargin=0.5in, xrightmargin=0.5in,captionpos=b, % |
16 |
|
|
abovecaptionskip=0.5cm, belowcaptionskip=0.5cm} |
17 |
mmeineke |
1044 |
|
18 |
mmeineke |
1051 |
\section{\label{oopseSec:foreword}Foreword} |
19 |
mmeineke |
1044 |
|
20 |
mmeineke |
1051 |
In this chapter, I present and detail the capabilities of the open |
21 |
|
|
source simulation package {\sc oopse}. It is important to note, that a |
22 |
|
|
simulation package of this size and scope would not have been possible |
23 |
|
|
without the collaborative efforts of my colleagues: Charles |
24 |
|
|
F.~Vardeman II, Teng Lin, Christopher J.~Fennell and J.~Daniel |
25 |
mmeineke |
1068 |
Gezelter. Although my contributions to {\sc oopse} are major, |
26 |
|
|
consideration of my work apart from the others would not give a |
27 |
mmeineke |
1051 |
complete description to the package's capabilities. As such, all |
28 |
|
|
contributions to {\sc oopse} to date are presented in this chapter. |
29 |
mmeineke |
1044 |
|
30 |
mmeineke |
1068 |
Charles Vardeman is responsible for the parallelization of the long |
31 |
|
|
range forces in {\sc oopse} (Sec.~\ref{oopseSec:parallelization}) as |
32 |
|
|
well as the inclusion of the embedded-atom potential for transition |
33 |
|
|
metals (Sec.~\ref{oopseSec:eam}). Teng Lin's contributions include |
34 |
|
|
refinement of the periodic boundary conditions |
35 |
mmeineke |
1054 |
(Sec.~\ref{oopseSec:pbc}), the z-constraint method |
36 |
|
|
(Sec.~\ref{oopseSec:zcons}), refinement of the property analysis |
37 |
|
|
programs (Sec.~\ref{oopseSec:props}), and development in the extended |
38 |
mmeineke |
1068 |
system integrators (Sec.~\ref{oopseSec:noseHooverThermo}). Christopher |
39 |
mmeineke |
1054 |
Fennell worked on the symplectic integrator |
40 |
|
|
(Sec.~\ref{oopseSec:integrate}) and the refinement of the {\sc ssd} |
41 |
|
|
water model (Sec.~\ref{oopseSec:SSD}). Daniel Gezelter lent his |
42 |
|
|
talents in the development of the extended system integrators |
43 |
|
|
(Sec.~\ref{oopseSec:noseHooverThermo}) as well as giving general |
44 |
|
|
direction and oversight to the entire project. My responsibilities |
45 |
|
|
covered the creation and specification of {\sc bass} |
46 |
|
|
(Sec.~\ref{oopseSec:IOfiles}), the original development of the single |
47 |
|
|
processor version of {\sc oopse}, contributions to the extended state |
48 |
|
|
integrators (Sec.~\ref{oopseSec:noseHooverThermo}), the implementation |
49 |
|
|
of the Lennard-Jones (Sec.~\ref{sec:LJPot}) and {\sc duff} |
50 |
|
|
(Sec.~\ref{oopseSec:DUFF}) force fields, and initial implementation of |
51 |
|
|
the property analysis (Sec.~\ref{oopseSec:props}) and system |
52 |
mmeineke |
1068 |
initialization (Sec.~\ref{oopseSec:initCoords}) utility programs. {\sc |
53 |
|
|
oopse}, like many other Molecular Dynamics programs, is a work in |
54 |
|
|
progress, and will continue to be so for many graduate student |
55 |
|
|
lifetimes. |
56 |
mmeineke |
1044 |
|
57 |
mmeineke |
1051 |
\section{\label{sec:intro}Introduction} |
58 |
mmeineke |
1044 |
|
59 |
mmeineke |
1051 |
When choosing to simulate a chemical system with molecular dynamics, |
60 |
|
|
there are a variety of options available. For simple systems, one |
61 |
|
|
might consider writing one's own programming code. However, as systems |
62 |
|
|
grow larger and more complex, building and maintaining code for the |
63 |
|
|
simulations becomes a time consuming task. In such cases it is usually |
64 |
mmeineke |
1054 |
more convenient for a researcher to turn to pre-existing simulation |
65 |
mmeineke |
1051 |
packages. These packages, such as {\sc amber}\cite{pearlman:1995} and |
66 |
|
|
{\sc charmm}\cite{Brooks83}, provide powerful tools for researchers to |
67 |
|
|
conduct simulations of their systems without spending their time |
68 |
|
|
developing a code base to conduct their research. This then frees them |
69 |
mmeineke |
1054 |
to perhaps explore experimental analogues to their models. |
70 |
mmeineke |
1044 |
|
71 |
mmeineke |
1051 |
Despite their utility, problems with these packages arise when |
72 |
|
|
researchers try to develop techniques or energetic models that the |
73 |
mmeineke |
1068 |
code was not originally designed to simulate. Examples of uncommonly |
74 |
mmeineke |
1051 |
implemented techniques and energetics include; dipole-dipole |
75 |
mmeineke |
1054 |
interactions, rigid body dynamics, and metallic embedded |
76 |
mmeineke |
1051 |
potentials. When faced with these obstacles, a researcher must either |
77 |
|
|
develop their own code or license and extend one of the commercial |
78 |
|
|
packages. What we have elected to do, is develop a package of |
79 |
|
|
simulation code capable of implementing the types of models upon which |
80 |
|
|
our research is based. |
81 |
mmeineke |
1044 |
|
82 |
mmeineke |
1068 |
In developing {\sc oopse}, we have adhered to the precepts of Open |
83 |
|
|
Source development, and are releasing our source code with a |
84 |
|
|
permissive license. It is our intent that by doing so, other |
85 |
|
|
researchers might benefit from our work, and add their own |
86 |
|
|
contributions to the package. The license under which {\sc oopse} is |
87 |
|
|
distributed allows any researcher to download and modify the source |
88 |
|
|
code for their own use. In this way further development of {\sc oopse} |
89 |
|
|
is not limited to only the models of interest to ourselves, but also |
90 |
|
|
those of the community of scientists who contribute back to the |
91 |
|
|
project. |
92 |
mmeineke |
1044 |
|
93 |
mmeineke |
1054 |
We have structured this chapter to first discuss the empirical energy |
94 |
mmeineke |
1051 |
functions that {\sc oopse } implements in |
95 |
mmeineke |
1054 |
Sec.~\ref{oopseSec:empiricalEnergy}. Following that is a discussion of |
96 |
mmeineke |
1051 |
the various input and output files associated with the package |
97 |
mmeineke |
1068 |
(Sec.~\ref{oopseSec:IOfiles}). Sec.~\ref{oopseSec:mechanics} |
98 |
mmeineke |
1051 |
elucidates the various Molecular Dynamics algorithms {\sc oopse} |
99 |
mmeineke |
1054 |
implements in the integration of the Newtonian equations of |
100 |
mmeineke |
1051 |
motion. Basic analysis of the trajectories obtained from the |
101 |
|
|
simulation is discussed in Sec.~\ref{oopseSec:props}. Program design |
102 |
mmeineke |
1068 |
considerations are presented in Sec.~\ref{oopseSec:design}. And |
103 |
|
|
lastly, Sec.~\ref{oopseSec:conclusion} concludes the chapter. |
104 |
mmeineke |
1044 |
|
105 |
mmeineke |
1051 |
\section{\label{oopseSec:empiricalEnergy}The Empirical Energy Functions} |
106 |
|
|
|
107 |
|
|
\subsection{\label{oopseSec:atomsMolecules}Atoms, Molecules and Rigid Bodies} |
108 |
|
|
|
109 |
mmeineke |
1044 |
The basic unit of an {\sc oopse} simulation is the atom. The |
110 |
|
|
parameters describing the atom are generalized to make the atom as |
111 |
|
|
flexible a representation as possible. They may represent specific |
112 |
|
|
atoms of an element, or be used for collections of atoms such as |
113 |
|
|
methyl and carbonyl groups. The atoms are also capable of having |
114 |
|
|
directional components associated with them (\emph{e.g.}~permanent |
115 |
mmeineke |
1054 |
dipoles). Charges, permanent dipoles, and Lennard-Jones parameters for |
116 |
mmeineke |
1068 |
a given atom type are set in the force field parameter files. |
117 |
mmeineke |
1044 |
|
118 |
mmeineke |
1054 |
\begin{lstlisting}[float,caption={[Specifier for molecules and atoms] A sample specification of an Ar molecule},label=sch:AtmMole] |
119 |
mmeineke |
1044 |
molecule{ |
120 |
|
|
name = "Ar"; |
121 |
|
|
nAtoms = 1; |
122 |
|
|
atom[0]{ |
123 |
|
|
type="Ar"; |
124 |
|
|
position( 0.0, 0.0, 0.0 ); |
125 |
|
|
} |
126 |
|
|
} |
127 |
|
|
\end{lstlisting} |
128 |
|
|
|
129 |
mmeineke |
1045 |
|
130 |
mmeineke |
1054 |
Atoms can be collected into secondary structures such as rigid bodies |
131 |
mmeineke |
1044 |
or molecules. The molecule is a way for {\sc oopse} to keep track of |
132 |
|
|
the atoms in a simulation in logical manner. Molecular units store the |
133 |
mmeineke |
1054 |
identities of all the atoms and rigid bodies associated with |
134 |
|
|
themselves, and are responsible for the evaluation of their own |
135 |
|
|
internal interactions (\emph{i.e.}~bonds, bends, and torsions). Scheme |
136 |
mmeineke |
1068 |
\ref{sch:AtmMole} shows how one creates a molecule in a ``model'' or |
137 |
mmeineke |
1054 |
\texttt{.mdl} file. The position of the atoms given in the |
138 |
|
|
declaration are relative to the origin of the molecule, and is used |
139 |
|
|
when creating a system containing the molecule. |
140 |
mmeineke |
1044 |
|
141 |
|
|
As stated previously, one of the features that sets {\sc oopse} apart |
142 |
|
|
from most of the current molecular simulation packages is the ability |
143 |
|
|
to handle rigid body dynamics. Rigid bodies are non-spherical |
144 |
|
|
particles or collections of particles that have a constant internal |
145 |
|
|
potential and move collectively.\cite{Goldstein01} They are not |
146 |
mmeineke |
1068 |
included in most simulation packages because of the algorithmic |
147 |
|
|
complexity involved in propagating orientational degrees of |
148 |
|
|
freedom. Until recently, integrators which propagate orientational |
149 |
|
|
motion have been much worse than those available for translational |
150 |
|
|
motion. |
151 |
mmeineke |
1044 |
|
152 |
|
|
Moving a rigid body involves determination of both the force and |
153 |
|
|
torque applied by the surroundings, which directly affect the |
154 |
|
|
translational and rotational motion in turn. In order to accumulate |
155 |
|
|
the total force on a rigid body, the external forces and torques must |
156 |
|
|
first be calculated for all the internal particles. The total force on |
157 |
|
|
the rigid body is simply the sum of these external forces. |
158 |
|
|
Accumulation of the total torque on the rigid body is more complex |
159 |
mmeineke |
1068 |
than the force because the torque is applied to the center of mass of |
160 |
|
|
the rigid body. The torque on rigid body $i$ is |
161 |
mmeineke |
1044 |
\begin{equation} |
162 |
|
|
\boldsymbol{\tau}_i= |
163 |
mmeineke |
1068 |
\sum_{a}\biggl[(\mathbf{r}_{ia}-\mathbf{r}_i)\times \mathbf{f}_{ia} |
164 |
|
|
+ \boldsymbol{\tau}_{ia}\biggr] |
165 |
mmeineke |
1044 |
\label{eq:torqueAccumulate} |
166 |
|
|
\end{equation} |
167 |
|
|
where $\boldsymbol{\tau}_i$ and $\mathbf{r}_i$ are the torque on and |
168 |
|
|
position of the center of mass respectively, while $\mathbf{f}_{ia}$, |
169 |
|
|
$\mathbf{r}_{ia}$, and $\boldsymbol{\tau}_{ia}$ are the force on, |
170 |
|
|
position of, and torque on the component particles of the rigid body. |
171 |
|
|
|
172 |
|
|
The summation of the total torque is done in the body fixed axis of |
173 |
mmeineke |
1068 |
each rigid body. In order to move between the space fixed and body |
174 |
mmeineke |
1044 |
fixed coordinate axes, parameters describing the orientation must be |
175 |
|
|
maintained for each rigid body. At a minimum, the rotation matrix |
176 |
|
|
(\textbf{A}) can be described by the three Euler angles ($\phi, |
177 |
|
|
\theta,$ and $\psi$), where the elements of \textbf{A} are composed of |
178 |
|
|
trigonometric operations involving $\phi, \theta,$ and |
179 |
|
|
$\psi$.\cite{Goldstein01} In order to avoid numerical instabilities |
180 |
|
|
inherent in using the Euler angles, the four parameter ``quaternion'' |
181 |
|
|
scheme is often used. The elements of \textbf{A} can be expressed as |
182 |
|
|
arithmetic operations involving the four quaternions ($q_0, q_1, q_2,$ |
183 |
|
|
and $q_3$).\cite{allen87:csl} Use of quaternions also leads to |
184 |
|
|
performance enhancements, particularly for very small |
185 |
|
|
systems.\cite{Evans77} |
186 |
|
|
|
187 |
|
|
{\sc oopse} utilizes a relatively new scheme that propagates the |
188 |
mmeineke |
1068 |
entire nine parameter rotation matrix. Further discussion |
189 |
mmeineke |
1054 |
on this choice can be found in Sec.~\ref{oopseSec:integrate}. An example |
190 |
|
|
definition of a rigid body can be seen in Scheme |
191 |
mmeineke |
1044 |
\ref{sch:rigidBody}. The positions in the atom definitions are the |
192 |
|
|
placements of the atoms relative to the origin of the rigid body, |
193 |
|
|
which itself has a position relative to the origin of the molecule. |
194 |
|
|
|
195 |
mmeineke |
1045 |
\begin{lstlisting}[float,caption={[Defining rigid bodies]A sample definition of a rigid body},label={sch:rigidBody}] |
196 |
mmeineke |
1044 |
molecule{ |
197 |
|
|
name = "TIP3P_water"; |
198 |
|
|
nRigidBodies = 1; |
199 |
|
|
rigidBody[0]{ |
200 |
|
|
nAtoms = 3; |
201 |
|
|
atom[0]{ |
202 |
|
|
type = "O_TIP3P"; |
203 |
|
|
position( 0.0, 0.0, -0.06556 ); |
204 |
|
|
} |
205 |
|
|
atom[1]{ |
206 |
|
|
type = "H_TIP3P"; |
207 |
|
|
position( 0.0, 0.75695, 0.52032 ); |
208 |
|
|
} |
209 |
|
|
atom[2]{ |
210 |
|
|
type = "H_TIP3P"; |
211 |
|
|
position( 0.0, -0.75695, 0.52032 ); |
212 |
|
|
} |
213 |
|
|
position( 0.0, 0.0, 0.0 ); |
214 |
|
|
orientation( 0.0, 0.0, 1.0 ); |
215 |
|
|
} |
216 |
|
|
} |
217 |
|
|
\end{lstlisting} |
218 |
|
|
|
219 |
mmeineke |
1054 |
\subsection{\label{sec:LJPot}The Lennard Jones Force Field} |
220 |
mmeineke |
1044 |
|
221 |
|
|
The most basic force field implemented in {\sc oopse} is the |
222 |
mmeineke |
1054 |
Lennard-Jones force field, which mimics the van der Waals interaction at |
223 |
mmeineke |
1044 |
long distances, and uses an empirical repulsion at short |
224 |
|
|
distances. The Lennard-Jones potential is given by: |
225 |
|
|
\begin{equation} |
226 |
|
|
V_{\text{LJ}}(r_{ij}) = |
227 |
|
|
4\epsilon_{ij} \biggl[ |
228 |
|
|
\biggl(\frac{\sigma_{ij}}{r_{ij}}\biggr)^{12} |
229 |
|
|
- \biggl(\frac{\sigma_{ij}}{r_{ij}}\biggr)^{6} |
230 |
|
|
\biggr] |
231 |
|
|
\label{eq:lennardJonesPot} |
232 |
|
|
\end{equation} |
233 |
|
|
Where $r_{ij}$ is the distance between particles $i$ and $j$, |
234 |
|
|
$\sigma_{ij}$ scales the length of the interaction, and |
235 |
|
|
$\epsilon_{ij}$ scales the well depth of the potential. Scheme |
236 |
mmeineke |
1054 |
\ref{sch:LJFF} gives and example \texttt{.bass} file that |
237 |
|
|
sets up a system of 108 Ar particles to be simulated using the |
238 |
|
|
Lennard-Jones force field. |
239 |
mmeineke |
1044 |
|
240 |
mmeineke |
1045 |
\begin{lstlisting}[float,caption={[Invocation of the Lennard-Jones force field] A sample system using the Lennard-Jones force field.},label={sch:LJFF}] |
241 |
mmeineke |
1044 |
|
242 |
|
|
#include "argon.mdl" |
243 |
|
|
|
244 |
|
|
nComponents = 1; |
245 |
|
|
component{ |
246 |
|
|
type = "Ar"; |
247 |
|
|
nMol = 108; |
248 |
|
|
} |
249 |
|
|
|
250 |
|
|
initialConfig = "./argon.init"; |
251 |
|
|
|
252 |
|
|
forceField = "LJ"; |
253 |
|
|
\end{lstlisting} |
254 |
|
|
|
255 |
|
|
Because this potential is calculated between all pairs, the force |
256 |
|
|
evaluation can become computationally expensive for large systems. To |
257 |
|
|
keep the pair evaluations to a manageable number, {\sc oopse} employs |
258 |
mmeineke |
1068 |
a cut-off radius.\cite{allen87:csl} The cutoff radius can either be |
259 |
|
|
specified in the \texttt{.bass} file, or left as its default value of |
260 |
mmeineke |
1044 |
$2.5\sigma_{ii}$, where $\sigma_{ii}$ is the largest Lennard-Jones |
261 |
|
|
length parameter present in the simulation. Truncating the calculation |
262 |
|
|
at $r_{\text{cut}}$ introduces a discontinuity into the potential |
263 |
mmeineke |
1068 |
energy and the force. To offset this discontinuity in the potential, |
264 |
|
|
the energy value at $r_{\text{cut}}$ is subtracted from the |
265 |
|
|
potential. This causes the potential to go to zero smoothly at the |
266 |
|
|
cut-off radius, and preserves conservation of energy in integrating |
267 |
|
|
the equations of motion. |
268 |
mmeineke |
1044 |
|
269 |
|
|
Interactions between dissimilar particles requires the generation of |
270 |
|
|
cross term parameters for $\sigma$ and $\epsilon$. These are |
271 |
|
|
calculated through the Lorentz-Berthelot mixing |
272 |
|
|
rules:\cite{allen87:csl} |
273 |
|
|
\begin{equation} |
274 |
|
|
\sigma_{ij} = \frac{1}{2}[\sigma_{ii} + \sigma_{jj}] |
275 |
|
|
\label{eq:sigmaMix} |
276 |
|
|
\end{equation} |
277 |
|
|
and |
278 |
|
|
\begin{equation} |
279 |
|
|
\epsilon_{ij} = \sqrt{\epsilon_{ii} \epsilon_{jj}} |
280 |
|
|
\label{eq:epsilonMix} |
281 |
|
|
\end{equation} |
282 |
|
|
|
283 |
mmeineke |
1051 |
\subsection{\label{oopseSec:DUFF}Dipolar Unified-Atom Force Field} |
284 |
mmeineke |
1044 |
|
285 |
|
|
The dipolar unified-atom force field ({\sc duff}) was developed to |
286 |
|
|
simulate lipid bilayers. The simulations require a model capable of |
287 |
|
|
forming bilayers, while still being sufficiently computationally |
288 |
mmeineke |
1054 |
efficient to allow large systems ($\sim$100's of phospholipids, |
289 |
|
|
$\sim$1000's of waters) to be simulated for long times |
290 |
|
|
($\sim$10's of nanoseconds). |
291 |
mmeineke |
1044 |
|
292 |
|
|
With this goal in mind, {\sc duff} has no point |
293 |
|
|
charges. Charge-neutral distributions were replaced with dipoles, |
294 |
|
|
while most atoms and groups of atoms were reduced to Lennard-Jones |
295 |
|
|
interaction sites. This simplification cuts the length scale of long |
296 |
mmeineke |
1068 |
range interactions from $\frac{1}{r}$ to $\frac{1}{r^3}$, and allows |
297 |
|
|
us to avoid the computationally expensive Ewald sum. Instead, we can |
298 |
|
|
use neighbor-lists and cutoff radii for the dipolar interactions, or |
299 |
|
|
include a reaction field to mimic larger range interactions. |
300 |
mmeineke |
1044 |
|
301 |
|
|
As an example, lipid head-groups in {\sc duff} are represented as |
302 |
mmeineke |
1068 |
point dipole interaction sites. By placing a dipole at the head group |
303 |
|
|
center of mass, our model mimics the charge separation found in common |
304 |
|
|
phospholipids such as phosphatidylcholine.\cite{Cevc87} Additionally, |
305 |
|
|
a large Lennard-Jones site is located at the pseudoatom's center of |
306 |
|
|
mass. The model is illustrated by the red atom in |
307 |
|
|
Fig.~\ref{oopseFig:lipidModel}. The water model we use to complement |
308 |
|
|
the dipoles of the lipids is our reparameterization of the soft sticky |
309 |
|
|
dipole (SSD) model of Ichiye |
310 |
mmeineke |
1044 |
\emph{et al.}\cite{liu96:new_model} |
311 |
|
|
|
312 |
|
|
\begin{figure} |
313 |
mmeineke |
1045 |
\centering |
314 |
|
|
\includegraphics[width=\linewidth]{lipidModel.eps} |
315 |
mmeineke |
1044 |
\caption{A representation of the lipid model. $\phi$ is the torsion angle, $\theta$ % |
316 |
|
|
is the bend angle, $\mu$ is the dipole moment of the head group, and n |
317 |
|
|
is the chain length.} |
318 |
mmeineke |
1045 |
\label{oopseFig:lipidModel} |
319 |
mmeineke |
1044 |
\end{figure} |
320 |
|
|
|
321 |
|
|
We have used a set of scalable parameters to model the alkyl groups |
322 |
|
|
with Lennard-Jones sites. For this, we have borrowed parameters from |
323 |
|
|
the TraPPE force field of Siepmann |
324 |
|
|
\emph{et al}.\cite{Siepmann1998} TraPPE is a unified-atom |
325 |
|
|
representation of n-alkanes, which is parametrized against phase |
326 |
|
|
equilibria using Gibbs ensemble Monte Carlo simulation |
327 |
|
|
techniques.\cite{Siepmann1998} One of the advantages of TraPPE is that |
328 |
|
|
it generalizes the types of atoms in an alkyl chain to keep the number |
329 |
mmeineke |
1068 |
of pseudoatoms to a minimum; the parameters for a unified atom such as |
330 |
mmeineke |
1044 |
$\text{CH}_2$ do not change depending on what species are bonded to |
331 |
|
|
it. |
332 |
|
|
|
333 |
|
|
TraPPE also constrains all bonds to be of fixed length. Typically, |
334 |
|
|
bond vibrations are the fastest motions in a molecular dynamic |
335 |
|
|
simulation. Small time steps between force evaluations must be used to |
336 |
mmeineke |
1068 |
ensure adequate energy conservation in the bond degrees of freedom. By |
337 |
|
|
constraining the bond lengths, larger time steps may be used when |
338 |
|
|
integrating the equations of motion. A simulation using {\sc duff} is |
339 |
|
|
illustrated in Scheme \ref{sch:DUFF}. |
340 |
mmeineke |
1044 |
|
341 |
mmeineke |
1045 |
\begin{lstlisting}[float,caption={[Invocation of {\sc duff}]Sample \texttt{.bass} file showing a simulation utilizing {\sc duff}},label={sch:DUFF}] |
342 |
mmeineke |
1044 |
|
343 |
|
|
#include "water.mdl" |
344 |
|
|
#include "lipid.mdl" |
345 |
|
|
|
346 |
|
|
nComponents = 2; |
347 |
|
|
component{ |
348 |
|
|
type = "simpleLipid_16"; |
349 |
|
|
nMol = 60; |
350 |
|
|
} |
351 |
|
|
|
352 |
|
|
component{ |
353 |
|
|
type = "SSD_water"; |
354 |
|
|
nMol = 1936; |
355 |
|
|
} |
356 |
|
|
|
357 |
|
|
initialConfig = "bilayer.init"; |
358 |
|
|
|
359 |
|
|
forceField = "DUFF"; |
360 |
|
|
|
361 |
|
|
\end{lstlisting} |
362 |
|
|
|
363 |
mmeineke |
1051 |
\subsection{\label{oopseSec:energyFunctions}{\sc duff} Energy Functions} |
364 |
mmeineke |
1044 |
|
365 |
|
|
The total potential energy function in {\sc duff} is |
366 |
|
|
\begin{equation} |
367 |
|
|
V = \sum^{N}_{I=1} V^{I}_{\text{Internal}} |
368 |
mmeineke |
1054 |
+ \sum^{N-1}_{I=1} \sum_{J>I} V^{IJ}_{\text{Cross}} |
369 |
mmeineke |
1044 |
\label{eq:totalPotential} |
370 |
|
|
\end{equation} |
371 |
|
|
Where $V^{I}_{\text{Internal}}$ is the internal potential of molecule $I$: |
372 |
|
|
\begin{equation} |
373 |
|
|
V^{I}_{\text{Internal}} = |
374 |
|
|
\sum_{\theta_{ijk} \in I} V_{\text{bend}}(\theta_{ijk}) |
375 |
|
|
+ \sum_{\phi_{ijkl} \in I} V_{\text{torsion}}(\phi_{ijkl}) |
376 |
|
|
+ \sum_{i \in I} \sum_{(j>i+4) \in I} |
377 |
|
|
\biggl[ V_{\text{LJ}}(r_{ij}) + V_{\text{dipole}} |
378 |
|
|
(\mathbf{r}_{ij},\boldsymbol{\Omega}_{i},\boldsymbol{\Omega}_{j}) |
379 |
|
|
\biggr] |
380 |
|
|
\label{eq:internalPotential} |
381 |
|
|
\end{equation} |
382 |
|
|
Here $V_{\text{bend}}$ is the bend potential for all 1, 3 bonded pairs |
383 |
|
|
within the molecule $I$, and $V_{\text{torsion}}$ is the torsion potential |
384 |
|
|
for all 1, 4 bonded pairs. The pairwise portions of the internal |
385 |
|
|
potential are excluded for pairs that are closer than three bonds, |
386 |
|
|
i.e.~atom pairs farther away than a torsion are included in the |
387 |
|
|
pair-wise loop. |
388 |
|
|
|
389 |
|
|
|
390 |
|
|
The bend potential of a molecule is represented by the following function: |
391 |
|
|
\begin{equation} |
392 |
|
|
V_{\text{bend}}(\theta_{ijk}) = k_{\theta}( \theta_{ijk} - \theta_0 )^2 \label{eq:bendPot} |
393 |
|
|
\end{equation} |
394 |
|
|
Where $\theta_{ijk}$ is the angle defined by atoms $i$, $j$, and $k$ |
395 |
mmeineke |
1054 |
(see Fig.~\ref{oopseFig:lipidModel}), $\theta_0$ is the equilibrium |
396 |
mmeineke |
1044 |
bond angle, and $k_{\theta}$ is the force constant which determines the |
397 |
|
|
strength of the harmonic bend. The parameters for $k_{\theta}$ and |
398 |
|
|
$\theta_0$ are borrowed from those in TraPPE.\cite{Siepmann1998} |
399 |
|
|
|
400 |
|
|
The torsion potential and parameters are also borrowed from TraPPE. It is |
401 |
|
|
of the form: |
402 |
|
|
\begin{equation} |
403 |
|
|
V_{\text{torsion}}(\phi) = c_1[1 + \cos \phi] |
404 |
|
|
+ c_2[1 + \cos(2\phi)] |
405 |
|
|
+ c_3[1 + \cos(3\phi)] |
406 |
|
|
\label{eq:origTorsionPot} |
407 |
|
|
\end{equation} |
408 |
mmeineke |
1068 |
Where: |
409 |
mmeineke |
1044 |
\begin{equation} |
410 |
mmeineke |
1068 |
\cos\phi = (\hat{\mathbf{r}}_{ij} \times \hat{\mathbf{r}}_{jk}) \cdot |
411 |
|
|
(\hat{\mathbf{r}}_{jk} \times \hat{\mathbf{r}}_{kl}) |
412 |
|
|
\label{eq:torsPhi} |
413 |
|
|
\end{equation} |
414 |
|
|
Here, $\hat{\mathbf{r}}_{\alpha\beta}$ are the set of unit bond |
415 |
|
|
vectors between atoms $i$, $j$, $k$, and $l$. For computational |
416 |
|
|
efficiency, the torsion potential has been recast after the method of |
417 |
|
|
{\sc charmm},\cite{Brooks83} in which the angle series is converted to |
418 |
|
|
a power series of the form: |
419 |
|
|
\begin{equation} |
420 |
mmeineke |
1044 |
V_{\text{torsion}}(\phi) = |
421 |
|
|
k_3 \cos^3 \phi + k_2 \cos^2 \phi + k_1 \cos \phi + k_0 |
422 |
|
|
\label{eq:torsionPot} |
423 |
|
|
\end{equation} |
424 |
|
|
Where: |
425 |
|
|
\begin{align*} |
426 |
|
|
k_0 &= c_1 + c_3 \\ |
427 |
|
|
k_1 &= c_1 - 3c_3 \\ |
428 |
|
|
k_2 &= 2 c_2 \\ |
429 |
|
|
k_3 &= 4c_3 |
430 |
|
|
\end{align*} |
431 |
|
|
By recasting the potential as a power series, repeated trigonometric |
432 |
|
|
evaluations are avoided during the calculation of the potential energy. |
433 |
|
|
|
434 |
|
|
|
435 |
|
|
The cross potential between molecules $I$ and $J$, $V^{IJ}_{\text{Cross}}$, is |
436 |
|
|
as follows: |
437 |
|
|
\begin{equation} |
438 |
|
|
V^{IJ}_{\text{Cross}} = |
439 |
|
|
\sum_{i \in I} \sum_{j \in J} |
440 |
|
|
\biggl[ V_{\text{LJ}}(r_{ij}) + V_{\text{dipole}} |
441 |
|
|
(\mathbf{r}_{ij},\boldsymbol{\Omega}_{i},\boldsymbol{\Omega}_{j}) |
442 |
|
|
+ V_{\text{sticky}} |
443 |
|
|
(\mathbf{r}_{ij},\boldsymbol{\Omega}_{i},\boldsymbol{\Omega}_{j}) |
444 |
|
|
\biggr] |
445 |
|
|
\label{eq:crossPotentail} |
446 |
|
|
\end{equation} |
447 |
|
|
Where $V_{\text{LJ}}$ is the Lennard Jones potential, |
448 |
|
|
$V_{\text{dipole}}$ is the dipole dipole potential, and |
449 |
|
|
$V_{\text{sticky}}$ is the sticky potential defined by the SSD model |
450 |
mmeineke |
1054 |
(Sec.~\ref{oopseSec:SSD}). Note that not all atom types include all |
451 |
mmeineke |
1044 |
interactions. |
452 |
|
|
|
453 |
|
|
The dipole-dipole potential has the following form: |
454 |
|
|
\begin{equation} |
455 |
|
|
V_{\text{dipole}}(\mathbf{r}_{ij},\boldsymbol{\Omega}_{i}, |
456 |
|
|
\boldsymbol{\Omega}_{j}) = \frac{|\mu_i||\mu_j|}{4\pi\epsilon_{0}r_{ij}^{3}} \biggl[ |
457 |
|
|
\boldsymbol{\hat{u}}_{i} \cdot \boldsymbol{\hat{u}}_{j} |
458 |
|
|
- |
459 |
mmeineke |
1068 |
3(\boldsymbol{\hat{u}}_i \cdot \hat{\mathbf{r}}_{ij}) % |
460 |
|
|
(\boldsymbol{\hat{u}}_j \cdot \hat{\mathbf{r}}_{ij}) \biggr] |
461 |
mmeineke |
1044 |
\label{eq:dipolePot} |
462 |
|
|
\end{equation} |
463 |
|
|
Here $\mathbf{r}_{ij}$ is the vector starting at atom $i$ pointing |
464 |
|
|
towards $j$, and $\boldsymbol{\Omega}_i$ and $\boldsymbol{\Omega}_j$ |
465 |
|
|
are the orientational degrees of freedom for atoms $i$ and $j$ |
466 |
|
|
respectively. $|\mu_i|$ is the magnitude of the dipole moment of atom |
467 |
mmeineke |
1054 |
$i$, $\boldsymbol{\hat{u}}_i$ is the standard unit orientation vector |
468 |
|
|
of $\boldsymbol{\Omega}_i$, and $\boldsymbol{\hat{r}}_{ij}$ is the |
469 |
|
|
unit vector pointing along $\mathbf{r}_{ij}$ |
470 |
|
|
($\boldsymbol{\hat{r}}_{ij}=\mathbf{r}_{ij}/|\mathbf{r}_{ij}|$). |
471 |
mmeineke |
1044 |
|
472 |
mmeineke |
1068 |
To improve computational efficiency of the dipole-dipole interactions, |
473 |
|
|
{\sc oopse} employs an electrostatic cutoff radius. This parameter can |
474 |
|
|
be set in the \texttt{.bass} file, and controls the length scale over |
475 |
|
|
which dipole interactions are felt. To compensate for the |
476 |
|
|
discontinuity in the potential and the forces at the cutoff radius, we |
477 |
|
|
have implemented a switching function to smoothly scale the |
478 |
|
|
dipole-dipole interaction at the cutoff. |
479 |
|
|
\begin{equation} |
480 |
|
|
S(r_{ij}) = |
481 |
|
|
\begin{cases} |
482 |
|
|
1 & \text{if $r_{ij} \le r_t$},\\ |
483 |
|
|
\frac{(r_{\text{cut}} + 2r_{ij} - 3r_t)(r_{\text{cut}} - r_{ij})^2} |
484 |
|
|
{(r_{\text{cut}} - r_t)^2} |
485 |
|
|
& \text{if $r_t < r_{ij} \le r_{\text{cut}}$}, \\ |
486 |
|
|
0 & \text{if $r_{ij} > r_{\text{cut}}$.} |
487 |
|
|
\end{cases} |
488 |
|
|
\label{eq:dipoleSwitching} |
489 |
|
|
\end{equation} |
490 |
|
|
Here $S(r_{ij})$ scales the potential at a given $r_{ij}$, and $r_t$ |
491 |
|
|
is the taper radius some given thickness less than the electrostatic |
492 |
|
|
cutoff. The switching thickness can be set in the \texttt{.bass} file. |
493 |
mmeineke |
1044 |
|
494 |
mmeineke |
1068 |
\subsection{\label{oopseSec:SSD}The {\sc duff} Water Models: SSD/E and SSD/RF} |
495 |
mmeineke |
1044 |
|
496 |
|
|
In the interest of computational efficiency, the default solvent used |
497 |
|
|
by {\sc oopse} is the extended Soft Sticky Dipole (SSD/E) water |
498 |
|
|
model.\cite{Gezelter04} The original SSD was developed by Ichiye |
499 |
|
|
\emph{et al.}\cite{liu96:new_model} as a modified form of the hard-sphere |
500 |
|
|
water model proposed by Bratko, Blum, and |
501 |
|
|
Luzar.\cite{Bratko85,Bratko95} It consists of a single point dipole |
502 |
|
|
with a Lennard-Jones core and a sticky potential that directs the |
503 |
|
|
particles to assume the proper hydrogen bond orientation in the first |
504 |
|
|
solvation shell. Thus, the interaction between two SSD water molecules |
505 |
|
|
\emph{i} and \emph{j} is given by the potential |
506 |
|
|
\begin{equation} |
507 |
|
|
V_{ij} = |
508 |
|
|
V_{ij}^{LJ} (r_{ij})\ + V_{ij}^{dp} |
509 |
|
|
(\mathbf{r}_{ij},\boldsymbol{\Omega}_i,\boldsymbol{\Omega}_j)\ + |
510 |
|
|
V_{ij}^{sp} |
511 |
|
|
(\mathbf{r}_{ij},\boldsymbol{\Omega}_i,\boldsymbol{\Omega}_j), |
512 |
|
|
\label{eq:ssdPot} |
513 |
|
|
\end{equation} |
514 |
|
|
where the $\mathbf{r}_{ij}$ is the position vector between molecules |
515 |
|
|
\emph{i} and \emph{j} with magnitude equal to the distance $r_{ij}$, and |
516 |
|
|
$\boldsymbol{\Omega}_i$ and $\boldsymbol{\Omega}_j$ represent the |
517 |
|
|
orientations of the respective molecules. The Lennard-Jones and dipole |
518 |
|
|
parts of the potential are given by equations \ref{eq:lennardJonesPot} |
519 |
|
|
and \ref{eq:dipolePot} respectively. The sticky part is described by |
520 |
|
|
the following, |
521 |
|
|
\begin{equation} |
522 |
|
|
u_{ij}^{sp}(\mathbf{r}_{ij},\boldsymbol{\Omega}_i,\boldsymbol{\Omega}_j)= |
523 |
|
|
\frac{\nu_0}{2}[s(r_{ij})w(\mathbf{r}_{ij}, |
524 |
|
|
\boldsymbol{\Omega}_i,\boldsymbol{\Omega}_j) + |
525 |
|
|
s^\prime(r_{ij})w^\prime(\mathbf{r}_{ij}, |
526 |
|
|
\boldsymbol{\Omega}_i,\boldsymbol{\Omega}_j)]\ , |
527 |
|
|
\label{eq:stickyPot} |
528 |
|
|
\end{equation} |
529 |
|
|
where $\nu_0$ is a strength parameter for the sticky potential, and |
530 |
|
|
$s$ and $s^\prime$ are cubic switching functions which turn off the |
531 |
|
|
sticky interaction beyond the first solvation shell. The $w$ function |
532 |
|
|
can be thought of as an attractive potential with tetrahedral |
533 |
|
|
geometry: |
534 |
|
|
\begin{equation} |
535 |
|
|
w({\bf r}_{ij},{\bf \Omega}_i,{\bf \Omega}_j)= |
536 |
|
|
\sin\theta_{ij}\sin2\theta_{ij}\cos2\phi_{ij}, |
537 |
|
|
\label{eq:stickyW} |
538 |
|
|
\end{equation} |
539 |
|
|
while the $w^\prime$ function counters the normal aligned and |
540 |
|
|
anti-aligned structures favored by point dipoles: |
541 |
|
|
\begin{equation} |
542 |
|
|
w^\prime({\bf r}_{ij},{\bf \Omega}_i,{\bf \Omega}_j)= |
543 |
|
|
(\cos\theta_{ij}-0.6)^2(\cos\theta_{ij}+0.8)^2-w^0, |
544 |
|
|
\label{eq:stickyWprime} |
545 |
|
|
\end{equation} |
546 |
|
|
It should be noted that $w$ is proportional to the sum of the $Y_3^2$ |
547 |
|
|
and $Y_3^{-2}$ spherical harmonics (a linear combination which |
548 |
|
|
enhances the tetrahedral geometry for hydrogen bonded structures), |
549 |
|
|
while $w^\prime$ is a purely empirical function. A more detailed |
550 |
|
|
description of the functional parts and variables in this potential |
551 |
|
|
can be found in the original SSD |
552 |
|
|
articles.\cite{liu96:new_model,liu96:monte_carlo,chandra99:ssd_md,Ichiye03} |
553 |
|
|
|
554 |
mmeineke |
1068 |
Since SSD/E is a single-point {\it dipolar} model, the force |
555 |
mmeineke |
1044 |
calculations are simplified significantly relative to the standard |
556 |
|
|
{\it charged} multi-point models. In the original Monte Carlo |
557 |
|
|
simulations using this model, Ichiye {\it et al.} reported that using |
558 |
|
|
SSD decreased computer time by a factor of 6-7 compared to other |
559 |
|
|
models.\cite{liu96:new_model} What is most impressive is that these savings |
560 |
|
|
did not come at the expense of accurate depiction of the liquid state |
561 |
mmeineke |
1068 |
properties. Indeed, SSD/E maintains reasonable agreement with the Head-Gordon |
562 |
mmeineke |
1044 |
diffraction data for the structural features of liquid |
563 |
mmeineke |
1068 |
water.\cite{hura00,liu96:new_model} Additionally, the dynamical properties |
564 |
|
|
exhibited by SSD/E agree with experiment better than those of more |
565 |
mmeineke |
1044 |
computationally expensive models (like TIP3P and |
566 |
|
|
SPC/E).\cite{chandra99:ssd_md} The combination of speed and accurate depiction |
567 |
mmeineke |
1068 |
of solvent properties makes SSD/E a very attractive model for the |
568 |
mmeineke |
1044 |
simulation of large scale biochemical simulations. |
569 |
|
|
|
570 |
|
|
Recent constant pressure simulations revealed issues in the original |
571 |
|
|
SSD model that led to lower than expected densities at all target |
572 |
|
|
pressures.\cite{Ichiye03,Gezelter04} The default model in {\sc oopse} |
573 |
|
|
is therefore SSD/E, a density corrected derivative of SSD that |
574 |
|
|
exhibits improved liquid structure and transport behavior. If the use |
575 |
|
|
of a reaction field long-range interaction correction is desired, it |
576 |
|
|
is recommended that the parameters be modified to those of the SSD/RF |
577 |
|
|
model. Solvent parameters can be easily modified in an accompanying |
578 |
mmeineke |
1068 |
\texttt{.bass} file as illustrated in the scheme below. A table of the |
579 |
mmeineke |
1044 |
parameter values and the drawbacks and benefits of the different |
580 |
mmeineke |
1054 |
density corrected SSD models can be found in |
581 |
|
|
reference~\cite{Gezelter04}. |
582 |
mmeineke |
1044 |
|
583 |
mmeineke |
1045 |
\begin{lstlisting}[float,caption={[A simulation of {\sc ssd} water]An example file showing a simulation including {\sc ssd} water.},label={sch:ssd}] |
584 |
mmeineke |
1044 |
|
585 |
|
|
#include "water.mdl" |
586 |
|
|
|
587 |
|
|
nComponents = 1; |
588 |
|
|
component{ |
589 |
|
|
type = "SSD_water"; |
590 |
|
|
nMol = 864; |
591 |
|
|
} |
592 |
|
|
|
593 |
|
|
initialConfig = "liquidWater.init"; |
594 |
|
|
|
595 |
|
|
forceField = "DUFF"; |
596 |
|
|
|
597 |
|
|
/* |
598 |
|
|
* The following two flags set the cutoff |
599 |
|
|
* radius for the electrostatic forces |
600 |
|
|
* as well as the skin thickness of the switching |
601 |
|
|
* function. |
602 |
|
|
*/ |
603 |
|
|
|
604 |
|
|
electrostaticCutoffRadius = 9.2; |
605 |
|
|
electrostaticSkinThickness = 1.38; |
606 |
|
|
|
607 |
|
|
\end{lstlisting} |
608 |
|
|
|
609 |
|
|
|
610 |
mmeineke |
1051 |
\subsection{\label{oopseSec:eam}Embedded Atom Method} |
611 |
mmeineke |
1044 |
|
612 |
mmeineke |
1054 |
There are Molecular Dynamics packages which have the |
613 |
mmeineke |
1044 |
capacity to simulate metallic systems, including some that have |
614 |
|
|
parallel computational abilities\cite{plimpton93}. Potentials that |
615 |
|
|
describe bonding transition metal |
616 |
mmeineke |
1068 |
systems\cite{Finnis84,Ercolessi88,Chen90,Qi99,Ercolessi02} have an |
617 |
mmeineke |
1044 |
attractive interaction which models ``Embedding'' |
618 |
|
|
a positively charged metal ion in the electron density due to the |
619 |
|
|
free valance ``sea'' of electrons created by the surrounding atoms in |
620 |
mmeineke |
1068 |
the system. A mostly-repulsive pairwise part of the potential |
621 |
mmeineke |
1044 |
describes the interaction of the positively charged metal core ions |
622 |
|
|
with one another. A particular potential description called the |
623 |
|
|
Embedded Atom Method\cite{Daw84,FBD86,johnson89,Lu97}({\sc eam}) that has |
624 |
|
|
particularly wide adoption has been selected for inclusion in {\sc oopse}. A |
625 |
mmeineke |
1068 |
good review of {\sc eam} and other metallic potential formulations was written |
626 |
mmeineke |
1044 |
by Voter.\cite{voter} |
627 |
|
|
|
628 |
|
|
The {\sc eam} potential has the form: |
629 |
|
|
\begin{eqnarray} |
630 |
|
|
V & = & \sum_{i} F_{i}\left[\rho_{i}\right] + \sum_{i} \sum_{j \neq i} |
631 |
|
|
\phi_{ij}({\bf r}_{ij}) \\ |
632 |
|
|
\rho_{i} & = & \sum_{j \neq i} f_{j}({\bf r}_{ij}) |
633 |
mmeineke |
1068 |
\end{eqnarray} |
634 |
mmeineke |
1044 |
where $F_{i} $ is the embedding function that equates the energy required to embed a |
635 |
|
|
positively-charged core ion $i$ into a linear superposition of |
636 |
|
|
spherically averaged atomic electron densities given by |
637 |
|
|
$\rho_{i}$. $\phi_{ij}$ is a primarily repulsive pairwise interaction |
638 |
|
|
between atoms $i$ and $j$. In the original formulation of |
639 |
mmeineke |
1054 |
{\sc eam}\cite{Daw84}, $\phi_{ij}$ was an entirely repulsive term, however |
640 |
mmeineke |
1044 |
in later refinements to EAM have shown that non-uniqueness between $F$ |
641 |
|
|
and $\phi$ allow for more general forms for $\phi$.\cite{Daw89} |
642 |
|
|
There is a cutoff distance, $r_{cut}$, which limits the |
643 |
|
|
summations in the {\sc eam} equation to the few dozen atoms |
644 |
|
|
surrounding atom $i$ for both the density $\rho$ and pairwise $\phi$ |
645 |
mmeineke |
1054 |
interactions. Foiles et al. fit EAM potentials for fcc metals Cu, Ag, Au, Ni, Pd, Pt and alloys of these metals\cite{FBD86}. These potential fits are in the DYNAMO 86 format and are included with {\sc oopse}. |
646 |
mmeineke |
1044 |
|
647 |
|
|
|
648 |
mmeineke |
1051 |
\subsection{\label{oopseSec:pbc}Periodic Boundary Conditions} |
649 |
mmeineke |
1044 |
|
650 |
|
|
\newcommand{\roundme}{\operatorname{round}} |
651 |
|
|
|
652 |
mmeineke |
1068 |
\textit{Periodic boundary conditions} are widely used to simulate bulk properties with a relatively small number of particles. The |
653 |
|
|
simulation box is replicated throughout space to form an infinite |
654 |
|
|
lattice. During the simulation, when a particle moves in the primary |
655 |
|
|
cell, its image in other cells move in exactly the same direction with |
656 |
|
|
exactly the same orientation. Thus, as a particle leaves the primary |
657 |
|
|
cell, one of its images will enter through the opposite face. If the |
658 |
|
|
simulation box is large enough to avoid ``feeling'' the symmetries of |
659 |
|
|
the periodic lattice, surface effects can be ignored. The available |
660 |
|
|
periodic cells in OOPSE are cubic, orthorhombic and parallelepiped. We |
661 |
|
|
use a $3 \times 3$ matrix, $\mathbf{H}$, to describe the shape and |
662 |
|
|
size of the simulation box. $\mathbf{H}$ is defined: |
663 |
mmeineke |
1044 |
\begin{equation} |
664 |
mmeineke |
1068 |
\mathbf{H} = ( \mathbf{h}_x, \mathbf{h}_y, \mathbf{h}_z ) |
665 |
mmeineke |
1044 |
\end{equation} |
666 |
mmeineke |
1068 |
Where $\mathbf{h}_j$ is the column vector of the $j$th axis of the |
667 |
|
|
box. During the course of the simulation both the size and shape of |
668 |
|
|
the box can be changed to allow volume fluctations when constraining |
669 |
|
|
the pressure. |
670 |
mmeineke |
1044 |
|
671 |
mmeineke |
1068 |
A real space vector, $\mathbf{r}$ can be transformed in to a box space |
672 |
|
|
vector, $\mathbf{s}$, and back through the following transformations: |
673 |
|
|
\begin{align} |
674 |
|
|
\mathbf{s} &= \mathbf{H}^{-1} \mathbf{r} \\ |
675 |
|
|
\mathbf{r} &= \mathbf{H} \mathbf{s} |
676 |
|
|
\end{align} |
677 |
|
|
The vector $\mathbf{s}$ is now a vector expressed as the number of box |
678 |
|
|
lengths in the $\mathbf{h}_x$, $\mathbf{h}_y$, and $\mathbf{h}_z$ |
679 |
|
|
directions. To find the minimum image of a vector $\mathbf{r}$, we |
680 |
|
|
first convert it to its corresponding vector in box space, and then, |
681 |
|
|
cast each element to lie on the in the range $[-0.5,0.5]$: |
682 |
mmeineke |
1044 |
\begin{equation} |
683 |
|
|
s_{i}^{\prime}=s_{i}-\roundme(s_{i}) |
684 |
|
|
\end{equation} |
685 |
mmeineke |
1068 |
Where $s_i$ is the $i$th element of $\mathbf{s}$, and |
686 |
|
|
$\roundme(s_i)$is given by |
687 |
mmeineke |
1044 |
\begin{equation} |
688 |
mmeineke |
1068 |
\roundme(x) = |
689 |
|
|
\begin{cases} |
690 |
|
|
\lfloor x+0.5 \rfloor & \text{if $x \ge 0$} \\ |
691 |
|
|
\lceil x-0.5 \rceil & \text{if $x < 0$ } |
692 |
|
|
\end{cases} |
693 |
mmeineke |
1044 |
\end{equation} |
694 |
mmeineke |
1068 |
Here $\lfloor x \rfloor$ is the floor operator, and gives the largest |
695 |
|
|
integer value that is not greater than $x$, and $\lceil x \rceil$ is |
696 |
|
|
the ceiling operator, and gives the smallest integer that is not less |
697 |
|
|
than $x$. For example, $\roundme(3.6)=4$, $\roundme(3.1)=3$, |
698 |
|
|
$\roundme(-3.6)=-4$, $\roundme(-3.1)=-3$. |
699 |
mmeineke |
1044 |
|
700 |
|
|
Finally, we obtain the minimum image coordinates $\mathbf{r}^{\prime}$ by |
701 |
mmeineke |
1068 |
transforming back to real space, |
702 |
mmeineke |
1044 |
\begin{equation} |
703 |
mmeineke |
1068 |
\mathbf{r}^{\prime}=\mathbf{H}^{-1}\mathbf{s}^{\prime}% |
704 |
mmeineke |
1044 |
\end{equation} |
705 |
mmeineke |
1068 |
In this way, particles are allowed to diffuse freely in $\mathbf{r}$, |
706 |
|
|
but their minimum images, $\mathbf{r}^{\prime}$ are used to compute |
707 |
|
|
the interatomic forces. |
708 |
mmeineke |
1044 |
|
709 |
|
|
|
710 |
mmeineke |
1051 |
\section{\label{oopseSec:IOfiles}Input and Output Files} |
711 |
mmeineke |
1044 |
|
712 |
|
|
\subsection{{\sc bass} and Model Files} |
713 |
|
|
|
714 |
mmeineke |
1071 |
Every {\sc oopse} simulation begins with a Bizarre Atom Simulation |
715 |
|
|
Syntax ({\sc bass}) file. {\sc bass} is a script syntax that is parsed |
716 |
|
|
by {\sc oopse} at runtime. The {\sc bass} file allows for the user to |
717 |
|
|
completely describe the system they wish to simulate, as well as tailor |
718 |
|
|
{\sc oopse}'s behavior during the simulation. {\sc bass} files are |
719 |
|
|
denoted with the extension |
720 |
mmeineke |
1044 |
\texttt{.bass}, an example file is shown in |
721 |
mmeineke |
1071 |
Scheme~\ref{sch:bassExample}. |
722 |
mmeineke |
1044 |
|
723 |
mmeineke |
1071 |
\begin{lstlisting}[float,caption={[An example of a complete {\sc bass} file] An example showing a complete {\sc bass} file.},label={sch:bassExample}] |
724 |
mmeineke |
1044 |
|
725 |
mmeineke |
1071 |
molecule{ |
726 |
|
|
name = "Ar"; |
727 |
|
|
nAtoms = 1; |
728 |
|
|
atom[0]{ |
729 |
|
|
type="Ar"; |
730 |
|
|
position( 0.0, 0.0, 0.0 ); |
731 |
|
|
} |
732 |
|
|
} |
733 |
|
|
|
734 |
|
|
nComponents = 1; |
735 |
|
|
component{ |
736 |
|
|
type = "Ar"; |
737 |
|
|
nMol = 108; |
738 |
|
|
} |
739 |
|
|
|
740 |
|
|
initialConfig = "./argon.init"; |
741 |
|
|
|
742 |
|
|
forceField = "LJ"; |
743 |
|
|
ensemble = "NVE"; // specify the simulation enesemble |
744 |
|
|
dt = 1.0; // the time step for integration |
745 |
|
|
runTime = 1e3; // the total simulation run time |
746 |
|
|
sampleTime = 100; // trajectory file frequency |
747 |
|
|
statusTime = 50; // statistics file frequency |
748 |
|
|
|
749 |
|
|
\end{lstlisting} |
750 |
|
|
|
751 |
mmeineke |
1054 |
Within the \texttt{.bass} file it is necessary to provide a complete |
752 |
mmeineke |
1044 |
description of the molecule before it is actually placed in the |
753 |
mmeineke |
1071 |
simulation. The {\sc bass} syntax was originally developed with this |
754 |
|
|
goal in mind, and allows for the specification of all the atoms in a |
755 |
|
|
molecular prototype, as well as any bonds, bends, or torsions. These |
756 |
mmeineke |
1044 |
descriptions can become lengthy for complex molecules, and it would be |
757 |
mmeineke |
1071 |
inconvenient to duplicate the simulation at the beginning of each {\sc |
758 |
|
|
bass} script. Addressing this issue {\sc bass} allows for the |
759 |
|
|
inclusion of model files at the top of a \texttt{.bass} file. These |
760 |
|
|
model files, denoted with the \texttt{.mdl} extension, allow the user |
761 |
|
|
to describe a molecular prototype once, then simply include it into |
762 |
|
|
each simulation containing that molecule. Returning to the example in |
763 |
|
|
Scheme~\ref{sch:bassExample}, the \texttt{.mdl} file's contents would |
764 |
|
|
be Scheme~\ref{sch:mdlExample}, and the new \texttt{.bass} file would |
765 |
|
|
become Scheme~\ref{sch:bassExPrime}. |
766 |
mmeineke |
1044 |
|
767 |
mmeineke |
1071 |
\begin{lstlisting}[float,caption={An example \texttt{.mdl} file.},label={sch:mdlExample}] |
768 |
|
|
|
769 |
|
|
molecule{ |
770 |
|
|
name = "Ar"; |
771 |
|
|
nAtoms = 1; |
772 |
|
|
atom[0]{ |
773 |
|
|
type="Ar"; |
774 |
|
|
position( 0.0, 0.0, 0.0 ); |
775 |
|
|
} |
776 |
|
|
} |
777 |
|
|
|
778 |
|
|
\end{lstlisting} |
779 |
|
|
|
780 |
|
|
\begin{lstlisting}[float,caption={Revised {\sc bass} example.},label={sch:bassExPrime}] |
781 |
|
|
|
782 |
|
|
#include "argon.mdl" |
783 |
|
|
|
784 |
|
|
molecule{ |
785 |
|
|
name = "Ar"; |
786 |
|
|
nAtoms = 1; |
787 |
|
|
atom[0]{ |
788 |
|
|
type="Ar"; |
789 |
|
|
position( 0.0, 0.0, 0.0 ); |
790 |
|
|
} |
791 |
|
|
} |
792 |
|
|
|
793 |
|
|
nComponents = 1; |
794 |
|
|
component{ |
795 |
|
|
type = "Ar"; |
796 |
|
|
nMol = 108; |
797 |
|
|
} |
798 |
|
|
|
799 |
|
|
initialConfig = "./argon.init"; |
800 |
|
|
|
801 |
|
|
forceField = "LJ"; |
802 |
|
|
ensemble = "NVE"; |
803 |
|
|
dt = 1.0; |
804 |
|
|
runTime = 1e3; |
805 |
|
|
sampleTime = 100; |
806 |
|
|
statusTime = 50; |
807 |
|
|
|
808 |
|
|
\end{lstlisting} |
809 |
|
|
|
810 |
mmeineke |
1051 |
\subsection{\label{oopseSec:coordFiles}Coordinate Files} |
811 |
mmeineke |
1044 |
|
812 |
|
|
The standard format for storage of a systems coordinates is a modified |
813 |
|
|
xyz-file syntax, the exact details of which can be seen in |
814 |
mmeineke |
1071 |
Scheme~\ref{sch:dumpFormat}. As all bonding and molecular information |
815 |
|
|
is stored in the \texttt{.bass} and \texttt{.mdl} files, the |
816 |
|
|
coordinate files are simply the complete set of coordinates for each |
817 |
|
|
atom at a given simulation time. One important note, although the |
818 |
|
|
simulation propagates the complete rotation matrix, directional |
819 |
|
|
entities are written out using quanternions, to save space in the |
820 |
|
|
output files. |
821 |
mmeineke |
1044 |
|
822 |
mmeineke |
1071 |
\begin{lstlisting}[float,caption={[The format of the coordinate files]Shows the format of the coordinate files. The fist line is the number of atoms. The second line begins with the time stamp followed by the three $\mathbf{H}$ column vectors. The next lines are the atomic coordinates for all atoms in the system. First is the name followed by position, velocity, quanternions, and lastly angular momentum.},label=sch:dumpFormat] |
823 |
|
|
|
824 |
|
|
nAtoms |
825 |
|
|
time; Hxx Hyx Hzx; Hxy Hyy Hzy; Hxz Hyz Hzz; |
826 |
|
|
Name1 x y z vx vy vz q0 q1 q2 q3 jx jy jz |
827 |
|
|
Name2 x y z vx vy vz q0 q1 q2 q3 jx jy jz |
828 |
|
|
etc... |
829 |
|
|
|
830 |
|
|
\end{lstlisting} |
831 |
|
|
|
832 |
|
|
|
833 |
|
|
There are three major files used by {\sc oopse} written in the |
834 |
|
|
coordinate format, they are as follows: the initialization file |
835 |
|
|
(\texttt{.init}), the simulation trajectory file (\texttt{.dump}), and |
836 |
|
|
the final coordinates of the simulation. The initialization file is |
837 |
|
|
necessary for {\sc oopse} to start the simulation with the proper |
838 |
|
|
coordinates, and is generated before the simulation run. The |
839 |
|
|
trajectory file is created at the beginning of the simulation, and is |
840 |
|
|
used to store snapshots of the simulation at regular intervals. The |
841 |
|
|
first frame is a duplication of the |
842 |
|
|
\texttt{.init} file, and each subsequent frame is appended to the file |
843 |
|
|
at an interval specified in the \texttt{.bass} file with the |
844 |
|
|
\texttt{sampleTime} flag. The final coordinate file is the end of run file. The |
845 |
mmeineke |
1054 |
\texttt{.eor} file stores the final configuration of the system for a |
846 |
mmeineke |
1044 |
given simulation. The file is updated at the same time as the |
847 |
mmeineke |
1071 |
\texttt{.dump} file, however, it only contains the most recent |
848 |
mmeineke |
1044 |
frame. In this way, an \texttt{.eor} file may be used as the |
849 |
mmeineke |
1071 |
initialization file to a second simulation in order to continue a |
850 |
|
|
simulation or recover one from a processor that has crashed during the |
851 |
|
|
course of the run. |
852 |
mmeineke |
1044 |
|
853 |
mmeineke |
1054 |
\subsection{\label{oopseSec:initCoords}Generation of Initial Coordinates} |
854 |
mmeineke |
1044 |
|
855 |
mmeineke |
1071 |
As was stated in Sec.~\ref{oopseSec:coordFiles}, an initialization |
856 |
|
|
file is needed to provide the starting coordinates for a |
857 |
|
|
simulation. The {\sc oopse} package provides a program called |
858 |
|
|
\texttt{sysBuilder} to aid in the creation of the \texttt{.init} |
859 |
|
|
file. \texttt{sysBuilder} uses {\sc bass}, and will recognize |
860 |
|
|
arguments and parameters in the \texttt{.bass} file that would |
861 |
|
|
otherwise be ignored by the simulation. |
862 |
mmeineke |
1044 |
|
863 |
|
|
\subsection{The Statistics File} |
864 |
|
|
|
865 |
mmeineke |
1071 |
The last output file generated by {\sc oopse} is the statistics |
866 |
|
|
file. This file records such statistical quantities as the |
867 |
|
|
instantaneous temperature, volume, pressure, etc. It is written out |
868 |
|
|
with the frequency specified in the \texttt{.bass} file with the |
869 |
|
|
\texttt{statusTime} keyword. The file allows the user to observe the |
870 |
|
|
system variables as a function of simulation time while the simulation |
871 |
|
|
is in progress. One useful function the statistics file serves is to |
872 |
|
|
monitor the conserved quantity of a given simulation ensemble, this |
873 |
|
|
allows the user to observe the stability of the integrator. The |
874 |
|
|
statistics file is denoted with the \texttt{.stat} file extension. |
875 |
mmeineke |
1044 |
|
876 |
mmeineke |
1051 |
\section{\label{oopseSec:mechanics}Mechanics} |
877 |
mmeineke |
1044 |
|
878 |
mmeineke |
1054 |
\subsection{\label{oopseSec:integrate}Integrating the Equations of Motion: the Symplectic Step Integrator} |
879 |
mmeineke |
1044 |
|
880 |
|
|
Integration of the equations of motion was carried out using the |
881 |
|
|
symplectic splitting method proposed by Dullweber \emph{et |
882 |
mmeineke |
1071 |
al.}.\cite{Dullweber1997} The reason for the selection of this |
883 |
|
|
integrator, is the poor energy conservation of rigid body systems |
884 |
|
|
using quaternion dynamics. While quaternions work well for |
885 |
|
|
orientational motion in alternate ensembles, the microcanonical |
886 |
|
|
ensemble has a constant energy requirement that is quite sensitive to |
887 |
|
|
errors in the equations of motion. The original implementation of {\sc |
888 |
|
|
oopse} utilized quaternions for rotational motion propagation; |
889 |
|
|
however, a detailed investigation showed that they resulted in a |
890 |
|
|
steady drift in the total energy, something that has been observed by |
891 |
|
|
others.\cite{Laird97} |
892 |
mmeineke |
1044 |
|
893 |
|
|
The key difference in the integration method proposed by Dullweber |
894 |
mmeineke |
1071 |
\emph{et al}.~({\sc dlm}) is that the entire rotation matrix is propagated from |
895 |
|
|
one time step to the next. In the past, this would not have been a |
896 |
|
|
feasible option, since the rotation matrix for a single body is nine |
897 |
|
|
elements long as opposed to three or four elements for Euler angles |
898 |
|
|
and quaternions respectively. System memory has become much less of an |
899 |
|
|
issue in recent times, and the {\sc dlm} method has used memory in |
900 |
|
|
exchange for substantial benefits in energy conservation. |
901 |
mmeineke |
1044 |
|
902 |
mmeineke |
1071 |
The {\sc dlm} method allows for Verlet style integration of both |
903 |
|
|
linear and angular motion of rigid bodies. In the integration method, |
904 |
|
|
the orientational propagation involves a sequence of matrix |
905 |
mmeineke |
1044 |
evaluations to update the rotation matrix.\cite{Dullweber1997} These |
906 |
mmeineke |
1071 |
matrix rotations are more costly computationally than the simpler |
907 |
|
|
arithmetic quaternion propagation. With the same time step, a 1000 SSD |
908 |
|
|
particle simulation shows an average 7\% increase in computation time |
909 |
|
|
using the {\sc dlm} method in place of quaternions. This cost is more |
910 |
|
|
than justified when comparing the energy conservation of the two |
911 |
|
|
methods as illustrated in Fig.~\ref{timestep}. |
912 |
mmeineke |
1044 |
|
913 |
|
|
\begin{figure} |
914 |
mmeineke |
1045 |
\centering |
915 |
|
|
\includegraphics[width=\linewidth]{timeStep.eps} |
916 |
mmeineke |
1071 |
\caption[Energy conservation for quaternion versus {\sc dlm} dynamics]{Energy conservation using quaternion based integration versus |
917 |
|
|
the {\sc dlm} method with |
918 |
mmeineke |
1044 |
increasing time step. For each time step, the dotted line is total |
919 |
mmeineke |
1071 |
energy using the {\sc dlm} integrator, and the solid line comes |
920 |
mmeineke |
1044 |
from the quaternion integrator. The larger time step plots are shifted |
921 |
|
|
up from the true energy baseline for clarity.} |
922 |
|
|
\label{timestep} |
923 |
|
|
\end{figure} |
924 |
|
|
|
925 |
mmeineke |
1071 |
In Fig.~\ref{timestep}, the resulting energy drift at various time |
926 |
|
|
steps for both the {\sc dlm} and quaternion integration schemes |
927 |
mmeineke |
1044 |
is compared. All of the 1000 SSD particle simulations started with the |
928 |
|
|
same configuration, and the only difference was the method for |
929 |
|
|
handling rotational motion. At time steps of 0.1 and 0.5 fs, both |
930 |
|
|
methods for propagating particle rotation conserve energy fairly well, |
931 |
|
|
with the quaternion method showing a slight energy drift over time in |
932 |
|
|
the 0.5 fs time step simulation. At time steps of 1 and 2 fs, the |
933 |
mmeineke |
1071 |
energy conservation benefits of the {\sc dlm} method are clearly |
934 |
mmeineke |
1044 |
demonstrated. Thus, while maintaining the same degree of energy |
935 |
|
|
conservation, one can take considerably longer time steps, leading to |
936 |
|
|
an overall reduction in computation time. |
937 |
|
|
|
938 |
|
|
Energy drift in these SSD particle simulations was unnoticeable for |
939 |
|
|
time steps up to three femtoseconds. A slight energy drift on the |
940 |
|
|
order of 0.012 kcal/mol per nanosecond was observed at a time step of |
941 |
|
|
four femtoseconds, and as expected, this drift increases dramatically |
942 |
mmeineke |
1071 |
with increasing time step. |
943 |
mmeineke |
1044 |
|
944 |
|
|
|
945 |
|
|
\subsection{\label{sec:extended}Extended Systems for other Ensembles} |
946 |
|
|
|
947 |
|
|
|
948 |
|
|
{\sc oopse} implements a |
949 |
|
|
|
950 |
|
|
|
951 |
mmeineke |
1071 |
\subsection{\label{oopseSec:noseHooverThermo}Nose-Hoover Thermostatting} |
952 |
mmeineke |
1044 |
|
953 |
|
|
To mimic the effects of being in a constant temperature ({\sc nvt}) |
954 |
|
|
ensemble, {\sc oopse} uses the Nose-Hoover extended system |
955 |
|
|
approach.\cite{Hoover85} In this method, the equations of motion for |
956 |
|
|
the particle positions and velocities are |
957 |
|
|
\begin{eqnarray} |
958 |
|
|
\dot{{\bf r}} & = & {\bf v} \\ |
959 |
|
|
\dot{{\bf v}} & = & \frac{{\bf f}}{m} - \chi {\bf v} |
960 |
|
|
\label{eq:nosehoovereom} |
961 |
|
|
\end{eqnarray} |
962 |
|
|
|
963 |
|
|
$\chi$ is an ``extra'' variable included in the extended system, and |
964 |
|
|
it is propagated using the first order equation of motion |
965 |
|
|
\begin{equation} |
966 |
|
|
\dot{\chi} = \frac{1}{\tau_{T}} \left( \frac{T}{T_{target}} - 1 \right) |
967 |
|
|
\label{eq:nosehooverext} |
968 |
|
|
\end{equation} |
969 |
|
|
where $T_{target}$ is the target temperature for the simulation, and |
970 |
|
|
$\tau_{T}$ is a time constant for the thermostat. |
971 |
|
|
|
972 |
|
|
To select the Nose-Hoover {\sc nvt} ensemble, the {\tt ensemble = NVT;} |
973 |
|
|
command would be used in the simulation's {\sc bass} file. There is |
974 |
|
|
some subtlety in choosing values for $\tau_{T}$, and it is usually set |
975 |
|
|
to values of a few ps. Within a {\sc bass} file, $\tau_{T}$ could be |
976 |
|
|
set to 1 ps using the {\tt tauThermostat = 1000; } command. |
977 |
|
|
|
978 |
mmeineke |
1071 |
\subsection{\label{oopseSec:rattle}The {\sc rattle} Method for Bond |
979 |
|
|
Constraints} |
980 |
mmeineke |
1044 |
|
981 |
mmeineke |
1071 |
In order to satisfy the constraints of fixed bond lengths within {\sc |
982 |
|
|
oopse}, we have implemented the {\sc rattle} algorithm of |
983 |
|
|
Andersen.\cite{andersen83} The algorithm is a velocity verlet |
984 |
|
|
formulation of the {\sc shake} method\cite{ryckaert77} of iteratively |
985 |
|
|
solving the Lagrange multipliers of constraint. The system of lagrange |
986 |
|
|
multipliers allows one to reformulate the equations of motion with |
987 |
|
|
explicit constraint forces on the equations of |
988 |
|
|
motion.\cite{fowles99:lagrange} |
989 |
|
|
|
990 |
|
|
Consider a system described by qoordinates $q_1$ and $q_2$ subject to an |
991 |
|
|
equation of constraint: |
992 |
|
|
\begin{equation} |
993 |
|
|
\sigma(q_1, q_2,t) = 0 |
994 |
|
|
\label{oopseEq:lm1} |
995 |
|
|
\end{equation} |
996 |
|
|
The Lagrange formulation of the equations of motion can be written: |
997 |
|
|
\begin{equation} |
998 |
|
|
\delta\int_{t_1}^{t_2}L\, dt = |
999 |
|
|
\int_{t_1}^{t_2} \sum_i \biggl [ \frac{\partial L}{\partial q_i} |
1000 |
|
|
- \frac{d}{dt}\biggl(\frac{\partial L}{\partial \dot{q}_i} |
1001 |
|
|
\biggr ) \biggr] \delta q_i \, dt = 0 |
1002 |
|
|
\label{oopseEq:lm2} |
1003 |
|
|
\end{equation} |
1004 |
|
|
Here, $\delta q_i$ is not independent for each $q$, as $q_1$ and $q_2$ |
1005 |
|
|
are linked by $\sigma$. However, $\sigma$ is fixed at any given |
1006 |
|
|
instant of time, giving: |
1007 |
|
|
\begin{align} |
1008 |
|
|
\delta\sigma &= \biggl( \frac{\partial\sigma}{\partial q_1} \delta q_1 % |
1009 |
|
|
+ \frac{\partial\sigma}{\partial q_2} \delta q_2 \biggr) = 0 \\ |
1010 |
|
|
% |
1011 |
|
|
\frac{\partial\sigma}{\partial q_1} \delta q_1 &= % |
1012 |
|
|
- \frac{\partial\sigma}{\partial q_2} \delta q_2 \\ |
1013 |
|
|
% |
1014 |
|
|
\delta q_2 &= - \biggl(\frac{\partial\sigma}{\partial q_1} \bigg / % |
1015 |
|
|
\frac{\partial\sigma}{\partial q_2} \biggr) \delta q_1 |
1016 |
|
|
\end{align} |
1017 |
|
|
Substituted back into Eq.~\ref{oopseEq:lm2}, |
1018 |
|
|
\begin{equation} |
1019 |
|
|
\int_{t_1}^{t_2}\biggl [ \biggl(\frac{\partial L}{\partial q_1} |
1020 |
|
|
- \frac{d}{dt}\,\frac{\partial L}{\partial \dot{q}_1} |
1021 |
|
|
\biggr) |
1022 |
|
|
- \biggl( \frac{\partial L}{\partial q_1} |
1023 |
|
|
- \frac{d}{dt}\,\frac{\partial L}{\partial \dot{q}_1} |
1024 |
|
|
\biggr) \biggl(\frac{\partial\sigma}{\partial q_1} \bigg / % |
1025 |
|
|
\frac{\partial\sigma}{\partial q_2} \biggr)\biggr] \delta q_1 \, dt = 0 |
1026 |
|
|
\label{oopseEq:lm3} |
1027 |
|
|
\end{equation} |
1028 |
|
|
Leading to, |
1029 |
|
|
\begin{equation} |
1030 |
|
|
\frac{\biggl(\frac{\partial L}{\partial q_1} |
1031 |
|
|
- \frac{d}{dt}\,\frac{\partial L}{\partial \dot{q}_1} |
1032 |
|
|
\biggr)}{\frac{\partial\sigma}{\partial q_1}} = |
1033 |
|
|
\frac{\biggl(\frac{\partial L}{\partial q_2} |
1034 |
|
|
- \frac{d}{dt}\,\frac{\partial L}{\partial \dot{q}_2} |
1035 |
|
|
\biggr)}{\frac{\partial\sigma}{\partial q_2}} |
1036 |
|
|
\label{oopseEq:lm4} |
1037 |
|
|
\end{equation} |
1038 |
|
|
This relation can only be statisfied, if both are equal to a single |
1039 |
|
|
function $-\lambda(t)$, |
1040 |
|
|
\begin{align} |
1041 |
|
|
\frac{\biggl(\frac{\partial L}{\partial q_1} |
1042 |
|
|
- \frac{d}{dt}\,\frac{\partial L}{\partial \dot{q}_1} |
1043 |
|
|
\biggr)}{\frac{\partial\sigma}{\partial q_1}} &= -\lambda(t) \\ |
1044 |
|
|
% |
1045 |
|
|
\frac{\partial L}{\partial q_1} |
1046 |
|
|
- \frac{d}{dt}\,\frac{\partial L}{\partial \dot{q}_1} &= |
1047 |
|
|
-\lambda(t)\,\frac{\partial\sigma}{\partial q_1} \\ |
1048 |
|
|
% |
1049 |
|
|
\frac{\partial L}{\partial q_1} |
1050 |
|
|
- \frac{d}{dt}\,\frac{\partial L}{\partial \dot{q}_1} |
1051 |
|
|
+ \mathcal{G}_i &= 0 |
1052 |
|
|
\end{align} |
1053 |
|
|
Where $\mathcal{G}_i$, the force of constraint on $i$, is: |
1054 |
|
|
\begin{equation} |
1055 |
|
|
\mathcal{G}_i = \lambda(t)\,\frac{\partial\sigma}{\partial q_1} |
1056 |
|
|
\label{oopseEq:lm5} |
1057 |
|
|
\end{equation} |
1058 |
|
|
|
1059 |
|
|
In a simulation, this would involve the solution of a set of $(m + n)$ |
1060 |
|
|
number of equations. Where $m$ is the number of constraints, and $n$ |
1061 |
|
|
is the number of constrained coordinates. In practice, this is not |
1062 |
|
|
done, as the matrix inversion neccassary to solve the system of |
1063 |
|
|
equations would be very time consuming to solve. Additionally, the |
1064 |
|
|
numerical error in the solution of the set of $\lambda$'s would be |
1065 |
|
|
compounded by the error inherent in propagating by the Velocity Verlet |
1066 |
|
|
algorithm ($\Delta t^4$). The verlet propagation error is negligible |
1067 |
|
|
in an unconstrained system, as one is interested in the statisitics of |
1068 |
|
|
the run, and not that the run be numerically exact to the ``true'' |
1069 |
|
|
integration. This relates back to the ergodic hypothesis that a time |
1070 |
|
|
integral of a valid trajectory will still give the correct enesemble |
1071 |
|
|
average. However, in the case of constraints, if the equations of |
1072 |
|
|
motion leave the ``true'' trajectory, they are departing from the |
1073 |
|
|
constrained surface. The method that is used, is to iteratively solve |
1074 |
|
|
for $\lambda(t)$ at each time step. |
1075 |
|
|
|
1076 |
|
|
In {\sc rattle} the equations of motion are modified subject to the |
1077 |
|
|
following two constraints: |
1078 |
|
|
\begin{align} |
1079 |
|
|
\sigma_{ij}[\mathbf{r}(t)] \equiv |
1080 |
|
|
[ \mathbf{r}_i(t) - \mathbf{r}_j(t)]^2 - d_{ij}^2 &= 0 % |
1081 |
|
|
\label{oopseEq:c1} \\ |
1082 |
|
|
% |
1083 |
|
|
[\mathbf{\dot{r}}_i(t) - \mathbf{\dot{r}}_j(t)] \cdot |
1084 |
|
|
[\mathbf{r}_i(t) - \mathbf{r}_j(t)] &= 0 \label{oopseEq:c2} |
1085 |
|
|
\end{align} |
1086 |
|
|
Eq.~\ref{oopseEq:c1} is the set of bond constraints, where $d_{ij}$ is |
1087 |
|
|
the constrained distance between atom $i$ and |
1088 |
|
|
$j$. Eq.~\ref{oopseEq:c2} constrains the velocities of $i$ and $j$ to |
1089 |
|
|
be perpindicular to the bond vector, so that the bond can neither grow |
1090 |
|
|
nor shrink. The constrained dynamics equations become: |
1091 |
|
|
\begin{equation} |
1092 |
|
|
m_i \mathbf{\ddot{r}}_i = \mathbf{F}_i + \mathbf{\mathcal{G}}_i |
1093 |
|
|
\label{oopseEq:r1} |
1094 |
|
|
\end{equation} |
1095 |
|
|
Where, |
1096 |
|
|
\begin{equation} |
1097 |
|
|
\mathbf{\mathcal{G}}_i = - \sum_j \lambda_{ij}(t)\,\nabla \sigma_{ij} |
1098 |
|
|
\label{oopseEq:r2} |
1099 |
|
|
\end{equation} |
1100 |
|
|
|
1101 |
|
|
In Velocity Verlet, if $\Delta t = h$, the propagation can be written: |
1102 |
|
|
\begin{align} |
1103 |
|
|
\mathbf{r}_i(t+h) &= |
1104 |
|
|
\mathbf{r}_i(t) + h\mathbf{\dot{r}}(t) + |
1105 |
|
|
\frac{h^2}{2m_i}\,\Bigl[ \mathbf{F}_i(t) + |
1106 |
|
|
\mathbf{\mathcal{G}}_{Ri}(t) \Bigr] \label{oopseEq:vv1} \\ |
1107 |
|
|
% |
1108 |
|
|
\mathbf{\dot{r}}_i(t+h) &= |
1109 |
|
|
\mathbf{\dot{r}}_i(t) + \frac{h}{2m_i} |
1110 |
|
|
\Bigl[ \mathbf{F}_i(t) + \mathbf{\mathcal{G}}_{Ri}(t) + |
1111 |
|
|
\mathbf{F}_i(t+h) + \mathbf{\mathcal{G}}_{Vi}(t+h) \Bigr] % |
1112 |
|
|
\label{oopseEq:vv2} |
1113 |
|
|
\end{align} |
1114 |
|
|
|
1115 |
|
|
|
1116 |
|
|
|
1117 |
mmeineke |
1054 |
\subsection{\label{oopseSec:zcons}Z-Constraint Method} |
1118 |
mmeineke |
1044 |
|
1119 |
mmeineke |
1071 |
Based on fluctuation-dissipation theorem, a force auto-correlation |
1120 |
mmeineke |
1044 |
method was developed to investigate the dynamics of ions inside the ion |
1121 |
|
|
channels.\cite{Roux91} Time-dependent friction coefficient can be calculated |
1122 |
mmeineke |
1054 |
from the deviation of the instantaneous force from its mean force. |
1123 |
mmeineke |
1044 |
|
1124 |
|
|
% |
1125 |
|
|
|
1126 |
|
|
\begin{equation} |
1127 |
|
|
\xi(z,t)=\langle\delta F(z,t)\delta F(z,0)\rangle/k_{B}T |
1128 |
|
|
\end{equation} |
1129 |
|
|
where% |
1130 |
|
|
\begin{equation} |
1131 |
|
|
\delta F(z,t)=F(z,t)-\langle F(z,t)\rangle |
1132 |
|
|
\end{equation} |
1133 |
|
|
|
1134 |
|
|
|
1135 |
|
|
If the time-dependent friction decay rapidly, static friction coefficient can |
1136 |
|
|
be approximated by% |
1137 |
|
|
|
1138 |
|
|
\begin{equation} |
1139 |
|
|
\xi^{static}(z)=\int_{0}^{\infty}\langle\delta F(z,t)\delta F(z,0)\rangle dt |
1140 |
|
|
\end{equation} |
1141 |
|
|
|
1142 |
|
|
|
1143 |
|
|
Hence, diffusion constant can be estimated by |
1144 |
|
|
\begin{equation} |
1145 |
|
|
D(z)=\frac{k_{B}T}{\xi^{static}(z)}=\frac{(k_{B}T)^{2}}{\int_{0}^{\infty |
1146 |
|
|
}\langle\delta F(z,t)\delta F(z,0)\rangle dt}% |
1147 |
|
|
\end{equation} |
1148 |
|
|
|
1149 |
|
|
|
1150 |
|
|
\bigskip Z-Constraint method, which fixed the z coordinates of the molecules |
1151 |
|
|
with respect to the center of the mass of the system, was proposed to obtain |
1152 |
|
|
the forces required in force auto-correlation method.\cite{Marrink94} However, |
1153 |
|
|
simply resetting the coordinate will move the center of the mass of the whole |
1154 |
|
|
system. To avoid this problem, a new method was used at {\sc oopse}. Instead of |
1155 |
|
|
resetting the coordinate, we reset the forces of z-constraint molecules as |
1156 |
|
|
well as subtract the total constraint forces from the rest of the system after |
1157 |
|
|
force calculation at each time step. |
1158 |
mmeineke |
1054 |
\begin{align} |
1159 |
|
|
F_{\alpha i}&=0\\ |
1160 |
|
|
V_{\alpha i}&=V_{\alpha i}-\frac{\sum\limits_{i}M_{_{\alpha i}}V_{\alpha i}}{\sum\limits_{i}M_{_{\alpha i}}}\\ |
1161 |
|
|
F_{\alpha i}&=F_{\alpha i}-\frac{M_{_{\alpha i}}}{\sum\limits_{\alpha}\sum\limits_{i}M_{_{\alpha i}}}\sum\limits_{\beta}F_{\beta}\\ |
1162 |
|
|
V_{\alpha i}&=V_{\alpha i}-\frac{\sum\limits_{\alpha}\sum\limits_{i}M_{_{\alpha i}}V_{\alpha i}}{\sum\limits_{\alpha}\sum\limits_{i}M_{_{\alpha i}}} |
1163 |
|
|
\end{align} |
1164 |
mmeineke |
1044 |
|
1165 |
|
|
At the very beginning of the simulation, the molecules may not be at its |
1166 |
|
|
constraint position. To move the z-constraint molecule to the specified |
1167 |
|
|
position, a simple harmonic potential is used% |
1168 |
|
|
|
1169 |
|
|
\begin{equation} |
1170 |
|
|
U(t)=\frac{1}{2}k_{Harmonic}(z(t)-z_{cons})^{2}% |
1171 |
|
|
\end{equation} |
1172 |
|
|
where $k_{Harmonic}$\bigskip\ is the harmonic force constant, $z(t)$ is |
1173 |
|
|
current z coordinate of the center of mass of the z-constraint molecule, and |
1174 |
|
|
$z_{cons}$ is the restraint position. Therefore, the harmonic force operated |
1175 |
|
|
on the z-constraint molecule at time $t$ can be calculated by% |
1176 |
|
|
\begin{equation} |
1177 |
|
|
F_{z_{Harmonic}}(t)=-\frac{\partial U(t)}{\partial z(t)}=-k_{Harmonic}% |
1178 |
|
|
(z(t)-z_{cons}) |
1179 |
|
|
\end{equation} |
1180 |
|
|
Worthy of mention, other kinds of potential functions can also be used to |
1181 |
|
|
drive the z-constraint molecule. |
1182 |
|
|
|
1183 |
mmeineke |
1051 |
\section{\label{oopseSec:props}Trajectory Analysis} |
1184 |
mmeineke |
1044 |
|
1185 |
mmeineke |
1051 |
\subsection{\label{oopseSec:staticProps}Static Property Analysis} |
1186 |
mmeineke |
1044 |
|
1187 |
|
|
The static properties of the trajectories are analyzed with the |
1188 |
|
|
program \texttt{staticProps}. The code is capable of calculating the following |
1189 |
|
|
pair correlations between species A and B: |
1190 |
|
|
\begin{itemize} |
1191 |
|
|
\item $g_{\text{AB}}(r)$: Eq.~\ref{eq:gofr} |
1192 |
|
|
\item $g_{\text{AB}}(r, \cos \theta)$: Eq.~\ref{eq:gofrCosTheta} |
1193 |
|
|
\item $g_{\text{AB}}(r, \cos \omega)$: Eq.~\ref{eq:gofrCosOmega} |
1194 |
|
|
\item $g_{\text{AB}}(x, y, z)$: Eq.~\ref{eq:gofrXYZ} |
1195 |
|
|
\item $\langle \cos \omega \rangle_{\text{AB}}(r)$: |
1196 |
|
|
Eq.~\ref{eq:cosOmegaOfR} |
1197 |
|
|
\end{itemize} |
1198 |
|
|
|
1199 |
|
|
The first pair correlation, $g_{\text{AB}}(r)$, is defined as follows: |
1200 |
|
|
\begin{equation} |
1201 |
|
|
g_{\text{AB}}(r) = \frac{V}{N_{\text{A}}N_{\text{B}}}\langle %% |
1202 |
|
|
\sum_{i \in \text{A}} \sum_{j \in \text{B}} %% |
1203 |
|
|
\delta( r - |\mathbf{r}_{ij}|) \rangle \label{eq:gofr} |
1204 |
|
|
\end{equation} |
1205 |
|
|
Where $\mathbf{r}_{ij}$ is the vector |
1206 |
|
|
\begin{equation*} |
1207 |
|
|
\mathbf{r}_{ij} = \mathbf{r}_j - \mathbf{r}_i \notag |
1208 |
|
|
\end{equation*} |
1209 |
|
|
and $\frac{V}{N_{\text{A}}N_{\text{B}}}$ normalizes the average over |
1210 |
|
|
the expected pair density at a given $r$. |
1211 |
|
|
|
1212 |
|
|
The next two pair correlations, $g_{\text{AB}}(r, \cos \theta)$ and |
1213 |
|
|
$g_{\text{AB}}(r, \cos \omega)$, are similar in that they are both two |
1214 |
|
|
dimensional histograms. Both use $r$ for the primary axis then a |
1215 |
|
|
$\cos$ for the secondary axis ($\cos \theta$ for |
1216 |
|
|
Eq.~\ref{eq:gofrCosTheta} and $\cos \omega$ for |
1217 |
|
|
Eq.~\ref{eq:gofrCosOmega}). This allows for the investigator to |
1218 |
|
|
correlate alignment on directional entities. $g_{\text{AB}}(r, \cos |
1219 |
|
|
\theta)$ is defined as follows: |
1220 |
|
|
\begin{equation} |
1221 |
|
|
g_{\text{AB}}(r, \cos \theta) = \frac{V}{N_{\text{A}}N_{\text{B}}}\langle |
1222 |
|
|
\sum_{i \in \text{A}} \sum_{j \in \text{B}} |
1223 |
|
|
\delta( \cos \theta - \cos \theta_{ij}) |
1224 |
|
|
\delta( r - |\mathbf{r}_{ij}|) \rangle |
1225 |
|
|
\label{eq:gofrCosTheta} |
1226 |
|
|
\end{equation} |
1227 |
|
|
Where |
1228 |
|
|
\begin{equation*} |
1229 |
|
|
\cos \theta_{ij} = \mathbf{\hat{i}} \cdot \mathbf{\hat{r}}_{ij} |
1230 |
|
|
\end{equation*} |
1231 |
|
|
Here $\mathbf{\hat{i}}$ is the unit directional vector of species $i$ |
1232 |
|
|
and $\mathbf{\hat{r}}_{ij}$ is the unit vector associated with vector |
1233 |
|
|
$\mathbf{r}_{ij}$. |
1234 |
|
|
|
1235 |
|
|
The second two dimensional histogram is of the form: |
1236 |
|
|
\begin{equation} |
1237 |
|
|
g_{\text{AB}}(r, \cos \omega) = |
1238 |
|
|
\frac{V}{N_{\text{A}}N_{\text{B}}}\langle |
1239 |
|
|
\sum_{i \in \text{A}} \sum_{j \in \text{B}} |
1240 |
|
|
\delta( \cos \omega - \cos \omega_{ij}) |
1241 |
|
|
\delta( r - |\mathbf{r}_{ij}|) \rangle \label{eq:gofrCosOmega} |
1242 |
|
|
\end{equation} |
1243 |
|
|
Here |
1244 |
|
|
\begin{equation*} |
1245 |
|
|
\cos \omega_{ij} = \mathbf{\hat{i}} \cdot \mathbf{\hat{j}} |
1246 |
|
|
\end{equation*} |
1247 |
|
|
Again, $\mathbf{\hat{i}}$ and $\mathbf{\hat{j}}$ are the unit |
1248 |
|
|
directional vectors of species $i$ and $j$. |
1249 |
|
|
|
1250 |
|
|
The static analysis code is also cable of calculating a three |
1251 |
|
|
dimensional pair correlation of the form: |
1252 |
|
|
\begin{equation}\label{eq:gofrXYZ} |
1253 |
|
|
g_{\text{AB}}(x, y, z) = |
1254 |
|
|
\frac{V}{N_{\text{A}}N_{\text{B}}}\langle |
1255 |
|
|
\sum_{i \in \text{A}} \sum_{j \in \text{B}} |
1256 |
|
|
\delta( x - x_{ij}) |
1257 |
|
|
\delta( y - y_{ij}) |
1258 |
|
|
\delta( z - z_{ij}) \rangle |
1259 |
|
|
\end{equation} |
1260 |
|
|
Where $x_{ij}$, $y_{ij}$, and $z_{ij}$ are the $x$, $y$, and $z$ |
1261 |
|
|
components respectively of vector $\mathbf{r}_{ij}$. |
1262 |
|
|
|
1263 |
|
|
The final pair correlation is similar to |
1264 |
|
|
Eq.~\ref{eq:gofrCosOmega}. $\langle \cos \omega |
1265 |
|
|
\rangle_{\text{AB}}(r)$ is calculated in the following way: |
1266 |
|
|
\begin{equation}\label{eq:cosOmegaOfR} |
1267 |
|
|
\langle \cos \omega \rangle_{\text{AB}}(r) = |
1268 |
|
|
\langle \sum_{i \in \text{A}} \sum_{j \in \text{B}} |
1269 |
|
|
(\cos \omega_{ij}) \delta( r - |\mathbf{r}_{ij}|) \rangle |
1270 |
|
|
\end{equation} |
1271 |
|
|
Here $\cos \omega_{ij}$ is defined in the same way as in |
1272 |
|
|
Eq.~\ref{eq:gofrCosOmega}. This equation is a single dimensional pair |
1273 |
|
|
correlation that gives the average correlation of two directional |
1274 |
|
|
entities as a function of their distance from each other. |
1275 |
|
|
|
1276 |
|
|
All static properties are calculated on a frame by frame basis. The |
1277 |
|
|
trajectory is read a single frame at a time, and the appropriate |
1278 |
|
|
calculations are done on each frame. Once one frame is finished, the |
1279 |
|
|
next frame is read in, and a running average of the property being |
1280 |
|
|
calculated is accumulated in each frame. The program allows for the |
1281 |
|
|
user to specify more than one property be calculated in single run, |
1282 |
|
|
preventing the need to read a file multiple times. |
1283 |
|
|
|
1284 |
|
|
\subsection{\label{dynamicProps}Dynamic Property Analysis} |
1285 |
|
|
|
1286 |
|
|
The dynamic properties of a trajectory are calculated with the program |
1287 |
|
|
\texttt{dynamicProps}. The program will calculate the following properties: |
1288 |
|
|
\begin{gather} |
1289 |
|
|
\langle | \mathbf{r}(t) - \mathbf{r}(0) |^2 \rangle \label{eq:rms}\\ |
1290 |
|
|
\langle \mathbf{v}(t) \cdot \mathbf{v}(0) \rangle \label{eq:velCorr} \\ |
1291 |
|
|
\langle \mathbf{j}(t) \cdot \mathbf{j}(0) \rangle \label{eq:angularVelCorr} |
1292 |
|
|
\end{gather} |
1293 |
|
|
|
1294 |
|
|
Eq.~\ref{eq:rms} is the root mean square displacement |
1295 |
|
|
function. Eq.~\ref{eq:velCorr} and Eq.~\ref{eq:angularVelCorr} are the |
1296 |
|
|
velocity and angular velocity correlation functions respectively. The |
1297 |
|
|
latter is only applicable to directional species in the simulation. |
1298 |
|
|
|
1299 |
|
|
The \texttt{dynamicProps} program handles he file in a manner different from |
1300 |
|
|
\texttt{staticProps}. As the properties calculated by this program are time |
1301 |
|
|
dependent, multiple frames must be read in simultaneously by the |
1302 |
|
|
program. For small trajectories this is no problem, and the entire |
1303 |
|
|
trajectory is read into memory. However, for long trajectories of |
1304 |
|
|
large systems, the files can be quite large. In order to accommodate |
1305 |
|
|
large files, \texttt{dynamicProps} adopts a scheme whereby two blocks of memory |
1306 |
|
|
are allocated to read in several frames each. |
1307 |
|
|
|
1308 |
|
|
In this two block scheme, the correlation functions are first |
1309 |
|
|
calculated within each memory block, then the cross correlations |
1310 |
|
|
between the frames contained within the two blocks are |
1311 |
|
|
calculated. Once completed, the memory blocks are incremented, and the |
1312 |
|
|
process is repeated. A diagram illustrating the process is shown in |
1313 |
mmeineke |
1054 |
Fig.~\ref{oopseFig:dynamicPropsMemory}. As was the case with |
1314 |
|
|
\texttt{staticProps}, multiple properties may be calculated in a |
1315 |
|
|
single run to avoid multiple reads on the same file. |
1316 |
mmeineke |
1044 |
|
1317 |
|
|
|
1318 |
mmeineke |
1054 |
|
1319 |
mmeineke |
1051 |
\section{\label{oopseSec:design}Program Design} |
1320 |
mmeineke |
1044 |
|
1321 |
mmeineke |
1054 |
\subsection{\label{sec:architecture} {\sc oopse} Architecture} |
1322 |
mmeineke |
1044 |
|
1323 |
mmeineke |
1054 |
The core of OOPSE is divided into two main object libraries: |
1324 |
|
|
\texttt{libBASS} and \texttt{libmdtools}. \texttt{libBASS} is the |
1325 |
|
|
library developed around the parsing engine and \texttt{libmdtools} |
1326 |
|
|
is the software library developed around the simulation engine. These |
1327 |
|
|
two libraries are designed to encompass all the basic functions and |
1328 |
|
|
tools that {\sc oopse} provides. Utility programs, such as the |
1329 |
|
|
property analyzers, need only link against the software libraries to |
1330 |
|
|
gain access to parsing, force evaluation, and input / output |
1331 |
|
|
routines. |
1332 |
mmeineke |
1044 |
|
1333 |
mmeineke |
1054 |
Contained in \texttt{libBASS} are all the routines associated with |
1334 |
|
|
reading and parsing the \texttt{.bass} input files. Given a |
1335 |
|
|
\texttt{.bass} file, \texttt{libBASS} will open it and any associated |
1336 |
|
|
\texttt{.mdl} files; then create structures in memory that are |
1337 |
|
|
templates of all the molecules specified in the input files. In |
1338 |
|
|
addition, any simulation parameters set in the \texttt{.bass} file |
1339 |
|
|
will be placed in a structure for later query by the controlling |
1340 |
|
|
program. |
1341 |
mmeineke |
1044 |
|
1342 |
mmeineke |
1054 |
Located in \texttt{libmdtools} are all other routines necessary to a |
1343 |
|
|
Molecular Dynamics simulation. The library uses the main data |
1344 |
|
|
structures returned by \texttt{libBASS} to initialize the various |
1345 |
|
|
parts of the simulation: the atom structures and positions, the force |
1346 |
|
|
field, the integrator, \emph{et cetera}. After initialization, the |
1347 |
|
|
library can be used to perform a variety of tasks: integrate a |
1348 |
|
|
Molecular Dynamics trajectory, query phase space information from a |
1349 |
|
|
specific frame of a completed trajectory, or even recalculate force or |
1350 |
|
|
energetic information about specific frames from a completed |
1351 |
|
|
trajectory. |
1352 |
mmeineke |
1044 |
|
1353 |
mmeineke |
1054 |
With these core libraries in place, several programs have been |
1354 |
|
|
developed to utilize the routines provided by \texttt{libBASS} and |
1355 |
|
|
\texttt{libmdtools}. The main program of the package is \texttt{oopse} |
1356 |
|
|
and the corresponding parallel version \texttt{oopse\_MPI}. These two |
1357 |
|
|
programs will take the \texttt{.bass} file, and create then integrate |
1358 |
|
|
the simulation specified in the script. The two analysis programs |
1359 |
|
|
\texttt{staticProps} and \texttt{dynamicProps} utilize the core |
1360 |
|
|
libraries to initialize and read in trajectories from previously |
1361 |
|
|
completed simulations, in addition to the ability to use functionality |
1362 |
|
|
from \texttt{libmdtools} to recalculate forces and energies at key |
1363 |
|
|
frames in the trajectories. Lastly, the family of system building |
1364 |
|
|
programs (Sec.~\ref{oopseSec:initCoords}) also use the libraries to |
1365 |
|
|
store and output the system configurations they create. |
1366 |
|
|
|
1367 |
|
|
\subsection{\label{oopseSec:parallelization} Parallelization of {\sc oopse}} |
1368 |
|
|
|
1369 |
|
|
Although processor power is continually growing month by month, it is |
1370 |
|
|
still unreasonable to simulate systems of more then a 1000 atoms on a |
1371 |
|
|
single processor. To facilitate study of larger system sizes or |
1372 |
|
|
smaller systems on long time scales in a reasonable period of time, |
1373 |
|
|
parallel methods were developed allowing multiple CPU's to share the |
1374 |
|
|
simulation workload. Three general categories of parallel |
1375 |
|
|
decomposition method's have been developed including atomic, spatial |
1376 |
|
|
and force decomposition methods. |
1377 |
|
|
|
1378 |
mmeineke |
1044 |
Algorithmically simplest of the three method's is atomic decomposition |
1379 |
|
|
where N particles in a simulation are split among P processors for the |
1380 |
|
|
duration of the simulation. Computational cost scales as an optimal |
1381 |
|
|
$O(N/P)$ for atomic decomposition. Unfortunately all processors must |
1382 |
|
|
communicate positions and forces with all other processors leading |
1383 |
|
|
communication to scale as an unfavorable $O(N)$ independent of the |
1384 |
|
|
number of processors. This communication bottleneck led to the |
1385 |
|
|
development of spatial and force decomposition methods in which |
1386 |
|
|
communication among processors scales much more favorably. Spatial or |
1387 |
|
|
domain decomposition divides the physical spatial domain into 3D boxes |
1388 |
|
|
in which each processor is responsible for calculation of forces and |
1389 |
|
|
positions of particles located in its box. Particles are reassigned to |
1390 |
|
|
different processors as they move through simulation space. To |
1391 |
|
|
calculate forces on a given particle, a processor must know the |
1392 |
|
|
positions of particles within some cutoff radius located on nearby |
1393 |
|
|
processors instead of the positions of particles on all |
1394 |
|
|
processors. Both communication between processors and computation |
1395 |
|
|
scale as $O(N/P)$ in the spatial method. However, spatial |
1396 |
|
|
decomposition adds algorithmic complexity to the simulation code and |
1397 |
|
|
is not very efficient for small N since the overall communication |
1398 |
|
|
scales as the surface to volume ratio $(N/P)^{2/3}$ in three |
1399 |
|
|
dimensions. |
1400 |
|
|
|
1401 |
|
|
Force decomposition assigns particles to processors based on a block |
1402 |
|
|
decomposition of the force matrix. Processors are split into a |
1403 |
|
|
optimally square grid forming row and column processor groups. Forces |
1404 |
|
|
are calculated on particles in a given row by particles located in |
1405 |
|
|
that processors column assignment. Force decomposition is less complex |
1406 |
|
|
to implement then the spatial method but still scales computationally |
1407 |
|
|
as $O(N/P)$ and scales as $(N/\sqrt{p})$ in communication |
1408 |
|
|
cost. Plimpton also found that force decompositions scales more |
1409 |
|
|
favorably then spatial decomposition up to 10,000 atoms and favorably |
1410 |
|
|
competes with spatial methods for up to 100,000 atoms. |
1411 |
|
|
|
1412 |
mmeineke |
1054 |
\subsection{\label{oopseSec:memAlloc}Memory Issues in Trajectory Analysis} |
1413 |
mmeineke |
1044 |
|
1414 |
mmeineke |
1054 |
For large simulations, the trajectory files can sometimes reach sizes |
1415 |
|
|
in excess of several gigabytes. In order to effectively analyze that |
1416 |
|
|
amount of data+, two memory management schemes have been devised for |
1417 |
|
|
\texttt{staticProps} and for \texttt{dynamicProps}. The first scheme, |
1418 |
|
|
developed for \texttt{staticProps}, is the simplest. As each frame's |
1419 |
|
|
statistics are calculated independent of each other, memory is |
1420 |
|
|
allocated for each frame, then freed once correlation calculations are |
1421 |
|
|
complete for the snapshot. To prevent multiple passes through a |
1422 |
|
|
potentially large file, \texttt{staticProps} is capable of calculating |
1423 |
|
|
all requested correlations per frame with only a single pair loop in |
1424 |
|
|
each frame and a single read through of the file. |
1425 |
mmeineke |
1044 |
|
1426 |
mmeineke |
1054 |
The second, more advanced memory scheme, is used by |
1427 |
|
|
\texttt{dynamicProps}. Here, the program must have multiple frames in |
1428 |
|
|
memory to calculate time dependent correlations. In order to prevent a |
1429 |
|
|
situation where the program runs out of memory due to large |
1430 |
|
|
trajectories, the user is able to specify that the trajectory be read |
1431 |
|
|
in blocks. The number of frames in each block is specified by the |
1432 |
|
|
user, and upon reading a block of the trajectory, |
1433 |
|
|
\texttt{dynamicProps} will calculate all of the time correlation frame |
1434 |
|
|
pairs within the block. After in block correlations are complete, a |
1435 |
|
|
second block of the trajectory is read, and the cross correlations are |
1436 |
|
|
calculated between the two blocks. this second block is then freed and |
1437 |
|
|
then incremented and the process repeated until the end of the |
1438 |
|
|
trajectory. Once the end is reached, the first block is freed then |
1439 |
|
|
incremented, and the again the internal time correlations are |
1440 |
|
|
calculated. The algorithm with the second block is then repeated with |
1441 |
|
|
the new origin block, until all frame pairs have been correlated in |
1442 |
|
|
time. This process is illustrated in |
1443 |
|
|
Fig.~\ref{oopseFig:dynamicPropsMemory}. |
1444 |
mmeineke |
1044 |
|
1445 |
mmeineke |
1054 |
\begin{figure} |
1446 |
|
|
\centering |
1447 |
|
|
\includegraphics[width=\linewidth]{dynamicPropsMem.eps} |
1448 |
|
|
\caption[A representation of the block correlations in \texttt{dynamicProps}]{This diagram illustrates the memory management used by \texttt{dynamicProps}, which follows the scheme: $\sum^{N_{\text{memory blocks}}}_{i=1}[ \operatorname{self}(i) + \sum^{N_{\text{memory blocks}}}_{j>i} \operatorname{cross}(i,j)]$. The shaded region represents the self correlation of the memory block, and the open blocks are read one at a time and the cross correlations between blocks are calculated.} |
1449 |
|
|
\label{oopseFig:dynamicPropsMemory} |
1450 |
|
|
\end{figure} |
1451 |
mmeineke |
1044 |
|
1452 |
mmeineke |
1054 |
\subsection{\label{openSource}Open Source and Distribution License} |
1453 |
mmeineke |
1044 |
|
1454 |
mmeineke |
1054 |
\section{\label{oopseSec:conclusion}Conclusion} |
1455 |
mmeineke |
1044 |
|
1456 |
mmeineke |
1054 |
We have presented the design and implementation of our open source |
1457 |
|
|
simulation package {\sc oopse}. The package offers novel |
1458 |
|
|
capabilities to the field of Molecular Dynamics simulation packages in |
1459 |
|
|
the form of dipolar force fields, and symplectic integration of rigid |
1460 |
|
|
body dynamics. It is capable of scaling across multiple processors |
1461 |
|
|
through the use of MPI. It also implements several integration |
1462 |
|
|
ensembles allowing the end user control over temperature and |
1463 |
|
|
pressure. In addition, it is capable of integrating constrained |
1464 |
|
|
dynamics through both the {\sc rattle} algorithm and the z-constraint |
1465 |
|
|
method. |
1466 |
mmeineke |
1044 |
|
1467 |
mmeineke |
1054 |
These features are all brought together in a single open-source |
1468 |
|
|
development package. This allows researchers to not only benefit from |
1469 |
|
|
{\sc oopse}, but also contribute to {\sc oopse}'s development as |
1470 |
|
|
well.Documentation and source code for {\sc oopse} can be downloaded |
1471 |
|
|
from \texttt{http://www.openscience.org/oopse/}. |
1472 |
mmeineke |
1044 |
|