-
Notifications
You must be signed in to change notification settings - Fork 1
Expand file tree
/
Copy pathopenmdDoc.tex
More file actions
6575 lines (5996 loc) · 300 KB
/
Copy pathopenmdDoc.tex
File metadata and controls
6575 lines (5996 loc) · 300 KB
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
\documentclass[letterpaper]{report}
\usepackage[text={6.5in,9in},centering]{geometry}
\usepackage[T1]{fontenc}
\usepackage{amssymb}
\usepackage{amsmath}
\usepackage{newpxtext,eulerpx}
\usepackage{sourcecodepro}
\usepackage{listings}
\usepackage{graphicx}
\usepackage{setspace}
\usepackage{tabularx}
\usepackage{longtable}
\pagestyle{plain}
\pagenumbering{arabic}
\usepackage{floatrow}
\usepackage{url}
\usepackage[margin=0.5cm,font=small,format=hang]{caption}
\usepackage[version=4]{mhchem}
\usepackage{color}
\newcommand{\angstrom}{\mbox{\normalfont\AA}}
\let\textDJ\DJ
\let\DJ\relax
\DeclareRobustCommand{\DJ}{%
\ifmmode
\mathDJ
\else
\textDJ
\fi
}
\makeatletter
\newcommand{\mathDJ}{\text{\raisebox{0.25ex}{-}\kern-0.4em$\m@th D$}}
\renewcommand{\l@section}{\@dottedtocline{1}{2em}{3em}}
\renewcommand{\l@subsection}{\@dottedtocline{2}{4.0em}{4em}}
\renewcommand{\l@subsubsection}{\@dottedtocline{3}{7.4em}{5em}}
\makeatother
\renewcommand{\baselinestretch}{1.2}
\usepackage[super,comma,sort&compress]{natbib}
\bibpunct{}{}{,}{s}{}{;}
%\bibpunct{[}{]}{,}{n}{}{;}
\DeclareFloatFont{tiny}{\scriptsize}% "scriptsize" is defined by floatrow, "tiny" not
\floatsetup[table]{font=small}
\definecolor{mygreen}{rgb}{0,0.6,0}
\definecolor{mygray}{rgb}{0.5,0.5,0.5}
\definecolor{mymauve}{rgb}{0.58,0,0.82}
%\renewcommand\citemid{\ } % no comma in optional reference note
\lstset{
language=C,
frame=TB,
basicstyle=\footnotesize\ttfamily, %
xleftmargin=0.25in,
xrightmargin=0.25in,
captionpos=b, %
abovecaptionskip=0.5cm,
belowcaptionskip=0.5cm,
escapeinside={~}{~},
identifierstyle=\color{black},
commentstyle=\color{mygreen}\textit, % comment style
keywordstyle=\color{blue}, % keyword style
stringstyle=\color{mymauve} % string literal style
}
\renewcommand{\lstlistingname}{Example}
\lstnewenvironment{code}[1][]%
{\noindent\minipage{\linewidth}\vspace{0.5\baselineskip}
\lstset{
language=C,
basicstyle=\footnotesize\ttfamily,%
captionpos=b,
aboveskip=0.5cm,
belowskip=0.5cm,
abovecaptionskip=0.5cm,
belowcaptionskip=0.5cm,
escapeinside={~}{~},
identifierstyle=\color{black},
commentstyle=\color{mygreen}\textit, % comment style
keywordstyle=\color{blue}, % keyword style
stringstyle=\color{mymauve}, % string literal style
frame=single,#1}}
{\endminipage}
\begin{document}
\newcolumntype{A}{p{1.5in}}
\newcolumntype{B}{p{0.75in}}
\newcolumntype{C}{p{1.5in}}
\newcolumntype{D}{p{2in}}
\newcolumntype{E}{p{0.25in}}
\newcolumntype{F}{p{2in}}
\newcolumntype{G}{p{3.75in}}
\newcolumntype{H}{p{0.75in}}
\newcolumntype{I}{p{5in}}
\newcolumntype{J}{p{1.5in}}
\newcolumntype{K}{p{0.75in}}
\newcolumntype{L}{p{1.5in}}
\newcolumntype{M}{p{2in}}
\title{OpenMD-3.1: Molecular Dynamics in the Open}
\author{C.R. Drisko, Sydney A. Shavalier, Benjamin M. Harless, Veronica Freund, \\
Teng Lin, Charles F. Vardeman II, Christopher J. Fennell, Matthew A. Meineke, \\
Patrick B. Louden, Hemanta Bhattarai, Joseph R. Michalka, Kelsey M. Stocker, \\
James M. Marr, Anderson D.S. Duraes, Suzanne M. Neidhart, Shenyu Kuang, \\
Madan Lamichhane, Xiuquan Sun, Chunlei Li, Kyle Daily, \\
Alexander Mazanek, Yang Zheng, \\ and \\
J. Daniel Gezelter \\ \\
Department of Chemistry and Biochemistry\\
University of Notre Dame\\
Notre Dame, Indiana 46556}
\maketitle
\section*{Preface}
OpenMD is an open source molecular dynamics engine which is
capable of efficiently simulating liquids, proteins, nanoparticles,
interfaces, metals, and other complex systems using atom types with
orientational degrees of freedom (e.g. ``sticky'' atoms, point
multipoles, and coarse-grained assemblies). It is a test-bed for new
molecular simulation methodology, but is also efficient and easy to
use. Liquids, proteins, zeolites, lipids, inorganic nanomaterials,
transition metals (bulk, flat interfaces, and nanoparticles) and a
wide array of other systems have all been simulated using this
code. OpenMD works on parallel computers using the Message
Passing Interface (MPI), and comes with a number of trajectory
analysis and utility programs that are easy to use and modify. An
OpenMD simulation is specified using a very simple meta-data language
that is easy to learn.
\tableofcontents
%\listoffigures
%\listoftables
%\mainmatter
\chapter{\label{sec:intro}Introduction}
There are a number of excellent molecular dynamics packages available
to the chemical physics
community.\cite{Brooks83,MacKerell98,pearlman:1995,Gromacs,Gromacs3,DL_POLY,Tinker,Paradyn,namd,macromodel}
All of these packages are stable, polished programs which solve many
problems of interest. Most are now capable of performing molecular
dynamics simulations on parallel computers. Some have source code
which is freely available to the entire scientific community. Few,
however, are capable of efficiently integrating the equations of
motion for atom types with orientational degrees of freedom
(e.g. point multipoles, and ``sticky'' atoms). And only one of the
programs referenced can handle transition metal force fields like the
Embedded Atom Method (EAM). The direction our research program
has taken us now involves the use of atoms with orientational degrees
of freedom as well as transition metals. Since these simulation
methods may be of some use to other researchers, we have decided to
release our program (and all related source code) to the scientific
community.
This document communicates the algorithmic details of our program,
OpenMD. We have structured this document to first discuss the
underlying concepts in this simulation package (Chapter
\ref{section:IOfiles}). The empirical energy functions implemented
are discussed in Chapter~\ref{chapter:forceFields}.
Section~\ref{section:mechanics} describes the various Molecular Dynamics
algorithms OpenMD implements in the integration of Hamilton's
equations of motion. Program design considerations for parallel
computing are presented in Sec.~\ref{section:parallelization}.
Concluding remarks are presented in Sec.~\ref{section:conclusion}.
\chapter{\label{section:IOfiles}Concepts \& Files}
A simulation in OpenMD is built using a few fundamental
conceptual building blocks, most of which are chemically intuitive.
The basic unit of a simulation is an {\tt atom}. The parameters
describing an {\tt atom} have been generalized to make it as flexible
as possible; this means that in addition to translational degrees of
freedom, {\tt atoms} may also have {\it orientational} degrees of
freedom.
The fundamental (static) properties of {\tt atoms} are defined by the
{\tt forceField} chosen for the simulation. The atomic properties
specified by a {\tt forceField} might include (but are not limited to)
charge, $\sigma$ and $\epsilon$ values for Lennard-Jones interactions,
the strength of the dipole moment ($\mu$), the mass, and the moments
of inertia. Other more complicated properties of atoms might also be
specified by the {\tt forceField}.
{\tt Atoms} can be grouped together in many ways. A {\tt rigidBody}
contains atoms that exert no forces on one another and which move as a
single rigid unit. A {\tt cutoffGroup} may contain atoms which
function together as a (rigid {\it or} non-rigid) unit for potential
energy calculations,
\begin{equation}
V_{ab} = s(r_{ab}) \sum_{i \in a} \sum_{j \in b} V_{ij}(r_{ij})
\end{equation}
Here, $a$ and $b$ are two {\tt cutoffGroups} containing multiple atoms
($a = \left\{i\right\}$ and $b = \left\{j\right\}$). $s(r_{ab})$ is a
generalized switching function which insures that the atoms in the two
{\tt cutoffGroups} are treated identically as the two groups enter or
leave an interaction region.
{\tt Atoms} may also be grouped in more traditional ways into {\tt
bonds}, {\tt bends}, {\tt torsions}, and {\tt inversions}. These
groupings allow the correct choice of interaction parameters for
short-range interactions to be chosen from the definitions in the {\tt
forceField}.
All of these groups of {\tt atoms} are brought together in the {\tt
molecule}, which is the fundamental structure for setting up an
OpenMD simulation. {\tt Molecules} contain lists of {\tt atoms}
followed by listings of the other atomic groupings ({\tt bonds}, {\tt
bends}, {\tt torsions}, {\tt rigidBodies}, and {\tt cutoffGroups})
which relate the atoms to one another. Since a {\tt rigidBody} is a
collection of atoms that are propagated in fixed relationships to one
another, OpenMD uses an internal structure called a {\tt
StuntDouble} to store information about those objects that can change
position {\it independently} during a simulation. That is, an atom
that is part of a rigid body is not itself a StuntDouble. In this
case, the rigid body is the StuntDouble. However, an atom that is
free to move independently {\it is} its own StuntDouble.
Simulations often involve heterogeneous collections of molecules. To
specify a mixture of {\tt molecule} types, OpenMD uses {\tt
components}. Even simulations containing only one type of molecule
must specify a single {\tt component}.
Starting a simulation requires two types of information: {\it
meta-data}, which describes the types of objects present in the
simulation, and {\it configuration} information, which describes the
initial state of these objects. An OpenMD file is a single
combined file format that describes both of these kinds of data. An
OpenMD file contains one {\tt <MetaData>} block and {\it at least
one} {\tt <Snapshot>} block.
The language for the {\tt <MetaData>} block is a C-based syntax that
is parsed at the beginning of the simulation. Configuration
information is specified for all {\tt integrableObjects} in a {\tt
<Snapshot>} block. Both the {\tt <MetaData>} and {\tt <Snapshot>}
formats are described in the following sections.
\begin{code}[caption={[The structure of an OpenMD file]
The basic structure of an OpenMD file contains HTML-like tags to
define simulation meta-data and subsequent instantaneous configuration
information. A well-formed OpenMD file must contain one {\tt <MetaData>}
block and {\it at least one} {\tt <Snapshot>} block. Each
{\tt <Snapshot>} is further divided into {\tt <FrameData>} and
{\tt <StuntDoubles>} sections.},label={sch:mdFormat}]
<OpenMD>
<MetaData>
// see section ~\ref{sec:miscConcepts}~ for details on the formatting
// of information contained inside the <MetaData> tags
</MetaData>
<Snapshot> // An instantaneous configuration
<FrameData>
// FrameData contains information on the time
// stamp, the size of the simulation box, and
// the current state of extended system
// ensemble variables.
</FrameData>
<StuntDoubles>
// StuntDouble information comprises the
// positions, velocities, orientations, and
// angular velocities of anything that is
// capable of independent motion during
// the simulation.
</StuntDoubles>
<SiteData>
// SiteData displays atomic-level details
// of dynamic quantities, e.g. fluctuating charges,
// electric fields, etc., computed during a
// simulation. This block is not always present.
</SiteData>
</Snapshot>
<Snapshot> // Multiple <Snapshot> sections can be
</Snapshot> // present in a well-formed OpenMD file
<Snapshot> // Further information on <Snapshot> blocks
</Snapshot> // can be found in section ~\ref{section:coordFiles}~.
</OpenMD>
\end{code}
\section{OpenMD Files and {\tt <MetaData>} blocks}
OpenMD uses HTML-like delimiters to separate {\tt
<MetaData>} and {\tt <Snapshot>} blocks. A C-based syntax
is used to parse the {\tt <MetaData>} blocks at run time. These
blocks allow the user to completely describe the system they wish to
simulate, as well as tailor OpenMD's behavior during the
simulation. OpenMD files are typically denoted with the
extension {\tt .omd}. An overview of an OpenMD file is shown in
Example~\ref{sch:mdFormat} and example file is shown in
Example~\ref{sch:mdExample}.
\begin{code}[caption={[An example of a complete OpenMD
file] An example showing a complete OpenMD file.},
label={sch:mdExample}]
<OpenMD>
<MetaData>
molecule{
name = "Ar";
atom[0]{
type="Ar";
position( 0.0, 0.0, 0.0 );
}
}
component{
type = "Ar";
nMol = 3;
}
forceField = "LJ";
ensemble = "NVE"; // specify the simulation ensemble
dt = 1.0; // the time step for integration
runTime = 1e3; // the total simulation run time
sampleTime = 100; // trajectory file frequency
statusTime = 50; // statistics file frequency
</MetaData>
<Snapshot>
<FrameData>
Time: 0
Hmat: {{ 28.569, 0, 0 }, { 0, 28.569, 0 }, { 0, 0, 28.569 }}
Thermostat: 0 , 0
Barostat: {{ 0, 0, 0 }, { 0, 0, 0 }, { 0, 0, 0 }}
</FrameData>
<StuntDoubles>
0 pv 17.5 13.3 12.8 1.181e-03 -1.630e-03 -1.369e-03
1 pv -12.8 -14.9 -8.4 -4.440e-04 -2.048e-03 1.130e-03
2 pv -10.0 -15.2 -6.5 2.239e-03 -6.310e-03 1.810e-03
</StuntDoubles>
</Snapshot>
</OpenMD>
\end{code}
In the {\tt <MetaData>} block, it is necessary to provide a
complete description of the molecule before it is actually placed in
the simulation. OpenMD's meta-data syntax allows for the use of
{\it include files} to specify all atoms in a molecular prototype, as
well as any bonds, bends, torsions, or other structural groupings of
atoms. Include files allow the user to describe a molecular prototype
once, then simply include it into each simulation containing that
molecule. Returning to the example in Scheme~\ref{sch:mdExample}, the
include file's contents would be Scheme~\ref{sch:mdIncludeExample},
and the new OpenMD file would become Scheme~\ref{sch:mdExPrime}.
\begin{code}[caption={An example molecule definition in an
include file.},label={sch:mdIncludeExample}]
molecule{
name = "Ar";
atom[0]{
type="Ar";
position( 0.0, 0.0, 0.0 );
}
}
\end{code}
\begin{code}[caption={Revised OpenMD input file
example.},label={sch:mdExPrime}]
<OpenMD>
<MetaData>
#include "argon.inc"
component{
type = "Ar";
nMol = 3;
}
forceField = "LJ";
ensemble = "NVE";
dt = 1.0;
runTime = 1e3;
sampleTime = 100;
statusTime = 50;
</MetaData>
</MetaData>
<Snapshot>
<FrameData>
Time: 0
Hmat: {{ 28.569, 0, 0 }, { 0, 28.569, 0 }, { 0, 0, 28.569 }}
Thermostat: 0 , 0
Barostat: {{ 0, 0, 0 }, { 0, 0, 0 }, { 0, 0, 0 }}
</FrameData>
<StuntDoubles>
0 pv 17.5 13.3 12.8 1.181e-03 -1.630e-03 -1.369e-03
1 pv -12.8 -14.9 -8.4 -4.440e-04 -2.048e-03 1.130e-03
2 pv -10.0 -15.2 -6.5 2.239e-03 -6.310e-03 1.810e-03
</StuntDoubles>
</Snapshot>
</OpenMD>
\end{code}
\section{\label{section:atomsMolecules}Atoms, Molecules, and other
ways of grouping atoms}
As mentioned above, the fundamental unit for an OpenMD
simulation is the {\tt atom}. Atoms can be collected into secondary
structures such as {\tt rigidBodies}, {\tt cutoffGroups}, or {\tt
molecules}. The {\tt molecule} is a way for OpenMD to keep
track of the atoms in a simulation in logical manner. Molecular units
store the identities of all the atoms and rigid bodies associated with
themselves, and they are responsible for the evaluation of their own
internal interactions (\emph{i.e.}~bonds, bends, torsions, and
inversions). Scheme \ref{sch:mdIncludeExample} shows how one creates a
molecule in an included meta-data file. The positions of the atoms
given in the declaration are relative to the origin of the molecule,
and the origin is used when creating a system containing the molecule.
One of the features that sets OpenMD apart from most of the
current molecular simulation packages is the ability to handle rigid
body dynamics. Rigid bodies are non-spherical particles or collections
of particles (e.g. a phenyl ring) that have a constant internal
potential and move collectively.\cite{Goldstein01} They are not
included in most simulation packages because of the algorithmic
complexity involved in propagating orientational degrees of freedom.
Integrators which propagate orientational motion with an acceptable
level of energy conservation for molecular dynamics are relatively new
inventions.
Moving a rigid body involves determination of both the force and
torque applied by the surroundings, which directly affect the
translational and rotational motion in turn. In order to accumulate
the total force on a rigid body, the external forces and torques must
first be calculated for all the internal particles. The total force on
the rigid body is simply the sum of these external forces.
Accumulation of the total torque on the rigid body is more complex
than the force because the torque is applied to the center of mass of
the rigid body. The space-fixed torque on rigid body $i$ is
\begin{equation}
\boldsymbol{\tau}_i=
\sum_{a}\biggl[(\mathbf{r}_{ia}-\mathbf{r}_i)\times \mathbf{f}_{ia}
+ \boldsymbol{\tau}_{ia}\biggr],
\label{eq:torqueAccumulate}
\end{equation}
where $\boldsymbol{\tau}_i$ and $\mathbf{r}_i$ are the torque on and
position of the center of mass respectively, while $\mathbf{f}_{ia}$,
$\mathbf{r}_{ia}$, and $\boldsymbol{\tau}_{ia}$ are the force on,
position of, and torque on the component particles of the rigid body.
The summation of the total torque is done in the body fixed axis of
each rigid body. In order to move between the space fixed and body
fixed coordinate axes, parameters describing the orientation must be
maintained for each rigid body. At a minimum, the rotation matrix
($\mathsf{A}$) can be described by the three Euler angles ($\phi,
\theta,$ and $\psi$), where the elements of $\mathsf{A}$ are composed of
trigonometric operations involving $\phi, \theta,$ and
$\psi$.\cite{Goldstein01} In order to avoid numerical instabilities
inherent in using the Euler angles, the four parameter ``quaternion''
scheme is often used. The elements of $\mathsf{A}$ can be expressed as
arithmetic operations involving the four quaternions ($q_w, q_x, q_y,$
and $q_z$).\cite{Allen87} Use of quaternions also leads to
performance enhancements, particularly for very small
systems.\cite{Evans77}
Rather than use one of the previously stated methods, OpenMD
utilizes a relatively new scheme that propagates the entire nine
parameter rotation matrix. Further discussion on this choice can be
found in Sec.~\ref{section:integrate}. An example definition of a
rigid body can be seen in Scheme
\ref{sch:rigidBody}.
\begin{code}[caption={[Defining rigid bodies]A sample
definition of a molecule containing a rigid body and a cutoff
group},label={sch:rigidBody}]
molecule{
name = "TIP3P";
atom[0]{
type = "O_TIP3P";
position( 0.0, 0.0, -0.06556 );
}
atom[1]{
type = "H_TIP3P";
position( 0.0, 0.75695, 0.52032 );
}
atom[2]{
type = "H_TIP3P";
position( 0.0, -0.75695, 0.52032 );
}
rigidBody[0]{
members(0, 1, 2);
}
cutoffGroup{
members(0, 1, 2);
}
}
\end{code}
\section{\label{sec:miscConcepts}Creating a {\tt <MetaData>} block}
The actual creation of a {\tt <MetaData>} block requires several key
components. The first part of the file needs to be the declaration of
all of the molecule prototypes used in the simulation. This is
typically done through included prototype files. Only the molecules
actually present in the simulation need to be declared; however,
OpenMD allows for the declaration of more molecules than are
needed. This gives the user the ability to build up a library of
commonly used molecules into a single include file.
Once all prototypes are declared, the ordering of the rest of the
block is less stringent. The molecular composition of the simulation
is specified with {\tt component} statements. Each different type of
molecule present in the simulation is considered a separate
component (an example is shown in
Sch.~\ref{sch:mdExPrime}). The component blocks tell OpenMD the
number of molecules that will be in the simulation, and the order in
which the components blocks are declared sets the ordering of the real
atoms in the {\tt <Snapshot>} block as well as in the output files. The
remainder of the script then sets the various simulation parameters
for the system of interest.
The required set of parameters that must be present in all simulations
is given in Table~\ref{table:reqParams}. Since the user can use
OpenMD to perform energy minimizations as well as molecular
dynamics simulations, either a {\tt minimizer} block or the {\tt
ensemble} keyword must be present. The {\tt ensemble} keyword is
responsible for selecting the integration method used for the
calculation of the equations of motion. An in depth discussion of the
various methods available in OpenMD can be found in
Sec.~\ref{section:mechanics}. The {\tt minimizer} block selects which
minimization method to use, and more details on the choices of
minimizer parameters can be found in Sec.~\ref{section:minimizer}. The
{\tt forceField} statement is important for the selection of which
forces will be used in the course of the simulation. OpenMD
supports several force fields, and allows the user to create their own
using a range of pre-defined empirical energy functions. The format
of force field files is outlined
Chapter~\ref{chapter:forceFields}. The force fields are
interchangeable between simulations, with the only requirement being
that all atoms needed by the simulation are defined within the
selected force field.
For molecular dynamics simulations, the time step between force
evaluations is set with the {\tt dt} parameter, and {\tt runTime} will
set the time length of the simulation. Note, that {\tt runTime} is an
absolute time, meaning if the simulation is started at t = 10.0~ns
with a {\tt runTime} of 25.0~ns, the simulation will only run for an
additional 15.0~ns.
For energy minimizations, it is not necessary to specify {\tt dt} or
{\tt runTime}.
To set the initial positions and velocities of all the integrable
objects in the simulation, OpenMD will use the last good {\tt
<Snapshot>} block that was found in the startup file that it was
called with. If the {\tt useInitalTime} flag is set to {\tt true},
the time stamp from this snapshot will also set the initial time stamp
for the simulation. Additional parameters are summarized in
Table~\ref{table:genParams}.
It is important to note the fundamental units in all files which are
read and written by OpenMD. Energies are in $\mbox{kcal
mol}^{-1}$, distances are in $\mbox{\AA}$, times are in $\mbox{fs}$,
translational velocities are in $\mbox{\AA~fs}^{-1}$, and masses are
in $\mbox{amu}$. Orientational degrees of freedom are described using
quaternions (unitless, but $q_w^2 + q_x^2 + q_y^2 + q_z^2 = 1$),
body-fixed angular momenta ($\mbox{amu \AA}^{2} \mbox{radians
fs}^{-1}$), and body-fixed moments of inertia ($\mbox{amu \AA}^{2}$).
\begin{longtable}[c]{ABCD}
\caption{Meta-data Keywords: Required Parameters}
\\
{\bf keyword} & {\bf units} & {\bf use} & {\bf remarks} \\ \hline
\endhead
\hline
\endfoot
{\tt forceField} & string & Sets the base name for the force field file &
OpenMD appends a {\tt .frc} to the end of this to look for a force
field file.\\
{\tt component} & & Defines the molecular components of the system &
Every {\tt <MetaData>} block must have a component statement. \\
{\tt minimizer} & block & Sets parameters for the minimizer & Either {\tt ensemble} or {\tt minimizer} must be specified. \\
{\tt ensemble} & string & Sets the ensemble. & Some possible ensembles are
{\tt NVE, NVT, NPTi, NPAT, NPTf, NPTxyz, LD} and {\tt LangevinHull}. Either {\tt ensemble}
or {\tt minimizer} must be specified. \\
{\tt dt} & fs & Sets the time step. & Selection of {\tt dt} should be
small enough to sample the fastest motion of the simulation. ({\tt
dt} is required for molecular dynamics simulations)\\
{\tt runTime} & fs & Sets the time at which the simulation should
end. & This is an absolute time, and will end the simulation when the
current time meets or exceeds the {\tt runTime}. ({\tt runTime} is
required for molecular dynamics simulations)
\label{table:reqParams}
\end{longtable}
\begin{longtable}[c]{ABCD}
\caption{Meta-data Keywords: Optional Parameters}
\\
{\bf keyword} & {\bf units} & {\bf use} & {\bf remarks} \\ \hline
\endhead
\hline
\endfoot
{\tt forceFieldVariant} & string & Sets the name of the variant of the
force field. & EAM has three variants: {\tt u3}, {\tt u6}, and
{\tt VC}. \\
{\tt forceFieldFileName} & string & Overrides the default force field
file name & Each force field has a default file name, and this
parameter can override the default file name for the chosen force
field. \\
{\tt usePeriodicBoundaryConditions} & & & \\
& logical & Turns periodic boundary conditions on/off. & Default is true. \\
{\tt orthoBoxTolerance} & double & & decides how orthogonal the periodic
box must be before we can use cheaper box calculations \\
{\tt cutoffRadius} & $\mbox{\AA}$ & Manually sets the cutoff radius &
default is set from {\tt cutoffPolicy} \\
{\tt skinThickness} & \AA & thickness of the skin for the Verlet
neighbor lists & defaults to 1 \AA \\
{\tt switchingRadius} & $\mbox{\AA}$ & Manually sets the inner radius
for the switching function. & Defaults to 85~\% of the {\tt
cutoffRadius}. \\
{\tt switchingFunctionType} & & & \\
& string & {\tt cubic} or
{\tt fifth\_order\_polynomial} & Default is {\tt cubic} \\
{\tt useInitialTime} & logical & Sets whether the initial time is
taken from the last {\tt <Snapshot>} in the startup file. & Useful when recovering a simulation from a crashed processor. Default is {\tt false}. \\
{\tt useInitialExtendedSystemState} & & & \\
& logical & keep the extended
system variables? & Should the extended
variables (the thermostat and barostat) be kept from the {\tt <Snapshot>} block? \\
{\tt sampleTime} & fs & Sets the frequency at which the {\tt .dump} file is written. & Default is set to {\tt runTime}. \\
{\tt resetTime} & fs & Sets the frequency at which the extended system
variables are reset to zero & The default is to never reset these
variables. \\
{\tt statusTime} & fs & Sets the frequency at which the {\tt .stat} file is written. & Default is set to {\tt sampleTime}. \\
{\tt finalConfig} & string & Sets the name of the final output file. & Useful when stringing simulations together. Defaults to the root name of the initial meta-data file but with an {\tt .eor} extension. \\
{\tt compressDumpFile} & logical & & should the {\tt .dump} file be
compressed on the fly? \\
{\tt statFileFormat} & string & columns to print in the {\tt .stat}
file where each column is separated by a pipe ($\mid$) symbol. & (The
default is the first eight of these columns in order.) \\
& \multicolumn{3}{p{4.25in}}{Allowed
column names are: {\tt TIME, TOTAL\_ENERGY, POTENTIAL\_ENERGY, KINETIC\_ENERGY,
TEMPERATURE, PRESSURE, VOLUME, CONSERVED\_QUANTITY, HULLVOLUME, GYRVOLUME,
TRANSLATIONAL\_KINETIC, ROTATIONAL\_KINETIC, LONG\_RANGE\_POTENTIAL,
SHORT\_RANGE\_POTENTIAL, VANDERWAALS\_POTENTIAL,
ELECTROSTATIC\_POTENTIAL, METALLIC\_POTENTIAL,
HYDROGEN\_BONDING\_POTENTIAL, RECIPROCAL\_POTENTIAL,
SURFACE\_POTENTIAL, BOND\_POTENTIAL, BEND\_POTENTIAL,
DIHEDRAL\_POTENTIAL, INVERSION\_POTENTIAL, RAW\_POTENTIAL, RESTRAINT\_POTENTIAL,
PRESSURE\_TENSOR, SYSTEM\_DIPOLE, SYSTEM\_QUADRUPOLE, HEATFLUX,
ELECTRONIC\_TEMPERATURE, COM, COM\_VELOCITY, ANGULAR\_MOMENTUM, POTENTIAL\_SELECTION }} \\
{\tt printPressureTensor} & logical & sets whether OpenMD will print
out the pressure tensor & can be useful for calculations of the bulk
modulus \\
{\tt electrostaticSummationMethod} & & & \\
& string & {\tt shifted\_force,
shifted\_potential, hard, switched, taylor\_shifted,} or {\tt reaction\_field} &
default is {\tt shifted\_force} \\
{\tt electrostaticScreeningMethod} & & & \\
& string & {\tt undamped} or {\tt damped} & default is {\tt damped} \\
{\tt dielectric} & unitless & Sets the dielectric constant for
reaction field. & If {\tt electrostaticSummationMethod} is set to {\tt
reaction\_field}, then {\tt dielectric} must be set. \\
{\tt dampingAlpha} & $\mbox{\AA}^{-1}$ & governs strength of
electrostatic damping & default = max(0, 0.425 - {\tt cutoffRadius} * 0.02) \\
{\tt tempSet} & logical & resample velocities from a Maxwell-Boltzmann
distribution set to {\tt targetTemp} & default is {\tt false} \\
{\tt thermalTime} & fs & how often to perform a {\tt tempSet} &
default is never \\
{\tt targetTemp} & K & sets the target temperature & no default value \\
{\tt tauThermostat} & fs & time constant for Nos\'{e}-Hoover
thermostat & times from 100-10,000 fs are reasonable \\
{\tt targetPressure} & atm & sets the target pressure & no default value\\
{\tt surfaceTension} & & sets the target surface tension in the x-y
plane & no default value \\
{\tt tauBarostat} & fs & time constant for the
Nos\'{e}-Hoover-Andersen barostat & times from 1000 to 100,000 fs
are reasonable \\
{\tt seed } & integer & Sets the seed value for the random number generator. & The seed needs to be at least 9 digits long. The default is to take the seed from the CPU clock. \\
\label{table:genParams}
\end{longtable}
\section{\label{section:coordFiles}{\tt <Snapshot>} Blocks}
The standard format for storage of a system's coordinates is the {\tt
<Snapshot>} block , the exact details of which can be seen in
Scheme~\ref{sch:dumpFormat}. As all bonding and molecular information
is stored in the {\tt <MetaData>} blocks, the {\tt <Snapshot>} blocks
contain only the coordinates of the objects which move independently
during the simulation. It is important to note that {\it not all
atoms} are capable of independent motion. Atoms which are part of
rigid bodies are not ``integrable objects'' in the equations of
motion; the rigid bodies themselves are the integrable objects.
Therefore, the coordinate file contains coordinates of all the {\tt
integrableObjects} in the system. For systems without rigid bodies,
this is simply the coordinates of all the atoms.
It is important to note that although the simulation propagates the
complete rotation matrix, directional entities are written out using
quaternions to save space in the output files.
{\tt <Snapshot>} blocks contain an initial sub-block called {\tt
<FrameData>} which contains the time stamp, the three column vectors
of $\mathsf{H}$, and optional extra information for the extended sytem
ensembles. Following this, the lines in the {\tt <StuntDoubles>} sub-block provide
information about the instantaneous configuration of each integrable
object. For each integrable object, the global index is followed by a
short string describing what additional information is present on the
line. Atoms with only position and velocity information use the {\tt
pv} string which must then be followed by the position and velocity
vectors for that atom. Directional atoms and Rigid Bodies typically
use the {\tt pvqj} string which is followed by position, velocity,
quaternions, and lastly, body fixed angular momentum for that
integrable object.
Some {\tt <Snapshot>} blocks may also contain a {\tt <SiteData>} block
that can provide information about fluctuating charge ({\tt c}),
charge velocity ({\tt w}), electric field ({\tt e}) or other data at
atomic sites (with two indices) or at the centers of mass of rigid
bodies (with one index). The {\tt <SiteData>} block is usually not
needed to initiate a simulation, but are often necessary for obtaining
detailed atomic-level information following a simulation.
\begin{code}[caption={[The format of the {\tt <Snapshot>} block]
An example of the format of the {\tt <Snapshot>} block. There is an
initial sub-block called {\tt <FrameData>} which contains the time
stamp, the three column vectors of $\mathsf{H}$, and optional extra
information for the extended sytem ensembles. The lines in the {\tt
<StuntDoubles>} sub-block provide information about the instantaneous
configuration of each integrable object. For each integrable object,
the global index is followed by a short string describing what
additional information is present on the line. Atoms with only
position and velocity information use the {\tt pv} string which must
then be followed by the position and velocity vectors for that atom.
Directional atoms and Rigid Bodies typically use the {\tt pvqj} string
which is followed by position, velocity, quaternions, and
lastly, body fixed angular momentum for that integrable object.
Some files will also contain a {\tt <SiteData>} block that can provide
additional atomic-level information.},label={sch:dumpFormat}]
<Snapshot>
<FrameData>
Time: 0
Hmat: {{ Hxx, Hyx, Hzx }, { Hxy, Hyy, Hzy }, { Hxz, Hyz, Hzz }}
Thermostat: 0 , 0
Barostat: {{ 0, 0, 0 }, { 0, 0, 0 }, { 0, 0, 0 }}
</FrameData>
<StuntDoubles>
0 pv x y z vx vy vz
1 pv x y z vx vy vz
2 pvqj x y z vx vy vz qw qx qy qz jx jy jz
3 pvqj x y z vx vy vz qw qx qy qz jx jy jz
</StuntDoubles>
<SiteData>
0 0 cwe c w ex ey ez
1 0 cwe c w ex ey ez
2 cwe c w ex ey ez
2 0 cwe c w ex ey ez
2 1 cwe c w ex ey ez
2 2 cwe c w ex ey ez
3 cwe c w ex ey ez
3 0 cwe c w ex ey ez
3 1 cwe c w ex ey ez
3 2 cwe c w ex ey ez
</SiteData>
</Snapshot>
\end{code}
There are three OpenMD files that are written using the combined
format. They are: the initial startup file (\texttt{.omd}), the
simulation trajectory file (\texttt{.dump}), and the final coordinates
or ``end-of-run'' for the simulation (\texttt{.eor}). The initial
startup file is necessary for OpenMD to start the simulation with
the proper coordinates, and this file must be generated by the user
before the simulation run. The trajectory (or ``dump'') file is
updated during simulation and is used to store snapshots of the
coordinates at regular intervals. The first frame is a duplication of
the initial configuration (the last good {\tt <Snapshot>} in the
startup file), and each subsequent frame is appended to the dump file
at an interval specified in the meta-data file with the
\texttt{sampleTime} flag. The final coordinate file is the
``end-of-run'' file. The \texttt{.eor} file stores the final
configuration of the system for a given simulation. The file is
updated at the same time as the \texttt{.dump} file, but it only
contains the most recent frame. In this way, an \texttt{.eor} file may
be used to initialize a second simulation should it be necessary to
recover from a crash or power outage. The coordinate files generated
by OpenMD (both \texttt{.dump} and \texttt{.eor}) all contain the
same {\tt <MetaData>} block as the startup file, so they may be
used to start up a new simulation if desired.
\section{\label{section:initCoords}Generation of Initial Coordinates}
As was stated in Sec.~\ref{section:coordFiles}, a meaningful {\tt
<Snapshot>} block is necessary for specifying for the starting
coordinates for a simulation. Since each simulation is different,
system creation is left to the end user; however, we have included a
few sample programs which make some specialized structures. The {\tt
<Snapshot>} block must index the integrable objects in the correct
order. The ordering of the integrable objects relies on the ordering
of molecules within the {\tt <MetaData>} block. OpenMD
expects the order to comply with the following guidelines:
\begin{enumerate}
\item All of the molecules of the first declared component are given
before proceeding to the molecules of the second component, and so on
for all subsequently declared components.
\item The ordering of the atoms for each molecule follows the order
declared in the molecule's declaration within the model file.
\item Only atoms which are not members of a {\tt rigidBody} are
included.
\item Rigid Body coordinates for a molecule are listed immediately
after the the other atoms in a molecule. Some molecules may be
entirely rigid, in which case, only the rigid body coordinates are
given.
\end{enumerate}
An example is given in the OpenMD file in Scheme~\ref{sch:initEx1}.
\begin{code}[caption={Example declaration of the
$\text{I}_2$ molecule and the HCl molecule in {\tt <MetaData>} and
{\tt <Snapshot>} blocks. Note that even though $\text{I}_2$ is
declared before HCl, the {\tt <Snapshot>} block follows the order {\it in
which the components were included}.}, label=sch:initEx1]
<OpenMD>
<MetaData>
molecule{
name = "I2";
atom[0]{ type = "I"; }
atom[1]{ type = "I"; }
bond{ members( 0, 1); }
}
molecule{
name = "HCl"
atom[0]{ type = "H";}
atom[1]{ type = "Cl";}
bond{ members( 0, 1); }
}
component{
type = "HCl";
nMol = 4;
}
component{
type = "I2";
nMol = 1;
}
</MetaData>
<Snapshot>
<FrameData>
Time: 0
Hmat: {{ 10.0, 0.0, 0.0 }, { 0.0, 10.0, 0.0 }, { 0.0, 0.0, 10.0 }}
</FrameData>
<StuntDoubles>
0 pv x y z vx vy vz // H from first HCl molecule
1 pv x y z vx vy vz // Cl from first HCl molecule
2 pv x y z vx vy vz // H from second HCl molecule
3 pv x y z vx vy vz // Cl from second HCl molecule
4 pv x y z vx vy vz // H from third HCl molecule
5 pv x y z vx vy vz // Cl from third HCl molecule
6 pv x y z vx vy vz // H from fourth HCl molecule
7 pv x y z vx vy vz // Cl from fourth HCl molecule
8 pv x y z vx vy vz // First I from I2 molecule
9 pv x y z vx vy vz // Second I from I2 molecule
</StuntDoubles>
</Snapshot>
</OpenMD>
\end{code}
\section{The Statistics ({\tt .stat}) and Report ({\tt .report}) Files}
An important output file generated by OpenMD is the statistics
file. This file records such statistical quantities as the
instantaneous temperature (in $K$), volume (in $\mbox{\AA}^{3}$),
pressure (in $\mbox{atm}$), etc. It is written out with the frequency
specified in the meta-data file with the
\texttt{statusTime} keyword. The file allows the user to observe the
system variables as a function of simulation time while the simulation
is in progress. One useful function the statistics file serves is to
monitor the conserved quantity of a given simulation ensemble,
allowing the user to gauge the stability of the integrator. The
statistics file is denoted with the \texttt{.stat} file extension.
The fields printed in the {\tt .stat} file are determined by the {\tt
statFileFormat} keyword that the user sets in the {\tt <MetaData>} block.
At the end of a simulation, OpenMD will also produce a {\tt
.report} file. This file contain statistical means and 95\%
confidence intervals for the fields printed in the {\tt .stat} file.
Report files can be very useful in determining how to alter the
simulation box to set an average volume after a constant pressure
simulation, and for verifying that the energy or conserved quantity
has been conserved reasonably well during a simulation.
\chapter{\label{chapter:forceFields}Force Fields}
Like many molecular simulation packages, OpenMD splits the
potential energy into the short-ranged (bonded) portion and a
long-range (non-bonded) potential,
\begin{equation}
V = V_{\mathrm{short-range}} + V_{\mathrm{long-range}}.
\end{equation}
The short-ranged portion includes the bonds, bends, torsions, and
inversions which have been defined in the meta-data file for the
molecules. The functional forms and parameters for these interactions
are defined by the force field which is selected in the MetaData
section.
\section{\label{section:divisionOfLabor}Separation into Internal and
Cross interactions}
The classical potential energy function for a system composed of $N$
molecules is traditionally written
\begin{equation}
V = \sum^{N}_{I=1} V^{I}_{\text{Internal}}
+ \sum^{N-1}_{I=1} \sum_{J>I} V^{IJ}_{\text{Cross}},
\label{eq:totalPotential}
\end{equation}
where $V^{I}_{\text{Internal}}$ contains all of the terms internal to
molecule $I$ (e.g. bonding, bending, torsional, and inversion terms)
and $V^{IJ}_{\text{Cross}}$ contains all intermolecular interactions
between molecules $I$ and $J$. For large molecules, the internal
potential may also include some non-bonded terms like electrostatic or
van der Waals interactions.
The types of atoms being simulated, as well as the specific functional
forms and parameters of the intra-molecular functions and the
long-range potentials are defined by the force field. In the following
sections we discuss the stucture of an OpenMD force field file and the
specification of blocks that may be present within these files.
\section{\label{section:frcFile}Force Field Files}
Force field files have a number of ``Blocks'' to delineate different
types of information. The blocks contain AtomType data, which provide
properties belonging to a single AtomType, as well as interaction
information which provides information about bonded or non-bonded
interactions that cannot be deduced from AtomType information alone.
A simple example of a forceField file is shown in scheme
\ref{sch:frcExample}.
\begin{code}[caption={[An example of a complete OpenMD
force field file for straight-chain united-atom alkanes.] An example
showing a complete OpenMD force field for straight-chain united-atom
alkanes.}, label={sch:frcExample}]
begin Options
Name = "alkane"
end Options
begin BaseAtomTypes
//name mass
C 12.0107
end BaseAtomTypes
begin AtomTypes
//name base mass
CH4 C 16.05
CH3 C 15.04
CH2 C 14.03
end AtomTypes
begin LennardJonesAtomTypes
//name epsilon sigma
CH4 0.2941 3.73
CH3 0.1947 3.75
CH2 0.09140 3.95
end LennardJonesAtomTypes
begin BondTypes
//AT1 AT2 Type r0 k
CH3 CH3 Harmonic 1.526 260
CH3 CH2 Harmonic 1.526 260
CH2 CH2 Harmonic 1.526 260
end BondTypes
begin BendTypes
//AT1 AT2 AT3 Type theta0 k
CH3 CH2 CH3 Harmonic 114.0 124.19
CH3 CH2 CH2 Harmonic 114.0 124.19
CH2 CH2 CH2 Harmonic 114.0 124.19
end BendTypes
begin TorsionTypes
//AT1 AT2 AT3 AT4 Type
CH3 CH2 CH2 CH3 Trappe 0.0 0.70544 -0.13549 1.5723
CH3 CH2 CH2 CH2 Trappe 0.0 0.70544 -0.13549 1.5723
CH2 CH2 CH2 CH2 Trappe 0.0 0.70544 -0.13549 1.5723
end TorsionTypes
\end{code}
\section{\label{section:ffOptions}The {\tt Options} block}
The Options block defines properties governing how the force field
interactions are carried out. This block is delineated with the text
tags {\tt begin Options} and {\tt end Options}. Most options don't
need to be set as they come with fairly sensible default values, but
the various keywords and their possible values are given in Scheme
\ref{sch:optionsBlock}.
\begin{code}[caption={[A force field Options block showing default values
for many force field options.] A force field Options block showing default values
for many force field options. Most of these options do not need to be
specified if the default values are working.},
label={sch:optionsBlock}]
begin Options
Name "alkane" // any string