17 |
|
\textwidth 6.5in |
18 |
|
\brokenpenalty=10000 |
19 |
|
\renewcommand{\baselinestretch}{1.2} |
20 |
+ |
\usepackage[square, comma, sort&compress]{natbib} |
21 |
+ |
\bibpunct{[}{]}{,}{n}{}{;} |
22 |
|
|
23 |
+ |
|
24 |
|
%\renewcommand\citemid{\ } % no comma in optional reference note |
25 |
< |
\lstset{language=C,frame=TB,basicstyle=\tiny,basicstyle=\ttfamily, % |
25 |
> |
\lstset{language=C,frame=TB,basicstyle=\footnotesize,basicstyle=\ttfamily, % |
26 |
|
xleftmargin=0.25in, xrightmargin=0.25in,captionpos=b, % |
27 |
|
abovecaptionskip=0.5cm, belowcaptionskip=0.5cm, escapeinside={~}{~}} |
28 |
|
\renewcommand{\lstlistingname}{Scheme} |
41 |
|
\newcolumntype{H}{p{0.75in}} |
42 |
|
\newcolumntype{I}{p{5in}} |
43 |
|
|
44 |
+ |
\newcolumntype{J}{p{1.5in}} |
45 |
+ |
\newcolumntype{K}{p{1.2in}} |
46 |
+ |
\newcolumntype{L}{p{1.5in}} |
47 |
+ |
\newcolumntype{M}{p{1.55in}} |
48 |
|
|
42 |
– |
\title{{\sc OpenMD}: Molecular Dynamics in the Open} |
49 |
|
|
50 |
< |
\author{Kelsey M. Stocker, Shenyu Kuang, Charles F. Vardeman II, \\ |
51 |
< |
Teng Lin, Christopher J. Fennell, Xiuquan Sun, \\ |
50 |
> |
\title{{\sc OpenMD-2.1}: Molecular Dynamics in the Open} |
51 |
> |
|
52 |
> |
\author{Joseph Michalka, James Marr, Kelsey Stocker, Madan Lamichhane, |
53 |
> |
Patrick Louden, \\ |
54 |
> |
Teng Lin, Charles F. Vardeman II, Christopher J. Fennell, Shenyu |
55 |
> |
Kuang, Xiuquan Sun, \\ |
56 |
|
Chunlei Li, Kyle Daily, Yang Zheng, Matthew A. Meineke, and \\ |
57 |
|
J. Daniel Gezelter \\ |
58 |
|
Department of Chemistry and Biochemistry\\ |
143 |
|
leave an interaction region. |
144 |
|
|
145 |
|
{\tt Atoms} may also be grouped in more traditional ways into {\tt |
146 |
< |
bonds}, {\tt bends}, and {\tt torsions}. These groupings allow the |
147 |
< |
correct choice of interaction parameters for short-range interactions |
148 |
< |
to be chosen from the definitions in the {\tt forceField}. |
146 |
> |
bonds}, {\tt bends}, {\tt torsions}, and {\tt inversions}. These |
147 |
> |
groupings allow the correct choice of interaction parameters for |
148 |
> |
short-range interactions to be chosen from the definitions in the {\tt |
149 |
> |
forceField}. |
150 |
|
|
151 |
|
All of these groups of {\tt atoms} are brought together in the {\tt |
152 |
|
molecule}, which is the fundamental structure for setting up and {\sc |
502 |
|
\endhead |
503 |
|
\hline |
504 |
|
\endfoot |
505 |
< |
{\tt forceField} & string & Sets the force field. & Possible force |
506 |
< |
fields are DUFF, WATER, LJ, EAM, SC, and CLAY. \\ |
505 |
> |
{\tt forceField} & string & Sets the base name for the force field file & |
506 |
> |
OpenMD appends a {\tt .frc} to the end of this to look for a force |
507 |
> |
field file.\\ |
508 |
|
{\tt component} & & Defines the molecular components of the system & |
509 |
|
Every {\tt $<$MetaData$>$} block must have a component statement. \\ |
510 |
|
{\tt minimizer} & string & Chooses a minimizer & Possible minimizers |
574 |
|
default is the first eight of these columns in order.) \\ |
575 |
|
& & \multicolumn{2}{p{3.5in}}{Allowed |
576 |
|
column names are: {\sc time, total\_energy, potential\_energy, kinetic\_energy, |
577 |
< |
temperature, pressure, volume, conserved\_quantity, |
577 |
> |
temperature, pressure, volume, conserved\_quantity, hullvolume, gyrvolume, |
578 |
|
translational\_kinetic, rotational\_kinetic, long\_range\_potential, |
579 |
|
short\_range\_potential, vanderwaals\_potential, |
580 |
< |
electrostatic\_potential, bond\_potential, bend\_potential, |
581 |
< |
dihedral\_potential, improper\_potential, vraw, vharm, |
582 |
< |
pressure\_tensor\_x, pressure\_tensor\_y, pressure\_tensor\_z}} \\ |
580 |
> |
electrostatic\_potential, metallic\_potential, |
581 |
> |
hydrogen\_bonding\_potential, bond\_potential, bend\_potential, |
582 |
> |
dihedral\_potential, inversion\_potential, raw\_potential, restraint\_potential, |
583 |
> |
pressure\_tensor, system\_dipole, heatflux, electronic\_temperature}} \\ |
584 |
|
{\tt printPressureTensor} & logical & sets whether {\sc OpenMD} will print |
585 |
|
out the pressure tensor & can be useful for calculations of the bulk |
586 |
|
modulus \\ |
775 |
|
allowing the user to gauge the stability of the integrator. The |
776 |
|
statistics file is denoted with the \texttt{.stat} file extension. |
777 |
|
|
778 |
< |
\chapter{\label{section:empiricalEnergy}The Empirical Energy |
766 |
< |
Functions} |
778 |
> |
\chapter{\label{section:forceFields}Force Fields} |
779 |
|
|
780 |
< |
Like many simulation packages, {\sc OpenMD} splits the potential energy |
781 |
< |
into the short-ranged (bonded) portion and a long-range (non-bonded) |
782 |
< |
potential, |
780 |
> |
Like many molecular simulation packages, {\sc OpenMD} splits the |
781 |
> |
potential energy into the short-ranged (bonded) portion and a |
782 |
> |
long-range (non-bonded) potential, |
783 |
|
\begin{equation} |
784 |
|
V = V_{\mathrm{short-range}} + V_{\mathrm{long-range}}. |
785 |
|
\end{equation} |
786 |
< |
The short-ranged portion includes the explicit bonds, bends, and |
787 |
< |
torsions which have been defined in the meta-data file for the |
788 |
< |
molecules which are present in the simulation. The functional forms and |
789 |
< |
parameters for these interactions are defined by the force field which |
790 |
< |
is chosen. |
786 |
> |
The short-ranged portion includes the bonds, bends, torsions, and |
787 |
> |
inversions which have been defined in the meta-data file for the |
788 |
> |
molecules. The functional forms and parameters for these interactions |
789 |
> |
are defined by the force field which is selected in the MetaData |
790 |
> |
section. |
791 |
|
|
792 |
< |
Calculating the long-range (non-bonded) potential involves a sum over |
793 |
< |
all pairs of atoms (except for those atoms which are involved in a |
794 |
< |
bond, bend, or torsion with each other). If done poorly, calculating |
795 |
< |
the the long-range interactions for $N$ atoms would involve $N(N-1)/2$ |
784 |
< |
evaluations of atomic distances. To reduce the number of distance |
785 |
< |
evaluations between pairs of atoms, {\sc OpenMD} uses a switched cutoff |
786 |
< |
with Verlet neighbor lists.\cite{Allen87} It is well known that |
787 |
< |
neutral groups which contain charges will exhibit pathological forces |
788 |
< |
unless the cutoff is applied to the neutral groups evenly instead of |
789 |
< |
to the individual atoms.\cite{leach01:mm} {\sc OpenMD} allows users to |
790 |
< |
specify cutoff groups which may contain an arbitrary number of atoms |
791 |
< |
in the molecule. Atoms in a cutoff group are treated as a single unit |
792 |
< |
for the evaluation of the switching function: |
792 |
> |
\section{\label{section:shortRange}The basic interactions} |
793 |
> |
|
794 |
> |
The energy function for a system composed of $N$ molecules is |
795 |
> |
traditionally written |
796 |
|
\begin{equation} |
797 |
< |
V_{\mathrm{long-range}} = \sum_{a} \sum_{b>a} s(r_{ab}) \sum_{i \in a} \sum_{j \in b} V_{ij}(r_{ij}), |
798 |
< |
\end{equation} |
799 |
< |
where $r_{ab}$ is the distance between the centers of mass of the two |
797 |
< |
cutoff groups ($a$ and $b$). |
798 |
< |
|
799 |
< |
The sums over $a$ and $b$ are over the cutoff groups that are present |
800 |
< |
in the simulation. Atoms which are not explicitly defined as members |
801 |
< |
of a {\tt cutoffGroup} are treated as a group consisting of only one |
802 |
< |
atom. The switching function, $s(r)$ is the standard cubic switching |
803 |
< |
function, |
804 |
< |
\begin{equation} |
805 |
< |
S(r) = |
806 |
< |
\begin{cases} |
807 |
< |
1 & \text{if $r \le r_{\text{sw}}$},\\ |
808 |
< |
\frac{(r_{\text{cut}} + 2r - 3r_{\text{sw}})(r_{\text{cut}} - r)^2} |
809 |
< |
{(r_{\text{cut}} - r_{\text{sw}})^3} |
810 |
< |
& \text{if $r_{\text{sw}} < r \le r_{\text{cut}}$}, \\ |
811 |
< |
0 & \text{if $r > r_{\text{cut}}$.} |
812 |
< |
\end{cases} |
813 |
< |
\label{eq:dipoleSwitching} |
797 |
> |
V = \sum^{N}_{I=1} V^{I}_{\text{Internal}} |
798 |
> |
+ \sum^{N-1}_{I=1} \sum_{J>I} V^{IJ}_{\text{Cross}}, |
799 |
> |
\label{eq:totalPotential} |
800 |
|
\end{equation} |
801 |
< |
Here, $r_{\text{sw}}$ is the {\tt switchingRadius}, or the distance |
802 |
< |
beyond which interactions are reduced, and $r_{\text{cut}}$ is the |
803 |
< |
{\tt cutoffRadius}, or the distance at which interactions are |
804 |
< |
truncated. |
801 |
> |
where $V^{IJ}_{\text{Cross}}$ contains all intermolecular interactions |
802 |
> |
between molecules $I$ and $J$, and $V^{I}_{\text{Internal}}$ is the |
803 |
> |
internal potential of molecule $I$: |
804 |
> |
\begin{align*} |
805 |
> |
V^{I}_{\text{Internal}} = & |
806 |
> |
\sum_{r_{ij} \in I} V_{\text{bond}}(r_{ij}) |
807 |
> |
+ \sum_{\theta_{ijk} \in I} V_{\text{bend}}(\theta_{ijk}) |
808 |
> |
+ \sum_{\phi_{ijkl} \in I} V_{\text{torsion}}(\phi_{ijkl}) |
809 |
> |
+ \sum_{\omega_{ijkl} \in I} V_{\text{inversion}}(\omega_{ijkl}) \\ |
810 |
> |
& + \sum_{i \in I} \sum_{(j>i+4) \in I} |
811 |
> |
\biggl[ V_{\text{dispersion}}(r_{ij}) + V_{\text{electrostatic}} |
812 |
> |
(\mathbf{r}_{ij},\boldsymbol{\Omega}_{i},\boldsymbol{\Omega}_{j}) |
813 |
> |
\biggr]. |
814 |
> |
\label{eq:internalPotential} |
815 |
> |
\end{align*} |
816 |
> |
Here $V_{\text{bond}}, V_{\text{bend}}, |
817 |
> |
V_{\text{torsion}},\mathrm{~and~} V_{\text{inversion}}$ represent the |
818 |
> |
bond, bend, torsion, and inversion potentials for all |
819 |
> |
topologically-connected sets of atoms within the molecule. Bonds are |
820 |
> |
the primary way of specifying how the atoms are connected together to |
821 |
> |
form the molecule (i.e. they define the molecular topology). The |
822 |
> |
other short-range interactions may be specified explicitly in the |
823 |
> |
molecule definition, or they may be deduced from bonding information. |
824 |
> |
For example, bends can be implicitly deduced from two bonds which |
825 |
> |
share an atom. Torsions can be deduced from two bends that share a |
826 |
> |
bond. Inversion potentials are utilized primarily to enforce |
827 |
> |
planarity around $sp^2$-hybridized sites, and these are specified with |
828 |
> |
central atoms and satellites (or an atom with bonds to exactly three |
829 |
> |
satellites). The pairwise portions of the non-bonded interactions are |
830 |
> |
usually excluded for atom pairs that are involved in the same bond, |
831 |
> |
bend, or torsion. All other atom pairs within a molecule are subject |
832 |
> |
to non-bonded pair potentials. |
833 |
|
|
834 |
< |
Users of {\sc OpenMD} do not need to specify the {\tt cutoffRadius} or |
835 |
< |
{\tt switchingRadius}. In simulations containing only Lennard-Jones |
836 |
< |
atoms, the cutoff radius has a default value of $2.5\sigma_{ii}$, |
837 |
< |
where $\sigma_{ii}$ is the largest Lennard-Jones length parameter |
838 |
< |
present in the simulation. In simulations containing charged or |
825 |
< |
dipolar atoms, the default cutoff radius is $15 \mbox{\AA}$. |
834 |
> |
The types of atoms being simulated, as well as the specific functional |
835 |
> |
forms and parameters of the intra-molecular functions and the |
836 |
> |
long-range potentials are defined by the force field. In the following |
837 |
> |
sections we discuss the stucture of an OpenMD force field file and the |
838 |
> |
specification of blocks that may be present within these files. |
839 |
|
|
840 |
< |
The {\tt switchingRadius} is set to a default value of 95\% of the |
828 |
< |
{\tt cutoffRadius}. In the special case of a simulation containing |
829 |
< |
{\it only} Lennard-Jones atoms, the default switching radius takes the |
830 |
< |
same value as the cutoff radius, and {\sc OpenMD} will use a shifted |
831 |
< |
potential to remove discontinuities in the potential at the cutoff. |
832 |
< |
Both radii may be specified in the meta-data file. |
840 |
> |
\section{\label{section:frcFile}Force Field Files} |
841 |
|
|
842 |
< |
Force fields can be added to {\sc OpenMD}, although it comes with a few |
843 |
< |
simple examples (Lennard-Jones, {\sc duff}, {\sc water}, and {\sc |
844 |
< |
eam}) which are explained in the following sections. |
842 |
> |
Force field files have a number of ``Blocks'' to demarkate different |
843 |
> |
types of information. The blocks contain AtomType data, which provide |
844 |
> |
properties belonging to a single AtomType, as well as interaction |
845 |
> |
information which provides information about bonded or non-bonded |
846 |
> |
interactions that cannot be deduced from AtomType information alone. |
847 |
> |
A simple example of a forceField file is shown in scheme |
848 |
> |
\ref{sch:frcExample}. |
849 |
|
|
850 |
< |
\section{\label{sec:LJPot}The Lennard Jones Force Field} |
850 |
> |
\begin{lstlisting}[float,caption={[An example of a complete OpenMD |
851 |
> |
force field file for straight-chain united-atom alkanes.] An example |
852 |
> |
showing a complete OpenMD force field for straight-chain united-atom |
853 |
> |
alkanes.}, label={sch:frcExample}] |
854 |
> |
begin Options |
855 |
> |
Name = "alkane" end |
856 |
> |
Options |
857 |
|
|
858 |
< |
The most basic force field implemented in {\sc OpenMD} is the |
859 |
< |
Lennard-Jones force field, which mimics the van der Waals interaction |
860 |
< |
at long distances and uses an empirical repulsion at short |
858 |
> |
begin BaseAtomTypes |
859 |
> |
//name mass |
860 |
> |
C 12.0107 |
861 |
> |
end BaseAtomTypes |
862 |
> |
|
863 |
> |
begin AtomTypes |
864 |
> |
//name base mass |
865 |
> |
CH4 C 16.05 |
866 |
> |
CH3 C 15.04 |
867 |
> |
CH2 C 14.03 |
868 |
> |
end AtomTypes |
869 |
> |
|
870 |
> |
begin LennardJonesAtomTypes |
871 |
> |
//name epsilon sigma |
872 |
> |
CH4 0.2941 3.73 |
873 |
> |
CH3 0.1947 3.75 |
874 |
> |
CH2 0.09140 3.95 |
875 |
> |
end LennardJonesAtomTypes |
876 |
> |
|
877 |
> |
begin BondTypes |
878 |
> |
//AT1 AT2 Type r0 k |
879 |
> |
CH3 CH3 Harmonic 1.526 260 |
880 |
> |
CH3 CH2 Harmonic 1.526 260 |
881 |
> |
CH2 CH2 Harmonic 1.526 260 |
882 |
> |
end BondTypes |
883 |
> |
|
884 |
> |
begin BendTypes |
885 |
> |
//AT1 AT2 AT3 Type theta0 k |
886 |
> |
CH3 CH2 CH3 Harmonic 114.0 124.19 |
887 |
> |
CH3 CH2 CH2 Harmonic 114.0 124.19 |
888 |
> |
CH2 CH2 CH2 Harmonic 114.0 124.19 |
889 |
> |
end BendTypes |
890 |
> |
|
891 |
> |
begin TorsionTypes |
892 |
> |
//AT1 AT2 AT3 AT4 Type |
893 |
> |
CH3 CH2 CH2 CH3 Trappe 0.0 0.70544 -0.13549 1.5723 |
894 |
> |
CH3 CH2 CH2 CH2 Trappe 0.0 0.70544 -0.13549 1.5723 |
895 |
> |
CH2 CH2 CH2 CH2 Trappe 0.0 0.70544 -0.13549 1.5723 |
896 |
> |
end TorsionTypes |
897 |
> |
\end{lstlisting} |
898 |
> |
|
899 |
> |
\subsection{\label{section:ffOptions}The Options block} |
900 |
> |
|
901 |
> |
The Options block defines properties governing how the force field |
902 |
> |
interactions are carried out. This block is delineated with the text |
903 |
> |
tags {\tt begin Options} and {\tt end Options}. Most options don't |
904 |
> |
need to be set as they come with fairly sensible default values, but |
905 |
> |
the various keywords and their possible values are given in Scheme |
906 |
> |
\ref{sch:optionsBlock}. |
907 |
> |
|
908 |
> |
\begin{lstlisting}[caption={[A force field Options block showing default values |
909 |
> |
for many force field options.] A force field Options block showing default values |
910 |
> |
for many force field options. Most of these options do not need to be |
911 |
> |
specified if the default values are working.}, |
912 |
> |
label={sch:optionsBlock}] |
913 |
> |
begin Options |
914 |
> |
Name = "alkane" // any string |
915 |
> |
vdWtype = "Lennard-Jones" |
916 |
> |
DistanceMixingRule = "arithmetic" // can also be "geometric" or "cubic" |
917 |
> |
DistanceType = "sigma" // can also be Rmin |
918 |
> |
EnergyMixingRule = "geometric" // can also be "arithmetic" or "hhg" |
919 |
> |
EnergyUnitScaling = 1.0 |
920 |
> |
MetallicEnergyUnitScaling = 1.0 |
921 |
> |
DistanceUnitScaling = 1.0 |
922 |
> |
AngleUnitScaling = 1.0 |
923 |
> |
TorsionAngleConvention = "180_is_trans" // can also be "0_is_trans" |
924 |
> |
vdW-12-scale = 0.0 |
925 |
> |
vdW-13-scale = 0.0 |
926 |
> |
vdW-14-scale = 0.0 |
927 |
> |
electrostatic-12-scale = 0.0 |
928 |
> |
electrostatic-13-scale = 0.0 |
929 |
> |
electrostatic-14-scale = 0.0 |
930 |
> |
GayBerneMu = 2.0 |
931 |
> |
GayBerneNu = 1.0 |
932 |
> |
EAMMixingMethod = "Johnson" // can also be "Daw" |
933 |
> |
end Options |
934 |
> |
\end{lstlisting} |
935 |
> |
|
936 |
> |
\subsection{\label{section:ffBase}The BaseAtomTypes block} |
937 |
> |
|
938 |
> |
An AtomType the primary data structure that OpenMD uses to store |
939 |
> |
static data about an atom. Things that belong to AtomType are |
940 |
> |
universal properties (i.e. does this atom have a fixed charge? What |
941 |
> |
is its mass?) Dynamic properties of an atom are not intended to be |
942 |
> |
properties of an atom type. A BaseAtomType can be used to build |
943 |
> |
extended sets of related atom types that all fall back to one |
944 |
> |
particular type. For example, one might want a series of atomTypes |
945 |
> |
that inherit from more basic types: |
946 |
> |
\begin{displaymath} |
947 |
> |
\mathtt{ALA-CA} \rightarrow \mathtt{CT} \rightarrow \mathtt{CSP3} \rightarrow \mathtt{C} |
948 |
> |
\end{displaymath} |
949 |
> |
where for each step to the right, the atomType falls back to more |
950 |
> |
primitive data. That is, the mass could be a property of the {\tt C} |
951 |
> |
type, while Lennard-Jones parameters could be properties of the {\tt |
952 |
> |
CSP3} type. {\tt CT} could have charge information and its own set |
953 |
> |
of Lennard-Jones parameter that override the CSP3 parameters. And the |
954 |
> |
{\tt ALA-CA} type might have specific torsion or charge information |
955 |
> |
that override the lower level types. A BaseAtomType contains only |
956 |
> |
information a primitive name and the mass of this atom type. |
957 |
> |
BaseAtomTypes can also be useful in creating files that can be easily |
958 |
> |
viewed in visualization programs. The {\tt Dump2XYZ} utility has the |
959 |
> |
ability to print out the names of the base atom types for displaying |
960 |
> |
simulations in Jmol or VMD. |
961 |
> |
|
962 |
> |
\begin{lstlisting}[caption={[A simple example of a BaseAtomType |
963 |
> |
block.] A simple example of a BaseAtomType block.}, |
964 |
> |
label={sch:baseAtomTypesBlock}] |
965 |
> |
begin BaseAtomTypes |
966 |
> |
//Name mass (amu) |
967 |
> |
H 1.0079 |
968 |
> |
O 15.9994 |
969 |
> |
Si 28.0855 |
970 |
> |
Al 26.981538 |
971 |
> |
Mg 24.3050 |
972 |
> |
Ca 40.078 |
973 |
> |
Fe 55.845 |
974 |
> |
Li 6.941 |
975 |
> |
Na 22.98977 |
976 |
> |
K 39.0983 |
977 |
> |
Cs 132.90545 |
978 |
> |
Ca 40.078 |
979 |
> |
Ba 137.327 |
980 |
> |
Cl 35.453 |
981 |
> |
end BaseAtomTypes |
982 |
> |
\end{lstlisting} |
983 |
> |
|
984 |
> |
\subsection{\label{section:ffAtom}The AtomTypes block} |
985 |
> |
|
986 |
> |
AtomTypes inherit most properties from BaseAtomTypes, but can override |
987 |
> |
their lower-level properties as well. Scheme \ref{sch:atomTypesBlock} |
988 |
> |
shows an example where multiple types of oxygen atoms can inherit mass |
989 |
> |
from the oxygen base type. |
990 |
> |
|
991 |
> |
\begin{lstlisting}[caption={[An example of a AtomTypes block.] A |
992 |
> |
simple example of an AtomType block which |
993 |
> |
shows how multiple types can inherit from the same base type.}, |
994 |
> |
label={sch:atomTypesBlock}] |
995 |
> |
begin AtomTypes |
996 |
> |
//Name baseatomtype |
997 |
> |
h* H |
998 |
> |
ho H |
999 |
> |
o* O |
1000 |
> |
oh O |
1001 |
> |
ob O |
1002 |
> |
obos O |
1003 |
> |
obts O |
1004 |
> |
obss O |
1005 |
> |
ohs O |
1006 |
> |
st Si |
1007 |
> |
ao Al |
1008 |
> |
at Al |
1009 |
> |
mgo Mg |
1010 |
> |
mgh Mg |
1011 |
> |
cao Ca |
1012 |
> |
cah Ca |
1013 |
> |
feo Fe |
1014 |
> |
lio Li |
1015 |
> |
end AtomTypes |
1016 |
> |
\end{lstlisting} |
1017 |
> |
|
1018 |
> |
\subsection{\label{section:ffDirectionalAtom}The DirectionalAtomTypes |
1019 |
> |
block} |
1020 |
> |
DirectionalAtoms have orientational degrees of freedom as well as |
1021 |
> |
translation, so they have moment of inertia tensors. |
1022 |
> |
|
1023 |
> |
\begin{lstlisting}[caption={[An example of a DirectionalAtomTypes block.] A |
1024 |
> |
simple example of a DirectionalAtomTypes block.}, |
1025 |
> |
label={sch:datomTypesBlock}] |
1026 |
> |
begin DirectionalAtomTypes |
1027 |
> |
//Name I_xx I_yy I_zz (All moments in (amu*Ang^2) |
1028 |
> |
SSD 1.7696 0.6145 1.1550 |
1029 |
> |
SSD_E 1.7696 0.6145 1.1550 |
1030 |
> |
GBC6H6 88.781 88.781 177.561 |
1031 |
> |
GBCH3OH 4.056 20.258 20.999 |
1032 |
> |
GBH2O 1.777 0.581 1.196 |
1033 |
> |
end DirectionalAtomTypes |
1034 |
> |
|
1035 |
> |
\end{lstlisting} |
1036 |
> |
|
1037 |
> |
|
1038 |
> |
\subsection{\label{section::ffAtomProperties}AtomType properties} |
1039 |
> |
\subsubsection{\label{section:ffLJ}The LennardJonesAtomTypes block} |
1040 |
> |
The most basic interatomic interaction implemented in {\sc OpenMD} is |
1041 |
> |
the Lennard-Jones potential, which mimics the van der Waals |
1042 |
> |
interaction at long distances and uses an empirical repulsion at short |
1043 |
|
distances. The Lennard-Jones potential is given by: |
1044 |
|
\begin{equation} |
1045 |
|
V_{\text{LJ}}(r_{ij}) = |
1051 |
|
\end{equation} |
1052 |
|
where $r_{ij}$ is the distance between particles $i$ and $j$, |
1053 |
|
$\sigma_{ij}$ scales the length of the interaction, and |
1054 |
< |
$\epsilon_{ij}$ scales the well depth of the potential. Scheme |
855 |
< |
\ref{sch:LJFF} gives an example meta-data file that |
856 |
< |
sets up a system of 108 Ar particles to be simulated using the |
857 |
< |
Lennard-Jones force field. |
1054 |
> |
$\epsilon_{ij}$ scales the well depth of the potential. |
1055 |
|
|
859 |
– |
\begin{lstlisting}[float,caption={[Invocation of the Lennard-Jones |
860 |
– |
force field] A sample startup file for a small Lennard-Jones |
861 |
– |
simulation.},label={sch:LJFF}] |
862 |
– |
<OpenMD> |
863 |
– |
<MetaData> |
864 |
– |
#include "argon.md" |
865 |
– |
|
866 |
– |
component{ |
867 |
– |
type = "Ar"; |
868 |
– |
nMol = 108; |
869 |
– |
} |
870 |
– |
|
871 |
– |
forceField = "LJ"; |
872 |
– |
</MetaData> |
873 |
– |
<Snapshot> // not shown in this scheme |
874 |
– |
</Snapshot> |
875 |
– |
</OpenMD> |
876 |
– |
\end{lstlisting} |
877 |
– |
|
1056 |
|
Interactions between dissimilar particles requires the generation of |
1057 |
|
cross term parameters for $\sigma$ and $\epsilon$. These parameters |
1058 |
|
are determined using the Lorentz-Berthelot mixing |
1067 |
|
\label{eq:epsilonMix} |
1068 |
|
\end{equation} |
1069 |
|
|
1070 |
< |
\section{\label{section:DUFF}Dipolar Unified-Atom Force Field} |
1071 |
< |
|
894 |
< |
The dipolar unified-atom force field ({\sc duff}) was developed to |
895 |
< |
simulate lipid bilayers. These types of simulations require a model |
896 |
< |
capable of forming bilayers, while still being sufficiently |
897 |
< |
computationally efficient to allow large systems ($\sim$100's of |
898 |
< |
phospholipids, $\sim$1000's of waters) to be simulated for long times |
899 |
< |
($\sim$10's of nanoseconds). With this goal in mind, {\sc duff} has no |
900 |
< |
point charges. Charge-neutral distributions are replaced with dipoles, |
901 |
< |
while most atoms and groups of atoms are reduced to Lennard-Jones |
902 |
< |
interaction sites. This simplification reduces the length scale of |
903 |
< |
long range interactions from $\frac{1}{r}$ to $\frac{1}{r^3}$, |
904 |
< |
removing the need for the computationally expensive Ewald |
905 |
< |
sum. Instead, Verlet neighbor-lists and cutoff radii are used for the |
906 |
< |
dipolar interactions, and, if desired, a reaction field may be added |
907 |
< |
to mimic longer range interactions. |
908 |
< |
|
909 |
< |
As an example, lipid head-groups in {\sc duff} are represented as |
910 |
< |
point dipole interaction sites. Placing a dipole at the head group's |
911 |
< |
center of mass mimics the charge separation found in common |
912 |
< |
phospholipid head groups such as phosphatidylcholine.\cite{Cevc87} |
913 |
< |
Additionally, a large Lennard-Jones site is located at the |
914 |
< |
pseudoatom's center of mass. The model is illustrated by the red atom |
915 |
< |
in Fig.~\ref{fig:lipidModel}. The water model we use to |
916 |
< |
complement the dipoles of the lipids is a |
917 |
< |
reparameterization\cite{fennell04} of the soft sticky dipole (SSD) |
918 |
< |
model of Ichiye |
919 |
< |
\emph{et al.}\cite{liu96:new_model} |
920 |
< |
|
921 |
< |
\begin{figure} |
922 |
< |
\centering |
923 |
< |
\includegraphics[width=\linewidth]{lipidModel.pdf} |
924 |
< |
\caption[A representation of a lipid model in {\sc duff}]{A |
925 |
< |
representation of the lipid model. $\phi$ is the torsion angle, |
926 |
< |
$\theta$ is the bend angle, and $\mu$ is the dipole moment of the head |
927 |
< |
group.} |
928 |
< |
\label{fig:lipidModel} |
929 |
< |
\end{figure} |
930 |
< |
|
931 |
< |
A set of scalable parameters has been used to model the alkyl groups |
932 |
< |
with Lennard-Jones sites. For this, parameters from the TraPPE force |
933 |
< |
field of Siepmann \emph{et al.}\cite{Siepmann1998} have been |
934 |
< |
utilized. TraPPE is a unified-atom representation of n-alkanes which |
935 |
< |
is parametrized against phase equilibria using Gibbs ensemble Monte |
936 |
< |
Carlo simulation techniques.\cite{Siepmann1998} One of the advantages |
937 |
< |
of TraPPE is that it generalizes the types of atoms in an alkyl chain |
938 |
< |
to keep the number of pseudoatoms to a minimum; thus, the parameters |
939 |
< |
for a unified atom such as $\text{CH}_2$ do not change depending on |
940 |
< |
what species are bonded to it. |
941 |
< |
|
942 |
< |
As is required by TraPPE, {\sc duff} also constrains all bonds to be |
943 |
< |
of fixed length. Typically, bond vibrations are the fastest motions in |
944 |
< |
a molecular dynamic simulation. With these vibrations present, small |
945 |
< |
time steps between force evaluations must be used to ensure adequate |
946 |
< |
energy conservation in the bond degrees of freedom. By constraining |
947 |
< |
the bond lengths, larger time steps may be used when integrating the |
948 |
< |
equations of motion. A simulation using {\sc duff} is illustrated in |
949 |
< |
Scheme \ref{sch:DUFF}. |
950 |
< |
|
951 |
< |
\begin{lstlisting}[float,caption={[Invocation of {\sc duff}]A portion |
952 |
< |
of a startup file showing a simulation utilizing {\sc |
953 |
< |
duff}},label={sch:DUFF}] |
954 |
< |
<OpenMD> |
955 |
< |
<MetaData> |
956 |
< |
#include "water.md" |
957 |
< |
#include "lipid.md" |
958 |
< |
|
959 |
< |
component{ |
960 |
< |
type = "simpleLipid_16"; |
961 |
< |
nMol = 60; |
962 |
< |
} |
963 |
< |
|
964 |
< |
component{ |
965 |
< |
type = "SSD_water"; |
966 |
< |
nMol = 1936; |
967 |
< |
} |
968 |
< |
|
969 |
< |
forceField = "DUFF"; |
970 |
< |
</MetaData> |
971 |
< |
<Snapshot> // not shown in this scheme |
972 |
< |
</Snapshot> |
973 |
< |
</OpenMD> |
974 |
< |
\end{lstlisting} |
975 |
< |
|
976 |
< |
\subsection{\label{section:energyFunctions}{\sc duff} Energy Functions} |
977 |
< |
|
978 |
< |
The total potential energy function in {\sc duff} is |
979 |
< |
\begin{equation} |
980 |
< |
V = \sum^{N}_{I=1} V^{I}_{\text{Internal}} |
981 |
< |
+ \sum^{N-1}_{I=1} \sum_{J>I} V^{IJ}_{\text{Cross}}, |
982 |
< |
\label{eq:totalPotential} |
983 |
< |
\end{equation} |
984 |
< |
where $V^{I}_{\text{Internal}}$ is the internal potential of molecule $I$: |
985 |
< |
\begin{equation} |
986 |
< |
V^{I}_{\text{Internal}} = |
987 |
< |
\sum_{\theta_{ijk} \in I} V_{\text{bend}}(\theta_{ijk}) |
988 |
< |
+ \sum_{\phi_{ijkl} \in I} V_{\text{torsion}}(\phi_{ijkl}) |
989 |
< |
+ \sum_{i \in I} \sum_{(j>i+4) \in I} |
990 |
< |
\biggl[ V_{\text{LJ}}(r_{ij}) + V_{\text{dipole}} |
991 |
< |
(\mathbf{r}_{ij},\boldsymbol{\Omega}_{i},\boldsymbol{\Omega}_{j}) |
992 |
< |
\biggr]. |
993 |
< |
\label{eq:internalPotential} |
994 |
< |
\end{equation} |
995 |
< |
Here $V_{\text{bend}}$ is the bend potential for all 1, 3 bonded pairs |
996 |
< |
within the molecule $I$, and $V_{\text{torsion}}$ is the torsion |
997 |
< |
potential for all 1, 4 bonded pairs. The pairwise portions of the |
998 |
< |
non-bonded interactions are excluded for atom pairs that are involved |
999 |
< |
in the smae bond, bend, or torsion. All other atom pairs within a |
1000 |
< |
molecule are subject to the LJ pair potential. |
1001 |
< |
|
1002 |
< |
The bend potential of a molecule is represented by the following function: |
1003 |
< |
\begin{equation} |
1004 |
< |
V_{\text{bend}}(\theta_{ijk}) = k_{\theta}( \theta_{ijk} - \theta_0 |
1005 |
< |
)^2, \label{eq:bendPot} |
1006 |
< |
\end{equation} |
1007 |
< |
where $\theta_{ijk}$ is the angle defined by atoms $i$, $j$, and $k$ |
1008 |
< |
(see Fig.~\ref{fig:lipidModel}), $\theta_0$ is the equilibrium |
1009 |
< |
bond angle, and $k_{\theta}$ is the force constant which determines the |
1010 |
< |
strength of the harmonic bend. The parameters for $k_{\theta}$ and |
1011 |
< |
$\theta_0$ are borrowed from those in TraPPE.\cite{Siepmann1998} |
1012 |
< |
|
1013 |
< |
The torsion potential and parameters are also borrowed from TraPPE. It is |
1014 |
< |
of the form: |
1015 |
< |
\begin{equation} |
1016 |
< |
V_{\text{torsion}}(\phi) = c_1[1 + \cos \phi] |
1017 |
< |
+ c_2[1 + \cos(2\phi)] |
1018 |
< |
+ c_3[1 + \cos(3\phi)], |
1019 |
< |
\label{eq:origTorsionPot} |
1020 |
< |
\end{equation} |
1021 |
< |
where: |
1022 |
< |
\begin{equation} |
1023 |
< |
\cos\phi = (\hat{\mathbf{r}}_{ij} \times \hat{\mathbf{r}}_{jk}) \cdot |
1024 |
< |
(\hat{\mathbf{r}}_{jk} \times \hat{\mathbf{r}}_{kl}). |
1025 |
< |
\label{eq:torsPhi} |
1026 |
< |
\end{equation} |
1027 |
< |
Here, $\hat{\mathbf{r}}_{\alpha\beta}$ are the set of unit bond |
1028 |
< |
vectors between atoms $i$, $j$, $k$, and $l$. For computational |
1029 |
< |
efficiency, the torsion potential has been recast after the method of |
1030 |
< |
{\sc charmm},\cite{Brooks83} in which the angle series is converted to |
1031 |
< |
a power series of the form: |
1032 |
< |
\begin{equation} |
1033 |
< |
V_{\text{torsion}}(\phi) = |
1034 |
< |
k_3 \cos^3 \phi + k_2 \cos^2 \phi + k_1 \cos \phi + k_0, |
1035 |
< |
\label{eq:torsionPot} |
1036 |
< |
\end{equation} |
1037 |
< |
where: |
1038 |
< |
\begin{align*} |
1039 |
< |
k_0 &= c_1 + c_3, \\ |
1040 |
< |
k_1 &= c_1 - 3c_3, \\ |
1041 |
< |
k_2 &= 2 c_2, \\ |
1042 |
< |
k_3 &= 4c_3. |
1043 |
< |
\end{align*} |
1044 |
< |
By recasting the potential as a power series, repeated trigonometric |
1045 |
< |
evaluations are avoided during the calculation of the potential |
1046 |
< |
energy. |
1047 |
< |
|
1048 |
< |
|
1049 |
< |
The cross potential between molecules $I$ and $J$, |
1050 |
< |
$V^{IJ}_{\text{Cross}}$, is as follows: |
1051 |
< |
\begin{equation} |
1052 |
< |
V^{IJ}_{\text{Cross}} = |
1053 |
< |
\sum_{i \in I} \sum_{j \in J} |
1054 |
< |
\biggl[ V_{\text{LJ}}(r_{ij}) + V_{\text{dipole}} |
1055 |
< |
(\mathbf{r}_{ij},\boldsymbol{\Omega}_{i},\boldsymbol{\Omega}_{j}) |
1056 |
< |
+ V_{\text{sticky}} |
1057 |
< |
(\mathbf{r}_{ij},\boldsymbol{\Omega}_{i},\boldsymbol{\Omega}_{j}) |
1058 |
< |
\biggr], |
1059 |
< |
\label{eq:crossPotentail} |
1060 |
< |
\end{equation} |
1061 |
< |
where $V_{\text{LJ}}$ is the Lennard Jones potential, |
1062 |
< |
$V_{\text{dipole}}$ is the dipole dipole potential, and |
1063 |
< |
$V_{\text{sticky}}$ is the sticky potential defined by the SSD model |
1064 |
< |
(Sec.~\ref{section:SSD}). Note that not all atom types include all |
1065 |
< |
interactions. |
1066 |
< |
|
1070 |
> |
\subsubsection{\label{section:ffCharge}The ChargeAtomTypes block} |
1071 |
> |
\subsubsection{\label{section:ffMultipole}The MultipoleAtomTypes block} |
1072 |
|
The dipole-dipole potential has the following form: |
1073 |
|
\begin{equation} |
1074 |
|
V_{\text{dipole}}(\mathbf{r}_{ij},\boldsymbol{\Omega}_{i}, |
1088 |
|
the unit vector pointing along $\mathbf{r}_{ij}$ |
1089 |
|
($\boldsymbol{\hat{r}}_{ij}=\mathbf{r}_{ij}/|\mathbf{r}_{ij}|$). |
1090 |
|
|
1091 |
< |
\subsection{\label{section:SSD}The {\sc duff} Water Models: SSD/E |
1092 |
< |
and SSD/RF} |
1091 |
> |
\subsubsection{\label{section:ffGB}The FluctuatingChargeAtomTypes block} |
1092 |
> |
\subsubsection{\label{section:ffPol}The PolarizableAtomTypes block} |
1093 |
> |
\subsubsection{\label{section:ffGB}The GayBerneAtomTypes block} |
1094 |
> |
\subsubsection{\label{section:ffSticky}The StickyAtomTypes block} |
1095 |
|
|
1096 |
< |
In the interest of computational efficiency, the default solvent used |
1097 |
< |
by {\sc OpenMD} is the extended Soft Sticky Dipole (SSD/E) water |
1098 |
< |
model.\cite{fennell04} The original SSD was developed by Ichiye |
1099 |
< |
\emph{et al.}\cite{liu96:new_model} as a modified form of the hard-sphere |
1093 |
< |
water model proposed by Bratko, Blum, and |
1096 |
> |
One of the solvents used by {\sc OpenMD} is the extended Soft Sticky |
1097 |
> |
Dipole (SSD/E) water model.\cite{fennell04} The original SSD was |
1098 |
> |
developed by Ichiye \emph{et al.}\cite{liu96:new_model} as a modified |
1099 |
> |
form of the hard-sphere water model proposed by Bratko, Blum, and |
1100 |
|
Luzar.\cite{Bratko85,Bratko95} It consists of a single point dipole |
1101 |
|
with a Lennard-Jones core and a sticky potential that directs the |
1102 |
|
particles to assume the proper hydrogen bond orientation in the first |
1161 |
|
\label{fig:ssd} |
1162 |
|
\end{figure} |
1163 |
|
|
1158 |
– |
|
1164 |
|
Since SSD/E is a single-point {\it dipolar} model, the force |
1165 |
|
calculations are simplified significantly relative to the standard |
1166 |
|
{\it charged} multi-point models. In the original Monte Carlo |
1190 |
|
and the drawbacks and benefits of the different density corrected SSD |
1191 |
|
models can be found in reference~\cite{fennell04}. |
1192 |
|
|
1193 |
< |
\section{\label{section:WATER}The {\sc water} Force Field} |
1194 |
< |
|
1190 |
< |
In addition to the {\sc duff} force field's solvent description, a |
1191 |
< |
separate {\sc water} force field has been included for simulating most |
1192 |
< |
of the common rigid-body water models. This force field includes the |
1193 |
< |
simple and point-dipolar models (SSD, SSD1, SSD/E, SSD/RF, and DPD |
1194 |
< |
water), as well as the common charge-based models (SPC, SPC/E, TIP3P, |
1195 |
< |
TIP4P, and |
1196 |
< |
TIP5P).\cite{liu96:new_model,Ichiye03,fennell04,Marrink01,Berendsen81,Berendsen87,Jorgensen83,Mahoney00} |
1197 |
< |
In order to handle these models, charge-charge interactions were |
1198 |
< |
included in the force-loop: |
1199 |
< |
\begin{equation} |
1200 |
< |
V_{\text{charge}}(r_{ij}) = \sum_{ij}\frac{q_iq_je^2}{r_{ij}}, |
1201 |
< |
\end{equation} |
1202 |
< |
where $q$ represents the charge on particle $i$ or $j$, and $e$ is the |
1203 |
< |
charge of an electron in Coulombs. The charge-charge interaction |
1204 |
< |
support is rudimentary in the current version of {\sc OpenMD}. As with |
1205 |
< |
the other pair interactions, charges can be simulated with a pure |
1206 |
< |
cutoff or a reaction field. The various methods for performing the |
1207 |
< |
Ewald summation have not yet been included. The {\sc water} force |
1208 |
< |
field can be easily expanded through modification of the {\sc water} |
1209 |
< |
force field file ({\tt WATER.frc}). By adding atom types and inserting |
1210 |
< |
the appropriate parameters, it is possible to extend the force field |
1211 |
< |
to handle rigid molecules other than water. |
1212 |
< |
|
1213 |
< |
\section{\label{section:eam}Embedded Atom Method} |
1214 |
< |
|
1193 |
> |
\subsection{\label{section::ffMetals}Metallic Atom Types} |
1194 |
> |
\subsubsection{\label{section:ffEAM}The EAMAtomTypes block} |
1195 |
|
{\sc OpenMD} implements a potential that describes bonding in |
1196 |
|
transition metal |
1197 |
|
systems.~\cite{Finnis84,Ercolessi88,Chen90,Qi99,Ercolessi02} This |
1269 |
|
$\mbox{kcal mol}^{-1}$ as in the rest of the {\sc OpenMD} force field |
1270 |
|
files. |
1271 |
|
|
1272 |
< |
\section{\label{section:sc}The Sutton-Chen Force Field} |
1272 |
> |
\subsubsection{\label{section:ffSC}The SuttonChenAtomTypes block} |
1273 |
|
|
1274 |
|
The Sutton-Chen ({\sc sc})~\cite{Chen90} potential has been used to |
1275 |
|
study a wide range of phenomena in metals. Although it is similar in |
1305 |
|
quantum-adapted variant of the {\sc sc} potential, the user would add |
1306 |
|
the {\tt forceFieldVariant = "QSC";} line to the meta-data file. |
1307 |
|
|
1308 |
+ |
\subsection{\label{section::ffShortRange}Short Range Interactions} |
1309 |
+ |
\subsubsection{\label{section:ffBond}The BondTypes block} |
1310 |
+ |
\subsubsection{\label{section:ffBend}The BendTypes block} |
1311 |
+ |
A harmonic bend potential is represented by the following function: |
1312 |
+ |
\begin{equation} |
1313 |
+ |
V_{\text{bend}}(\theta_{ijk}) = k_{\theta}( \theta_{ijk} - \theta_0 |
1314 |
+ |
)^2, \label{eq:bendPot} |
1315 |
+ |
\end{equation} |
1316 |
+ |
where $\theta_{ijk}$ is the angle defined by atoms $i$, $j$, and $k$, |
1317 |
+ |
$\theta_0$ is the equilibrium bond angle, and $k_{\theta}$ is the |
1318 |
+ |
force constant which determines the strength of the harmonic bend. |
1319 |
+ |
|
1320 |
+ |
\subsubsection{\label{section:ffTorsion}The TorsionTypes block} |
1321 |
+ |
The torsion potential is often represented as a cosine series of the |
1322 |
+ |
form: |
1323 |
+ |
\begin{equation} |
1324 |
+ |
V_{\text{torsion}}(\phi) = c_1[1 + \cos \phi] |
1325 |
+ |
+ c_2[1 + \cos(2\phi)] |
1326 |
+ |
+ c_3[1 + \cos(3\phi)], |
1327 |
+ |
\label{eq:origTorsionPot} |
1328 |
+ |
\end{equation} |
1329 |
+ |
where: |
1330 |
+ |
\begin{equation} |
1331 |
+ |
\cos\phi = (\hat{\mathbf{r}}_{ij} \times \hat{\mathbf{r}}_{jk}) \cdot |
1332 |
+ |
(\hat{\mathbf{r}}_{jk} \times \hat{\mathbf{r}}_{kl}). |
1333 |
+ |
\label{eq:torsPhi} |
1334 |
+ |
\end{equation} |
1335 |
+ |
Here, $\hat{\mathbf{r}}_{\alpha\beta}$ are the set of unit bond |
1336 |
+ |
vectors between atoms $i$, $j$, $k$, and $l$. For computational |
1337 |
+ |
efficiency, the torsion potential has been recast after the method of |
1338 |
+ |
{\sc charmm},\cite{Brooks83} in which the angle series is converted to |
1339 |
+ |
a power series of the form: |
1340 |
+ |
\begin{equation} |
1341 |
+ |
V_{\text{torsion}}(\phi) = |
1342 |
+ |
k_3 \cos^3 \phi + k_2 \cos^2 \phi + k_1 \cos \phi + k_0, |
1343 |
+ |
\label{eq:torsionPot} |
1344 |
+ |
\end{equation} |
1345 |
+ |
where: |
1346 |
+ |
\begin{align*} |
1347 |
+ |
k_0 &= c_1 + c_3, \\ |
1348 |
+ |
k_1 &= c_1 - 3c_3, \\ |
1349 |
+ |
k_2 &= 2 c_2, \\ |
1350 |
+ |
k_3 &= 4c_3. |
1351 |
+ |
\end{align*} |
1352 |
+ |
By recasting the potential as a power series, repeated trigonometric |
1353 |
+ |
evaluations are avoided during the calculation of the potential |
1354 |
+ |
energy. |
1355 |
+ |
|
1356 |
+ |
\subsubsection{\label{section:ffInversion}The InversionTypes block} |
1357 |
+ |
\subsection{\label{section::ffLongRange}Long Range Interactions} |
1358 |
+ |
\subsubsection{\label{section:ffInversion}The NonBondedInteraction block} |
1359 |
+ |
|
1360 |
+ |
|
1361 |
+ |
|
1362 |
+ |
(see Fig.~\ref{fig:lipidModel}), The parameters for $k_{\theta}$ and |
1363 |
+ |
$\theta_0$ are borrowed from those in TraPPE.\cite{Siepmann1998} |
1364 |
+ |
|
1365 |
+ |
Calculating the long-range (non-bonded) potential involves a sum over |
1366 |
+ |
all pairs of atoms (except for those atoms which are involved in a |
1367 |
+ |
bond, bend, or torsion with each other). If done poorly, calculating |
1368 |
+ |
the the long-range interactions for $N$ atoms would involve $N(N-1)/2$ |
1369 |
+ |
evaluations of atomic distances. To reduce the number of distance |
1370 |
+ |
evaluations between pairs of atoms, {\sc OpenMD} allows the use of |
1371 |
+ |
switched cutoffs with Verlet neighbor lists.\cite{Allen87} Neutral |
1372 |
+ |
groups which contain charges will exhibit pathological forces unless |
1373 |
+ |
the cutoff is applied to the neutral groups evenly instead of to the |
1374 |
+ |
individual atoms.\cite{leach01:mm} {\sc OpenMD} allows users to |
1375 |
+ |
specify cutoff groups which may contain an arbitrary number of atoms |
1376 |
+ |
in the molecule. Atoms in a cutoff group are treated as a single unit |
1377 |
+ |
for the evaluation of the switching function: |
1378 |
+ |
\begin{equation} |
1379 |
+ |
V_{\mathrm{long-range}} = \sum_{a} \sum_{b>a} s(r_{ab}) \sum_{i \in a} \sum_{j \in b} V_{ij}(r_{ij}), |
1380 |
+ |
\end{equation} |
1381 |
+ |
where $r_{ab}$ is the distance between the centers of mass of the two |
1382 |
+ |
cutoff groups ($a$ and $b$). |
1383 |
+ |
|
1384 |
+ |
The sums over $a$ and $b$ are over the cutoff groups that are present |
1385 |
+ |
in the simulation. Atoms which are not explicitly defined as members |
1386 |
+ |
of a {\tt cutoffGroup} are treated as a group consisting of only one |
1387 |
+ |
atom. The switching function, $s(r)$ is the standard cubic switching |
1388 |
+ |
function, |
1389 |
+ |
\begin{equation} |
1390 |
+ |
S(r) = |
1391 |
+ |
\begin{cases} |
1392 |
+ |
1 & \text{if $r \le r_{\text{sw}}$},\\ |
1393 |
+ |
\frac{(r_{\text{cut}} + 2r - 3r_{\text{sw}})(r_{\text{cut}} - r)^2} |
1394 |
+ |
{(r_{\text{cut}} - r_{\text{sw}})^3} |
1395 |
+ |
& \text{if $r_{\text{sw}} < r \le r_{\text{cut}}$}, \\ |
1396 |
+ |
0 & \text{if $r > r_{\text{cut}}$.} |
1397 |
+ |
\end{cases} |
1398 |
+ |
\label{eq:dipoleSwitching} |
1399 |
+ |
\end{equation} |
1400 |
+ |
Here, $r_{\text{sw}}$ is the {\tt switchingRadius}, or the distance |
1401 |
+ |
beyond which interactions are reduced, and $r_{\text{cut}}$ is the |
1402 |
+ |
{\tt cutoffRadius}, or the distance at which interactions are |
1403 |
+ |
truncated. |
1404 |
+ |
|
1405 |
+ |
Users of {\sc OpenMD} do not need to specify the {\tt cutoffRadius} or |
1406 |
+ |
{\tt switchingRadius}. In simulations containing only Lennard-Jones |
1407 |
+ |
atoms, the cutoff radius has a default value of $2.5\sigma_{ii}$, |
1408 |
+ |
where $\sigma_{ii}$ is the largest Lennard-Jones length parameter |
1409 |
+ |
present in the simulation. In simulations containing charged or |
1410 |
+ |
dipolar atoms, the default cutoff radius is $15 \mbox{\AA}$. |
1411 |
+ |
|
1412 |
+ |
The {\tt switchingRadius} is set to a default value of 95\% of the |
1413 |
+ |
{\tt cutoffRadius}. In the special case of a simulation containing |
1414 |
+ |
{\it only} Lennard-Jones atoms, the default switching radius takes the |
1415 |
+ |
same value as the cutoff radius, and {\sc OpenMD} will use a shifted |
1416 |
+ |
potential to remove discontinuities in the potential at the cutoff. |
1417 |
+ |
Both radii may be specified in the meta-data file. |
1418 |
+ |
|
1419 |
+ |
Force fields can be added to {\sc OpenMD}, although it comes with a few |
1420 |
+ |
simple examples (Lennard-Jones, {\sc duff}, {\sc water}, and {\sc |
1421 |
+ |
eam}) which are explained in the following sections. |
1422 |
+ |
|
1423 |
+ |
\section{\label{sec:LJPot}The Lennard Jones Force Field} |
1424 |
+ |
|
1425 |
+ |
Scheme |
1426 |
+ |
\ref{sch:LJFF} gives an example meta-data file that |
1427 |
+ |
sets up a system of 108 Ar particles to be simulated using the |
1428 |
+ |
Lennard-Jones force field. |
1429 |
+ |
|
1430 |
+ |
\begin{lstlisting}[float,caption={[Invocation of the Lennard-Jones |
1431 |
+ |
force field] A sample startup file for a small Lennard-Jones |
1432 |
+ |
simulation.},label={sch:LJFF}] |
1433 |
+ |
<OpenMD> |
1434 |
+ |
<MetaData> |
1435 |
+ |
#include "argon.md" |
1436 |
+ |
|
1437 |
+ |
component{ |
1438 |
+ |
type = "Ar"; |
1439 |
+ |
nMol = 108; |
1440 |
+ |
} |
1441 |
+ |
|
1442 |
+ |
forceField = "LJ"; |
1443 |
+ |
</MetaData> |
1444 |
+ |
<Snapshot> // not shown in this scheme |
1445 |
+ |
</Snapshot> |
1446 |
+ |
</OpenMD> |
1447 |
+ |
\end{lstlisting} |
1448 |
+ |
|
1449 |
+ |
|
1450 |
+ |
\section{\label{section:DUFF}Dipolar Unified-Atom Force Field} |
1451 |
+ |
|
1452 |
+ |
The dipolar unified-atom force field ({\sc duff}) was developed to |
1453 |
+ |
simulate lipid bilayers. These types of simulations require a model |
1454 |
+ |
capable of forming bilayers, while still being sufficiently |
1455 |
+ |
computationally efficient to allow large systems ($\sim$100's of |
1456 |
+ |
phospholipids, $\sim$1000's of waters) to be simulated for long times |
1457 |
+ |
($\sim$10's of nanoseconds). With this goal in mind, {\sc duff} has no |
1458 |
+ |
point charges. Charge-neutral distributions are replaced with dipoles, |
1459 |
+ |
while most atoms and groups of atoms are reduced to Lennard-Jones |
1460 |
+ |
interaction sites. This simplification reduces the length scale of |
1461 |
+ |
long range interactions from $\frac{1}{r}$ to $\frac{1}{r^3}$, |
1462 |
+ |
removing the need for the computationally expensive Ewald |
1463 |
+ |
sum. Instead, Verlet neighbor-lists and cutoff radii are used for the |
1464 |
+ |
dipolar interactions, and, if desired, a reaction field may be added |
1465 |
+ |
to mimic longer range interactions. |
1466 |
+ |
|
1467 |
+ |
As an example, lipid head-groups in {\sc duff} are represented as |
1468 |
+ |
point dipole interaction sites. Placing a dipole at the head group's |
1469 |
+ |
center of mass mimics the charge separation found in common |
1470 |
+ |
phospholipid head groups such as phosphatidylcholine.\cite{Cevc87} |
1471 |
+ |
Additionally, a large Lennard-Jones site is located at the |
1472 |
+ |
pseudoatom's center of mass. The model is illustrated by the red atom |
1473 |
+ |
in Fig.~\ref{fig:lipidModel}. The water model we use to |
1474 |
+ |
complement the dipoles of the lipids is a |
1475 |
+ |
reparameterization\cite{fennell04} of the soft sticky dipole (SSD) |
1476 |
+ |
model of Ichiye |
1477 |
+ |
\emph{et al.}\cite{liu96:new_model} |
1478 |
+ |
|
1479 |
+ |
\begin{figure} |
1480 |
+ |
\centering |
1481 |
+ |
\includegraphics[width=\linewidth]{lipidModel.pdf} |
1482 |
+ |
\caption[A representation of a lipid model in {\sc duff}]{A |
1483 |
+ |
representation of the lipid model. $\phi$ is the torsion angle, |
1484 |
+ |
$\theta$ is the bend angle, and $\mu$ is the dipole moment of the head |
1485 |
+ |
group.} |
1486 |
+ |
\label{fig:lipidModel} |
1487 |
+ |
\end{figure} |
1488 |
+ |
|
1489 |
+ |
A set of scalable parameters has been used to model the alkyl groups |
1490 |
+ |
with Lennard-Jones sites. For this, parameters from the TraPPE force |
1491 |
+ |
field of Siepmann \emph{et al.}\cite{Siepmann1998} have been |
1492 |
+ |
utilized. TraPPE is a unified-atom representation of n-alkanes which |
1493 |
+ |
is parametrized against phase equilibria using Gibbs ensemble Monte |
1494 |
+ |
Carlo simulation techniques.\cite{Siepmann1998} One of the advantages |
1495 |
+ |
of TraPPE is that it generalizes the types of atoms in an alkyl chain |
1496 |
+ |
to keep the number of pseudoatoms to a minimum; thus, the parameters |
1497 |
+ |
for a unified atom such as $\text{CH}_2$ do not change depending on |
1498 |
+ |
what species are bonded to it. |
1499 |
+ |
|
1500 |
+ |
As is required by TraPPE, {\sc duff} also constrains all bonds to be |
1501 |
+ |
of fixed length. Typically, bond vibrations are the fastest motions in |
1502 |
+ |
a molecular dynamic simulation. With these vibrations present, small |
1503 |
+ |
time steps between force evaluations must be used to ensure adequate |
1504 |
+ |
energy conservation in the bond degrees of freedom. By constraining |
1505 |
+ |
the bond lengths, larger time steps may be used when integrating the |
1506 |
+ |
equations of motion. A simulation using {\sc duff} is illustrated in |
1507 |
+ |
Scheme \ref{sch:DUFF}. |
1508 |
+ |
|
1509 |
+ |
\begin{lstlisting}[float,caption={[Invocation of {\sc duff}]A portion |
1510 |
+ |
of a startup file showing a simulation utilizing {\sc |
1511 |
+ |
duff}},label={sch:DUFF}] |
1512 |
+ |
<OpenMD> |
1513 |
+ |
<MetaData> |
1514 |
+ |
#include "water.md" |
1515 |
+ |
#include "lipid.md" |
1516 |
+ |
|
1517 |
+ |
component{ |
1518 |
+ |
type = "simpleLipid_16"; |
1519 |
+ |
nMol = 60; |
1520 |
+ |
} |
1521 |
+ |
|
1522 |
+ |
component{ |
1523 |
+ |
type = "SSD_water"; |
1524 |
+ |
nMol = 1936; |
1525 |
+ |
} |
1526 |
+ |
|
1527 |
+ |
forceField = "DUFF"; |
1528 |
+ |
</MetaData> |
1529 |
+ |
<Snapshot> // not shown in this scheme |
1530 |
+ |
</Snapshot> |
1531 |
+ |
</OpenMD> |
1532 |
+ |
\end{lstlisting} |
1533 |
+ |
|
1534 |
+ |
|
1535 |
+ |
|
1536 |
+ |
The cross potential between molecules $I$ and $J$, |
1537 |
+ |
$V^{IJ}_{\text{Cross}}$, is as follows: |
1538 |
+ |
\begin{equation} |
1539 |
+ |
V^{IJ}_{\text{Cross}} = |
1540 |
+ |
\sum_{i \in I} \sum_{j \in J} |
1541 |
+ |
\biggl[ V_{\text{LJ}}(r_{ij}) + V_{\text{dipole}} |
1542 |
+ |
(\mathbf{r}_{ij},\boldsymbol{\Omega}_{i},\boldsymbol{\Omega}_{j}) |
1543 |
+ |
+ V_{\text{sticky}} |
1544 |
+ |
(\mathbf{r}_{ij},\boldsymbol{\Omega}_{i},\boldsymbol{\Omega}_{j}) |
1545 |
+ |
\biggr], |
1546 |
+ |
\label{eq:crossPotentail} |
1547 |
+ |
\end{equation} |
1548 |
+ |
where $V_{\text{LJ}}$ is the Lennard Jones potential, |
1549 |
+ |
$V_{\text{dipole}}$ is the dipole dipole potential, and |
1550 |
+ |
$V_{\text{sticky}}$ is the sticky potential defined by the SSD model |
1551 |
+ |
(Sec.~\ref{section:SSD}). Note that not all atom types include all |
1552 |
+ |
interactions. |
1553 |
+ |
|
1554 |
+ |
|
1555 |
+ |
\section{\label{section:WATER}The {\sc water} Force Field} |
1556 |
+ |
|
1557 |
+ |
In addition to the {\sc duff} force field's solvent description, a |
1558 |
+ |
separate {\sc water} force field has been included for simulating most |
1559 |
+ |
of the common rigid-body water models. This force field includes the |
1560 |
+ |
simple and point-dipolar models (SSD, SSD1, SSD/E, SSD/RF, and DPD |
1561 |
+ |
water), as well as the common charge-based models (SPC, SPC/E, TIP3P, |
1562 |
+ |
TIP4P, and |
1563 |
+ |
TIP5P).\cite{liu96:new_model,Ichiye03,fennell04,Marrink01,Berendsen81,Berendsen87,Jorgensen83,Mahoney00} |
1564 |
+ |
In order to handle these models, charge-charge interactions were |
1565 |
+ |
included in the force-loop: |
1566 |
+ |
\begin{equation} |
1567 |
+ |
V_{\text{charge}}(r_{ij}) = \sum_{ij}\frac{q_iq_je^2}{r_{ij}}, |
1568 |
+ |
\end{equation} |
1569 |
+ |
where $q$ represents the charge on particle $i$ or $j$, and $e$ is the |
1570 |
+ |
charge of an electron in Coulombs. The charge-charge interaction |
1571 |
+ |
support is rudimentary in the current version of {\sc OpenMD}. As with |
1572 |
+ |
the other pair interactions, charges can be simulated with a pure |
1573 |
+ |
cutoff or a reaction field. The various methods for performing the |
1574 |
+ |
Ewald summation have not yet been included. The {\sc water} force |
1575 |
+ |
field can be easily expanded through modification of the {\sc water} |
1576 |
+ |
force field file ({\tt WATER.frc}). By adding atom types and inserting |
1577 |
+ |
the appropriate parameters, it is possible to extend the force field |
1578 |
+ |
to handle rigid molecules other than water. |
1579 |
+ |
|
1580 |
+ |
|
1581 |
+ |
\section{\label{section:sc}The Sutton-Chen Force Field} |
1582 |
+ |
|
1583 |
+ |
|
1584 |
|
\section{\label{section:clay}The CLAY force field} |
1585 |
|
|
1586 |
|
The {\sc clay} force field is based on an ionic (nonbonded) |
2877 |
|
|
2878 |
|
\section{Constant Pressure without periodic boundary conditions (The LangevinHull)} |
2879 |
|
|
2880 |
< |
The Langevin Hull uses an external bath at a fixed constant pressure |
2880 |
> |
The Langevin Hull\cite{Vardeman2011} uses an external bath at a fixed constant pressure |
2881 |
|
($P$) and temperature ($T$) with an effective solvent viscosity |
2882 |
|
($\eta$). This bath interacts only with the objects on the exterior |
2883 |
|
hull of the system. Defining the hull of the atoms in a simulation is |
3189 |
|
\label{table:zconParams} |
3190 |
|
\end{longtable} |
3191 |
|
|
3192 |
< |
\chapter{\label{section:restraints}Restraints} |
3193 |
< |
Restraints are external potentials that are added to a system to keep |
3194 |
< |
particular molecules or collections of particles close to some |
3195 |
< |
reference structure. A restraint can be a collective |
3192 |
> |
% \chapter{\label{section:restraints}Restraints} |
3193 |
> |
% Restraints are external potentials that are added to a system to |
3194 |
> |
% keep particular molecules or collections of particles close to some |
3195 |
> |
% reference structure. A restraint can be a collective |
3196 |
|
|
3197 |
|
\chapter{\label{section:thermInt}Thermodynamic Integration} |
3198 |
|
|
3332 |
|
\label{table:thermIntParams} |
3333 |
|
\end{longtable} |
3334 |
|
|
3335 |
+ |
\chapter{\label{section:rnemd}Reverse Non-Equilibrium Molecular Dynamics (RNEMD)} |
3336 |
|
|
3337 |
+ |
There are many ways to compute transport properties from molecular |
3338 |
+ |
dynamics simulations. Equilibrium Molecular Dynamics (EMD) |
3339 |
+ |
simulations can be used by computing relevant time correlation |
3340 |
+ |
functions and assuming linear response theory holds. For some transport properties (notably thermal conductivity), EMD approaches |
3341 |
+ |
are subject to noise and poor convergence of the relevant |
3342 |
+ |
correlation functions. Traditional Non-equilibrium Molecular Dynamics |
3343 |
+ |
(NEMD) methods impose a gradient (e.g. thermal or momentum) on a |
3344 |
+ |
simulation. However, the resulting flux is often difficult to |
3345 |
+ |
measure. Furthermore, problems arise for NEMD simulations of |
3346 |
+ |
heterogeneous systems, such as phase-phase boundaries or interfaces, |
3347 |
+ |
where the type of gradient to enforce at the boundary between |
3348 |
+ |
materials is unclear. |
3349 |
+ |
|
3350 |
+ |
{\it Reverse} Non-Equilibrium Molecular Dynamics (RNEMD) methods adopt |
3351 |
+ |
a different approach in that an unphysical {\it flux} is imposed |
3352 |
+ |
between different regions or ``slabs'' of the simulation box. The |
3353 |
+ |
response of the system is to develop a temperature or momentum {\it |
3354 |
+ |
gradient} between the two regions. Since the amount of the applied |
3355 |
+ |
flux is known exactly, and the measurement of gradient is generally |
3356 |
+ |
less complicated, imposed-flux methods typically take shorter |
3357 |
+ |
simulation times to obtain converged results for transport properties. |
3358 |
+ |
|
3359 |
+ |
\begin{figure} |
3360 |
+ |
\includegraphics[width=\linewidth]{rnemdDemo} |
3361 |
+ |
\caption{The (VSS) RNEMD approach imposes unphysical transfer of both |
3362 |
+ |
linear momentum and kinetic energy between a ``hot'' slab and a |
3363 |
+ |
``cold'' slab in the simulation box. The system responds to this |
3364 |
+ |
imposed flux by generating both momentum and temperature gradients. |
3365 |
+ |
The slope of the gradients can then be used to compute transport |
3366 |
+ |
properties (e.g. shear viscosity and thermal conductivity).} |
3367 |
+ |
\label{rnemdDemo} |
3368 |
+ |
\end{figure} |
3369 |
+ |
|
3370 |
+ |
\section{\label{section:algo}Three algorithms for carrying out RNEMD simulations} |
3371 |
+ |
\subsection{\label{subsection:swapping}The swapping algorithm} |
3372 |
+ |
The original ``swapping'' approaches by M\"{u}ller-Plathe {\it et |
3373 |
+ |
al.}\cite{ISI:000080382700030,MullerPlathe:1997xw} can be understood |
3374 |
+ |
as a sequence of imaginary elastic collisions between particles in |
3375 |
+ |
opposite slabs. In each collision, the entire momentum vectors of |
3376 |
+ |
both particles may be exchanged to generate a thermal |
3377 |
+ |
flux. Alternatively, a single component of the momentum vectors may be |
3378 |
+ |
exchanged to generate a shear flux. This algorithm turns out to be |
3379 |
+ |
quite useful in many simulations. However, the M\"{u}ller-Plathe |
3380 |
+ |
swapping approach perturbs the system away from ideal |
3381 |
+ |
Maxwell-Boltzmann distributions, and this may leads to undesirable |
3382 |
+ |
side-effects when the applied flux becomes large.\cite{Maginn:2010} |
3383 |
+ |
This limits the applicability of the swapping algorithm, so in OpenMD, |
3384 |
+ |
we have implemented two additional algorithms for RNEMD in addition to the |
3385 |
+ |
original swapping approach. |
3386 |
+ |
|
3387 |
+ |
\subsection{\label{subsection:nivs}Non-Isotropic Velocity Scaling (NIVS)} |
3388 |
+ |
Instead of having momentum exchange between {\it individual particles} |
3389 |
+ |
in each slab, the NIVS algorithm applies velocity scaling to all of |
3390 |
+ |
the selected particles in both slabs.\cite{kuang:164101} A combination of linear |
3391 |
+ |
momentum, kinetic energy, and flux constraint equations governs the |
3392 |
+ |
amount of velocity scaling performed at each step. Interested readers |
3393 |
+ |
should consult ref. \citealp{kuang:164101} for further details on the |
3394 |
+ |
methodology. |
3395 |
+ |
|
3396 |
+ |
NIVS has been shown to be very effective at producing thermal |
3397 |
+ |
gradients and for computing thermal conductivities, particularly for |
3398 |
+ |
heterogeneous interfaces. Although the NIVS algorithm can also be |
3399 |
+ |
applied to impose a directional momentum flux, thermal anisotropy was |
3400 |
+ |
observed in relatively high flux simulations, and the method is not |
3401 |
+ |
suitable for imposing a shear flux or for computing shear viscosities. |
3402 |
+ |
|
3403 |
+ |
\subsection{\label{subsection:vss}Velocity Shearing and Scaling (VSS)} |
3404 |
+ |
The third RNEMD algorithm implemented in OpenMD utilizes a series of |
3405 |
+ |
simultaneous velocity shearing and scaling exchanges between the two |
3406 |
+ |
slabs.\cite{2012MolPh.110..691K} This method results in a set of simpler equations to satisfy |
3407 |
+ |
the conservation constraints while creating a desired flux between the |
3408 |
+ |
two slabs. |
3409 |
+ |
|
3410 |
+ |
The VSS approach is versatile in that it may be used to implement both |
3411 |
+ |
thermal and shear transport either separately or simultaneously. |
3412 |
+ |
Perturbations of velocities away from the ideal Maxwell-Boltzmann |
3413 |
+ |
distributions are minimal, and thermal anisotropy is kept to a |
3414 |
+ |
minimum. This ability to generate simultaneous thermal and shear |
3415 |
+ |
fluxes has been utilized to map out the shear viscosity of SPC/E water |
3416 |
+ |
over a wide range of temperatures (90~K) just with a single simulation. |
3417 |
+ |
VSS-RNEMD also allows the directional momentum flux to have |
3418 |
+ |
arbitary directions, which could aid in the study of anisotropic solid |
3419 |
+ |
surfaces in contact with liquid environments. |
3420 |
+ |
|
3421 |
+ |
\section{\label{section:usingRNEMD}Using OpenMD to perform a RNEMD simulation} |
3422 |
+ |
\subsection{\label{section:rnemdParams} What the user needs to specify} |
3423 |
+ |
To carry out a RNEMD simulation, |
3424 |
+ |
a user must specify a number of parameters in the MetaData (.md) file. |
3425 |
+ |
Because the RNEMD methods have a large number of parameters, these |
3426 |
+ |
must be enclosed in a {\it separate} {\tt RNEMD\{...\}} block. The most important |
3427 |
+ |
parameters to specify are the {\tt useRNEMD}, {\tt fluxType} and flux |
3428 |
+ |
parameters. Most other parameters (summarized in table |
3429 |
+ |
\ref{table:rnemd}) have reasonable default values. {\tt fluxType} |
3430 |
+ |
sets up the kind of exchange that will be carried out between the two |
3431 |
+ |
slabs (either Kinetic Energy ({\tt KE}) or momentum ({\tt Px, Py, Pz, |
3432 |
+ |
Pvector}), or some combination of these). The flux is specified |
3433 |
+ |
with the use of three possible parameters: {\tt kineticFlux} for |
3434 |
+ |
kinetic energy exchange, as well as {\tt momentumFlux} or {\tt |
3435 |
+ |
momentumFluxVector} for simulations with directional exchange. |
3436 |
+ |
|
3437 |
+ |
\subsection{\label{section:rnemdResults} Processing the results} |
3438 |
+ |
OpenMD will generate a {\tt .rnemd} |
3439 |
+ |
file with the same prefix as the original {\tt .md} file. This file |
3440 |
+ |
contains a running average of properties of interest computed within a |
3441 |
+ |
set of bins that divide the simulation cell along the $z$-axis. The |
3442 |
+ |
first column of the {\tt .rnemd} file is the $z$ coordinate of the |
3443 |
+ |
center of each bin, while following columns may contain the average |
3444 |
+ |
temperature, velocity, or density within each bin. The output format |
3445 |
+ |
in the {\tt .rnemd} file can be altered with the {\tt outputFields}, |
3446 |
+ |
{\tt outputBins}, and {\tt outputFileName} parameters. A report at |
3447 |
+ |
the top of the {\tt .rnemd} file contains the current exchange totals |
3448 |
+ |
as well as the average flux applied during the simulation. Using the |
3449 |
+ |
slope of the temperature or velocity gradient obtaine from the {\tt |
3450 |
+ |
.rnemd} file along with the applied flux, the user can very simply |
3451 |
+ |
arrive at estimates of thermal conductivities ($\lambda$), |
3452 |
+ |
\begin{equation} |
3453 |
+ |
J_z = -\lambda \frac{\partial T}{\partial z}, |
3454 |
+ |
\end{equation} |
3455 |
+ |
and shear viscosities ($\eta$), |
3456 |
+ |
\begin{equation} |
3457 |
+ |
j_z(p_x) = -\eta \frac{\partial \langle v_x \rangle}{\partial z}. |
3458 |
+ |
\end{equation} |
3459 |
+ |
Here, the quantities on the left hand side are the actual flux values |
3460 |
+ |
(in the header of the {\tt .rnemd} file), while the slopes are |
3461 |
+ |
obtained from linear fits to the gradients observed in the {\tt |
3462 |
+ |
.rnemd} file. |
3463 |
+ |
|
3464 |
+ |
More complicated simulations (including interfaces) require a bit more |
3465 |
+ |
care. Here the second derivative may be required to compute the |
3466 |
+ |
interfacial thermal conductance, |
3467 |
+ |
\begin{align} |
3468 |
+ |
G^\prime &= \left(\nabla\lambda \cdot \mathbf{\hat{n}}\right)_{z_0} \\ |
3469 |
+ |
&= \frac{\partial}{\partial z}\left(-\frac{J_z}{ |
3470 |
+ |
\left(\frac{\partial T}{\partial z}\right)}\right)_{z_0} \\ |
3471 |
+ |
&= J_z\left(\frac{\partial^2 T}{\partial z^2}\right)_{z_0} \Big/ |
3472 |
+ |
\left(\frac{\partial T}{\partial z}\right)_{z_0}^2. |
3473 |
+ |
\label{derivativeG} |
3474 |
+ |
\end{align} |
3475 |
+ |
where $z_0$ is the location of the interface between two materials and |
3476 |
+ |
$\mathbf{\hat{n}}$ is a unit vector normal to the interface. We |
3477 |
+ |
suggest that users interested in interfacial conductance consult |
3478 |
+ |
reference \citealp{kuang:AuThl} for other approaches to computing $G$. |
3479 |
+ |
Users interested in {\it friction coefficients} at heterogeneous |
3480 |
+ |
interfaces may also find reference \citealp{2012MolPh.110..691K} |
3481 |
+ |
useful. |
3482 |
+ |
|
3483 |
+ |
\newpage |
3484 |
+ |
|
3485 |
+ |
\begin{longtable}[c]{JKLM} |
3486 |
+ |
\caption{Meta-data Keywords: Parameters for RNEMD simulations}\\ |
3487 |
+ |
\multicolumn{4}{c}{The following keywords must be enclosed inside a {\tt RNEMD\{...\}} block.} |
3488 |
+ |
\\ \hline |
3489 |
+ |
{\bf keyword} & {\bf units} & {\bf use} & {\bf remarks} \\ \hline |
3490 |
+ |
\endhead |
3491 |
+ |
\hline |
3492 |
+ |
\endfoot |
3493 |
+ |
{\tt useRNEMD} & logical & perform RNEMD? & default is ``false'' \\ |
3494 |
+ |
{\tt objectSelection} & string & see section \ref{section:syntax} |
3495 |
+ |
for selection syntax & default is ``select all'' \\ |
3496 |
+ |
{\tt method} & string & exchange method & one of the following: |
3497 |
+ |
{\tt Swap, NIVS,} or {\tt VSS} (default is {\tt VSS}) \\ |
3498 |
+ |
{\tt fluxType} & string & what is being exchanged between slabs? & one |
3499 |
+ |
of the following: {\tt KE, Px, Py, Pz, Pvector, KE+Px, KE+Py, KE+Pvector} \\ |
3500 |
+ |
{\tt kineticFlux} & kcal mol$^{-1}$ \AA$^{-2}$ fs$^{-1}$ & specify the kinetic energy flux & \\ |
3501 |
+ |
{\tt momentumFlux} & amu \AA$^{-1}$ fs$^{-2}$ & specify the momentum flux & \\ |
3502 |
+ |
{\tt momentumFluxVector} & amu \AA$^{-1}$ fs$^{-2}$ & specify the momentum flux when |
3503 |
+ |
{\tt Pvector} is part of the exchange & Vector3d input\\ |
3504 |
+ |
{\tt exchangeTime} & fs & how often to perform the exchange & default is 100 fs\\ |
3505 |
+ |
|
3506 |
+ |
{\tt slabWidth} & $\mbox{\AA}$ & width of the two exchange slabs & default is $\mathsf{H}_{zz} / 10.0$ \\ |
3507 |
+ |
{\tt slabAcenter} & $\mbox{\AA}$ & center of the end slab & default is 0 \\ |
3508 |
+ |
{\tt slabBcenter} & $\mbox{\AA}$ & center of the middle slab & default is $\mathsf{H}_{zz} / 2$ \\ |
3509 |
+ |
{\tt outputFileName} & string & file name for output histograms & default is the same prefix as the |
3510 |
+ |
.md file, but with the {\tt .rnemd} extension \\ |
3511 |
+ |
{\tt outputBins} & int & number of $z$-bins in the output histogram & |
3512 |
+ |
default is 20 \\ |
3513 |
+ |
{\tt outputFields} & string & columns to print in the {\tt .rnemd} |
3514 |
+ |
file where each column is separated by a pipe ($\mid$) symbol. & Allowed column names are: {\sc z, temperature, velocity, density} \\ |
3515 |
+ |
\label{table:rnemd} |
3516 |
+ |
\end{longtable} |
3517 |
+ |
|
3518 |
|
\chapter{\label{section:minimizer}Energy Minimization} |
3519 |
|
|
3520 |
< |
As one of the basic procedures of molecular modeling, energy |
3083 |
< |
minimization is used to identify local configurations that are stable |
3520 |
> |
Energy minimization is used to identify local configurations that are stable |
3521 |
|
points on the potential energy surface. There is a vast literature on |
3522 |
|
energy minimization algorithms have been developed to search for the |
3523 |
|
global energy minimum as well as to find local structures which are |
3644 |
|
\begin{figure} |
3645 |
|
\centering |
3646 |
|
\includegraphics[width=3in]{heirarchy.pdf} |
3647 |
< |
\caption[Class heirarchy for StuntDoubles in {\sc OpenMD}-4]{ \\ The |
3648 |
< |
class heirarchy of StuntDoubles in {\sc OpenMD}-4. The selection |
3647 |
> |
\caption[Class heirarchy for StuntDoubles in {\sc OpenMD}]{ \\ The |
3648 |
> |
class heirarchy of StuntDoubles in {\sc OpenMD}. The selection |
3649 |
|
syntax allows the user to select any of the objects that are descended |
3650 |
|
from a StuntDouble.} |
3651 |
|
\label{fig:heirarchy} |
3825 |
|
-z & {\tt -{}-zconstraint} & replace the atom types of zconstraint molecules (default=off) \\ |
3826 |
|
-r & {\tt -{}-rigidbody} & add a pseudo COM atom to rigidbody (default=off) \\ |
3827 |
|
-t & {\tt -{}-watertype} & replace the atom type of water model (default=on) \\ |
3828 |
< |
-b & {\tt -{}-basetype} & using base atom type (default=off) \\ |
3828 |
> |
-b & {\tt -{}-basetype} & using base atom type |
3829 |
> |
(default=off) \\ |
3830 |
> |
-v& {\tt -{}-velocities} & Print velocities in xyz file (default=off)\\ |
3831 |
> |
-f& {\tt -{}-forces} & Print forces xyz file (default=off)\\ |
3832 |
> |
-u& {\tt -{}-vectors} & Print vectors (dipoles, etc) in xyz file |
3833 |
> |
(default=off)\\ |
3834 |
> |
-c& {\tt -{}-charges} & Print charges in xyz file (default=off)\\ |
3835 |
> |
-e& {\tt -{}-efield} & Print electric field vector in xyz file |
3836 |
> |
(default=off)\\ |
3837 |
|
& {\tt -{}-repeatX=INT} & The number of images to repeat in the x direction (default=`0') \\ |
3838 |
|
& {\tt -{}-repeatY=INT} & The number of images to repeat in the y direction (default=`0') \\ |
3839 |
|
& {\tt -{}-repeatZ=INT} & The number of images to repeat in the z direction (default=`0') \\ |
3925 |
|
& {\tt -{}-sele1=selection script} & select the first StuntDouble set \\ |
3926 |
|
& {\tt -{}-sele2=selection script} & select the second StuntDouble set \\ |
3927 |
|
& {\tt -{}-sele3=selection script} & select the third StuntDouble set \\ |
3928 |
< |
& {\tt -{}-refsele=selection script} & select reference (can only be used with {\tt -{}-gxyz}) \\ |
3928 |
> |
& {\tt -{}-refsele=selection script} & select reference (can only |
3929 |
> |
be used with {\tt -{}-gxyz}) \\ |
3930 |
> |
& {\tt -{}-comsele=selection script} |
3931 |
> |
& select stunt doubles for center-of-mass |
3932 |
> |
reference point\\ |
3933 |
> |
& {\tt -{}-seleoffset=INT} & global index offset for a second object (used |
3934 |
> |
to define a vector between sites in molecule)\\ |
3935 |
> |
|
3936 |
|
& {\tt -{}-molname=STRING} & molecule name \\ |
3937 |
|
& {\tt -{}-begin=INT} & begin internal index \\ |
3938 |
|
& {\tt -{}-end=INT} & end internal index \\ |
3939 |
+ |
& {\tt -{}-radius=DOUBLE} & nanoparticle radius\\ |
3940 |
|
\hline |
3941 |
|
\multicolumn{3}{|l|}{One option from the following group of options is required:} \\ |
3942 |
|
\hline |
3943 |
< |
& {\tt -{}-gofr} & $g(r)$ \\ |
3944 |
< |
& {\tt -{}-r\_theta} & $g(r, \cos(\theta))$ \\ |
3945 |
< |
& {\tt -{}-r\_omega} & $g(r, \cos(\omega))$ \\ |
3946 |
< |
& {\tt -{}-theta\_omega} & $g(\cos(\theta), \cos(\omega))$ \\ |
3943 |
> |
& {\tt -{}-bo} & bond order parameter ({\tt -{}-rcut} must be specified) \\ |
3944 |
> |
& {\tt -{}-bor} & bond order parameter as a function of |
3945 |
> |
radius ({\tt -{}-rcut} must be specified) \\ |
3946 |
> |
& {\tt -{}-bad} & $N(\theta)$ bond angle density within ({\tt -{}-rcut} must be specified) \\ |
3947 |
> |
& {\tt -{}-count} & count of molecules matching selection |
3948 |
> |
criteria (and associated statistics) \\ |
3949 |
> |
-g& {\tt -{}-gofr} & $g(r)$ \\ |
3950 |
> |
& {\tt -{}-gofz} & $g(z)$ \\ |
3951 |
> |
& {\tt -{}-r\_theta} & $g(r, \cos(\theta))$ \\ |
3952 |
> |
& {\tt -{}-r\_omega} & $g(r, \cos(\omega))$ \\ |
3953 |
> |
& {\tt -{}-r\_z} & $g(r, z)$ \\ |
3954 |
> |
& {\tt -{}-theta\_omega} & $g(\cos(\theta), \cos(\omega))$ \\ |
3955 |
|
& {\tt -{}-gxyz} & $g(x, y, z)$ \\ |
3956 |
< |
& {\tt -{}-p2} & $P_2$ order parameter ({\tt -{}-sele1} and {\tt -{}-sele2} must be specified) \\ |
3956 |
> |
& {\tt -{}-twodgofr} & 2D $g(r)$ (Slab width {\tt -{}-dz} must be specified)\\ |
3957 |
> |
-p& {\tt -{}-p2} & $P_2$ order parameter ({\tt -{}-sele1} must be specified, {\tt -{}-sele2} is optional) \\ |
3958 |
> |
& {\tt -{}-rp2} & Ripple order parameter ({\tt -{}-sele1} and {\tt -{}-sele2} must be specified) \\ |
3959 |
|
& {\tt -{}-scd} & $S_{CD}$ order parameter(either {\tt -{}-sele1}, {\tt -{}-sele2}, {\tt -{}-sele3} are specified or {\tt -{}-molname}, {\tt -{}-begin}, {\tt -{}-end} are specified) \\ |
3960 |
< |
& {\tt -{}-density} & density plot ({\tt -{}-sele1} must be specified) \\ |
3961 |
< |
& {\tt -{}-slab\_density} & slab density ({\tt -{}-sele1} must be specified) |
3960 |
> |
-d& {\tt -{}-density} & density plot \\ |
3961 |
> |
& {\tt -{}-slab\_density} & slab density \\ |
3962 |
> |
& {\tt -{}-p\_angle} & $p(\cos(\theta))$ ($\theta$ |
3963 |
> |
is the angle between molecular axis and radial vector from origin\\ |
3964 |
> |
& {\tt -{}-hxy} & Calculates the undulation spectrum, $h(x,y)$, of an interface \\ |
3965 |
> |
& {\tt -{}-rho\_r} & $\rho(r)$\\ |
3966 |
> |
& {\tt -{}-angle\_r} & $\theta(r)$ (spatially resolves the |
3967 |
> |
angle between the molecular axis and the radial vector from the |
3968 |
> |
origin)\\ |
3969 |
> |
& {\tt -{}-hullvol} & hull volume of nanoparticle\\ |
3970 |
> |
& {\tt -{}-rodlength} & length of nanorod\\ |
3971 |
> |
-Q& {\tt -{}-tet\_param} & tetrahedrality order parameter ($Q$)\\ |
3972 |
> |
& {\tt -{}-tet\_param\_z} & spatially-resolved tetrahedrality order |
3973 |
> |
parameter $Q(z)$\\ |
3974 |
> |
& {\tt -{}-rnemdz} & slab-resolved RNEMD statistics (temperature, |
3975 |
> |
density, velocity)\\ |
3976 |
> |
& {\tt -{}-rnemdr} & shell-resolved RNEMD statistics (temperature, |
3977 |
> |
density, angular velocity) |
3978 |
|
\end{longtable} |
3979 |
|
|
3980 |
|
\subsection{\label{section:DynamicProps}DynamicProps} |
4015 |
|
-o& {\tt -{}-output=filename} & output file name \\ |
4016 |
|
& {\tt -{}-sele1=selection script} & select first StuntDouble set \\ |
4017 |
|
& {\tt -{}-sele2=selection script} & select second StuntDouble set (if sele2 is not set, use script from sele1) \\ |
4018 |
+ |
& {\tt -{}-order=INT} & Lengendre Polynomial Order\\ |
4019 |
+ |
-z& {\tt -{}-nzbins=INT} & Number of $z$ bins (default=`100`)\\ |
4020 |
+ |
-m& {\tt -{}-memory=memory specification} |
4021 |
+ |
&Available memory |
4022 |
+ |
(default=`2G`)\\ |
4023 |
|
\hline |
4024 |
|
\multicolumn{3}{|l|}{One option from the following group of options is required:} \\ |
4025 |
|
\hline |
4026 |
< |
-r& {\tt -{}-rcorr} & compute mean square displacement \\ |
4027 |
< |
-v& {\tt -{}-vcorr} & compute velocity correlation function \\ |
4028 |
< |
-d& {\tt -{}-dcorr} & compute dipole correlation function |
4026 |
> |
-s& {\tt -{}-selecorr} & selection correlation function \\ |
4027 |
> |
-r& {\tt -{}-rcorr} & compute mean squared displacement \\ |
4028 |
> |
-v& {\tt -{}-vcorr} & velocity autocorrelation function \\ |
4029 |
> |
-d& {\tt -{}-dcorr} & dipole correlation function \\ |
4030 |
> |
-l& {\tt -{}-lcorr} & Lengendre correlation function \\ |
4031 |
> |
& {\tt -{}-lcorrZ} & Lengendre correlation function binned by $z$ \\ |
4032 |
> |
& {\tt -{}-cohZ} & Lengendre correlation function for OH bond vectors binned by $z$\\ |
4033 |
> |
-M& {\tt -{}-sdcorr} & System dipole correlation function\\ |
4034 |
> |
& {\tt -{}-r\_rcorr} & Radial mean squared displacement\\ |
4035 |
> |
& {\tt -{}-thetacorr} & Angular mean squared displacement\\ |
4036 |
> |
& {\tt -{}-drcorr} & Directional mean squared displacement for particles with unit vectors\\ |
4037 |
> |
& {\tt -{}-helfandEcorr} & Helfand moment for thermal conductvity\\ |
4038 |
> |
-p& {\tt -{}-momentum} & Helfand momentum for viscosity\\ |
4039 |
> |
& {\tt -{}-stresscorr} & Stress tensor correlation function |
4040 |
|
\end{longtable} |
4041 |
|
|
4042 |
|
\chapter{\label{section:PreparingInput} Preparing Input Configurations} |
4103 |
|
to {\tt atom2md}, but they use a specific input format and do not |
4104 |
|
expect the the input specifier on the command line. |
4105 |
|
|
4106 |
+ |
|
4107 |
|
\section{\label{section:SimpleBuilder}SimpleBuilder} |
4108 |
|
|
4109 |
|
{\tt SimpleBuilder} creates simple lattice structures. It requires an |
4128 |
|
& {\tt -{}-nx=INT} & number of unit cells in x\\ |
4129 |
|
& {\tt -{}-ny=INT} & number of unit cells in y\\ |
4130 |
|
& {\tt -{}-nz=INT} & number of unit cells in z |
4131 |
+ |
\end{longtable} |
4132 |
+ |
|
4133 |
+ |
\section{\label{section:icosahedralBuilder}icosahedralBuilder} |
4134 |
+ |
|
4135 |
+ |
{\tt icosahedralBuilder} creates single-component geometric solids |
4136 |
+ |
that can be useful in simulating nanostructures. Like the other |
4137 |
+ |
builders, it requires an initial, but skeletal {\sc OpenMD} file to |
4138 |
+ |
specify the component that is to be placed on the lattice. The total |
4139 |
+ |
number of placed molecules will be shown at the top of the |
4140 |
+ |
configuration file that is generated, and that number may not match |
4141 |
+ |
the original meta-data file, so a new meta-data file is also generated |
4142 |
+ |
which matches the lattice structure. |
4143 |
+ |
|
4144 |
+ |
The options available for icosahedralBuilder are as follows: |
4145 |
+ |
\begin{longtable}[c]{|EFG|} |
4146 |
+ |
\caption{icosahedralBuilder Command-line Options} |
4147 |
+ |
\\ \hline |
4148 |
+ |
{\bf option} & {\bf verbose option} & {\bf behavior} \\ \hline |
4149 |
+ |
\endhead |
4150 |
+ |
\hline |
4151 |
+ |
\endfoot |
4152 |
+ |
-h& {\tt -{}-help} & Print help and exit\\ |
4153 |
+ |
-V& {\tt -{}-version} & Print version and exit\\ |
4154 |
+ |
-o& {\tt -{}-output=STRING} & Output file name\\ |
4155 |
+ |
-n& {\tt -{}-shells=INT} & Nanoparticle shells\\ |
4156 |
+ |
-d& {\tt -{}-latticeConstant=DOUBLE} & Lattice spacing in Angstroms for cubic lattice.\\ |
4157 |
+ |
-c& {\tt -{}-columnAtoms=INT} & Number of atoms along central |
4158 |
+ |
column (Decahedron only)\\ |
4159 |
+ |
-t& {\tt -{}-twinAtoms=INT} & Number of atoms along twin |
4160 |
+ |
boundary (Decahedron only) \\ |
4161 |
+ |
-p& {\tt -{}-truncatedPlanes=INT} & Number of truncated planes (Curling-stone Decahedron only)\\ |
4162 |
+ |
\hline |
4163 |
+ |
\multicolumn{3}{|l|}{One option from the following group of options is required:} \\ |
4164 |
+ |
\hline |
4165 |
+ |
& {\tt -{}-ico} & Create an Icosahedral cluster \\ |
4166 |
+ |
& {\tt -{}-deca} & Create a regualar Decahedral cluster\\ |
4167 |
+ |
& {\tt -{}-ino} & Create an Ino Decahedral cluster\\ |
4168 |
+ |
& {\tt -{}-marks} & Create a Marks Decahedral cluster\\ |
4169 |
+ |
& {\tt -{}-stone} & Create a Curling-stone Decahedral cluster |
4170 |
|
\end{longtable} |
4171 |
|
|
4172 |
+ |
|
4173 |
|
\section{\label{section:Hydro}Hydro} |
4174 |
|
{\tt Hydro} generates resistance tensor ({\tt .diff}) files which are |
4175 |
|
required when using the Langevin integrator using complex rigid |
4203 |
|
\end{longtable} |
4204 |
|
|
4205 |
|
|
4206 |
+ |
|
4207 |
+ |
|
4208 |
+ |
|
4209 |
|
\chapter{\label{section:parallelization} Parallel Simulation Implementation} |
4210 |
|
|
4211 |
|
Although processor power is continually improving, it is still |
4289 |
|
DMR-0079647. |
4290 |
|
|
4291 |
|
|
4292 |
< |
\bibliographystyle{jcc} |
4292 |
> |
\bibliographystyle{aip} |
4293 |
|
\bibliography{openmdDoc} |
4294 |
|
|
4295 |
|
\end{document} |