| 779 |
|
allowing the user to gauge the stability of the integrator. The |
| 780 |
|
statistics file is denoted with the \texttt{.stat} file extension. |
| 781 |
|
|
| 782 |
< |
\chapter{\label{section:forceFields}Force Fields} |
| 782 |
> |
\chapter{\label{chapter:forceFields}Force Fields} |
| 783 |
|
|
| 784 |
|
Like many molecular simulation packages, {\sc OpenMD} splits the |
| 785 |
|
potential energy into the short-ranged (bonded) portion and a |
| 793 |
|
are defined by the force field which is selected in the MetaData |
| 794 |
|
section. |
| 795 |
|
|
| 796 |
< |
\section{\label{section:shortRange}The basic interactions} |
| 796 |
> |
\section{\label{section:divisionOfLabor}Separation into Internal and |
| 797 |
> |
Cross interactions} |
| 798 |
|
|
| 799 |
< |
The energy function for a system composed of $N$ molecules is |
| 800 |
< |
traditionally written |
| 799 |
> |
The classical potential energy function for a system composed of $N$ |
| 800 |
> |
molecules is traditionally written |
| 801 |
|
\begin{equation} |
| 802 |
|
V = \sum^{N}_{I=1} V^{I}_{\text{Internal}} |
| 803 |
|
+ \sum^{N-1}_{I=1} \sum_{J>I} V^{IJ}_{\text{Cross}}, |
| 804 |
|
\label{eq:totalPotential} |
| 805 |
|
\end{equation} |
| 806 |
< |
where $V^{IJ}_{\text{Cross}}$ contains all intermolecular interactions |
| 807 |
< |
between molecules $I$ and $J$, and $V^{I}_{\text{Internal}}$ is the |
| 808 |
< |
internal potential of molecule $I$: |
| 809 |
< |
\begin{align*} |
| 810 |
< |
V^{I}_{\text{Internal}} = & |
| 811 |
< |
\sum_{r_{ij} \in I} V_{\text{bond}}(r_{ij}) |
| 811 |
< |
+ \sum_{\theta_{ijk} \in I} V_{\text{bend}}(\theta_{ijk}) |
| 812 |
< |
+ \sum_{\phi_{ijkl} \in I} V_{\text{torsion}}(\phi_{ijkl}) |
| 813 |
< |
+ \sum_{\omega_{ijkl} \in I} V_{\text{inversion}}(\omega_{ijkl}) \\ |
| 814 |
< |
& + \sum_{i \in I} \sum_{(j>i+4) \in I} |
| 815 |
< |
\biggl[ V_{\text{dispersion}}(r_{ij}) + V_{\text{electrostatic}} |
| 816 |
< |
(\mathbf{r}_{ij},\boldsymbol{\Omega}_{i},\boldsymbol{\Omega}_{j}) |
| 817 |
< |
\biggr]. |
| 818 |
< |
\label{eq:internalPotential} |
| 819 |
< |
\end{align*} |
| 820 |
< |
Here $V_{\text{bond}}, V_{\text{bend}}, |
| 821 |
< |
V_{\text{torsion}},\mathrm{~and~} V_{\text{inversion}}$ represent the |
| 822 |
< |
bond, bend, torsion, and inversion potentials for all |
| 823 |
< |
topologically-connected sets of atoms within the molecule. Bonds are |
| 824 |
< |
the primary way of specifying how the atoms are connected together to |
| 825 |
< |
form the molecule (i.e. they define the molecular topology). The |
| 826 |
< |
other short-range interactions may be specified explicitly in the |
| 827 |
< |
molecule definition, or they may be deduced from bonding information. |
| 828 |
< |
For example, bends can be implicitly deduced from two bonds which |
| 829 |
< |
share an atom. Torsions can be deduced from two bends that share a |
| 830 |
< |
bond. Inversion potentials are utilized primarily to enforce |
| 831 |
< |
planarity around $sp^2$-hybridized sites, and these are specified with |
| 832 |
< |
central atoms and satellites (or an atom with bonds to exactly three |
| 833 |
< |
satellites). The pairwise portions of the non-bonded interactions are |
| 834 |
< |
usually excluded for atom pairs that are involved in the same bond, |
| 835 |
< |
bend, or torsion. All other atom pairs within a molecule are subject |
| 836 |
< |
to non-bonded pair potentials. |
| 806 |
> |
where $V^{I}_{\text{Internal}}$ contains all of the terms internal to |
| 807 |
> |
molecule $I$ (e.g. bonding, bending, torsional, and inversion terms) |
| 808 |
> |
and $V^{IJ}_{\text{Cross}}$ contains all intermolecular interactions |
| 809 |
> |
between molecules $I$ and $J$. For large molecules, the internal |
| 810 |
> |
potential may also include some non-bonded terms like electrostatic or |
| 811 |
> |
van der Waals interactions. |
| 812 |
|
|
| 813 |
|
The types of atoms being simulated, as well as the specific functional |
| 814 |
|
forms and parameters of the intra-molecular functions and the |
| 875 |
|
end TorsionTypes |
| 876 |
|
\end{lstlisting} |
| 877 |
|
|
| 878 |
< |
\subsection{\label{section:ffOptions}The Options block} |
| 878 |
> |
\section{\label{section:ffOptions}The Options block} |
| 879 |
|
|
| 880 |
|
The Options block defines properties governing how the force field |
| 881 |
|
interactions are carried out. This block is delineated with the text |
| 893 |
|
Name = "alkane" // any string |
| 894 |
|
vdWtype = "Lennard-Jones" |
| 895 |
|
DistanceMixingRule = "arithmetic" // can also be "geometric" or "cubic" |
| 896 |
< |
DistanceType = "sigma" // can also be Rmin |
| 896 |
> |
DistanceType = "sigma" // can also be "Rmin" |
| 897 |
|
EnergyMixingRule = "geometric" // can also be "arithmetic" or "hhg" |
| 898 |
|
EnergyUnitScaling = 1.0 |
| 899 |
|
MetallicEnergyUnitScaling = 1.0 |
| 912 |
|
end Options |
| 913 |
|
\end{lstlisting} |
| 914 |
|
|
| 915 |
< |
\subsection{\label{section:ffBase}The BaseAtomTypes block} |
| 915 |
> |
\section{\label{section:ffBase}The BaseAtomTypes block} |
| 916 |
|
|
| 917 |
|
An AtomType the primary data structure that OpenMD uses to store |
| 918 |
|
static data about an atom. Things that belong to AtomType are |
| 960 |
|
end BaseAtomTypes |
| 961 |
|
\end{lstlisting} |
| 962 |
|
|
| 963 |
< |
\subsection{\label{section:ffAtom}The AtomTypes block} |
| 963 |
> |
\section{\label{section:ffAtom}The AtomTypes block} |
| 964 |
|
|
| 965 |
|
AtomTypes inherit most properties from BaseAtomTypes, but can override |
| 966 |
|
their lower-level properties as well. Scheme \ref{sch:atomTypesBlock} |
| 994 |
|
end AtomTypes |
| 995 |
|
\end{lstlisting} |
| 996 |
|
|
| 997 |
< |
\subsection{\label{section:ffDirectionalAtom}The DirectionalAtomTypes |
| 997 |
> |
\section{\label{section:ffDirectionalAtom}The DirectionalAtomTypes |
| 998 |
|
block} |
| 999 |
|
DirectionalAtoms have orientational degrees of freedom as well as |
| 1000 |
|
translation, so moving these atoms requires information about the |
| 1029 |
|
of freedom (three translations and two rotations), and will alter |
| 1030 |
|
calculation of temperatures to reflect this. |
| 1031 |
|
|
| 1032 |
< |
\subsection{\label{section::ffAtomProperties}AtomType properties} |
| 1033 |
< |
\subsubsection{\label{section:ffLJ}The LennardJonesAtomTypes block} |
| 1032 |
> |
\section{\label{section::ffAtomProperties}AtomType properties} |
| 1033 |
> |
\subsection{\label{section:ffLJ}The LennardJonesAtomTypes block} |
| 1034 |
|
One of the most basic interatomic interactions implemented in {\sc |
| 1035 |
|
OpenMD} is the Lennard-Jones potential, which mimics the van der |
| 1036 |
|
Waals interaction at long distances and uses an empirical repulsion at |
| 1088 |
|
end LennardJonesAtomTypes |
| 1089 |
|
\end{lstlisting} |
| 1090 |
|
|
| 1091 |
< |
\subsubsection{\label{section:ffCharge}The ChargeAtomTypes block} |
| 1091 |
> |
\subsection{\label{section:ffCharge}The ChargeAtomTypes block} |
| 1092 |
|
|
| 1093 |
|
In molecular simulations, proper accumulation of the electrostatic |
| 1094 |
|
interactions is essential and is one of the most |
| 1128 |
|
end ChargeAtomTypes |
| 1129 |
|
\end{lstlisting} |
| 1130 |
|
|
| 1131 |
< |
\subsubsection{\label{section:ffMultipole}The MultipoleAtomTypes |
| 1131 |
> |
\subsection{\label{section:ffMultipole}The MultipoleAtomTypes |
| 1132 |
|
block} |
| 1133 |
|
For complex charge distributions that are centered on single sites, it |
| 1134 |
|
is convenient to write the total electrostatic potential in terms of |
| 1144 |
|
\end{equation} |
| 1145 |
|
Here, the point charge, dipole, and quadrupole for site $\bf a$ are |
| 1146 |
|
given by $C_{\bf a}$, $D_{{\bf a}\alpha}$, and $Q_{{\bf |
| 1147 |
< |
a}\alpha\beta}$, respectively. These are the primitive |
| 1147 |
> |
a}\alpha\beta}$, respectively. These are the {\it primitive} |
| 1148 |
|
multipoles. If the site is representing a distribution of charges, |
| 1149 |
|
these can be expressed as, |
| 1150 |
|
\begin{align} |
| 1209 |
|
in units of degrees. Dipole moments are entered in units of Debye, |
| 1210 |
|
and Quadrupole moments in units of Debye \AA. |
| 1211 |
|
|
| 1212 |
< |
\subsubsection{\label{section:ffGB}The FluctuatingChargeAtomTypes block} |
| 1213 |
< |
\subsubsection{\label{section:ffPol}The PolarizableAtomTypes block} |
| 1239 |
< |
\subsubsection{\label{section:ffGB}The GayBerneAtomTypes block} |
| 1212 |
> |
\subsection{\label{section:ffGB}The FluctuatingChargeAtomTypes block} |
| 1213 |
> |
%\subsubsection{\label{section:ffPol}The PolarizableAtomTypes block} |
| 1214 |
|
|
| 1215 |
+ |
\subsection{\label{section:ffGB}The GayBerneAtomTypes block} |
| 1216 |
+ |
|
| 1217 |
|
The Gay-Berne potential has been widely used in the liquid crystal |
| 1218 |
< |
community to describe this anisotropic phase |
| 1218 |
> |
community to describe anisotropic phase |
| 1219 |
|
behavior.~\cite{Gay:1981yu,Berne:1972pb,Kushick:1976xy,Luckhurst:1990fy,Cleaver:1996rt} |
| 1220 |
|
The form of the Gay-Berne potential implemented in OpenMD was |
| 1221 |
|
generalized by Cleaver {\it et al.} and is appropriate for dissimilar |
| 1222 |
< |
uniaxial ellipsoids.\cite{Cleaver:1996rt} The potential is constructed in the |
| 1223 |
< |
familiar form of the Lennard-Jones function using |
| 1222 |
> |
uniaxial ellipsoids.\cite{Cleaver:1996rt} The potential is constructed |
| 1223 |
> |
in the familiar form of the Lennard-Jones function using |
| 1224 |
|
orientation-dependent $\sigma$ and $\epsilon$ parameters, |
| 1225 |
|
\begin{equation*} |
| 1226 |
|
V_{ij}({{\bf \hat u}_i}, {{\bf \hat u}_j}, {{\bf \hat |
| 1278 |
|
end GayBerneAtomTypes |
| 1279 |
|
\end{lstlisting} |
| 1280 |
|
|
| 1281 |
< |
\subsubsection{\label{section:ffSticky}The StickyAtomTypes block} |
| 1281 |
> |
\subsection{\label{section:ffSticky}The StickyAtomTypes block} |
| 1282 |
|
|
| 1283 |
|
One of the solvents that can be simulated by {\sc OpenMD} is the |
| 1284 |
|
extended Soft Sticky Dipole (SSD/E) water model.\cite{fennell04} The |
| 1388 |
|
end StickyAtomTypes |
| 1389 |
|
\end{lstlisting} |
| 1390 |
|
|
| 1391 |
< |
\subsection{\label{section::ffMetals}Metallic Atom Types} |
| 1391 |
> |
\section{\label{section::ffMetals}Metallic Atom Types} |
| 1392 |
|
|
| 1393 |
|
{\sc OpenMD} implements a number of related potentials that describe |
| 1394 |
|
bonding in transition metals. These potentials have an attractive |
| 1417 |
|
The pairwise portion of the potential, $\phi_{ij}$, is usually a |
| 1418 |
|
repulsive interaction between atoms $i$ and $j$. |
| 1419 |
|
|
| 1420 |
< |
\subsubsection{\label{section:ffEAM}The EAMAtomTypes block} |
| 1420 |
> |
\subsection{\label{section:ffEAM}The EAMAtomTypes block} |
| 1421 |
|
The Embedded Atom Method ({\sc eam}) is one of the most widely-used |
| 1422 |
|
potentials for transition |
| 1423 |
|
metals.~\cite{Finnis84,Ercolessi88,Chen90,Qi99,Ercolessi02,Daw84,FBD86,johnson89,Lu97} |
| 1483 |
|
end EAMAtomTypes |
| 1484 |
|
\end{lstlisting} |
| 1485 |
|
|
| 1486 |
+ |
\subsection{\label{section:ffSC}The SuttonChenAtomTypes block} |
| 1487 |
|
|
| 1511 |
– |
\subsubsection{\label{section:ffSC}The SuttonChenAtomTypes block} |
| 1512 |
– |
|
| 1488 |
|
The Sutton-Chen ({\sc sc})~\cite{Chen90} potential has been used to |
| 1489 |
|
study a wide range of phenomena in metals. Although it has the same |
| 1490 |
< |
basic form as the {\sc eam} potential, the Sutton-Chen model takes on |
| 1491 |
< |
a simpler form, |
| 1490 |
> |
basic form as the {\sc eam} potential, the Sutton-Chen model requires |
| 1491 |
> |
a simpler set of parameters, |
| 1492 |
|
\begin{equation} |
| 1493 |
|
\label{eq:SCP1} |
| 1494 |
|
U_{tot}=\sum _{i}\left[ \frac{1}{2}\sum _{j\neq |
| 1538 |
|
end SCAtomTypes |
| 1539 |
|
\end{lstlisting} |
| 1540 |
|
|
| 1541 |
< |
\subsection{\label{section::ffShortRange}Short Range Interactions} |
| 1542 |
< |
\subsubsection{\label{section:ffBond}The BondTypes block} |
| 1543 |
< |
\subsubsection{\label{section:ffBend}The BendTypes block} |
| 1544 |
< |
A harmonic bend potential is represented by the following function: |
| 1570 |
< |
\begin{equation} |
| 1571 |
< |
V_{\text{bend}}(\theta_{ijk}) = k_{\theta}( \theta_{ijk} - \theta_0 |
| 1572 |
< |
)^2, \label{eq:bendPot} |
| 1573 |
< |
\end{equation} |
| 1574 |
< |
where $\theta_{ijk}$ is the angle defined by atoms $i$, $j$, and $k$, |
| 1575 |
< |
$\theta_0$ is the equilibrium bond angle, and $k_{\theta}$ is the |
| 1576 |
< |
force constant which determines the strength of the harmonic bend. |
| 1577 |
< |
|
| 1578 |
< |
\subsubsection{\label{section:ffTorsion}The TorsionTypes block} |
| 1579 |
< |
The torsion potential is often represented as a cosine series of the |
| 1580 |
< |
form: |
| 1581 |
< |
\begin{equation} |
| 1582 |
< |
V_{\text{torsion}}(\phi) = c_1[1 + \cos \phi] |
| 1583 |
< |
+ c_2[1 + \cos(2\phi)] |
| 1584 |
< |
+ c_3[1 + \cos(3\phi)], |
| 1585 |
< |
\label{eq:origTorsionPot} |
| 1586 |
< |
\end{equation} |
| 1587 |
< |
where: |
| 1588 |
< |
\begin{equation} |
| 1589 |
< |
\cos\phi = (\hat{\mathbf{r}}_{ij} \times \hat{\mathbf{r}}_{jk}) \cdot |
| 1590 |
< |
(\hat{\mathbf{r}}_{jk} \times \hat{\mathbf{r}}_{kl}). |
| 1591 |
< |
\label{eq:torsPhi} |
| 1592 |
< |
\end{equation} |
| 1593 |
< |
Here, $\hat{\mathbf{r}}_{\alpha\beta}$ are the set of unit bond |
| 1594 |
< |
vectors between atoms $i$, $j$, $k$, and $l$. For computational |
| 1595 |
< |
efficiency, the torsion potential has been recast after the method of |
| 1596 |
< |
{\sc charmm},\cite{Brooks83} in which the angle series is converted to |
| 1597 |
< |
a power series of the form: |
| 1598 |
< |
\begin{equation} |
| 1599 |
< |
V_{\text{torsion}}(\phi) = |
| 1600 |
< |
k_3 \cos^3 \phi + k_2 \cos^2 \phi + k_1 \cos \phi + k_0, |
| 1601 |
< |
\label{eq:torsionPot} |
| 1602 |
< |
\end{equation} |
| 1603 |
< |
where: |
| 1541 |
> |
\section{\label{section::ffShortRange}Short Range Interactions} |
| 1542 |
> |
The internal structure of a molecule is usually specified in terms of |
| 1543 |
> |
a set of ``bonded'' terms in the potential energy function for |
| 1544 |
> |
molecule $I$, |
| 1545 |
|
\begin{align*} |
| 1546 |
< |
k_0 &= c_1 + c_3, \\ |
| 1547 |
< |
k_1 &= c_1 - 3c_3, \\ |
| 1548 |
< |
k_2 &= 2 c_2, \\ |
| 1549 |
< |
k_3 &= 4c_3. |
| 1546 |
> |
V^{I}_{\text{Internal}} = & |
| 1547 |
> |
\sum_{r_{ij} \in I} V_{\text{bond}}(r_{ij}) |
| 1548 |
> |
+ \sum_{\theta_{ijk} \in I} V_{\text{bend}}(\theta_{ijk}) |
| 1549 |
> |
+ \sum_{\phi_{ijkl} \in I} V_{\text{torsion}}(\phi_{ijkl}) |
| 1550 |
> |
+ \sum_{\omega_{ijkl} \in I} V_{\text{inversion}}(\omega_{ijkl}) \\ |
| 1551 |
> |
& + \sum_{i \in I} \sum_{(j>i+4) \in I} |
| 1552 |
> |
\biggl[ V_{\text{dispersion}}(r_{ij}) + V_{\text{electrostatic}} |
| 1553 |
> |
(\mathbf{r}_{ij},\boldsymbol{\Omega}_{i},\boldsymbol{\Omega}_{j}) |
| 1554 |
> |
\biggr]. |
| 1555 |
> |
\label{eq:internalPotential} |
| 1556 |
|
\end{align*} |
| 1557 |
< |
By recasting the potential as a power series, repeated trigonometric |
| 1558 |
< |
evaluations are avoided during the calculation of the potential |
| 1559 |
< |
energy. |
| 1557 |
> |
Here $V_{\text{bond}}, V_{\text{bend}}, |
| 1558 |
> |
V_{\text{torsion}},\mathrm{~and~} V_{\text{inversion}}$ represent the |
| 1559 |
> |
bond, bend, torsion, and inversion potentials for all |
| 1560 |
> |
topologically-connected sets of atoms within the molecule. Bonds are |
| 1561 |
> |
the primary way of specifying how the atoms are connected together to |
| 1562 |
> |
form the molecule (i.e. they define the molecular topology). The |
| 1563 |
> |
other short-range interactions may be specified explicitly in the |
| 1564 |
> |
molecule definition, or they may be deduced from bonding information. |
| 1565 |
> |
For example, bends can be implicitly deduced from two bonds which |
| 1566 |
> |
share an atom. Torsions can be deduced from two bends that share a |
| 1567 |
> |
bond. Inversion potentials are utilized primarily to enforce |
| 1568 |
> |
planarity around $sp^2$-hybridized sites, and these are specified with |
| 1569 |
> |
central atoms and satellites (or an atom with bonds to exactly three |
| 1570 |
> |
satellites). Non-bonded interactions are usually excluded for atom |
| 1571 |
> |
pairs that are involved in the same bond, bend, or torsion, but all |
| 1572 |
> |
other atom pairs are included in the standard non-bonded interactions. |
| 1573 |
|
|
| 1574 |
< |
\subsubsection{\label{section:ffInversion}The InversionTypes block} |
| 1575 |
< |
\subsection{\label{section::ffLongRange}Long Range Interactions} |
| 1576 |
< |
\subsubsection{\label{section:ffNBinteraction}The NonBondedInteraction block} |
| 1574 |
> |
Bond lengths, angles, and torsions (dihedrals) are ``natural'' |
| 1575 |
> |
coordinates to treat molecular motion, as it is usually in these |
| 1576 |
> |
coordinates that most chemists understand the behavior of molecules. |
| 1577 |
> |
The bond lengths and angles are often considered ``hard'' degrees of |
| 1578 |
> |
freedom. That is, we can't deform them very much without a |
| 1579 |
> |
significant energetic penalty. On the other hand, dihedral angles or |
| 1580 |
> |
torsions are ``soft'' and typically undergo significant deformation |
| 1581 |
> |
under normal conditions. |
| 1582 |
|
|
| 1583 |
+ |
\subsection{\label{section:ffBond}The BondTypes block} |
| 1584 |
|
|
| 1585 |
+ |
Bonds are the primary way to specify how the atoms are connected |
| 1586 |
+ |
together to form the molecule. In general, bonds exert strong |
| 1587 |
+ |
restoring forces to keep the molecule compact. The bond energy |
| 1588 |
+ |
functions are relatively simple functions of the distance between two |
| 1589 |
+ |
atomic sites, |
| 1590 |
+ |
\begin{equation} |
| 1591 |
+ |
b = \left| \vec{r}_{ij} \right| = \sqrt{(x_j - x_i)^2 + (y_j - y_i)^2 |
| 1592 |
+ |
+ (z_j - z_i)^2}. |
| 1593 |
+ |
\end{equation} |
| 1594 |
+ |
All BondTypes must specify two AtomType names ($i$ and $j$) that |
| 1595 |
+ |
describe when that bond should be applied, as well as an equilibrium |
| 1596 |
+ |
bond length, $b_{ij}^0$, in units of \AA. The most common forms for |
| 1597 |
+ |
bonding potentials are {\tt Harmonic} bonds, |
| 1598 |
+ |
\begin{equation} |
| 1599 |
+ |
V_{\text{bond}}(b) = \frac{k_{ij}}{2} \left(b - |
| 1600 |
+ |
b_{ij}^0 \right)^2 |
| 1601 |
+ |
\end{equation} |
| 1602 |
+ |
and {\tt Morse} bonds, |
| 1603 |
+ |
\begin{equation} |
| 1604 |
+ |
V_{\text{bond}}(b) = D_{ij} \left[ 1 - e^{-\beta_{ij} (b - b_{ij}^0)} \right]^2 |
| 1605 |
+ |
\end{equation} |
| 1606 |
|
|
| 1607 |
< |
(see Fig.~\ref{fig:lipidModel}), The parameters for $k_{\theta}$ and |
| 1608 |
< |
$\theta_0$ are borrowed from those in TraPPE.\cite{Siepmann1998} |
| 1607 |
> |
OpenMD can also simulate some less common types of bond potentials, |
| 1608 |
> |
including {\tt Fixed} bonds (which are constrained to be at a |
| 1609 |
> |
specified bond length), |
| 1610 |
> |
\begin{equation} |
| 1611 |
> |
V_{\text{bond}}(b) = 0. |
| 1612 |
> |
\end{equation} |
| 1613 |
> |
{\tt Cubic} bonds include terms to model anharmonicity, |
| 1614 |
> |
\begin{equation} |
| 1615 |
> |
V_{\text{bond}}(b) = K_3 (b - b_{ij}^0)^3 + K_2 (b - b_{ij}^0)^2 + K_1 (b - b_{ij}^0) + K_0, |
| 1616 |
> |
\end{equation} |
| 1617 |
> |
and {\tt Quartic} bonds, include another term in the Taylor |
| 1618 |
> |
expansion around $b_{ij}^0$, |
| 1619 |
> |
\begin{equation} |
| 1620 |
> |
V_{\text{bond}}(b) = K_4 (b - b_{ij}^0)^4 + K_3 (b - b_{ij}^0)^3 + |
| 1621 |
> |
K_2 (b - b_{ij}^0)^2 + K_1 (b - b_{ij}^0) + K_0, |
| 1622 |
> |
\end{equation} |
| 1623 |
> |
can also be simulated. Note that the factor of $1/2$ that is included |
| 1624 |
> |
in the {\tt Harmonic} bond type force constant is {\it not} present in |
| 1625 |
> |
either the {\tt Cubic} or {\tt Quartic} bond types. |
| 1626 |
|
|
| 1627 |
< |
Calculating the long-range (non-bonded) potential involves a sum over |
| 1624 |
< |
all pairs of atoms (except for those atoms which are involved in a |
| 1625 |
< |
bond, bend, or torsion with each other). If done poorly, calculating |
| 1626 |
< |
the the long-range interactions for $N$ atoms would involve $N(N-1)/2$ |
| 1627 |
< |
evaluations of atomic distances. To reduce the number of distance |
| 1628 |
< |
evaluations between pairs of atoms, {\sc OpenMD} allows the use of |
| 1629 |
< |
switched cutoffs with Verlet neighbor lists.\cite{Allen87} Neutral |
| 1630 |
< |
groups which contain charges will exhibit pathological forces unless |
| 1631 |
< |
the cutoff is applied to the neutral groups evenly instead of to the |
| 1632 |
< |
individual atoms.\cite{leach01:mm} {\sc OpenMD} allows users to |
| 1633 |
< |
specify cutoff groups which may contain an arbitrary number of atoms |
| 1634 |
< |
in the molecule. Atoms in a cutoff group are treated as a single unit |
| 1635 |
< |
for the evaluation of the switching function: |
| 1627 |
> |
{\tt Polynomial} bonds which can have any number of terms, |
| 1628 |
|
\begin{equation} |
| 1629 |
< |
V_{\mathrm{long-range}} = \sum_{a} \sum_{b>a} s(r_{ab}) \sum_{i \in a} \sum_{j \in b} V_{ij}(r_{ij}), |
| 1629 |
> |
V_{\text{bond}}(b) = \sum_n K_n (b - b_{ij}^0)^n. |
| 1630 |
|
\end{equation} |
| 1631 |
< |
where $r_{ab}$ is the distance between the centers of mass of the two |
| 1632 |
< |
cutoff groups ($a$ and $b$). |
| 1631 |
> |
can also be specified by giving a sequence of integer ($n$) and force |
| 1632 |
> |
constant ($K_n$) pairs. |
| 1633 |
|
|
| 1634 |
< |
The sums over $a$ and $b$ are over the cutoff groups that are present |
| 1635 |
< |
in the simulation. Atoms which are not explicitly defined as members |
| 1636 |
< |
of a {\tt cutoffGroup} are treated as a group consisting of only one |
| 1637 |
< |
atom. The switching function, $s(r)$ is the standard cubic switching |
| 1638 |
< |
function, |
| 1634 |
> |
The order of terms in the BondTypes block is: |
| 1635 |
> |
\begin{itemize} |
| 1636 |
> |
\item {\tt AtomType} 1 |
| 1637 |
> |
\item {\tt AtomType} 2 |
| 1638 |
> |
\item {\tt BondType} (one of {\tt Harmonic}, {\tt Morse}, {\tt Fixed}, {\tt |
| 1639 |
> |
Cubic}, {\tt Quartic}, or {\tt Polynomial}) |
| 1640 |
> |
\item $b_{ij}^0$, the equilibrium bond length in \AA |
| 1641 |
> |
\item any other parameters required by the {\tt BondType} |
| 1642 |
> |
\end{itemize} |
| 1643 |
> |
|
| 1644 |
> |
\begin{lstlisting}[caption={[An example of a BondTypes block.] A |
| 1645 |
> |
simple example of a BondTypes block. Distances ($b_0$) |
| 1646 |
> |
are given in \AA\ and force constants are given in |
| 1647 |
> |
units so that when multiplied by the correct power of distance they |
| 1648 |
> |
return energies in kcal/mol. For example $k$ for a Harmonic bond is |
| 1649 |
> |
in units of kcal/mol/\AA$^2$.}, |
| 1650 |
> |
label={sch:BondTypes}] |
| 1651 |
> |
begin BondTypes |
| 1652 |
> |
//Atom1 Atom2 Harmonic b0 k (kcal/mol/A^2) |
| 1653 |
> |
CH2 CH2 Harmonic 1.526 260 |
| 1654 |
> |
//Atom1 Atom2 Morse b0 D beta (A^-1) |
| 1655 |
> |
CN NC Morse 1.157437 212.95 2.5802 |
| 1656 |
> |
//Atom1 Atom2 Fixed b0 |
| 1657 |
> |
CT HC Fixed 1.09 |
| 1658 |
> |
//Atom1 Atom2 Cubic b0 K3 K2 K1 K0 |
| 1659 |
> |
//Atom1 Atom2 Quartic b0 K4 K3 K2 K1 K0 |
| 1660 |
> |
//Atom1 Atom2 Polynomial b0 n Kn [m Km] |
| 1661 |
> |
end BondTypes |
| 1662 |
> |
\end{lstlisting} |
| 1663 |
> |
|
| 1664 |
> |
There are advantages and disadvantages of all of the different types |
| 1665 |
> |
of bonds, but specific simulation tasks may call for specific |
| 1666 |
> |
behaviors. |
| 1667 |
> |
|
| 1668 |
> |
\subsection{\label{section:ffBend}The BendTypes block} |
| 1669 |
> |
The equilibrium geometries and energy functions for bending motions in |
| 1670 |
> |
a molecule are strongly dependent on the bonding environment of the |
| 1671 |
> |
central atomic site. For example, different types of hybridized |
| 1672 |
> |
carbon centers require different bending angles and force constants to |
| 1673 |
> |
describe the local geometry. |
| 1674 |
> |
|
| 1675 |
> |
The bending potential energy functions used in most force fields are |
| 1676 |
> |
often simple functions of the angle between two bonds, |
| 1677 |
|
\begin{equation} |
| 1678 |
< |
S(r) = |
| 1679 |
< |
\begin{cases} |
| 1680 |
< |
1 & \text{if $r \le r_{\text{sw}}$},\\ |
| 1681 |
< |
\frac{(r_{\text{cut}} + 2r - 3r_{\text{sw}})(r_{\text{cut}} - r)^2} |
| 1682 |
< |
{(r_{\text{cut}} - r_{\text{sw}})^3} |
| 1683 |
< |
& \text{if $r_{\text{sw}} < r \le r_{\text{cut}}$}, \\ |
| 1684 |
< |
0 & \text{if $r > r_{\text{cut}}$.} |
| 1685 |
< |
\end{cases} |
| 1686 |
< |
\label{eq:dipoleSwitching} |
| 1678 |
> |
\theta_{ijk} = \cos^{-1} \left(\frac{\vec{r}_{ij} \cdot |
| 1679 |
> |
\vec{r}_{jk}}{\left| \vec{r}_{ij} \right| \left| \vec{r}_{ij} |
| 1680 |
> |
\right|} \right) |
| 1681 |
> |
\end{equation} |
| 1682 |
> |
Here atom $j$ is the central atom that is bonded to two partners $i$ |
| 1683 |
> |
and $k$. |
| 1684 |
> |
|
| 1685 |
> |
All BendTypes must specify three AtomType names ($i$, $j$ and $k$) |
| 1686 |
> |
that describe when that bend potential should be applied, as well as |
| 1687 |
> |
an equilibrium bending angle, $\theta_{ijk}^0$, in units of |
| 1688 |
> |
degrees. The most common forms for bending potentials are {\tt |
| 1689 |
> |
Harmonic} bends, |
| 1690 |
> |
\begin{equation} |
| 1691 |
> |
V_{\text{bend}}(\theta_{ijk}) = \frac{k_{ijk}}{2}( \theta_{ijk} - \theta_{ijk}^0 |
| 1692 |
> |
)^2, \label{eq:bendPot} |
| 1693 |
|
\end{equation} |
| 1694 |
< |
Here, $r_{\text{sw}}$ is the {\tt switchingRadius}, or the distance |
| 1695 |
< |
beyond which interactions are reduced, and $r_{\text{cut}}$ is the |
| 1696 |
< |
{\tt cutoffRadius}, or the distance at which interactions are |
| 1697 |
< |
truncated. |
| 1694 |
> |
where $k_{ijk}$ is the force constant which determines the strength of |
| 1695 |
> |
the harmonic bend. {\tt UreyBradley} bends utilize an additional 1-3 |
| 1696 |
> |
bond-type interaction in addition to the harmonic bending potential: |
| 1697 |
> |
\begin{equation} |
| 1698 |
> |
V_{\text{bend}}(\vec{r}_i , \vec{r}_j, \vec{r}_k) |
| 1699 |
> |
=\frac{k_{ijk}}{2}( \theta_{ijk} - \theta_{ijk}^0)^2 |
| 1700 |
> |
+ \frac{k_{ub}}{2}( r_{ik} - s_0 )^2. \label{eq:ubBend} |
| 1701 |
> |
\end{equation} |
| 1702 |
|
|
| 1703 |
< |
Users of {\sc OpenMD} do not need to specify the {\tt cutoffRadius} or |
| 1704 |
< |
{\tt switchingRadius}. In simulations containing only Lennard-Jones |
| 1705 |
< |
atoms, the cutoff radius has a default value of $2.5\sigma_{ii}$, |
| 1706 |
< |
where $\sigma_{ii}$ is the largest Lennard-Jones length parameter |
| 1707 |
< |
present in the simulation. In simulations containing charged or |
| 1708 |
< |
dipolar atoms, the default cutoff radius is $15 \mbox{\AA}$. |
| 1703 |
> |
A {\tt Cosine} bend is a variant on the harmonic bend which utilizes |
| 1704 |
> |
the cosine of the angle instead of the angle itself, |
| 1705 |
> |
\begin{equation} |
| 1706 |
> |
V_{\text{bend}}(\theta_{ijk}) = \frac{k_{ijk}}{2}( \cos\theta_{ijk} - |
| 1707 |
> |
\cos \theta_{ijk}^0 )^2. \label{eq:cosBend} |
| 1708 |
> |
\end{equation} |
| 1709 |
|
|
| 1710 |
< |
The {\tt switchingRadius} is set to a default value of 95\% of the |
| 1711 |
< |
{\tt cutoffRadius}. In the special case of a simulation containing |
| 1712 |
< |
{\it only} Lennard-Jones atoms, the default switching radius takes the |
| 1713 |
< |
same value as the cutoff radius, and {\sc OpenMD} will use a shifted |
| 1714 |
< |
potential to remove discontinuities in the potential at the cutoff. |
| 1715 |
< |
Both radii may be specified in the meta-data file. |
| 1710 |
> |
OpenMD can also simulate some less common types of bend potentials, |
| 1711 |
> |
including {\tt Cubic} bends, which include terms to model anharmonicity, |
| 1712 |
> |
\begin{equation} |
| 1713 |
> |
V_{\text{bend}}(\theta_{ijk}) = K_3 (\theta_{ijk} - \theta_{ijk}^0)^3 + K_2 (\theta_{ijk} - \theta_{ijk}^0)^2 + K_1 (\theta_{ijk} - \theta_{ijk}^0) + K_0, |
| 1714 |
> |
\end{equation} |
| 1715 |
> |
and {\tt Quartic} bends, which include another term in the Taylor |
| 1716 |
> |
expansion around $\theta_{ijk}^0$, |
| 1717 |
> |
\begin{equation} |
| 1718 |
> |
V_{\text{bend}}(\theta_{ijk}) = K_4 (\theta_{ijk} - \theta_{ijk}^0)^4 + K_3 (\theta_{ijk} - \theta_{ijk}^0)^3 + |
| 1719 |
> |
K_2 (\theta_{ijk} - \theta_{ijk}^0)^2 + K_1 (\theta_{ijk} - |
| 1720 |
> |
\theta_{ijk}^0) + K_0, |
| 1721 |
> |
\end{equation} |
| 1722 |
> |
can also be simulated. Note that the factor of $1/2$ that is included |
| 1723 |
> |
in the {\tt Harmonic} bend type force constant is {\it not} present in |
| 1724 |
> |
either the {\tt Cubic} or {\tt Quartic} bend types. |
| 1725 |
|
|
| 1726 |
< |
Force fields can be added to {\sc OpenMD}, although it comes with a few |
| 1727 |
< |
simple examples (Lennard-Jones, {\sc duff}, {\sc water}, and {\sc |
| 1728 |
< |
eam}) which are explained in the following sections. |
| 1726 |
> |
{\tt Polynomial} bends which can have any number of terms, |
| 1727 |
> |
\begin{equation} |
| 1728 |
> |
V_{\text{bend}}(\theta_{ijk}) = \sum_n K_n (\theta_{ijk} - \theta_{ijk}^0)^n. |
| 1729 |
> |
\end{equation} |
| 1730 |
> |
can also be specified by giving a sequence of integer ($n$) and force |
| 1731 |
> |
constant ($K_n$) pairs. |
| 1732 |
|
|
| 1733 |
< |
\section{\label{sec:LJPot}The Lennard Jones Force Field} |
| 1733 |
> |
The order of terms in the BendTypes block is: |
| 1734 |
> |
\begin{itemize} |
| 1735 |
> |
\item {\tt AtomType} 1 |
| 1736 |
> |
\item {\tt AtomType} 2 (this is the central atom) |
| 1737 |
> |
\item {\tt AtomType} 3 |
| 1738 |
> |
\item {\tt BendType} (one of {\tt Harmonic}, {\tt UreyBradley}, {\tt |
| 1739 |
> |
Cosine}, {\tt Cubic}, {\tt Quartic}, or {\tt Polynomial}) |
| 1740 |
> |
\item $\theta_{ijk}^0$, the equilibrium bending angle in degrees. |
| 1741 |
> |
\item any other parameters required by the {\tt BendType} |
| 1742 |
> |
\end{itemize} |
| 1743 |
|
|
| 1744 |
< |
Scheme |
| 1745 |
< |
\ref{sch:LJFF} gives an example meta-data file that |
| 1746 |
< |
sets up a system of 108 Ar particles to be simulated using the |
| 1747 |
< |
Lennard-Jones force field. |
| 1744 |
> |
\begin{lstlisting}[caption={[An example of a BendTypes block.] A |
| 1745 |
> |
simple example of a BendTypes block. By convention, equilibrium angles |
| 1746 |
> |
($\theta_0$) are given in degrees but force constants are given in |
| 1747 |
> |
units so that when multiplied by the correct power of angle (in |
| 1748 |
> |
radians) they return energies in kcal/mol. For example $k$ for a |
| 1749 |
> |
Harmonic bend is in units of kcal/mol/radians$^2$.}, |
| 1750 |
> |
label={sch:BendTypes}] |
| 1751 |
> |
begin BendTypes |
| 1752 |
> |
//Atom1 Atom2 Atom3 Harmonic theta0(deg) Ktheta(kcal/mol/radians^2) |
| 1753 |
> |
CT CT CT Harmonic 109.5 80.000000 |
| 1754 |
> |
CH2 CH CH2 Harmonic 112.0 117.68 |
| 1755 |
> |
CH3 CH2 SH Harmonic 96.0 67.220 |
| 1756 |
> |
//UreyBradley |
| 1757 |
> |
//Atom1 Atom2 Atom3 UreyBradley theta0 Ktheta s0 Kub |
| 1758 |
> |
//Cosine |
| 1759 |
> |
//Atom1 Atom2 Atom3 Cosine theta0 Ktheta(kcal/mol) |
| 1760 |
> |
//Cubic |
| 1761 |
> |
//Atom1 Atom2 Atom3 Cubic theta0 K3 K2 K1 K0 |
| 1762 |
> |
//Quartic |
| 1763 |
> |
//Atom1 Atom2 Atom3 Quartic theta0 K4 K3 K2 K1 K0 |
| 1764 |
> |
//Polynomial |
| 1765 |
> |
//Atom1 Atom2 Atom3 Polynomial theta0 n Kn [m Km] |
| 1766 |
> |
end BendTypes |
| 1767 |
> |
\end{lstlisting} |
| 1768 |
|
|
| 1769 |
< |
\begin{lstlisting}[float,caption={[Invocation of the Lennard-Jones |
| 1770 |
< |
force field] A sample startup file for a small Lennard-Jones |
| 1771 |
< |
simulation.},label={sch:LJFF}] |
| 1772 |
< |
<OpenMD> |
| 1692 |
< |
<MetaData> |
| 1693 |
< |
#include "argon.md" |
| 1769 |
> |
Note that the parameters for a particular bend type are the same for |
| 1770 |
> |
any bending triplet of the same atomic types (in the same or reversed |
| 1771 |
> |
order). Changing the AtomType in the Atom2 position will change the |
| 1772 |
> |
matched bend types in the force field. |
| 1773 |
|
|
| 1774 |
< |
component{ |
| 1775 |
< |
type = "Ar"; |
| 1776 |
< |
nMol = 108; |
| 1777 |
< |
} |
| 1774 |
> |
\subsection{\label{section:ffTorsion}The TorsionTypes block} |
| 1775 |
> |
The torsion potential for rotation of groups around a central bond can |
| 1776 |
> |
often be represented with various cosine functions. For two |
| 1777 |
> |
tetrahedral ($sp^3$) carbons connected by a single bond, the torsion |
| 1778 |
> |
potential might be |
| 1779 |
> |
\begin{equation*} |
| 1780 |
> |
V_{\text{torsion}} = \frac{v}{2} \left[ 1 + \cos( 3 \phi ) \right] |
| 1781 |
> |
\end{equation*} |
| 1782 |
> |
where $v$ is the barrier for going from {\em staggered} $\rightarrow$ |
| 1783 |
> |
{\em eclipsed} conformations, while for $sp^2$ carbons connected by a |
| 1784 |
> |
double bond, the potential might be |
| 1785 |
> |
\begin{equation*} |
| 1786 |
> |
V_{\text{torsion}} = \frac{w}{2} \left[ 1 - \cos( 2 \phi ) \right] |
| 1787 |
> |
\end{equation*} |
| 1788 |
> |
where $w$ is the barrier for going from {\em cis} $\rightarrow$ {\em |
| 1789 |
> |
trans} conformations. |
| 1790 |
|
|
| 1791 |
< |
forceField = "LJ"; |
| 1792 |
< |
</MetaData> |
| 1793 |
< |
<Snapshot> // not shown in this scheme |
| 1794 |
< |
</Snapshot> |
| 1795 |
< |
</OpenMD> |
| 1796 |
< |
\end{lstlisting} |
| 1791 |
> |
A general torsion potentials can be represented as a cosine series of |
| 1792 |
> |
the form: |
| 1793 |
> |
\begin{equation} |
| 1794 |
> |
V_{\text{torsion}}(\phi_{ijkl}) = c_1[1 + \cos \phi_{ijkl}] |
| 1795 |
> |
+ c_2[1 - \cos(2\phi_{ijkl})] |
| 1796 |
> |
+ c_3[1 + \cos(3\phi_{ijkl})], |
| 1797 |
> |
\label{eq:origTorsionPot} |
| 1798 |
> |
\end{equation} |
| 1799 |
> |
where the angle $\phi_{ijkl}$ is defined |
| 1800 |
> |
\begin{equation} |
| 1801 |
> |
\cos\phi_{ijkl} = (\hat{\mathbf{r}}_{ij} \times \hat{\mathbf{r}}_{jk}) \cdot |
| 1802 |
> |
(\hat{\mathbf{r}}_{jk} \times \hat{\mathbf{r}}_{kl}). |
| 1803 |
> |
\label{eq:torsPhi} |
| 1804 |
> |
\end{equation} |
| 1805 |
> |
Here, $\hat{\mathbf{r}}_{\alpha\beta}$ are the set of unit bond |
| 1806 |
> |
vectors between atoms $i$, $j$, $k$, and $l$. |
| 1807 |
> |
|
| 1808 |
> |
For computational efficiency, OpenMD recasts torsion potential in the |
| 1809 |
> |
method of {\sc charmm},\cite{Brooks83} in which the angle series is |
| 1810 |
> |
converted to a power series of the form: |
| 1811 |
> |
\begin{equation} |
| 1812 |
> |
V_{\text{torsion}}(\phi_{ijkl}) = |
| 1813 |
> |
k_3 \cos^3 \phi_{ijkl} + k_2 \cos^2 \phi_{ijkl} + k_1 \cos \phi_{ijkl} + k_0, |
| 1814 |
> |
\label{eq:torsionPot} |
| 1815 |
> |
\end{equation} |
| 1816 |
> |
where: |
| 1817 |
> |
\begin{align*} |
| 1818 |
> |
k_0 &= c_1 + 2 c_2 + c_3, \\ |
| 1819 |
> |
k_1 &= c_1 - 3c_3, \\ |
| 1820 |
> |
k_2 &= - 2 c_2, \\ |
| 1821 |
> |
k_3 &= 4 c_3. |
| 1822 |
> |
\end{align*} |
| 1823 |
> |
By recasting the potential as a power series, repeated trigonometric |
| 1824 |
> |
evaluations are avoided during the calculation of the potential |
| 1825 |
> |
energy. |
| 1826 |
|
|
| 1827 |
+ |
Using this framework, OpenMD implements a variety of different |
| 1828 |
+ |
potential energy functions for torsions: |
| 1829 |
+ |
\begin{itemize} |
| 1830 |
+ |
\item {\tt Cubic}: |
| 1831 |
+ |
\begin{equation*} |
| 1832 |
+ |
V_{\text{torsion}}(\phi) = |
| 1833 |
+ |
k_3 \cos^3 \phi + k_2 \cos^2 \phi + k_1 \cos \phi + k_0, |
| 1834 |
+ |
\end{equation*} |
| 1835 |
+ |
\item {\tt Quartic}: |
| 1836 |
+ |
\begin{equation*} |
| 1837 |
+ |
V_{\text{torsion}}(\phi) = k_4 \cos^4 \phi + |
| 1838 |
+ |
k_3 \cos^3 \phi + k_2 \cos^2 \phi + k_1 \cos \phi + k_0, |
| 1839 |
+ |
\end{equation*} |
| 1840 |
+ |
\item {\tt Polynomial}: |
| 1841 |
+ |
\begin{equation*} |
| 1842 |
+ |
V_{\text{torsion}}(\phi) = \sum_n k_n \cos^n \phi , |
| 1843 |
+ |
\end{equation*} |
| 1844 |
+ |
\item {\tt Charmm}: |
| 1845 |
+ |
\begin{equation*} |
| 1846 |
+ |
V_{\text{torsion}}(\phi) = \sum_n K_n \left( 1 + cos(n |
| 1847 |
+ |
\phi - \delta_n) \right), |
| 1848 |
+ |
\end{equation*} |
| 1849 |
+ |
\item {\tt Opls}: |
| 1850 |
+ |
\begin{equation*} |
| 1851 |
+ |
V_{\text{torsion}}(\phi) = \frac{1}{2} \left(v_1 (1 + \cos \phi) \right) |
| 1852 |
+ |
+ v_2 (1 - \cos 2 \phi) + v_3 (1 + \cos 3 \phi), |
| 1853 |
+ |
\end{equation*} |
| 1854 |
+ |
\item {\tt Trappe}:\cite{Siepmann1998} |
| 1855 |
+ |
\begin{equation*} |
| 1856 |
+ |
V_{\text{torsion}}(\phi) = c_0 + c_1 (1 + \cos \phi) + c_2 (1 - \cos 2 \phi) + |
| 1857 |
+ |
c_3 (1 + \cos 3 \phi), |
| 1858 |
+ |
\end{equation*} |
| 1859 |
+ |
\item {\tt Harmonic}: |
| 1860 |
+ |
\begin{equation*} |
| 1861 |
+ |
V_{\text{torsion}}(\phi) = \frac{d_0}{2} \left(\phi - \phi^0\right). |
| 1862 |
+ |
\end{equation*} |
| 1863 |
+ |
\end{itemize} |
| 1864 |
|
|
| 1865 |
< |
\section{\label{section:DUFF}Dipolar Unified-Atom Force Field} |
| 1865 |
> |
Most torsion types don't require specific angle information in the |
| 1866 |
> |
parameters since they are typically expressed in cosine polynomials. |
| 1867 |
> |
{\tt Charmm} and {\tt Harmonic} torsions are a bit different. {\tt |
| 1868 |
> |
Charmm} torsion types require a set of phase angles, $\delta_n$ that |
| 1869 |
> |
are expressed in degrees, and associated periodicity numbers, $n$. |
| 1870 |
> |
{\tt Harmonic} torsions have an equilibrium torsion angle, $\phi_0$ |
| 1871 |
> |
that is measured in degrees, while $d_0$ has units of |
| 1872 |
> |
kcal/mol/degrees$^2$. All other torsion parameters are measured in |
| 1873 |
> |
units of kcal/mol. |
| 1874 |
|
|
| 1875 |
< |
The dipolar unified-atom force field ({\sc duff}) was developed to |
| 1876 |
< |
simulate lipid bilayers. These types of simulations require a model |
| 1877 |
< |
capable of forming bilayers, while still being sufficiently |
| 1878 |
< |
computationally efficient to allow large systems ($\sim$100's of |
| 1879 |
< |
phospholipids, $\sim$1000's of waters) to be simulated for long times |
| 1880 |
< |
($\sim$10's of nanoseconds). With this goal in mind, {\sc duff} has no |
| 1881 |
< |
point charges. Charge-neutral distributions are replaced with dipoles, |
| 1882 |
< |
while most atoms and groups of atoms are reduced to Lennard-Jones |
| 1883 |
< |
interaction sites. This simplification reduces the length scale of |
| 1884 |
< |
long range interactions from $\frac{1}{r}$ to $\frac{1}{r^3}$, |
| 1885 |
< |
removing the need for the computationally expensive Ewald |
| 1886 |
< |
sum. Instead, Verlet neighbor-lists and cutoff radii are used for the |
| 1887 |
< |
dipolar interactions, and, if desired, a reaction field may be added |
| 1888 |
< |
to mimic longer range interactions. |
| 1875 |
> |
\begin{lstlisting}[caption={[An example of a TorsionTypes block.] A |
| 1876 |
> |
simple example of a TorsionTypes block. Energy constants are given in |
| 1877 |
> |
kcal / mol, and when required by the form, $\delta$ angles are given |
| 1878 |
> |
in degrees.}, |
| 1879 |
> |
label={sch:TorsionTypes}] |
| 1880 |
> |
begin TorsionTypes |
| 1881 |
> |
//Cubic |
| 1882 |
> |
//Atom1 Atom2 Atom3 Atom4 Cubic k3 k2 k1 k0 |
| 1883 |
> |
CH2 CH2 CH2 CH2 Cubic 5.9602 -0.2568 -3.802 2.1586 |
| 1884 |
> |
CH2 CH CH CH2 Cubic 3.3254 -0.4215 -1.686 1.1661 |
| 1885 |
> |
//Trappe |
| 1886 |
> |
//Atom1 Atom2 Atom3 Atom4 Trappe c0 c1 c2 c3 |
| 1887 |
> |
CH3 CH2 CH2 SH Trappe 0.10507 -0.10342 0.03668 0.60874 |
| 1888 |
> |
//Charmm |
| 1889 |
> |
//Atom1 Atom2 Atom3 Atom4 Charmm Kchi n delta [Kchi n delta] |
| 1890 |
> |
CT CT CT C Charmm 0.156 3 0.00 |
| 1891 |
> |
OH CT CT OH Charmm 0.144 3 0.00 1.175 2 0 |
| 1892 |
> |
HC CT CM CM Charmm 1.150 1 0.00 0.38 3 180 |
| 1893 |
> |
//Quartic |
| 1894 |
> |
//Atom1 Atom2 Atom3 Atom4 Quartic k4 k3 k2 k1 k0 |
| 1895 |
> |
//Polynomial |
| 1896 |
> |
//Atom1 Atom2 Atom3 Atom4 Polynomial n Kn [m Km] |
| 1897 |
> |
S CH2 CH2 C Polynomial 0 2.218 1 2.905 2 -3.136 3 -0.7313 4 6.272 5 -7.528 |
| 1898 |
> |
end TorsionTypes |
| 1899 |
> |
\end{lstlisting} |
| 1900 |
|
|
| 1901 |
< |
As an example, lipid head-groups in {\sc duff} are represented as |
| 1902 |
< |
point dipole interaction sites. Placing a dipole at the head group's |
| 1903 |
< |
center of mass mimics the charge separation found in common |
| 1728 |
< |
phospholipid head groups such as phosphatidylcholine.\cite{Cevc87} |
| 1729 |
< |
Additionally, a large Lennard-Jones site is located at the |
| 1730 |
< |
pseudoatom's center of mass. The model is illustrated by the red atom |
| 1731 |
< |
in Fig.~\ref{fig:lipidModel}. The water model we use to |
| 1732 |
< |
complement the dipoles of the lipids is a |
| 1733 |
< |
reparameterization\cite{fennell04} of the soft sticky dipole (SSD) |
| 1734 |
< |
model of Ichiye |
| 1735 |
< |
\emph{et al.}\cite{liu96:new_model} |
| 1901 |
> |
Note that the parameters for a particular torsion type are the same |
| 1902 |
> |
for any torsional quartet of the same atomic types (in the same or |
| 1903 |
> |
reversed order). |
| 1904 |
|
|
| 1905 |
< |
\begin{figure} |
| 1738 |
< |
\centering |
| 1739 |
< |
\includegraphics[width=\linewidth]{lipidModel.pdf} |
| 1740 |
< |
\caption[A representation of a lipid model in {\sc duff}]{A |
| 1741 |
< |
representation of the lipid model. $\phi$ is the torsion angle, |
| 1742 |
< |
$\theta$ is the bend angle, and $\mu$ is the dipole moment of the head |
| 1743 |
< |
group.} |
| 1744 |
< |
\label{fig:lipidModel} |
| 1745 |
< |
\end{figure} |
| 1905 |
> |
\subsection{\label{section:ffInversion}The InversionTypes block} |
| 1906 |
|
|
| 1907 |
< |
A set of scalable parameters has been used to model the alkyl groups |
| 1908 |
< |
with Lennard-Jones sites. For this, parameters from the TraPPE force |
| 1909 |
< |
field of Siepmann \emph{et al.}\cite{Siepmann1998} have been |
| 1910 |
< |
utilized. TraPPE is a unified-atom representation of n-alkanes which |
| 1751 |
< |
is parametrized against phase equilibria using Gibbs ensemble Monte |
| 1752 |
< |
Carlo simulation techniques.\cite{Siepmann1998} One of the advantages |
| 1753 |
< |
of TraPPE is that it generalizes the types of atoms in an alkyl chain |
| 1754 |
< |
to keep the number of pseudoatoms to a minimum; thus, the parameters |
| 1755 |
< |
for a unified atom such as $\text{CH}_2$ do not change depending on |
| 1756 |
< |
what species are bonded to it. |
| 1907 |
> |
Inversion potentials are often added to force fields to enforce |
| 1908 |
> |
planarity around $sp^2$-hybridized carbons or to correct vibrational |
| 1909 |
> |
frequencies for umbrella-like vibrational modes for central atoms |
| 1910 |
> |
bonded to exactly three satellite groups. |
| 1911 |
|
|
| 1912 |
< |
As is required by TraPPE, {\sc duff} also constrains all bonds to be |
| 1913 |
< |
of fixed length. Typically, bond vibrations are the fastest motions in |
| 1914 |
< |
a molecular dynamic simulation. With these vibrations present, small |
| 1915 |
< |
time steps between force evaluations must be used to ensure adequate |
| 1916 |
< |
energy conservation in the bond degrees of freedom. By constraining |
| 1917 |
< |
the bond lengths, larger time steps may be used when integrating the |
| 1918 |
< |
equations of motion. A simulation using {\sc duff} is illustrated in |
| 1919 |
< |
Scheme \ref{sch:DUFF}. |
| 1912 |
> |
In OpenMD's version of an inversion, the central atom is entered {\it |
| 1913 |
> |
first} in each line in the {\tt InversionTypes} block. Note that |
| 1914 |
> |
this is quite different than how other codes treat Improper torsional |
| 1915 |
> |
potentials to mimic inversion behavior. In some other widely-used |
| 1916 |
> |
simulation packages, the central atom is treated as atom 3 in a |
| 1917 |
> |
standard torsion form: |
| 1918 |
> |
\begin{itemize} |
| 1919 |
> |
\item OpenMD: I - (J - K - L) (e.g. I is $sp^2$ hybridized carbon) |
| 1920 |
> |
\item AMBER: I - J - K - L (e.g. K is $sp^2$ hybridized carbon) |
| 1921 |
> |
\end{itemize} |
| 1922 |
|
|
| 1923 |
< |
\begin{lstlisting}[float,caption={[Invocation of {\sc duff}]A portion |
| 1924 |
< |
of a startup file showing a simulation utilizing {\sc |
| 1925 |
< |
duff}},label={sch:DUFF}] |
| 1926 |
< |
<OpenMD> |
| 1927 |
< |
<MetaData> |
| 1928 |
< |
#include "water.md" |
| 1929 |
< |
#include "lipid.md" |
| 1923 |
> |
The inversion angle itself is defined as: |
| 1924 |
> |
\begin{equation} |
| 1925 |
> |
\cos\omega_{i-jkl} = \left(\hat{\mathbf{r}}_{il} \times |
| 1926 |
> |
\hat{\mathbf{r}}_{ij}\right)\cdot\left( \hat{\mathbf{r}}_{il} \times |
| 1927 |
> |
\hat{\mathbf{r}}_{ik}\right) |
| 1928 |
> |
\end{equation} |
| 1929 |
> |
Here, $\hat{\mathbf{r}}_{\alpha\beta}$ are the set of unit bond |
| 1930 |
> |
vectors between the central atoms $i$, and the satellite atoms $j$, |
| 1931 |
> |
$k$, and $l$. Note that other definitions of inversion angles are |
| 1932 |
> |
possible, so users are encouraged to be particularly careful when |
| 1933 |
> |
converting other force field files for use with OpenMD. |
| 1934 |
|
|
| 1935 |
< |
component{ |
| 1936 |
< |
type = "simpleLipid_16"; |
| 1937 |
< |
nMol = 60; |
| 1938 |
< |
} |
| 1939 |
< |
|
| 1940 |
< |
component{ |
| 1941 |
< |
type = "SSD_water"; |
| 1942 |
< |
nMol = 1936; |
| 1943 |
< |
} |
| 1944 |
< |
|
| 1945 |
< |
forceField = "DUFF"; |
| 1946 |
< |
</MetaData> |
| 1947 |
< |
<Snapshot> // not shown in this scheme |
| 1948 |
< |
</Snapshot> |
| 1949 |
< |
</OpenMD> |
| 1935 |
> |
There are many common ways to create planarity or umbrella behavior in |
| 1936 |
> |
a potential energy function, and OpenMD implements a number of the |
| 1937 |
> |
more common functions: |
| 1938 |
> |
\begin{itemize} |
| 1939 |
> |
\item {\tt ImproperCosine}: |
| 1940 |
> |
\begin{equation*} |
| 1941 |
> |
V_{\text{torsion}}(\omega) = \sum_n \frac{K_n}{2} \left( 1 + cos(n |
| 1942 |
> |
\omega - \delta_n) \right), |
| 1943 |
> |
\end{equation*} |
| 1944 |
> |
\item {\tt AmberImproper}: |
| 1945 |
> |
\begin{equation*} |
| 1946 |
> |
V_{\text{torsion}}(\omega) = \frac{v}{2} (1 - \cos\left(2 \left(\omega - \omega_0\right)\right), |
| 1947 |
> |
\end{equation*} |
| 1948 |
> |
\item {\tt Harmonic}: |
| 1949 |
> |
\begin{equation*} |
| 1950 |
> |
V_{\text{torsion}}(\omega) = \frac{d}{2} \left(\omega - \omega_0\right). |
| 1951 |
> |
\end{equation*} |
| 1952 |
> |
\end{itemize} |
| 1953 |
> |
\begin{lstlisting}[caption={[An example of an InversionTypes block.] A |
| 1954 |
> |
simple example of a InversionTypes block. Angles ($\delta_n$ and |
| 1955 |
> |
$\omega_0$) angles are given in degrees, while energy parameters ($v, |
| 1956 |
> |
K_n$) are given in kcal / mol. The Harmonic Inversion type has a |
| 1957 |
> |
force constant that must be given in kcal/mol/degrees$^2$.}, |
| 1958 |
> |
label={sch:InversionTypes}] |
| 1959 |
> |
begin InversionTypes |
| 1960 |
> |
//Harmonic |
| 1961 |
> |
//Atom1 Atom2 Atom3 Atom4 Harmonic d(kcal/mol/deg^2) omega0 |
| 1962 |
> |
RCHar3 X X X Harmonic 1.21876e-2 180.0 |
| 1963 |
> |
//AmberImproper |
| 1964 |
> |
//Atom1 Atom2 Atom3 Atom4 AmberImproper v(kcal/mol) |
| 1965 |
> |
C CT N O AmberImproper 10.500000 |
| 1966 |
> |
CA CA CA CT AmberImproper 1.100000 |
| 1967 |
> |
//ImproperCosine |
| 1968 |
> |
//Atom1 Atom2 Atom3 Atom4 ImproperCosine Kn n delta_n [Kn n delta_n] |
| 1969 |
> |
end InversionTypes |
| 1970 |
|
\end{lstlisting} |
| 1971 |
|
|
| 1972 |
+ |
\section{\label{section::ffLongRange}Long Range Interactions} |
| 1973 |
|
|
| 1974 |
+ |
Calculating the long-range (non-bonded) potential involves a sum over |
| 1975 |
+ |
all pairs of atoms (except for those atoms which are involved in a |
| 1976 |
+ |
bond, bend, or torsion with each other). Many of these interactions |
| 1977 |
+ |
can be inferred from the AtomTypes, |
| 1978 |
|
|
| 1979 |
< |
The cross potential between molecules $I$ and $J$, |
| 1980 |
< |
$V^{IJ}_{\text{Cross}}$, is as follows: |
| 1796 |
< |
\begin{equation} |
| 1797 |
< |
V^{IJ}_{\text{Cross}} = |
| 1798 |
< |
\sum_{i \in I} \sum_{j \in J} |
| 1799 |
< |
\biggl[ V_{\text{LJ}}(r_{ij}) + V_{\text{dipole}} |
| 1800 |
< |
(\mathbf{r}_{ij},\boldsymbol{\Omega}_{i},\boldsymbol{\Omega}_{j}) |
| 1801 |
< |
+ V_{\text{sticky}} |
| 1802 |
< |
(\mathbf{r}_{ij},\boldsymbol{\Omega}_{i},\boldsymbol{\Omega}_{j}) |
| 1803 |
< |
\biggr], |
| 1804 |
< |
\label{eq:crossPotentail} |
| 1805 |
< |
\end{equation} |
| 1806 |
< |
where $V_{\text{LJ}}$ is the Lennard Jones potential, |
| 1807 |
< |
$V_{\text{dipole}}$ is the dipole dipole potential, and |
| 1808 |
< |
$V_{\text{sticky}}$ is the sticky potential defined by the SSD model |
| 1809 |
< |
(Sec.~\ref{section:SSD}). Note that not all atom types include all |
| 1810 |
< |
interactions. |
| 1979 |
> |
\subsection{\label{section:ffNBinteraction}The NonBondedInteractions |
| 1980 |
> |
block} |
| 1981 |
|
|
| 1982 |
+ |
The user might want like to specify explicit or special interactions |
| 1983 |
+ |
that override the default non-bodned interactions that are inferred |
| 1984 |
+ |
from the AtomTypes. To do this, OpenMD implements a |
| 1985 |
+ |
NonBondedInteractions block to allow the user to specify the following |
| 1986 |
+ |
(pair-wise) non-bonded interactions: |
| 1987 |
|
|
| 1988 |
< |
\section{\label{section:WATER}The {\sc water} Force Field} |
| 1988 |
> |
\begin{itemize} |
| 1989 |
> |
\item {\tt LennardJones}: |
| 1990 |
> |
\begin{equation*} |
| 1991 |
> |
V_{\text{NB}}(r) = 4 \epsilon_{ij} \left( |
| 1992 |
> |
\left(\frac{\sigma_{ij}}{r} \right)^{12} - |
| 1993 |
> |
\left(\frac{\sigma_{ij}}{r} \right)^{6} \right), |
| 1994 |
> |
\end{equation*} |
| 1995 |
> |
\item {\tt ShiftedMorse}: |
| 1996 |
> |
\begin{equation*} |
| 1997 |
> |
V_{\text{NB}}(r) = D_{ij} \left( e^{-2 \beta_{ij} (r - |
| 1998 |
> |
r^0)} - 2 e^{- \beta_{ij} (r - |
| 1999 |
> |
r^0)} \right), |
| 2000 |
> |
\end{equation*} |
| 2001 |
> |
\item {\tt RepulsiveMorse}: |
| 2002 |
> |
\begin{equation*} |
| 2003 |
> |
V_{\text{NB}}(r) = D_{ij} \left( e^{-2 \beta_{ij} (r - |
| 2004 |
> |
r^0)} \right), |
| 2005 |
> |
\end{equation*} |
| 2006 |
> |
\item {\tt RepulsivePower}: |
| 2007 |
> |
\begin{equation*} |
| 2008 |
> |
V_{\text{NB}}(r) = \epsilon_{ij} |
| 2009 |
> |
\left(\frac{\sigma_{ij}}{r} \right)^{n_{ij}}. |
| 2010 |
> |
\end{equation*} |
| 2011 |
> |
\end{itemize} |
| 2012 |
|
|
| 2013 |
< |
In addition to the {\sc duff} force field's solvent description, a |
| 2014 |
< |
separate {\sc water} force field has been included for simulating most |
| 2015 |
< |
of the common rigid-body water models. This force field includes the |
| 2016 |
< |
simple and point-dipolar models (SSD, SSD1, SSD/E, SSD/RF, and DPD |
| 2017 |
< |
water), as well as the common charge-based models (SPC, SPC/E, TIP3P, |
| 2018 |
< |
TIP4P, and |
| 2019 |
< |
TIP5P).\cite{liu96:new_model,Ichiye03,fennell04,Marrink01,Berendsen81,Berendsen87,Jorgensen83,Mahoney00} |
| 1822 |
< |
In order to handle these models, charge-charge interactions were |
| 1823 |
< |
included in the force-loop: |
| 1824 |
< |
\begin{equation} |
| 1825 |
< |
V_{\text{charge}}(r_{ij}) = \sum_{ij}\frac{q_iq_je^2}{r_{ij}}, |
| 1826 |
< |
\end{equation} |
| 1827 |
< |
where $q$ represents the charge on particle $i$ or $j$, and $e$ is the |
| 1828 |
< |
charge of an electron in Coulombs. The charge-charge interaction |
| 1829 |
< |
support is rudimentary in the current version of {\sc OpenMD}. As with |
| 1830 |
< |
the other pair interactions, charges can be simulated with a pure |
| 1831 |
< |
cutoff or a reaction field. The various methods for performing the |
| 1832 |
< |
Ewald summation have not yet been included. The {\sc water} force |
| 1833 |
< |
field can be easily expanded through modification of the {\sc water} |
| 1834 |
< |
force field file ({\tt WATER.frc}). By adding atom types and inserting |
| 1835 |
< |
the appropriate parameters, it is possible to extend the force field |
| 1836 |
< |
to handle rigid molecules other than water. |
| 2013 |
> |
\begin{lstlisting}[caption={[An example of a NonBondedInteractions block.] A |
| 2014 |
> |
simple example of a NonBondedInteractions block. Distances ($\sigma, |
| 2015 |
> |
r_0$) are given in \AA, while energies ($\epsilon, D0$) are in |
| 2016 |
> |
kcal/mol. The Morse potentials have an additional parameter $\beta_0$ |
| 2017 |
> |
which is in units of \AA$^{-1}$.}, |
| 2018 |
> |
label={sch:InversionTypes}] |
| 2019 |
> |
begin NonBondedInteractions |
| 2020 |
|
|
| 2021 |
+ |
//Lennard-Jones |
| 2022 |
+ |
//Atom1 Atom2 LennardJones sigma epsilon |
| 2023 |
+ |
Au CH3 LennardJones 3.54 0.2146 |
| 2024 |
+ |
Au CH2 LennardJones 3.54 0.1749 |
| 2025 |
+ |
Au CH LennardJones 3.54 0.1749 |
| 2026 |
+ |
Au S LennardJones 2.40 8.465 |
| 2027 |
|
|
| 2028 |
< |
\section{\label{section:sc}The Sutton-Chen Force Field} |
| 2028 |
> |
//Shifted Morse |
| 2029 |
> |
//Atom1 Atom2 ShiftedMorse r0 D0 beta0 |
| 2030 |
> |
Au O_SPCE ShiftedMorse 3.70 0.0424 0.769 |
| 2031 |
|
|
| 2032 |
+ |
//Repulsive Morse |
| 2033 |
+ |
//Atom1 Atom2 RepulsiveMorse r0 D0 beta0 |
| 2034 |
+ |
Au H_SPCE RepulsiveMorse -1.00 0.00850 0.769 |
| 2035 |
|
|
| 2036 |
< |
\section{\label{section:clay}The CLAY force field} |
| 2037 |
< |
|
| 2038 |
< |
The {\sc clay} force field is based on an ionic (nonbonded) |
| 2039 |
< |
description of the metal-oxygen interactions associated with hydrated |
| 2040 |
< |
phases. All atoms are represented as point charges and are allowed |
| 2041 |
< |
complete translational freedom. Metal-oxygen interactions are based on |
| 1848 |
< |
a simple Lennard-Jones potential combined with electrostatics. The |
| 1849 |
< |
empirical parameters were optimized by Cygan {\it et |
| 1850 |
< |
al.}\cite{Cygan04} on the basis of known mineral structures, and |
| 1851 |
< |
partial atomic charges were derived from periodic DFT quantum chemical |
| 1852 |
< |
calculations of simple oxide, hydroxide, and oxyhydroxide model |
| 1853 |
< |
compounds with well-defined structures. |
| 2036 |
> |
//Repulsive Power |
| 2037 |
> |
//Atom1 Atom2 RepulsivePower sigma epsilon n |
| 2038 |
> |
Au ON RepulsivePower 3.47005 0.186208 11 |
| 2039 |
> |
Au NO RepulsivePower 3.53955 0.168629 11 |
| 2040 |
> |
end NonBondedInteractions |
| 2041 |
> |
\end{lstlisting} |
| 2042 |
|
|
| 1855 |
– |
|
| 2043 |
|
\section{\label{section:electrostatics}Electrostatics} |
| 2044 |
|
|
| 2045 |
|
To aid in performing simulations in more traditional force fields, we |
| 2169 |
|
\end{equation} |
| 2170 |
|
the shifted potential (eq. (\ref{eq:SPPot})) becomes |
| 2171 |
|
\begin{equation} |
| 2172 |
< |
V_{\textrm{DSP}}(r) = q_iq_j\left(\frac{\textrm{erfc}\left(\alpha r\right)}{r}-\ |
| 1986 |
< |
frac{\textrm{erfc}\left(\alpha R_\textrm{c}\right)}{R_\textrm{c}}\right) \quad r |
| 2172 |
> |
V_{\textrm{DSP}}(r) = q_iq_j\left(\frac{\textrm{erfc}\left(\alpha r\right)}{r}-\frac{\textrm{erfc}\left(\alpha R_\textrm{c}\right)}{R_\textrm{c}}\right) \quad r |
| 2173 |
|
\leqslant R_\textrm{c}, |
| 2174 |
|
\label{eq:DSPPot} |
| 2175 |
|
\end{equation} |
| 2214 |
|
{\sc OpenMD} is the DSF (Eq. \ref{eq:DSFPot}) with a damping parameter |
| 2215 |
|
($\alpha$) that is set algorithmically from the cutoff radius. |
| 2216 |
|
|
| 2217 |
+ |
|
| 2218 |
+ |
\section{\label{section:cutoffGroups}Switching Functions} |
| 2219 |
+ |
|
| 2220 |
+ |
If done poorly, calculating the the long-range interactions for $N$ |
| 2221 |
+ |
atoms would involve $N(N-1)/2$ evaluations of atomic distances. To |
| 2222 |
+ |
reduce the number of distance evaluations between pairs of atoms, {\sc |
| 2223 |
+ |
OpenMD} allows the use of switched cutoffs with Verlet neighbor |
| 2224 |
+ |
lists.\cite{Allen87} Neutral groups which contain charges will exhibit |
| 2225 |
+ |
pathological forces unless the cutoff is applied to the neutral groups |
| 2226 |
+ |
evenly instead of to the individual atoms.\cite{leach01:mm} {\sc |
| 2227 |
+ |
OpenMD} allows users to specify cutoff groups which may contain an |
| 2228 |
+ |
arbitrary number of atoms in the molecule. Atoms in a cutoff group |
| 2229 |
+ |
are treated as a single unit for the evaluation of the switching |
| 2230 |
+ |
function: |
| 2231 |
+ |
\begin{equation} |
| 2232 |
+ |
V_{\mathrm{long-range}} = \sum_{a} \sum_{b>a} s(r_{ab}) \sum_{i \in a} \sum_{j \in b} V_{ij}(r_{ij}), |
| 2233 |
+ |
\end{equation} |
| 2234 |
+ |
where $r_{ab}$ is the distance between the centers of mass of the two |
| 2235 |
+ |
cutoff groups ($a$ and $b$). |
| 2236 |
+ |
|
| 2237 |
+ |
The sums over $a$ and $b$ are over the cutoff groups that are present |
| 2238 |
+ |
in the simulation. Atoms which are not explicitly defined as members |
| 2239 |
+ |
of a {\tt cutoffGroup} are treated as a group consisting of only one |
| 2240 |
+ |
atom. The switching function, $s(r)$ is the standard cubic switching |
| 2241 |
+ |
function, |
| 2242 |
+ |
\begin{equation} |
| 2243 |
+ |
S(r) = |
| 2244 |
+ |
\begin{cases} |
| 2245 |
+ |
1 & \text{if $r \le r_{\text{sw}}$},\\ |
| 2246 |
+ |
\frac{(r_{\text{cut}} + 2r - 3r_{\text{sw}})(r_{\text{cut}} - r)^2} |
| 2247 |
+ |
{(r_{\text{cut}} - r_{\text{sw}})^3} |
| 2248 |
+ |
& \text{if $r_{\text{sw}} < r \le r_{\text{cut}}$}, \\ |
| 2249 |
+ |
0 & \text{if $r > r_{\text{cut}}$.} |
| 2250 |
+ |
\end{cases} |
| 2251 |
+ |
\label{eq:dipoleSwitching} |
| 2252 |
+ |
\end{equation} |
| 2253 |
+ |
Here, $r_{\text{sw}}$ is the {\tt switchingRadius}, or the distance |
| 2254 |
+ |
beyond which interactions are reduced, and $r_{\text{cut}}$ is the |
| 2255 |
+ |
{\tt cutoffRadius}, or the distance at which interactions are |
| 2256 |
+ |
truncated. |
| 2257 |
+ |
|
| 2258 |
+ |
Users of {\sc OpenMD} do not need to specify the {\tt cutoffRadius} or |
| 2259 |
+ |
{\tt switchingRadius}. In simulations containing only Lennard-Jones |
| 2260 |
+ |
atoms, the cutoff radius has a default value of $2.5\sigma_{ii}$, |
| 2261 |
+ |
where $\sigma_{ii}$ is the largest Lennard-Jones length parameter |
| 2262 |
+ |
present in the simulation. In simulations containing charged or |
| 2263 |
+ |
dipolar atoms, the default cutoff radius is $15 \mbox{\AA}$. |
| 2264 |
+ |
|
| 2265 |
+ |
The {\tt switchingRadius} is set to a default value of 95\% of the |
| 2266 |
+ |
{\tt cutoffRadius}. In the special case of a simulation containing |
| 2267 |
+ |
{\it only} Lennard-Jones atoms, the default switching radius takes the |
| 2268 |
+ |
same value as the cutoff radius, and {\sc OpenMD} will use a shifted |
| 2269 |
+ |
potential to remove discontinuities in the potential at the cutoff. |
| 2270 |
+ |
Both radii may be specified in the meta-data file. |
| 2271 |
+ |
|
| 2272 |
+ |
|
| 2273 |
+ |
\section{\label{section:WATER}The {\sc water} Force Field} |
| 2274 |
+ |
|
| 2275 |
+ |
In addition to the {\sc duff} force field's solvent description, a |
| 2276 |
+ |
separate {\sc water} force field has been included for simulating most |
| 2277 |
+ |
of the common rigid-body water models. This force field includes the |
| 2278 |
+ |
simple and point-dipolar models (SSD, SSD1, SSD/E, SSD/RF, and DPD |
| 2279 |
+ |
water), as well as the common charge-based models (SPC, SPC/E, TIP3P, |
| 2280 |
+ |
TIP4P, and |
| 2281 |
+ |
TIP5P).\cite{liu96:new_model,Ichiye03,fennell04,Marrink01,Berendsen81,Berendsen87,Jorgensen83,Mahoney00} |
| 2282 |
+ |
In order to handle these models, charge-charge interactions were |
| 2283 |
+ |
included in the force-loop: |
| 2284 |
+ |
\begin{equation} |
| 2285 |
+ |
V_{\text{charge}}(r_{ij}) = \sum_{ij}\frac{q_iq_je^2}{r_{ij}}, |
| 2286 |
+ |
\end{equation} |
| 2287 |
+ |
where $q$ represents the charge on particle $i$ or $j$, and $e$ is the |
| 2288 |
+ |
charge of an electron in Coulombs. The charge-charge interaction |
| 2289 |
+ |
support is rudimentary in the current version of {\sc OpenMD}. As with |
| 2290 |
+ |
the other pair interactions, charges can be simulated with a pure |
| 2291 |
+ |
cutoff or a reaction field. The various methods for performing the |
| 2292 |
+ |
Ewald summation have not yet been included. The {\sc water} force |
| 2293 |
+ |
field can be easily expanded through modification of the {\sc water} |
| 2294 |
+ |
force field file ({\tt WATER.frc}). By adding atom types and inserting |
| 2295 |
+ |
the appropriate parameters, it is possible to extend the force field |
| 2296 |
+ |
to handle rigid molecules other than water. |
| 2297 |
+ |
|
| 2298 |
+ |
|
| 2299 |
+ |
|
| 2300 |
|
\section{\label{section:pbc}Periodic Boundary Conditions} |
| 2301 |
|
|
| 2302 |
|
\newcommand{\roundme}{\operatorname{round}} |
