| 9 |
|
\usepackage{longtable} |
| 10 |
|
\pagestyle{plain} |
| 11 |
|
\pagenumbering{arabic} |
| 12 |
+ |
\usepackage{floatrow} |
| 13 |
+ |
\usepackage[margin=0.5cm,font=small,format=hang]{caption} |
| 14 |
+ |
|
| 15 |
|
\oddsidemargin 0.0cm |
| 16 |
|
\evensidemargin 0.0cm |
| 17 |
|
\topmargin -21pt |
| 20 |
|
\textwidth 6.5in |
| 21 |
|
\brokenpenalty=10000 |
| 22 |
|
\renewcommand{\baselinestretch}{1.2} |
| 23 |
+ |
\usepackage[square, comma, sort&compress]{natbib} |
| 24 |
+ |
\bibpunct{[}{]}{,}{n}{}{;} |
| 25 |
|
|
| 26 |
+ |
\DeclareFloatFont{tiny}{\scriptsize}% "scriptsize" is defined by floatrow, "tiny" not |
| 27 |
+ |
\floatsetup[table]{font=tiny} |
| 28 |
+ |
|
| 29 |
+ |
|
| 30 |
|
%\renewcommand\citemid{\ } % no comma in optional reference note |
| 31 |
< |
\lstset{language=C,frame=TB,basicstyle=\tiny,basicstyle=\ttfamily, % |
| 31 |
> |
\lstset{language=C,frame=TB,basicstyle=\footnotesize\ttfamily, % |
| 32 |
|
xleftmargin=0.25in, xrightmargin=0.25in,captionpos=b, % |
| 33 |
|
abovecaptionskip=0.5cm, belowcaptionskip=0.5cm, escapeinside={~}{~}} |
| 34 |
< |
\renewcommand{\lstlistingname}{Scheme} |
| 34 |
> |
\renewcommand{\lstlistingname}{Example} |
| 35 |
|
|
| 36 |
+ |
\lstnewenvironment{code}[1][]% |
| 37 |
+ |
{\noindent\minipage{\linewidth}\vspace{0.5\baselineskip} |
| 38 |
+ |
\lstset{language=C,basicstyle=\footnotesize\ttfamily,% |
| 39 |
+ |
captionpos=b,aboveskip=0.5cm,belowskip=0.5cm,abovecaptionskip=0.5cm,% |
| 40 |
+ |
belowcaptionskip=0.5cm,% |
| 41 |
+ |
escapeinside={~}{~},frame=single,#1}} |
| 42 |
+ |
{\endminipage} |
| 43 |
+ |
|
| 44 |
+ |
|
| 45 |
+ |
|
| 46 |
|
\begin{document} |
| 47 |
|
|
| 48 |
|
\newcolumntype{A}{p{1.5in}} |
| 57 |
|
\newcolumntype{H}{p{0.75in}} |
| 58 |
|
\newcolumntype{I}{p{5in}} |
| 59 |
|
|
| 60 |
+ |
\newcolumntype{J}{p{1.5in}} |
| 61 |
+ |
\newcolumntype{K}{p{1.2in}} |
| 62 |
+ |
\newcolumntype{L}{p{1.5in}} |
| 63 |
+ |
\newcolumntype{M}{p{1.55in}} |
| 64 |
|
|
| 42 |
– |
\title{{\sc OpenMD}: Molecular Dynamics in the Open} |
| 65 |
|
|
| 66 |
< |
\author{Kelsey M. Stocker, Shenyu Kuang, Charles F. Vardeman II, \\ |
| 67 |
< |
Teng Lin, Christopher J. Fennell, Xiuquan Sun, \\ |
| 68 |
< |
Chunlei Li, Kyle Daily, Yang Zheng, Matthew A. Meineke, and \\ |
| 69 |
< |
J. Daniel Gezelter \\ |
| 66 |
> |
\title{{\sc OpenMD-2.3}: Molecular Dynamics in the Open} |
| 67 |
> |
|
| 68 |
> |
\author{Madan Lamichhane, Patrick Louden, Joseph Michalka, Suzanne |
| 69 |
> |
Niedhart, \\ |
| 70 |
> |
Teng Lin, Charles F. Vardeman II, Christopher J. Fennell, Matthew |
| 71 |
> |
A. Meineke, \\ Shenyu Kuang, |
| 72 |
> |
Kelsey Stocker, James Marr, \\ Xiuquan Sun, |
| 73 |
> |
Chunlei Li, Kyle Daily, Yang Zheng, \\ and \\ |
| 74 |
> |
J. Daniel Gezelter \\ \\ |
| 75 |
|
Department of Chemistry and Biochemistry\\ |
| 76 |
|
University of Notre Dame\\ |
| 77 |
|
Notre Dame, Indiana 46556} |
| 78 |
|
|
| 79 |
|
\maketitle |
| 80 |
|
|
| 81 |
< |
\section*{Preface} |
| 82 |
< |
{\sc OpenMD} is an open source molecular dynamics engine which is capable of |
| 83 |
< |
efficiently simulating liquids, proteins, nanoparticles, interfaces, |
| 84 |
< |
and other complex systems using atom types with orientational degrees |
| 85 |
< |
of freedom (e.g. ``sticky'' atoms, point dipoles, and coarse-grained |
| 86 |
< |
assemblies). Proteins, zeolites, lipids, transition metals (bulk, flat |
| 87 |
< |
interfaces, and nanoparticles) have all been simulated using force |
| 88 |
< |
fields included with the code. {\sc OpenMD} works on parallel computers |
| 89 |
< |
using the Message Passing Interface (MPI), and comes with a number of |
| 81 |
> |
\section*{Preface} |
| 82 |
> |
|
| 83 |
> |
{\sc OpenMD} is an open source molecular dynamics engine which is |
| 84 |
> |
capable of efficiently simulating liquids, proteins, nanoparticles, |
| 85 |
> |
interfaces, and other complex systems using atom types with |
| 86 |
> |
orientational degrees of freedom (e.g. ``sticky'' atoms, point |
| 87 |
> |
multipoles, and coarse-grained assemblies). It is a test-bed for new |
| 88 |
> |
molecular simulation methodology, but is also efficient and easy to |
| 89 |
> |
use. Liquids, proteins, zeolites, lipids, inorganic nanomaterials, |
| 90 |
> |
transition metals (bulk, flat interfaces, and nanoparticles) and a |
| 91 |
> |
wide array of other systems have all been simulated using this |
| 92 |
> |
code. {\sc OpenMD} works on parallel computers using the Message |
| 93 |
> |
Passing Interface (MPI), and comes with a number of trajectory |
| 94 |
|
analysis and utility programs that are easy to use and modify. An |
| 95 |
|
OpenMD simulation is specified using a very simple meta-data language |
| 96 |
|
that is easy to learn. |
| 112 |
|
which is freely available to the entire scientific community. Few, |
| 113 |
|
however, are capable of efficiently integrating the equations of |
| 114 |
|
motion for atom types with orientational degrees of freedom |
| 115 |
< |
(e.g. point dipoles, and ``sticky'' atoms). And only one of the |
| 115 |
> |
(e.g. point multipoles, and ``sticky'' atoms). And only one of the |
| 116 |
|
programs referenced can handle transition metal force fields like the |
| 117 |
|
Embedded Atom Method ({\sc eam}). The direction our research program |
| 118 |
|
has taken us now involves the use of atoms with orientational degrees |
| 123 |
|
|
| 124 |
|
This document communicates the algorithmic details of our program, |
| 125 |
|
{\sc OpenMD}. We have structured this document to first discuss the |
| 126 |
< |
underlying concepts in this simulation package (Sec. |
| 126 |
> |
underlying concepts in this simulation package (Chapter |
| 127 |
|
\ref{section:IOfiles}). The empirical energy functions implemented |
| 128 |
< |
are discussed in Sec.~\ref{section:empiricalEnergy}. |
| 129 |
< |
Sec.~\ref{section:mechanics} describes the various Molecular Dynamics |
| 128 |
> |
are discussed in Chapter~\ref{chapter:forceFields}. |
| 129 |
> |
Section~\ref{section:mechanics} describes the various Molecular Dynamics |
| 130 |
|
algorithms {\sc OpenMD} implements in the integration of Hamilton's |
| 131 |
|
equations of motion. Program design considerations for parallel |
| 132 |
|
computing are presented in Sec.~\ref{section:parallelization}. |
| 135 |
|
\chapter{\label{section:IOfiles}Concepts \& Files} |
| 136 |
|
|
| 137 |
|
A simulation in {\sc OpenMD} is built using a few fundamental |
| 138 |
< |
conceptual building blocks most of which are chemically intuitive. |
| 138 |
> |
conceptual building blocks, most of which are chemically intuitive. |
| 139 |
|
The basic unit of a simulation is an {\tt atom}. The parameters |
| 140 |
|
describing an {\tt atom} have been generalized to make it as flexible |
| 141 |
|
as possible; this means that in addition to translational degrees of |
| 142 |
< |
freedom, {\tt Atoms} may also have {\it orientational} degrees of freedom. |
| 142 |
> |
freedom, {\tt atoms} may also have {\it orientational} degrees of |
| 143 |
> |
freedom. |
| 144 |
|
|
| 145 |
|
The fundamental (static) properties of {\tt atoms} are defined by the |
| 146 |
|
{\tt forceField} chosen for the simulation. The atomic properties |
| 165 |
|
leave an interaction region. |
| 166 |
|
|
| 167 |
|
{\tt Atoms} may also be grouped in more traditional ways into {\tt |
| 168 |
< |
bonds}, {\tt bends}, and {\tt torsions}. These groupings allow the |
| 169 |
< |
correct choice of interaction parameters for short-range interactions |
| 170 |
< |
to be chosen from the definitions in the {\tt forceField}. |
| 168 |
> |
bonds}, {\tt bends}, {\tt torsions}, and {\tt inversions}. These |
| 169 |
> |
groupings allow the correct choice of interaction parameters for |
| 170 |
> |
short-range interactions to be chosen from the definitions in the {\tt |
| 171 |
> |
forceField}. |
| 172 |
|
|
| 173 |
|
All of these groups of {\tt atoms} are brought together in the {\tt |
| 174 |
|
molecule}, which is the fundamental structure for setting up and {\sc |
| 203 |
|
$<$Snapshot$>$} block. Both the {\tt $<$MetaData$>$} and {\tt $<$Snapshot$>$} |
| 204 |
|
formats are described in the following sections. |
| 205 |
|
|
| 206 |
< |
\begin{lstlisting}[float,caption={[The structure of an {\sc OpenMD} file] |
| 206 |
> |
\begin{code}[caption={[The structure of an {\sc OpenMD} file] |
| 207 |
|
The basic structure of an {\sc OpenMD} file contains HTML-like tags to |
| 208 |
|
define simulation meta-data and subsequent instantaneous configuration |
| 209 |
< |
information. A well-formed {\sc OpenMD} file must contain one $<$MetaData$>$ |
| 210 |
< |
block and {\it at least one} $<$Snapshot$>$ block. Each |
| 211 |
< |
$<$Snapshot$>$ is further divided into $<$FrameData$>$ and |
| 212 |
< |
$<$StuntDoubles$>$ sections.}, |
| 180 |
< |
label=sch:mdFormat] |
| 209 |
> |
information. A well-formed {\sc OpenMD} file must contain one {\tt <MetaData>} |
| 210 |
> |
block and {\it at least one} {\tt <Snapshot>} block. Each |
| 211 |
> |
{\tt <Snapshot>} is further divided into {\tt <FrameData>} and |
| 212 |
> |
{\tt <StuntDoubles>} sections.},label={sch:mdFormat}] |
| 213 |
|
<OpenMD> |
| 214 |
|
<MetaData> |
| 215 |
|
// see section ~\ref{sec:miscConcepts}~ for details on the formatting |
| 235 |
|
<Snapshot> // Further information on <Snapshot> blocks |
| 236 |
|
</Snapshot> // can be found in section ~\ref{section:coordFiles}~. |
| 237 |
|
</OpenMD> |
| 238 |
< |
\end{lstlisting} |
| 238 |
> |
\end{code} |
| 239 |
|
|
| 240 |
|
|
| 241 |
|
\section{OpenMD Files and $<$MetaData$>$ blocks} |
| 242 |
|
|
| 243 |
< |
{\sc OpenMD} uses a HTML-like syntax to separate {\tt $<$MetaData$>$} and |
| 243 |
> |
{\sc OpenMD} uses HTML-like delimiters to separate {\tt $<$MetaData$>$} and |
| 244 |
|
{\tt $<$Snapshot$>$} blocks. A C-based syntax is used to parse the {\tt |
| 245 |
|
$<$MetaData$>$} blocks at run time. These blocks allow the user to |
| 246 |
|
completely describe the system they wish to simulate, as well as |
| 248 |
|
files are typically denoted with the extension {\tt .md} (which can |
| 249 |
|
stand for Meta-Data or Molecular Dynamics or Molecule Definition |
| 250 |
|
depending on the user's mood). An overview of an {\sc OpenMD} file is |
| 251 |
< |
shown in Scheme~\ref{sch:mdFormat} and example file is shown in |
| 252 |
< |
Scheme~\ref{sch:mdExample}. |
| 251 |
> |
shown in Example~\ref{sch:mdFormat} and example file is shown in |
| 252 |
> |
Example~\ref{sch:mdExample}. |
| 253 |
|
|
| 254 |
< |
\begin{lstlisting}[float,caption={[An example of a complete OpenMD |
| 254 |
> |
\begin{code}[caption={[An example of a complete OpenMD |
| 255 |
|
file] An example showing a complete OpenMD file.}, |
| 256 |
|
label={sch:mdExample}] |
| 257 |
|
<OpenMD> |
| 290 |
|
</StuntDoubles> |
| 291 |
|
</Snapshot> |
| 292 |
|
</OpenMD> |
| 293 |
< |
\end{lstlisting} |
| 293 |
> |
\end{code} |
| 294 |
|
|
| 295 |
< |
Within the {\tt $<$MetaData$>$} block it is necessary to provide a |
| 295 |
> |
In the {\tt $<$MetaData$>$} block, it is necessary to provide a |
| 296 |
|
complete description of the molecule before it is actually placed in |
| 297 |
< |
the simulation. {\sc OpenMD}'s meta-data syntax was originally |
| 298 |
< |
developed with this goal in mind, and allows for the use of {\it |
| 299 |
< |
include files} to specify all atoms in a molecular prototype, as well |
| 300 |
< |
as any bonds, bends, or torsions. Include files allow the user to |
| 301 |
< |
describe a molecular prototype once, then simply include it into each |
| 302 |
< |
simulation containing that molecule. Returning to the example in |
| 303 |
< |
Scheme~\ref{sch:mdExample}, the include file's contents would be |
| 304 |
< |
Scheme~\ref{sch:mdIncludeExample}, and the new {\sc OpenMD} file would |
| 273 |
< |
become Scheme~\ref{sch:mdExPrime}. |
| 297 |
> |
the simulation. {\sc OpenMD}'s meta-data syntax allows for the use of |
| 298 |
> |
{\it include files} to specify all atoms in a molecular prototype, as |
| 299 |
> |
well as any bonds, bends, torsions, or other structural groupings of |
| 300 |
> |
atoms. Include files allow the user to describe a molecular prototype |
| 301 |
> |
once, then simply include it into each simulation containing that |
| 302 |
> |
molecule. Returning to the example in Scheme~\ref{sch:mdExample}, the |
| 303 |
> |
include file's contents would be Scheme~\ref{sch:mdIncludeExample}, |
| 304 |
> |
and the new {\sc OpenMD} file would become Scheme~\ref{sch:mdExPrime}. |
| 305 |
|
|
| 306 |
< |
\begin{lstlisting}[float,caption={An example molecule definition in an |
| 306 |
> |
\begin{code}[caption={An example molecule definition in an |
| 307 |
|
include file.},label={sch:mdIncludeExample}] |
| 308 |
|
molecule{ |
| 309 |
|
name = "Ar"; |
| 312 |
|
position( 0.0, 0.0, 0.0 ); |
| 313 |
|
} |
| 314 |
|
} |
| 315 |
< |
\end{lstlisting} |
| 315 |
> |
\end{code} |
| 316 |
|
|
| 317 |
< |
\begin{lstlisting}[float,caption={Revised OpenMD input file |
| 317 |
> |
\begin{code}[caption={Revised OpenMD input file |
| 318 |
|
example.},label={sch:mdExPrime}] |
| 319 |
|
<OpenMD> |
| 320 |
|
<MetaData> |
| 347 |
|
</StuntDoubles> |
| 348 |
|
</Snapshot> |
| 349 |
|
</OpenMD> |
| 350 |
< |
\end{lstlisting} |
| 350 |
> |
\end{code} |
| 351 |
|
|
| 352 |
|
\section{\label{section:atomsMolecules}Atoms, Molecules, and other |
| 353 |
|
ways of grouping atoms} |
| 354 |
|
|
| 355 |
< |
As mentioned above, the fundamental unit for an {\sc OpenMD} simulation |
| 356 |
< |
is the {\tt atom}. Atoms can be collected into secondary structures |
| 357 |
< |
such as {\tt rigidBodies}, {\tt cutoffGroups}, or {\tt molecules}. The |
| 358 |
< |
{\tt molecule} is a way for {\sc OpenMD} to keep track of the atoms in |
| 359 |
< |
a simulation in logical manner. Molecular units store the identities |
| 360 |
< |
of all the atoms and rigid bodies associated with themselves, and they |
| 361 |
< |
are responsible for the evaluation of their own internal interactions |
| 362 |
< |
(\emph{i.e.}~bonds, bends, and torsions). Scheme |
| 363 |
< |
\ref{sch:mdIncludeExample} shows how one creates a molecule in an |
| 364 |
< |
included meta-data file. The positions of the atoms given in the |
| 365 |
< |
declaration are relative to the origin of the molecule, and the origin |
| 366 |
< |
is used when creating a system containing the molecule. |
| 355 |
> |
As mentioned above, the fundamental unit for an {\sc OpenMD} |
| 356 |
> |
simulation is the {\tt atom}. Atoms can be collected into secondary |
| 357 |
> |
structures such as {\tt rigidBodies}, {\tt cutoffGroups}, or {\tt |
| 358 |
> |
molecules}. The {\tt molecule} is a way for {\sc OpenMD} to keep |
| 359 |
> |
track of the atoms in a simulation in logical manner. Molecular units |
| 360 |
> |
store the identities of all the atoms and rigid bodies associated with |
| 361 |
> |
themselves, and they are responsible for the evaluation of their own |
| 362 |
> |
internal interactions (\emph{i.e.}~bonds, bends, torsions, and |
| 363 |
> |
inversions). Scheme \ref{sch:mdIncludeExample} shows how one creates a |
| 364 |
> |
molecule in an included meta-data file. The positions of the atoms |
| 365 |
> |
given in the declaration are relative to the origin of the molecule, |
| 366 |
> |
and the origin is used when creating a system containing the molecule. |
| 367 |
|
|
| 368 |
|
One of the features that sets {\sc OpenMD} apart from most of the |
| 369 |
|
current molecular simulation packages is the ability to handle rigid |
| 370 |
|
body dynamics. Rigid bodies are non-spherical particles or collections |
| 371 |
< |
of particles (e.g. $\mbox{C}_{60}$) that have a constant internal |
| 371 |
> |
of particles (e.g. a phenyl ring) that have a constant internal |
| 372 |
|
potential and move collectively.\cite{Goldstein01} They are not |
| 373 |
|
included in most simulation packages because of the algorithmic |
| 374 |
|
complexity involved in propagating orientational degrees of freedom. |
| 375 |
|
Integrators which propagate orientational motion with an acceptable |
| 376 |
< |
level of energy conservation for molecular dynamics are relatively |
| 377 |
< |
new inventions. |
| 376 |
> |
level of energy conservation for molecular dynamics are relatively new |
| 377 |
> |
inventions. |
| 378 |
|
|
| 379 |
|
Moving a rigid body involves determination of both the force and |
| 380 |
|
torque applied by the surroundings, which directly affect the |
| 418 |
|
rigid body can be seen in Scheme |
| 419 |
|
\ref{sch:rigidBody}. |
| 420 |
|
|
| 421 |
< |
\begin{lstlisting}[float,caption={[Defining rigid bodies]A sample |
| 421 |
> |
\begin{code}[caption={[Defining rigid bodies]A sample |
| 422 |
|
definition of a molecule containing a rigid body and a cutoff |
| 423 |
|
group},label={sch:rigidBody}] |
| 424 |
|
molecule{ |
| 444 |
|
members(0, 1, 2); |
| 445 |
|
} |
| 446 |
|
} |
| 447 |
< |
\end{lstlisting} |
| 447 |
> |
\end{code} |
| 448 |
|
|
| 449 |
|
\section{\label{sec:miscConcepts}Creating a $<$MetaData$>$ block} |
| 450 |
|
|
| 482 |
|
minimizer parameters can be found in |
| 483 |
|
Sec.~\ref{section:minimizer}. The {\tt forceField} statement is |
| 484 |
|
important for the selection of which forces will be used in the course |
| 485 |
< |
of the simulation. {\sc OpenMD} supports several force fields, as |
| 486 |
< |
outlined in Sec.~\ref{section:empiricalEnergy}. The force fields are |
| 485 |
> |
of the simulation. {\sc OpenMD} supports several force fields, and |
| 486 |
> |
allows the user to create their own using a range of pre-defined |
| 487 |
> |
empirical energy functions. The format of force field files is |
| 488 |
> |
outlined Chapter~\ref{chapter:forceFields}. The force fields are |
| 489 |
|
interchangeable between simulations, with the only requirement being |
| 490 |
|
that all atoms needed by the simulation are defined within the |
| 491 |
|
selected force field. |
| 524 |
|
\endhead |
| 525 |
|
\hline |
| 526 |
|
\endfoot |
| 527 |
< |
{\tt forceField} & string & Sets the force field. & Possible force |
| 528 |
< |
fields are DUFF, WATER, LJ, EAM, SC, and CLAY. \\ |
| 527 |
> |
{\tt forceField} & string & Sets the base name for the force field file & |
| 528 |
> |
OpenMD appends a {\tt .frc} to the end of this to look for a force |
| 529 |
> |
field file.\\ |
| 530 |
|
{\tt component} & & Defines the molecular components of the system & |
| 531 |
|
Every {\tt $<$MetaData$>$} block must have a component statement. \\ |
| 532 |
|
{\tt minimizer} & string & Chooses a minimizer & Possible minimizers |
| 564 |
|
box must be before we can use cheaper box calculations \\ |
| 565 |
|
{\tt cutoffRadius} & $\mbox{\AA}$ & Manually sets the cutoff radius & |
| 566 |
|
the default value is set by the {\tt cutoffPolicy} \\ |
| 533 |
– |
{\tt cutoffPolicy} & string & one of mix, max, or |
| 534 |
– |
traditional & the traditional cutoff policy is to set the cutoff |
| 535 |
– |
radius for all atoms in the system to the same value (governed by the |
| 536 |
– |
largest atom). mix and max are pair-dependent cutoff |
| 537 |
– |
methods. \\ |
| 567 |
|
{\tt skinThickness} & \AA & thickness of the skin for the Verlet |
| 568 |
|
neighbor lists & defaults to 1 \AA \\ |
| 569 |
|
{\tt switchingRadius} & $\mbox{\AA}$ & Manually sets the inner radius |
| 591 |
|
default is the first eight of these columns in order.) \\ |
| 592 |
|
& & \multicolumn{2}{p{3.5in}}{Allowed |
| 593 |
|
column names are: {\sc time, total\_energy, potential\_energy, kinetic\_energy, |
| 594 |
< |
temperature, pressure, volume, conserved\_quantity, |
| 594 |
> |
temperature, pressure, volume, conserved\_quantity, hullvolume, gyrvolume, |
| 595 |
|
translational\_kinetic, rotational\_kinetic, long\_range\_potential, |
| 596 |
|
short\_range\_potential, vanderwaals\_potential, |
| 597 |
< |
electrostatic\_potential, bond\_potential, bend\_potential, |
| 598 |
< |
dihedral\_potential, improper\_potential, vraw, vharm, |
| 599 |
< |
pressure\_tensor\_x, pressure\_tensor\_y, pressure\_tensor\_z}} \\ |
| 597 |
> |
electrostatic\_potential, metallic\_potential, |
| 598 |
> |
hydrogen\_bonding\_potential, bond\_potential, bend\_potential, |
| 599 |
> |
dihedral\_potential, inversion\_potential, raw\_potential, restraint\_potential, |
| 600 |
> |
pressure\_tensor, system\_dipole, heatflux, electronic\_temperature}} \\ |
| 601 |
|
{\tt printPressureTensor} & logical & sets whether {\sc OpenMD} will print |
| 602 |
|
out the pressure tensor & can be useful for calculations of the bulk |
| 603 |
|
modulus \\ |
| 604 |
|
{\tt electrostaticSummationMethod} & & & \\ |
| 605 |
|
& string & shifted\_force, |
| 606 |
< |
shifted\_potential, shifted\_force, or reaction\_field & |
| 606 |
> |
shifted\_potential, hard, switched, taylor\_shifted, or reaction\_field & |
| 607 |
|
default is shifted\_force. \\ |
| 608 |
|
{\tt electrostaticScreeningMethod} & & & \\ |
| 609 |
|
& string & undamped or damped & default is damped \\ |
| 618 |
|
default is never \\ |
| 619 |
|
{\tt targetTemp} & K & sets the target temperature & no default value \\ |
| 620 |
|
{\tt tauThermostat} & fs & time constant for Nos\'{e}-Hoover |
| 621 |
< |
thermostat & times from 1000-10,000 fs are reasonable \\ |
| 621 |
> |
thermostat & times from 100-10,000 fs are reasonable \\ |
| 622 |
|
{\tt targetPressure} & atm & sets the target pressure & no default value\\ |
| 623 |
|
{\tt surfaceTension} & & sets the target surface tension in the x-y |
| 624 |
|
plane & no default value \\ |
| 625 |
|
{\tt tauBarostat} & fs & time constant for the |
| 626 |
< |
Nos\'{e}-Hoover-Andersen barostat & times from 10,000 to 100,000 fs |
| 626 |
> |
Nos\'{e}-Hoover-Andersen barostat & times from 1000 to 100,000 fs |
| 627 |
|
are reasonable \\ |
| 628 |
|
{\tt seed } & integer & Sets the seed value for the random number generator. & The seed needs to be at least 9 digits long. The default is to take the seed from the CPU clock. \\ |
| 629 |
|
\label{table:genParams} |
| 649 |
|
complete rotation matrix, directional entities are written out using |
| 650 |
|
quaternions to save space in the output files. |
| 651 |
|
|
| 652 |
< |
\begin{lstlisting}[float,caption={[The format of the {\tt $<$Snapshot$>$} block] |
| 652 |
> |
\begin{code}[caption={[The format of the {\tt $<$Snapshot$>$} block] |
| 653 |
|
An example of the format of the {\tt $<$Snapshot$>$} block. There is an |
| 654 |
|
initial sub-block called {\tt $<$FrameData$>$} which contains the time |
| 655 |
|
stamp, the three column vectors of $\mathsf{H}$, and optional extra |
| 658 |
|
configuration of each integrable object. For each integrable object, |
| 659 |
|
the global index is followed by a short string describing what |
| 660 |
|
additional information is present on the line. Atoms with only |
| 661 |
< |
position and velocity information use the ``pv'' string which must |
| 661 |
> |
position and velocity information use the {\tt pv} string which must |
| 662 |
|
then be followed by the position and velocity vectors for that atom. |
| 663 |
< |
Directional atoms and Rigid Bodies typically use the ``pvqj'' string |
| 663 |
> |
Directional atoms and Rigid Bodies typically use the {\tt pvqj} string |
| 664 |
|
which is followed by position, velocity, quaternions, and |
| 665 |
< |
lastly, body fixed angular momentum for that integrable object.}, |
| 636 |
< |
label=sch:dumpFormat] |
| 665 |
> |
lastly, body fixed angular momentum for that integrable object.},label={sch:dumpFormat}] |
| 666 |
|
<Snapshot> |
| 667 |
|
<FrameData> |
| 668 |
|
Time: 0 |
| 677 |
|
3 pvqj x y z vx vy vz qw qx qy qz jx jy jz |
| 678 |
|
</StuntDoubles> |
| 679 |
|
</Snapshot> |
| 680 |
< |
\end{lstlisting} |
| 680 |
> |
\end{code} |
| 681 |
|
|
| 682 |
|
There are three {\sc OpenMD} files that are written using the combined |
| 683 |
|
format. They are: the initial startup file (\texttt{.md}), the |
| 728 |
|
\end{enumerate} |
| 729 |
|
An example is given in the {\sc OpenMD} file in Scheme~\ref{sch:initEx1}. |
| 730 |
|
|
| 731 |
< |
\begin{lstlisting}[float,caption={Example declaration of the |
| 732 |
< |
$\text{I}_2$ molecule and the HCl molecule in $<$MetaData$>$ and |
| 733 |
< |
$<$Snapshot$>$ blocks. Note that even though $\text{I}_2$ is |
| 734 |
< |
declared before HCl, the $<$Snapshot$>$ block follows the order {\it in |
| 731 |
> |
\begin{code}[caption={Example declaration of the |
| 732 |
> |
$\text{I}_2$ molecule and the HCl molecule in {\tt <MetaData>} and |
| 733 |
> |
{\tt <Snapshot>} blocks. Note that even though $\text{I}_2$ is |
| 734 |
> |
declared before HCl, the {\tt <Snapshot>} block follows the order {\it in |
| 735 |
|
which the components were included}.}, label=sch:initEx1] |
| 736 |
|
<OpenMD> |
| 737 |
|
<MetaData> |
| 775 |
|
</StuntDoubles> |
| 776 |
|
</Snapshot> |
| 777 |
|
</OpenMD> |
| 778 |
< |
\end{lstlisting} |
| 778 |
> |
\end{code} |
| 779 |
|
|
| 780 |
|
\section{The Statistics File} |
| 781 |
|
|
| 791 |
|
allowing the user to gauge the stability of the integrator. The |
| 792 |
|
statistics file is denoted with the \texttt{.stat} file extension. |
| 793 |
|
|
| 794 |
< |
\chapter{\label{section:empiricalEnergy}The Empirical Energy |
| 766 |
< |
Functions} |
| 794 |
> |
\chapter{\label{chapter:forceFields}Force Fields} |
| 795 |
|
|
| 796 |
< |
Like many simulation packages, {\sc OpenMD} splits the potential energy |
| 797 |
< |
into the short-ranged (bonded) portion and a long-range (non-bonded) |
| 798 |
< |
potential, |
| 796 |
> |
Like many molecular simulation packages, {\sc OpenMD} splits the |
| 797 |
> |
potential energy into the short-ranged (bonded) portion and a |
| 798 |
> |
long-range (non-bonded) potential, |
| 799 |
|
\begin{equation} |
| 800 |
|
V = V_{\mathrm{short-range}} + V_{\mathrm{long-range}}. |
| 801 |
|
\end{equation} |
| 802 |
< |
The short-ranged portion includes the explicit bonds, bends, and |
| 803 |
< |
torsions which have been defined in the meta-data file for the |
| 804 |
< |
molecules which are present in the simulation. The functional forms and |
| 805 |
< |
parameters for these interactions are defined by the force field which |
| 806 |
< |
is chosen. |
| 802 |
> |
The short-ranged portion includes the bonds, bends, torsions, and |
| 803 |
> |
inversions which have been defined in the meta-data file for the |
| 804 |
> |
molecules. The functional forms and parameters for these interactions |
| 805 |
> |
are defined by the force field which is selected in the MetaData |
| 806 |
> |
section. |
| 807 |
|
|
| 808 |
< |
Calculating the long-range (non-bonded) potential involves a sum over |
| 809 |
< |
all pairs of atoms (except for those atoms which are involved in a |
| 782 |
< |
bond, bend, or torsion with each other). If done poorly, calculating |
| 783 |
< |
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: |
| 793 |
< |
\begin{equation} |
| 794 |
< |
V_{\mathrm{long-range}} = \sum_{a} \sum_{b>a} s(r_{ab}) \sum_{i \in a} \sum_{j \in b} V_{ij}(r_{ij}), |
| 795 |
< |
\end{equation} |
| 796 |
< |
where $r_{ab}$ is the distance between the centers of mass of the two |
| 797 |
< |
cutoff groups ($a$ and $b$). |
| 808 |
> |
\section{\label{section:divisionOfLabor}Separation into Internal and |
| 809 |
> |
Cross interactions} |
| 810 |
|
|
| 811 |
< |
The sums over $a$ and $b$ are over the cutoff groups that are present |
| 812 |
< |
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, |
| 811 |
> |
The classical potential energy function for a system composed of $N$ |
| 812 |
> |
molecules is traditionally written |
| 813 |
|
\begin{equation} |
| 814 |
< |
S(r) = |
| 815 |
< |
\begin{cases} |
| 816 |
< |
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} |
| 814 |
> |
V = \sum^{N}_{I=1} V^{I}_{\text{Internal}} |
| 815 |
> |
+ \sum^{N-1}_{I=1} \sum_{J>I} V^{IJ}_{\text{Cross}}, |
| 816 |
> |
\label{eq:totalPotential} |
| 817 |
|
\end{equation} |
| 818 |
< |
Here, $r_{\text{sw}}$ is the {\tt switchingRadius}, or the distance |
| 819 |
< |
beyond which interactions are reduced, and $r_{\text{cut}}$ is the |
| 820 |
< |
{\tt cutoffRadius}, or the distance at which interactions are |
| 821 |
< |
truncated. |
| 818 |
> |
where $V^{I}_{\text{Internal}}$ contains all of the terms internal to |
| 819 |
> |
molecule $I$ (e.g. bonding, bending, torsional, and inversion terms) |
| 820 |
> |
and $V^{IJ}_{\text{Cross}}$ contains all intermolecular interactions |
| 821 |
> |
between molecules $I$ and $J$. For large molecules, the internal |
| 822 |
> |
potential may also include some non-bonded terms like electrostatic or |
| 823 |
> |
van der Waals interactions. |
| 824 |
|
|
| 825 |
< |
Users of {\sc OpenMD} do not need to specify the {\tt cutoffRadius} or |
| 826 |
< |
{\tt switchingRadius}. In simulations containing only Lennard-Jones |
| 827 |
< |
atoms, the cutoff radius has a default value of $2.5\sigma_{ii}$, |
| 828 |
< |
where $\sigma_{ii}$ is the largest Lennard-Jones length parameter |
| 829 |
< |
present in the simulation. In simulations containing charged or |
| 825 |
< |
dipolar atoms, the default cutoff radius is $15 \mbox{\AA}$. |
| 825 |
> |
The types of atoms being simulated, as well as the specific functional |
| 826 |
> |
forms and parameters of the intra-molecular functions and the |
| 827 |
> |
long-range potentials are defined by the force field. In the following |
| 828 |
> |
sections we discuss the stucture of an OpenMD force field file and the |
| 829 |
> |
specification of blocks that may be present within these files. |
| 830 |
|
|
| 831 |
< |
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. |
| 831 |
> |
\section{\label{section:frcFile}Force Field Files} |
| 832 |
|
|
| 833 |
< |
Force fields can be added to {\sc OpenMD}, although it comes with a few |
| 834 |
< |
simple examples (Lennard-Jones, {\sc duff}, {\sc water}, and {\sc |
| 835 |
< |
eam}) which are explained in the following sections. |
| 833 |
> |
Force field files have a number of ``Blocks'' to delineate different |
| 834 |
> |
types of information. The blocks contain AtomType data, which provide |
| 835 |
> |
properties belonging to a single AtomType, as well as interaction |
| 836 |
> |
information which provides information about bonded or non-bonded |
| 837 |
> |
interactions that cannot be deduced from AtomType information alone. |
| 838 |
> |
A simple example of a forceField file is shown in scheme |
| 839 |
> |
\ref{sch:frcExample}. |
| 840 |
|
|
| 841 |
< |
\section{\label{sec:LJPot}The Lennard Jones Force Field} |
| 841 |
> |
\begin{code}[caption={[An example of a complete OpenMD |
| 842 |
> |
force field file for straight-chain united-atom alkanes.] An example |
| 843 |
> |
showing a complete OpenMD force field for straight-chain united-atom |
| 844 |
> |
alkanes.}, label={sch:frcExample}] |
| 845 |
> |
begin Options |
| 846 |
> |
Name = "alkane" |
| 847 |
> |
end Options |
| 848 |
|
|
| 849 |
< |
The most basic force field implemented in {\sc OpenMD} is the |
| 850 |
< |
Lennard-Jones force field, which mimics the van der Waals interaction |
| 851 |
< |
at long distances and uses an empirical repulsion at short |
| 852 |
< |
distances. The Lennard-Jones potential is given by: |
| 849 |
> |
begin BaseAtomTypes |
| 850 |
> |
//name mass |
| 851 |
> |
C 12.0107 |
| 852 |
> |
end BaseAtomTypes |
| 853 |
> |
|
| 854 |
> |
begin AtomTypes |
| 855 |
> |
//name base mass |
| 856 |
> |
CH4 C 16.05 |
| 857 |
> |
CH3 C 15.04 |
| 858 |
> |
CH2 C 14.03 |
| 859 |
> |
end AtomTypes |
| 860 |
> |
|
| 861 |
> |
begin LennardJonesAtomTypes |
| 862 |
> |
//name epsilon sigma |
| 863 |
> |
CH4 0.2941 3.73 |
| 864 |
> |
CH3 0.1947 3.75 |
| 865 |
> |
CH2 0.09140 3.95 |
| 866 |
> |
end LennardJonesAtomTypes |
| 867 |
> |
|
| 868 |
> |
begin BondTypes |
| 869 |
> |
//AT1 AT2 Type r0 k |
| 870 |
> |
CH3 CH3 Harmonic 1.526 260 |
| 871 |
> |
CH3 CH2 Harmonic 1.526 260 |
| 872 |
> |
CH2 CH2 Harmonic 1.526 260 |
| 873 |
> |
end BondTypes |
| 874 |
> |
|
| 875 |
> |
begin BendTypes |
| 876 |
> |
//AT1 AT2 AT3 Type theta0 k |
| 877 |
> |
CH3 CH2 CH3 Harmonic 114.0 124.19 |
| 878 |
> |
CH3 CH2 CH2 Harmonic 114.0 124.19 |
| 879 |
> |
CH2 CH2 CH2 Harmonic 114.0 124.19 |
| 880 |
> |
end BendTypes |
| 881 |
> |
|
| 882 |
> |
begin TorsionTypes |
| 883 |
> |
//AT1 AT2 AT3 AT4 Type |
| 884 |
> |
CH3 CH2 CH2 CH3 Trappe 0.0 0.70544 -0.13549 1.5723 |
| 885 |
> |
CH3 CH2 CH2 CH2 Trappe 0.0 0.70544 -0.13549 1.5723 |
| 886 |
> |
CH2 CH2 CH2 CH2 Trappe 0.0 0.70544 -0.13549 1.5723 |
| 887 |
> |
end TorsionTypes |
| 888 |
> |
\end{code} |
| 889 |
> |
|
| 890 |
> |
\section{\label{section:ffOptions}The Options block} |
| 891 |
> |
|
| 892 |
> |
The Options block defines properties governing how the force field |
| 893 |
> |
interactions are carried out. This block is delineated with the text |
| 894 |
> |
tags {\tt begin Options} and {\tt end Options}. Most options don't |
| 895 |
> |
need to be set as they come with fairly sensible default values, but |
| 896 |
> |
the various keywords and their possible values are given in Scheme |
| 897 |
> |
\ref{sch:optionsBlock}. |
| 898 |
> |
|
| 899 |
> |
\begin{code}[caption={[A force field Options block showing default values |
| 900 |
> |
for many force field options.] A force field Options block showing default values |
| 901 |
> |
for many force field options. Most of these options do not need to be |
| 902 |
> |
specified if the default values are working.}, |
| 903 |
> |
label={sch:optionsBlock}] |
| 904 |
> |
begin Options |
| 905 |
> |
Name = "alkane" // any string |
| 906 |
> |
vdWtype = "Lennard-Jones" |
| 907 |
> |
DistanceMixingRule = "arithmetic" // can also be "geometric" or "cubic" |
| 908 |
> |
DistanceType = "sigma" // can also be "Rmin" |
| 909 |
> |
EnergyMixingRule = "geometric" // can also be "arithmetic" or "hhg" |
| 910 |
> |
EnergyUnitScaling = 1.0 |
| 911 |
> |
MetallicEnergyUnitScaling = 1.0 |
| 912 |
> |
DistanceUnitScaling = 1.0 |
| 913 |
> |
AngleUnitScaling = 1.0 |
| 914 |
> |
TorsionAngleConvention = "180_is_trans" // can also be "0_is_trans" |
| 915 |
> |
vdW-12-scale = 0.0 |
| 916 |
> |
vdW-13-scale = 0.0 |
| 917 |
> |
vdW-14-scale = 0.0 |
| 918 |
> |
electrostatic-12-scale = 0.0 |
| 919 |
> |
electrostatic-13-scale = 0.0 |
| 920 |
> |
electrostatic-14-scale = 0.0 |
| 921 |
> |
GayBerneMu = 2.0 |
| 922 |
> |
GayBerneNu = 1.0 |
| 923 |
> |
EAMMixingMethod = "Johnson" // can also be "Daw" |
| 924 |
> |
end Options |
| 925 |
> |
\end{code} |
| 926 |
> |
|
| 927 |
> |
\section{\label{section:ffBase}The BaseAtomTypes block} |
| 928 |
> |
|
| 929 |
> |
An AtomType the primary data structure that OpenMD uses to store |
| 930 |
> |
static data about an atom. Things that belong to AtomType are |
| 931 |
> |
universal properties (i.e. does this atom have a fixed charge? What |
| 932 |
> |
is its mass?) Dynamic properties of an atom are not intended to be |
| 933 |
> |
properties of an atom type. A BaseAtomType can be used to build |
| 934 |
> |
extended sets of related atom types that all fall back to one |
| 935 |
> |
particular type. For example, one might want a series of atomTypes |
| 936 |
> |
that inherit from more basic types: |
| 937 |
> |
\begin{displaymath} |
| 938 |
> |
\mathtt{ALA-CA} \rightarrow \mathtt{CT} \rightarrow \mathtt{CSP3} \rightarrow \mathtt{C} |
| 939 |
> |
\end{displaymath} |
| 940 |
> |
where for each step to the right, the atomType falls back to more |
| 941 |
> |
primitive data. That is, the mass could be a property of the {\tt C} |
| 942 |
> |
type, while Lennard-Jones parameters could be properties of the {\tt |
| 943 |
> |
CSP3} type. {\tt CT} could have charge information and its own set |
| 944 |
> |
of Lennard-Jones parameter that override the CSP3 parameters. And the |
| 945 |
> |
{\tt ALA-CA} type might have specific torsion or charge information |
| 946 |
> |
that override the lower level types. A BaseAtomType contains only |
| 947 |
> |
information a primitive name and the mass of this atom type. |
| 948 |
> |
BaseAtomTypes can also be useful in creating files that can be easily |
| 949 |
> |
viewed in visualization programs. The {\tt Dump2XYZ} utility has the |
| 950 |
> |
ability to print out the names of the base atom types for displaying |
| 951 |
> |
simulations in Jmol or VMD. |
| 952 |
> |
|
| 953 |
> |
\begin{code}[caption={[A simple example of a BaseAtomTypes |
| 954 |
> |
block.] A simple example of a BaseAtomTypes block.}, |
| 955 |
> |
label={sch:baseAtomTypesBlock}] |
| 956 |
> |
begin BaseAtomTypes |
| 957 |
> |
//Name mass (amu) |
| 958 |
> |
H 1.0079 |
| 959 |
> |
O 15.9994 |
| 960 |
> |
Si 28.0855 |
| 961 |
> |
Al 26.981538 |
| 962 |
> |
Mg 24.3050 |
| 963 |
> |
Ca 40.078 |
| 964 |
> |
Fe 55.845 |
| 965 |
> |
Li 6.941 |
| 966 |
> |
Na 22.98977 |
| 967 |
> |
K 39.0983 |
| 968 |
> |
Cs 132.90545 |
| 969 |
> |
Ca 40.078 |
| 970 |
> |
Ba 137.327 |
| 971 |
> |
Cl 35.453 |
| 972 |
> |
end BaseAtomTypes |
| 973 |
> |
\end{code} |
| 974 |
> |
|
| 975 |
> |
\section{\label{section:ffAtom}The AtomTypes block} |
| 976 |
> |
|
| 977 |
> |
AtomTypes inherit most properties from BaseAtomTypes, but can override |
| 978 |
> |
their lower-level properties as well. Scheme \ref{sch:atomTypesBlock} |
| 979 |
> |
shows an example where multiple types of oxygen atoms can inherit mass |
| 980 |
> |
from the oxygen base type. |
| 981 |
> |
|
| 982 |
> |
\begin{code}[caption={[An example of a AtomTypes block.] A |
| 983 |
> |
simple example of an AtomTypes block which |
| 984 |
> |
shows how multiple types can inherit from the same base type.}, |
| 985 |
> |
label={sch:atomTypesBlock}] |
| 986 |
> |
begin AtomTypes |
| 987 |
> |
//Name baseatomtype |
| 988 |
> |
h* H |
| 989 |
> |
ho H |
| 990 |
> |
o* O |
| 991 |
> |
oh O |
| 992 |
> |
ob O |
| 993 |
> |
obos O |
| 994 |
> |
obts O |
| 995 |
> |
obss O |
| 996 |
> |
ohs O |
| 997 |
> |
st Si |
| 998 |
> |
ao Al |
| 999 |
> |
at Al |
| 1000 |
> |
mgo Mg |
| 1001 |
> |
mgh Mg |
| 1002 |
> |
cao Ca |
| 1003 |
> |
cah Ca |
| 1004 |
> |
feo Fe |
| 1005 |
> |
lio Li |
| 1006 |
> |
end AtomTypes |
| 1007 |
> |
\end{code} |
| 1008 |
> |
|
| 1009 |
> |
\section{\label{section:ffDirectionalAtom}The DirectionalAtomTypes |
| 1010 |
> |
block} |
| 1011 |
> |
DirectionalAtoms have orientational degrees of freedom as well as |
| 1012 |
> |
translation, so moving these atoms requires information about the |
| 1013 |
> |
moments of inertias in the same way that translational motion requires |
| 1014 |
> |
mass. For DirectionalAtoms, OpenMD treats the mass distribution with |
| 1015 |
> |
higher priority than electrostatic distributions; the moment of |
| 1016 |
> |
inertia tensor, $\overleftrightarrow{\mathsf I}$, should be |
| 1017 |
> |
diagonalized to obtain body-fixed axes, and the three diagonal moments |
| 1018 |
> |
should correspond to rotational motion \textit{around} each of these |
| 1019 |
> |
body-fixed axes. Charge distributions may then result in dipole |
| 1020 |
> |
vectors that are oriented along a linear combination of the body-axes, |
| 1021 |
> |
and in quadrupole tensors that are not necessarily diagonal in the |
| 1022 |
> |
body frame. |
| 1023 |
> |
|
| 1024 |
> |
\begin{code}[caption={[An example of a DirectionalAtomTypes block.] A |
| 1025 |
> |
simple example of a DirectionalAtomTypes block.}, |
| 1026 |
> |
label={sch:datomTypesBlock}] |
| 1027 |
> |
begin DirectionalAtomTypes |
| 1028 |
> |
//Name I_xx I_yy I_zz (All moments in (amu*Ang^2) |
| 1029 |
> |
SSD 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 |
> |
CO2 43.06 43.06 0.0 // single-site model for CO2 |
| 1034 |
> |
end DirectionalAtomTypes |
| 1035 |
> |
|
| 1036 |
> |
\end{code} |
| 1037 |
> |
|
| 1038 |
> |
For a DirectionalAtom that represents a linear object, it is |
| 1039 |
> |
appropriate for one of the moments of inertia to be zero. In this |
| 1040 |
> |
case, OpenMD identifies that DirectionalAtom as having only 5 degrees |
| 1041 |
> |
of freedom (three translations and two rotations), and will alter |
| 1042 |
> |
calculation of temperatures to reflect this. |
| 1043 |
> |
|
| 1044 |
> |
\section{\label{section::ffAtomProperties}AtomType properties} |
| 1045 |
> |
\subsection{\label{section:ffLJ}The LennardJonesAtomTypes block} |
| 1046 |
> |
One of the most basic interatomic interactions implemented in {\sc |
| 1047 |
> |
OpenMD} is the Lennard-Jones potential, which mimics the van der |
| 1048 |
> |
Waals interaction at long distances and uses an empirical repulsion at |
| 1049 |
> |
short distances. The Lennard-Jones potential is given by: |
| 1050 |
|
\begin{equation} |
| 1051 |
|
V_{\text{LJ}}(r_{ij}) = |
| 1052 |
|
4\epsilon_{ij} \biggl[ |
| 1057 |
|
\end{equation} |
| 1058 |
|
where $r_{ij}$ is the distance between particles $i$ and $j$, |
| 1059 |
|
$\sigma_{ij}$ scales the length of the interaction, and |
| 1060 |
< |
$\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. |
| 1060 |
> |
$\epsilon_{ij}$ scales the well depth of the potential. |
| 1061 |
|
|
| 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 |
– |
|
| 1062 |
|
Interactions between dissimilar particles requires the generation of |
| 1063 |
|
cross term parameters for $\sigma$ and $\epsilon$. These parameters |
| 1064 |
< |
are determined using the Lorentz-Berthelot mixing |
| 1064 |
> |
are usually determined using the Lorentz-Berthelot mixing |
| 1065 |
|
rules:\cite{Allen87} |
| 1066 |
|
\begin{equation} |
| 1067 |
|
\sigma_{ij} = \frac{1}{2}[\sigma_{ii} + \sigma_{jj}], |
| 1073 |
|
\label{eq:epsilonMix} |
| 1074 |
|
\end{equation} |
| 1075 |
|
|
| 1076 |
< |
\section{\label{section:DUFF}Dipolar Unified-Atom Force Field} |
| 1076 |
> |
The values of $\sigma_{ii}$ and $\epsilon_{ii}$ are properties of atom |
| 1077 |
> |
type $i$, and must be specified in a section of the force field file |
| 1078 |
> |
called the {\tt LennardJonesAtomTypes} block (see listing |
| 1079 |
> |
\ref{sch:LJatomTypesBlock}). Separate Lennard-Jones interactions |
| 1080 |
> |
which are not determined by the mixing rules can also be specified in |
| 1081 |
> |
the {\tt NonbondedInteractionTypes} block (see section |
| 1082 |
> |
\ref{section:ffNBinteraction}). |
| 1083 |
|
|
| 1084 |
< |
The dipolar unified-atom force field ({\sc duff}) was developed to |
| 1085 |
< |
simulate lipid bilayers. These types of simulations require a model |
| 1086 |
< |
capable of forming bilayers, while still being sufficiently |
| 1087 |
< |
computationally efficient to allow large systems ($\sim$100's of |
| 1088 |
< |
phospholipids, $\sim$1000's of waters) to be simulated for long times |
| 1089 |
< |
($\sim$10's of nanoseconds). With this goal in mind, {\sc duff} has no |
| 1090 |
< |
point charges. Charge-neutral distributions are replaced with dipoles, |
| 1091 |
< |
while most atoms and groups of atoms are reduced to Lennard-Jones |
| 1092 |
< |
interaction sites. This simplification reduces the length scale of |
| 1093 |
< |
long range interactions from $\frac{1}{r}$ to $\frac{1}{r^3}$, |
| 1094 |
< |
removing the need for the computationally expensive Ewald |
| 1095 |
< |
sum. Instead, Verlet neighbor-lists and cutoff radii are used for the |
| 1096 |
< |
dipolar interactions, and, if desired, a reaction field may be added |
| 1097 |
< |
to mimic longer range interactions. |
| 1084 |
> |
\begin{code}[caption={[An example of a LennardJonesAtomTypes block.] A |
| 1085 |
> |
simple example of a LennardJonesAtomTypee block. Units for |
| 1086 |
> |
$\epsilon$ are kcal / mol and for $\sigma$ are \AA\ .}, |
| 1087 |
> |
label={sch:LJatomTypesBlock}] |
| 1088 |
> |
begin LennardJonesAtomTypes |
| 1089 |
> |
//Name epsilon sigma |
| 1090 |
> |
O_TIP4P 0.1550 3.15365 |
| 1091 |
> |
O_TIP4P-Ew 0.16275 3.16435 |
| 1092 |
> |
O_TIP5P 0.16 3.12 |
| 1093 |
> |
O_TIP5P-E 0.178 3.097 |
| 1094 |
> |
O_SPCE 0.15532 3.16549 |
| 1095 |
> |
O_SPC 0.15532 3.16549 |
| 1096 |
> |
CH4 0.279 3.73 |
| 1097 |
> |
CH3 0.185 3.75 |
| 1098 |
> |
CH2 0.0866 3.95 |
| 1099 |
> |
CH 0.0189 4.68 |
| 1100 |
> |
end LennardJonesAtomTypes |
| 1101 |
> |
\end{code} |
| 1102 |
|
|
| 1103 |
< |
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} |
| 1103 |
> |
\subsection{\label{section:ffCharge}The ChargeAtomTypes block} |
| 1104 |
|
|
| 1105 |
< |
\begin{figure} |
| 1106 |
< |
\centering |
| 1107 |
< |
\includegraphics[width=\linewidth]{lipidModel.pdf} |
| 1108 |
< |
\caption[A representation of a lipid model in {\sc duff}]{A |
| 1109 |
< |
representation of the lipid model. $\phi$ is the torsion angle, |
| 1110 |
< |
$\theta$ is the bend angle, and $\mu$ is the dipole moment of the head |
| 1111 |
< |
group.} |
| 1112 |
< |
\label{fig:lipidModel} |
| 1113 |
< |
\end{figure} |
| 1114 |
< |
|
| 1115 |
< |
A set of scalable parameters has been used to model the alkyl groups |
| 1116 |
< |
with Lennard-Jones sites. For this, parameters from the TraPPE force |
| 1117 |
< |
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 |
| 1105 |
> |
In molecular simulations, proper accumulation of the electrostatic |
| 1106 |
> |
interactions is essential and is one of the most |
| 1107 |
> |
computationally-demanding tasks. Most common molecular mechanics |
| 1108 |
> |
force fields represent atomic sites with full or partial charges |
| 1109 |
> |
protected by Lennard-Jones (short range) interactions. Partial charge |
| 1110 |
> |
values, $q_i$ are empirical representations of the distribution of |
| 1111 |
> |
electronic charge in a molecule. This means that nearly every pair |
| 1112 |
> |
interaction involves a calculation of charge-charge forces. Coupled |
| 1113 |
> |
with relatively long-ranged $r^{-1}$ decay, the monopole interactions |
| 1114 |
> |
quickly become the most expensive part of molecular simulations. The |
| 1115 |
> |
interactions between point charges can be handled via a number of |
| 1116 |
> |
different algorithms, but Coulomb's law is the fundamental physical |
| 1117 |
> |
principle governing these interactions, |
| 1118 |
|
\begin{equation} |
| 1119 |
< |
V = \sum^{N}_{I=1} V^{I}_{\text{Internal}} |
| 1120 |
< |
+ \sum^{N-1}_{I=1} \sum_{J>I} V^{IJ}_{\text{Cross}}, |
| 982 |
< |
\label{eq:totalPotential} |
| 1119 |
> |
V_{\text{charge}}(r_{ij}) = \sum_{ij}\frac{q_iq_je^2}{4 \pi \epsilon_0 |
| 1120 |
> |
r_{ij}}, |
| 1121 |
|
\end{equation} |
| 1122 |
< |
where $V^{I}_{\text{Internal}}$ is the internal potential of molecule $I$: |
| 1123 |
< |
\begin{equation} |
| 1124 |
< |
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. |
| 1122 |
> |
where $q$ represents the charge on particle $i$ or $j$, and $e$ is the |
| 1123 |
> |
charge of an electron in Coulombs. $\epsilon_0$ is the permittivity |
| 1124 |
> |
of free space. |
| 1125 |
|
|
| 1126 |
< |
The bend potential of a molecule is represented by the following function: |
| 1127 |
< |
\begin{equation} |
| 1128 |
< |
V_{\text{bend}}(\theta_{ijk}) = k_{\theta}( \theta_{ijk} - \theta_0 |
| 1129 |
< |
)^2, \label{eq:bendPot} |
| 1130 |
< |
\end{equation} |
| 1131 |
< |
where $\theta_{ijk}$ is the angle defined by atoms $i$, $j$, and $k$ |
| 1132 |
< |
(see Fig.~\ref{fig:lipidModel}), $\theta_0$ is the equilibrium |
| 1133 |
< |
bond angle, and $k_{\theta}$ is the force constant which determines the |
| 1134 |
< |
strength of the harmonic bend. The parameters for $k_{\theta}$ and |
| 1135 |
< |
$\theta_0$ are borrowed from those in TraPPE.\cite{Siepmann1998} |
| 1126 |
> |
\begin{code}[caption={[An example of a ChargeAtomTypes block.] A |
| 1127 |
> |
simple example of a ChargeAtomTypes block. Units for |
| 1128 |
> |
charge are in multiples of electron charge.}, |
| 1129 |
> |
label={sch:ChargeAtomTypesBlock}] |
| 1130 |
> |
begin ChargeAtomTypes |
| 1131 |
> |
// Name charge |
| 1132 |
> |
O_TIP3P -0.834 |
| 1133 |
> |
O_SPCE -0.8476 |
| 1134 |
> |
H_TIP3P 0.417 |
| 1135 |
> |
H_TIP4P 0.520 |
| 1136 |
> |
H_SPCE 0.4238 |
| 1137 |
> |
EP_TIP4P -1.040 |
| 1138 |
> |
Na+ 1.0 |
| 1139 |
> |
Cl- -1.0 |
| 1140 |
> |
end ChargeAtomTypes |
| 1141 |
> |
\end{code} |
| 1142 |
|
|
| 1143 |
< |
The torsion potential and parameters are also borrowed from TraPPE. It is |
| 1144 |
< |
of the form: |
| 1143 |
> |
\subsection{\label{section:ffMultipole}The MultipoleAtomTypes |
| 1144 |
> |
block} |
| 1145 |
> |
For complex charge distributions that are centered on single sites, it |
| 1146 |
> |
is convenient to write the total electrostatic potential in terms of |
| 1147 |
> |
multipole moments, |
| 1148 |
|
\begin{equation} |
| 1149 |
< |
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} |
| 1149 |
> |
U_{\bf{ab}}(r)=\hat{M}_{\bf a} \hat{M}_{\bf b} \frac{1}{r} \label{kernel}. |
| 1150 |
|
\end{equation} |
| 1151 |
< |
where: |
| 1151 |
> |
where the multipole operator on site $\bf a$, |
| 1152 |
|
\begin{equation} |
| 1153 |
< |
\cos\phi = (\hat{\mathbf{r}}_{ij} \times \hat{\mathbf{r}}_{jk}) \cdot |
| 1154 |
< |
(\hat{\mathbf{r}}_{jk} \times \hat{\mathbf{r}}_{kl}). |
| 1155 |
< |
\label{eq:torsPhi} |
| 1153 |
> |
\hat{M}_{\bf a} = C_{\bf a} - D_{{\bf a}\alpha} \frac{\partial}{\partial r_{\alpha}} |
| 1154 |
> |
+ Q_{{\bf a}\alpha\beta} |
| 1155 |
> |
\frac{\partial^2}{\partial r_{\alpha} \partial r_{\beta}} + \dots |
| 1156 |
|
\end{equation} |
| 1157 |
< |
Here, $\hat{\mathbf{r}}_{\alpha\beta}$ are the set of unit bond |
| 1158 |
< |
vectors between atoms $i$, $j$, $k$, and $l$. For computational |
| 1159 |
< |
efficiency, the torsion potential has been recast after the method of |
| 1160 |
< |
{\sc charmm},\cite{Brooks83} in which the angle series is converted to |
| 1161 |
< |
a power series of the form: |
| 1162 |
< |
\begin{equation} |
| 1163 |
< |
V_{\text{torsion}}(\phi) = |
| 1164 |
< |
k_3 \cos^3 \phi + k_2 \cos^2 \phi + k_1 \cos \phi + k_0, |
| 1165 |
< |
\label{eq:torsionPot} |
| 1166 |
< |
\end{equation} |
| 1167 |
< |
where: |
| 1168 |
< |
\begin{align*} |
| 1169 |
< |
k_0 &= c_1 + c_3, \\ |
| 1170 |
< |
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 |
< |
|
| 1157 |
> |
Here, the point charge, dipole, and quadrupole for site $\bf a$ are |
| 1158 |
> |
given by $C_{\bf a}$, $D_{{\bf a}\alpha}$, and $Q_{{\bf |
| 1159 |
> |
a}\alpha\beta}$, respectively. These are the {\it primitive} |
| 1160 |
> |
multipoles. If the site is representing a distribution of charges, |
| 1161 |
> |
these can be expressed as, |
| 1162 |
> |
\begin{align} |
| 1163 |
> |
C_{\bf a} =&\sum_{k \, \text{in \bf a}} q_k , \label{eq:charge} \\ |
| 1164 |
> |
D_{{\bf a}\alpha} =&\sum_{k \, \text{in \bf a}} q_k r_{k\alpha}, \label{eq:dipole}\\ |
| 1165 |
> |
Q_{{\bf a}\alpha\beta} =& \frac{1}{2} \sum_{k \, \text{in \bf a}} q_k |
| 1166 |
> |
r_{k\alpha} r_{k\beta} . \label{eq:quadrupole} |
| 1167 |
> |
\end{align} |
| 1168 |
> |
Note that the definition of the primitive quadrupole here differs from |
| 1169 |
> |
the standard traceless form, and contains an additional Taylor-series |
| 1170 |
> |
based factor of $1/2$. |
| 1171 |
|
|
| 1172 |
< |
The cross potential between molecules $I$ and $J$, |
| 1173 |
< |
$V^{IJ}_{\text{Cross}}$, is as follows: |
| 1172 |
> |
The details of the multipolar interactions will be given later, but |
| 1173 |
> |
many readers are familiar with the dipole-dipole potential: |
| 1174 |
|
\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 |
– |
|
| 1067 |
– |
The dipole-dipole potential has the following form: |
| 1068 |
– |
\begin{equation} |
| 1175 |
|
V_{\text{dipole}}(\mathbf{r}_{ij},\boldsymbol{\Omega}_{i}, |
| 1176 |
< |
\boldsymbol{\Omega}_{j}) = \frac{|\mu_i||\mu_j|}{4\pi\epsilon_{0}r_{ij}^{3}} \biggl[ |
| 1176 |
> |
\boldsymbol{\Omega}_{j}) = \frac{|{\bf D}_i||{\bf D}_j|}{4\pi\epsilon_{0}r_{ij}^{3}} \biggl[ |
| 1177 |
|
\boldsymbol{\hat{u}}_{i} \cdot \boldsymbol{\hat{u}}_{j} |
| 1178 |
|
- |
| 1179 |
|
3(\boldsymbol{\hat{u}}_i \cdot \hat{\mathbf{r}}_{ij}) % |
| 1183 |
|
Here $\mathbf{r}_{ij}$ is the vector starting at atom $i$ pointing |
| 1184 |
|
towards $j$, and $\boldsymbol{\Omega}_i$ and $\boldsymbol{\Omega}_j$ |
| 1185 |
|
are the orientational degrees of freedom for atoms $i$ and $j$ |
| 1186 |
< |
respectively. The magnitude of the dipole moment of atom $i$ is |
| 1187 |
< |
$|\mu_i|$, $\boldsymbol{\hat{u}}_i$ is the standard unit orientation |
| 1186 |
> |
respectively. The magnitude of the dipole moment of atom $i$ is $|{\bf |
| 1187 |
> |
D}_i|$, $\boldsymbol{\hat{u}}_i$ is the standard unit orientation |
| 1188 |
|
vector of $\boldsymbol{\Omega}_i$, and $\boldsymbol{\hat{r}}_{ij}$ is |
| 1189 |
|
the unit vector pointing along $\mathbf{r}_{ij}$ |
| 1190 |
|
($\boldsymbol{\hat{r}}_{ij}=\mathbf{r}_{ij}/|\mathbf{r}_{ij}|$). |
| 1191 |
|
|
| 1086 |
– |
\subsection{\label{section:SSD}The {\sc duff} Water Models: SSD/E |
| 1087 |
– |
and SSD/RF} |
| 1192 |
|
|
| 1193 |
< |
In the interest of computational efficiency, the default solvent used |
| 1194 |
< |
by {\sc OpenMD} is the extended Soft Sticky Dipole (SSD/E) water |
| 1195 |
< |
model.\cite{fennell04} The original SSD was developed by Ichiye |
| 1196 |
< |
\emph{et al.}\cite{liu96:new_model} as a modified form of the hard-sphere |
| 1193 |
> |
\begin{code}[caption={[An example of a MultipoleAtomTypes block.] A |
| 1194 |
> |
simple example of a MultipoleAtomTypes block. Dipoles are given in |
| 1195 |
> |
units of Debyes, and Quadrupole moments are given in units of Debye |
| 1196 |
> |
\AA~(or $10^{-26} \mathrm{~esu~cm}^2$)}, |
| 1197 |
> |
label={sch:MultipoleAtomTypesBlock}] |
| 1198 |
> |
begin MultipoleAtomTypes |
| 1199 |
> |
// Euler angles are given in zxz convention in units of degrees. |
| 1200 |
> |
// |
| 1201 |
> |
// point dipoles: |
| 1202 |
> |
// name d phi theta psi dipole_moment |
| 1203 |
> |
DIP d 0.0 0.0 0.0 1.91 // dipole points along z-body axis |
| 1204 |
> |
// |
| 1205 |
> |
// point quadrupoles: |
| 1206 |
> |
// name q phi theta psi Qxx Qyy Qzz |
| 1207 |
> |
CO2 q 0.0 0.0 0.0 0.0 0.0 -0.430592 //quadrupole tensor has zz element |
| 1208 |
> |
// |
| 1209 |
> |
// Atoms with both dipole and quadrupole moments: |
| 1210 |
> |
// name dq phi theta psi dipole_moment Qxx Qyy Qzz |
| 1211 |
> |
SSD dq 0.0 0.0 0.0 2.35 -1.682 1.762 -0.08 |
| 1212 |
> |
end MultipoleAtomTypes |
| 1213 |
> |
\end{code} |
| 1214 |
> |
|
| 1215 |
> |
Specifying a MultipoleAtomType requires declaring how the |
| 1216 |
> |
electrostatic frame for the site is rotated relative to the body-fixed |
| 1217 |
> |
axes for that atom. The Euler angles $(\phi, \theta, \psi)$ for this |
| 1218 |
> |
rotation must be given, and then the dipole, quadrupole, or all of |
| 1219 |
> |
these moments are specified in the electrostatic frame. In OpenMD, |
| 1220 |
> |
the Euler angles are specified in the $zxz$ convention and are entered |
| 1221 |
> |
in units of degrees. Dipole moments are entered in units of Debye, |
| 1222 |
> |
and Quadrupole moments in units of Debye \AA. |
| 1223 |
> |
|
| 1224 |
> |
%\subsection{\label{section:ffGB}The FluctuatingChargeAtomTypes block} |
| 1225 |
> |
%\subsubsection{\label{section:ffPol}The PolarizableAtomTypes block} |
| 1226 |
> |
|
| 1227 |
> |
\subsection{\label{section:ffGB}The GayBerneAtomTypes block} |
| 1228 |
> |
|
| 1229 |
> |
The Gay-Berne potential has been widely used in the liquid crystal |
| 1230 |
> |
community to describe anisotropic phase |
| 1231 |
> |
behavior.~\cite{Gay:1981yu,Berne:1972pb,Kushick:1976xy,Luckhurst:1990fy,Cleaver:1996rt} |
| 1232 |
> |
The form of the Gay-Berne potential implemented in OpenMD was |
| 1233 |
> |
generalized by Cleaver {\it et al.} and is appropriate for dissimilar |
| 1234 |
> |
uniaxial ellipsoids.\cite{Cleaver:1996rt} The potential is constructed |
| 1235 |
> |
in the familiar form of the Lennard-Jones function using |
| 1236 |
> |
orientation-dependent $\sigma$ and $\epsilon$ parameters, |
| 1237 |
> |
\begin{equation*} |
| 1238 |
> |
V_{ij}({{\bf \hat u}_i}, {{\bf \hat u}_j}, {{\bf \hat |
| 1239 |
> |
r}_{ij}}) = 4\epsilon ({{\bf \hat u}_i}, {{\bf \hat u}_j}, |
| 1240 |
> |
{{\bf \hat r}_{ij}})\left[\left(\frac{\sigma_0}{r_{ij}-\sigma({{\bf \hat u |
| 1241 |
> |
}_i}, |
| 1242 |
> |
{{\bf \hat u}_j}, {{\bf \hat r}_{ij}})+\sigma_0}\right)^{12} |
| 1243 |
> |
-\left(\frac{\sigma_0}{r_{ij}-\sigma({{\bf \hat u}_i}, {{\bf \hat u}_j}, |
| 1244 |
> |
{{\bf \hat r}_{ij}})+\sigma_0}\right)^6\right] |
| 1245 |
> |
\label{eq:gb} |
| 1246 |
> |
\end{equation*} |
| 1247 |
> |
|
| 1248 |
> |
The range $(\sigma({\bf \hat{u}}_{i},{\bf \hat{u}}_{j},{\bf |
| 1249 |
> |
\hat{r}}_{ij}))$, and strength $(\epsilon({\bf \hat{u}}_{i},{\bf |
| 1250 |
> |
\hat{u}}_{j},{\bf \hat{r}}_{ij}))$ parameters |
| 1251 |
> |
are dependent on the relative orientations of the two ellipsoids (${\bf |
| 1252 |
> |
\hat{u}}_{i},{\bf \hat{u}}_{j}$) as well as the direction of the |
| 1253 |
> |
inter-ellipsoid separation (${\bf \hat{r}}_{ij}$). The shape and |
| 1254 |
> |
attractiveness of each ellipsoid is governed by a relatively small set |
| 1255 |
> |
of parameters: |
| 1256 |
> |
\begin{itemize} |
| 1257 |
> |
\item $d$: range parameter for the side-by-side (S) and cross (X) configurations |
| 1258 |
> |
\item $l$: range parameter for the end-to-end (E) configuration |
| 1259 |
> |
\item $\epsilon_X$: well-depth parameter for the cross (X) configuration |
| 1260 |
> |
\item $\epsilon_S$: well-depth parameter for the side-by-side (S) configuration |
| 1261 |
> |
\item $\epsilon_E$: well depth parameter for the end-to-end (E) configuration |
| 1262 |
> |
\item $dw$: The ``softness'' of the potential |
| 1263 |
> |
\end{itemize} |
| 1264 |
> |
Additionally, there are two universal paramters to govern the overall |
| 1265 |
> |
importance of the purely orientational ($\nu$) and the mixed |
| 1266 |
> |
orientational / translational ($\mu$) parts of strength of the |
| 1267 |
> |
interactions. These parameters have default or ``canonical'' values, |
| 1268 |
> |
but may be changed as a force field option: |
| 1269 |
> |
\begin{itemize} |
| 1270 |
> |
\item $\nu$: purely orientational part : defaults to 1 |
| 1271 |
> |
\item $\mu$: mixed orientational / translational part : defaults to |
| 1272 |
> |
2 |
| 1273 |
> |
\end{itemize} |
| 1274 |
> |
Further details of the potential are given |
| 1275 |
> |
elsewhere,\cite{Luckhurst:1990fy,Golubkov06,SunX._jp0762020} and an |
| 1276 |
> |
excellent overview of the computational methods that can be used to |
| 1277 |
> |
efficiently compute forces and torques for this potential can be found |
| 1278 |
> |
in Ref. \citealp{Golubkov06} |
| 1279 |
> |
|
| 1280 |
> |
\begin{code}[caption={[An example of a GayBerneAtomTypes block.] A |
| 1281 |
> |
simple example of a GayBerneAtomTypes block. Distances ($d$ and $l$) |
| 1282 |
> |
are given in \AA\ and energies ($\epsilon_X, \epsilon_S, \epsilon_E$) |
| 1283 |
> |
are in units of kcal/mol. $dw$ is unitless.}, |
| 1284 |
> |
label={sch:GayBerneAtomTypes}] |
| 1285 |
> |
begin GayBerneAtomTypes |
| 1286 |
> |
//Name d l eps_X eps_S eps_E dw |
| 1287 |
> |
GBlinear 2.8104 9.993 0.774729 0.774729 0.116839 1.0 |
| 1288 |
> |
GBC6H6 4.65 2.03 0.540 0.540 1.9818 0.6 |
| 1289 |
> |
GBCH3OH 2.55 3.18 0.542 0.542 0.55826 1.0 |
| 1290 |
> |
end GayBerneAtomTypes |
| 1291 |
> |
\end{code} |
| 1292 |
> |
|
| 1293 |
> |
\subsection{\label{section:ffSticky}The StickyAtomTypes block} |
| 1294 |
> |
|
| 1295 |
> |
One of the solvents that can be simulated by {\sc OpenMD} is the |
| 1296 |
> |
extended Soft Sticky Dipole (SSD/E) water model.\cite{fennell04} The |
| 1297 |
> |
original SSD was developed by Ichiye \emph{et |
| 1298 |
> |
al.}\cite{liu96:new_model} as a modified form of the hard-sphere |
| 1299 |
|
water model proposed by Bratko, Blum, and |
| 1300 |
|
Luzar.\cite{Bratko85,Bratko95} It consists of a single point dipole |
| 1301 |
|
with a Lennard-Jones core and a sticky potential that directs the |
| 1361 |
|
\label{fig:ssd} |
| 1362 |
|
\end{figure} |
| 1363 |
|
|
| 1158 |
– |
|
| 1364 |
|
Since SSD/E is a single-point {\it dipolar} model, the force |
| 1365 |
|
calculations are simplified significantly relative to the standard |
| 1366 |
|
{\it charged} multi-point models. In the original Monte Carlo |
| 1379 |
|
|
| 1380 |
|
Recent constant pressure simulations revealed issues in the original |
| 1381 |
|
SSD model that led to lower than expected densities at all target |
| 1382 |
< |
pressures.\cite{Ichiye03,fennell04} The default model in {\sc OpenMD} |
| 1383 |
< |
is therefore SSD/E, a density corrected derivative of SSD that |
| 1384 |
< |
exhibits improved liquid structure and transport behavior. If the use |
| 1385 |
< |
of a reaction field long-range interaction correction is desired, it |
| 1386 |
< |
is recommended that the parameters be modified to those of the SSD/RF |
| 1387 |
< |
model (an SSD variant parameterized for reaction field). These solvent |
| 1183 |
< |
parameters are listed and can be easily modified in the {\sc duff} |
| 1184 |
< |
force field file ({\tt DUFF.frc}). A table of the parameter values |
| 1185 |
< |
and the drawbacks and benefits of the different density corrected SSD |
| 1186 |
< |
models can be found in reference~\cite{fennell04}. |
| 1382 |
> |
pressures,\cite{Ichiye03,fennell04} so variants on the sticky |
| 1383 |
> |
potential can be specified by using one of a number of substitute atom |
| 1384 |
> |
types (see listing \ref{sch:StickyAtomTypes}). A table of the |
| 1385 |
> |
parameter values and the drawbacks and benefits of the different |
| 1386 |
> |
density corrected SSD models can be found in |
| 1387 |
> |
reference~\citealp{fennell04}. |
| 1388 |
|
|
| 1389 |
< |
\section{\label{section:WATER}The {\sc water} Force Field} |
| 1390 |
< |
|
| 1391 |
< |
In addition to the {\sc duff} force field's solvent description, a |
| 1392 |
< |
separate {\sc water} force field has been included for simulating most |
| 1393 |
< |
of the common rigid-body water models. This force field includes the |
| 1394 |
< |
simple and point-dipolar models (SSD, SSD1, SSD/E, SSD/RF, and DPD |
| 1395 |
< |
water), as well as the common charge-based models (SPC, SPC/E, TIP3P, |
| 1396 |
< |
TIP4P, and |
| 1397 |
< |
TIP5P).\cite{liu96:new_model,Ichiye03,fennell04,Marrink01,Berendsen81,Berendsen87,Jorgensen83,Mahoney00} |
| 1398 |
< |
In order to handle these models, charge-charge interactions were |
| 1399 |
< |
included in the force-loop: |
| 1400 |
< |
\begin{equation} |
| 1401 |
< |
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. |
| 1389 |
> |
\begin{code}[caption={[An example of a StickyAtomTypes block.] A |
| 1390 |
> |
simple example of a StickyAtomTypes block. Distances ($r_l$, $r_u$, |
| 1391 |
> |
$r_{l}'$ and $r_{u}'$) are given in \AA\ and energies ($v_0, v_{0}'$) |
| 1392 |
> |
are in units of kcal/mol. $w_0$ is unitless.}, |
| 1393 |
> |
label={sch:StickyAtomTypes}] |
| 1394 |
> |
begin StickyAtomTypes |
| 1395 |
> |
//name w0 v0 (kcal/mol) v0p rl (Ang) ru rlp rup |
| 1396 |
> |
SSD_E 0.07715 3.90 3.90 2.40 3.80 2.75 3.35 |
| 1397 |
> |
SSD_RF 0.07715 3.90 3.90 2.40 3.80 2.75 3.35 |
| 1398 |
> |
SSD 0.07715 3.7284 3.7284 2.75 3.35 2.75 4.0 |
| 1399 |
> |
SSD1 0.07715 3.6613 3.6613 2.75 3.35 2.75 4.0 |
| 1400 |
> |
end StickyAtomTypes |
| 1401 |
> |
\end{code} |
| 1402 |
|
|
| 1403 |
< |
\section{\label{section:eam}Embedded Atom Method} |
| 1403 |
> |
\section{\label{section::ffMetals}Metallic Atom Types} |
| 1404 |
|
|
| 1405 |
< |
{\sc OpenMD} implements a potential that describes bonding in |
| 1406 |
< |
transition metal |
| 1407 |
< |
systems.~\cite{Finnis84,Ercolessi88,Chen90,Qi99,Ercolessi02} This |
| 1408 |
< |
potential has an attractive interaction which models ``Embedding'' a |
| 1409 |
< |
positively charged pseudo-atom core in the electron density due to the |
| 1410 |
< |
free valance ``sea'' of electrons created by the surrounding atoms in |
| 1411 |
< |
the system. A pairwise part of the potential (which is primarily |
| 1412 |
< |
repulsive) describes the interaction of the positively charged metal |
| 1223 |
< |
core ions with one another. The Embedded Atom Method ({\sc |
| 1224 |
< |
eam})~\cite{Daw84,FBD86,johnson89,Lu97} has been widely adopted in the |
| 1225 |
< |
materials science community and has been included in {\sc OpenMD}. A |
| 1226 |
< |
good review of {\sc eam} and other formulations of metallic potentials |
| 1227 |
< |
was given by Voter.\cite{Voter:95} |
| 1228 |
< |
|
| 1229 |
< |
The {\sc eam} potential has the form: |
| 1405 |
> |
{\sc OpenMD} implements a number of related potentials that describe |
| 1406 |
> |
bonding in transition metals. These potentials have an attractive |
| 1407 |
> |
interaction which models ``Embedding'' a positively charged |
| 1408 |
> |
pseudo-atom core in the electron density due to the free valance |
| 1409 |
> |
``sea'' of electrons created by the surrounding atoms in the system. |
| 1410 |
> |
A pairwise part of the potential (which is primarily repulsive) |
| 1411 |
> |
describes the interaction of the positively charged metal core ions |
| 1412 |
> |
with one another. These potentials have the form: |
| 1413 |
|
\begin{equation} |
| 1414 |
|
V = \sum_{i} F_{i}\left[\rho_{i}\right] + \sum_{i} \sum_{j \neq i} |
| 1415 |
|
\phi_{ij}({\bf r}_{ij}) |
| 1426 |
|
transition metal potentials require two loops through the atom pairs |
| 1427 |
|
to compute the inter-atomic forces. |
| 1428 |
|
|
| 1429 |
< |
The pairwise portion of the potential, $\phi_{ij}$, is a primarily |
| 1430 |
< |
repulsive interaction between atoms $i$ and $j$. In the original |
| 1248 |
< |
formulation of {\sc eam}\cite{Daw84}, $\phi_{ij}$ was an entirely |
| 1249 |
< |
repulsive term; however later refinements to {\sc eam} allowed for |
| 1250 |
< |
more general forms for $\phi$.\cite{Daw89} The effective cutoff |
| 1251 |
< |
distance, $r_{{\text cut}}$ is the distance at which the values of |
| 1252 |
< |
$f(r)$ and $\phi(r)$ drop to zero for all atoms present in the |
| 1253 |
< |
simulation. In practice, this distance is fairly small, limiting the |
| 1254 |
< |
summations in the {\sc eam} equation to the few dozen atoms |
| 1255 |
< |
surrounding atom $i$ for both the density $\rho$ and pairwise $\phi$ |
| 1256 |
< |
interactions. |
| 1429 |
> |
The pairwise portion of the potential, $\phi_{ij}$, is usually a |
| 1430 |
> |
repulsive interaction between atoms $i$ and $j$. |
| 1431 |
|
|
| 1432 |
< |
In computing forces for alloys, mixing rules as outlined by |
| 1433 |
< |
Johnson~\cite{johnson89} are used to compute the heterogenous pair |
| 1434 |
< |
potential, |
| 1432 |
> |
\subsection{\label{section:ffEAM}The EAMAtomTypes block} |
| 1433 |
> |
The Embedded Atom Method ({\sc eam}) is one of the most widely-used |
| 1434 |
> |
potentials for transition |
| 1435 |
> |
metals.~\cite{Finnis84,Ercolessi88,Chen90,Qi99,Ercolessi02,Daw84,FBD86,johnson89,Lu97} |
| 1436 |
> |
It has been widely adopted in the materials science community and a |
| 1437 |
> |
good review of {\sc eam} and other formulations of metallic potentials |
| 1438 |
> |
was given by Voter.\cite{Voter:95} |
| 1439 |
> |
|
| 1440 |
> |
In the original formulation of {\sc eam}\cite{Daw84}, the pair |
| 1441 |
> |
potential, $\phi_{ij}$ was an entirely repulsive term; however later |
| 1442 |
> |
refinements to {\sc eam} allowed for more general forms for |
| 1443 |
> |
$\phi$.\cite{Daw89} The effective cutoff distance, $r_{{\text cut}}$ |
| 1444 |
> |
is the distance at which the values of $f(r)$ and $\phi(r)$ drop to |
| 1445 |
> |
zero for all atoms present in the simulation. In practice, this |
| 1446 |
> |
distance is fairly small, limiting the summations in the {\sc eam} |
| 1447 |
> |
equation to the few dozen atoms surrounding atom $i$ for both the |
| 1448 |
> |
density $\rho$ and pairwise $\phi$ interactions. |
| 1449 |
> |
|
| 1450 |
> |
In computing forces for alloys, OpenMD uses mixing rules outlined by |
| 1451 |
> |
Johnson~\cite{johnson89} to compute the heterogenous pair potential, |
| 1452 |
|
\begin{equation} |
| 1453 |
|
\label{eq:johnson} |
| 1454 |
|
\phi_{ab}(r)=\frac{1}{2}\left( |
| 1480 |
|
$\mbox{kcal mol}^{-1}$ as in the rest of the {\sc OpenMD} force field |
| 1481 |
|
files. |
| 1482 |
|
|
| 1483 |
< |
\section{\label{section:sc}The Sutton-Chen Force Field} |
| 1483 |
> |
\begin{code}[caption={[An example of a EAMAtomTypes block.] A |
| 1484 |
> |
simple example of a EAMAtomTypes block. Here the only data provided is |
| 1485 |
> |
the name of a {\tt funcfl} file which contains the raw data for spline |
| 1486 |
> |
interpolations for the density, functional, and pair potential.}, |
| 1487 |
> |
label={sch:EAMAtomTypes}] |
| 1488 |
> |
begin EAMAtomTypes |
| 1489 |
> |
Au Au.u3.funcfl |
| 1490 |
> |
Ag Ag.u3.funcfl |
| 1491 |
> |
Cu Cu.u3.funcfl |
| 1492 |
> |
Ni Ni.u3.funcfl |
| 1493 |
> |
Pd Pd.u3.funcfl |
| 1494 |
> |
Pt Pt.u3.funcfl |
| 1495 |
> |
end EAMAtomTypes |
| 1496 |
> |
\end{code} |
| 1497 |
|
|
| 1498 |
+ |
\subsection{\label{section:ffSC}The SuttonChenAtomTypes block} |
| 1499 |
+ |
|
| 1500 |
|
The Sutton-Chen ({\sc sc})~\cite{Chen90} potential has been used to |
| 1501 |
< |
study a wide range of phenomena in metals. Although it is similar in |
| 1502 |
< |
form to the {\sc eam} potential, the Sutton-Chen model takes on a |
| 1503 |
< |
simpler form, |
| 1501 |
> |
study a wide range of phenomena in metals. Although it has the same |
| 1502 |
> |
basic form as the {\sc eam} potential, the Sutton-Chen model requires |
| 1503 |
> |
a simpler set of parameters, |
| 1504 |
> |
\begin{equation} |
| 1505 |
> |
\label{eq:SCP1} |
| 1506 |
> |
U_{tot}=\sum _{i}\left[ \frac{1}{2}\sum _{j\neq |
| 1507 |
> |
i}\epsilon_{ij}V^{pair}_{ij}(r_{ij})-c_{i}\epsilon_{ii}\sqrt{\rho_{i}}\right] , |
| 1508 |
> |
\end{equation} |
| 1509 |
> |
where $V^{pair}_{ij}$ and $\rho_{i}$ are given by |
| 1510 |
> |
\begin{equation} |
| 1511 |
> |
\label{eq:SCP2} |
| 1512 |
> |
V^{pair}_{ij}(r)=\left( |
| 1513 |
> |
\frac{\alpha_{ij}}{r_{ij}}\right)^{n_{ij}} \hspace{1in} \rho_{i}=\sum_{j\neq i}\left( |
| 1514 |
> |
\frac{\alpha_{ij}}{r_{ij}}\right) ^{m_{ij}} |
| 1515 |
> |
\end{equation} |
| 1516 |
> |
|
| 1517 |
> |
$V^{pair}_{ij}$ is a repulsive pairwise potential that accounts for |
| 1518 |
> |
interactions of the pseudo-atom cores. The $\sqrt{\rho_i}$ term in |
| 1519 |
> |
Eq. (\ref{eq:SCP1}) is an attractive many-body potential that models |
| 1520 |
> |
the interactions between the valence electrons and the cores of the |
| 1521 |
> |
pseudo-atoms. $\epsilon_{ij}$, $\epsilon_{ii}$, $c_i$ and |
| 1522 |
> |
$\alpha_{ij}$ are parameters used to tune the potential for different |
| 1523 |
> |
transition metals. |
| 1524 |
> |
|
| 1525 |
> |
The {\sc sc} potential form has also been parameterized by Qi {\it et |
| 1526 |
> |
al.}\cite{Qi99} These parameters were obtained via empirical and {\it |
| 1527 |
> |
ab initio} calculations to match structural features of the FCC |
| 1528 |
> |
crystal. Interested readers are encouraged to consult reference |
| 1529 |
> |
\citealp{Qi99} for further details. |
| 1530 |
> |
|
| 1531 |
> |
\begin{code}[caption={[An example of a SCAtomTypes block.] A |
| 1532 |
> |
simple example of a SCAtomTypes block. Distances ($\alpha$) |
| 1533 |
> |
are given in \AA\ and energies ($\epsilon$) are (by convention) given in |
| 1534 |
> |
units of eV. These units must be specified in the {\tt Options} block |
| 1535 |
> |
using the keyword {\tt MetallicEnergyUnitScaling}. Without this {\tt |
| 1536 |
> |
Options} keyword, the default units for $\epsilon$ are kcal/mol. The |
| 1537 |
> |
other parameters, $m$, $n$, and $c$ are unitless.}, |
| 1538 |
> |
label={sch:SCAtomTypes}] |
| 1539 |
> |
begin SCAtomTypes |
| 1540 |
> |
// Name epsilon(eV) c m n alpha(angstroms) |
| 1541 |
> |
Ni 0.0073767 84.745 5.0 10.0 3.5157 |
| 1542 |
> |
Cu 0.0057921 84.843 5.0 10.0 3.6030 |
| 1543 |
> |
Rh 0.0024612 305.499 5.0 13.0 3.7984 |
| 1544 |
> |
Pd 0.0032864 148.205 6.0 12.0 3.8813 |
| 1545 |
> |
Ag 0.0039450 96.524 6.0 11.0 4.0691 |
| 1546 |
> |
Ir 0.0037674 224.815 6.0 13.0 3.8344 |
| 1547 |
> |
Pt 0.0097894 71.336 7.0 11.0 3.9163 |
| 1548 |
> |
Au 0.0078052 53.581 8.0 11.0 4.0651 |
| 1549 |
> |
Au2 0.0078052 53.581 8.0 11.0 4.0651 |
| 1550 |
> |
end SCAtomTypes |
| 1551 |
> |
\end{code} |
| 1552 |
> |
|
| 1553 |
> |
\section{\label{section::ffShortRange}Short Range Interactions} |
| 1554 |
> |
The internal structure of a molecule is usually specified in terms of |
| 1555 |
> |
a set of ``bonded'' terms in the potential energy function for |
| 1556 |
> |
molecule $I$, |
| 1557 |
> |
\begin{align*} |
| 1558 |
> |
V^{I}_{\text{Internal}} = & |
| 1559 |
> |
\sum_{r_{ij} \in I} V_{\text{bond}}(r_{ij}) |
| 1560 |
> |
+ \sum_{\theta_{ijk} \in I} V_{\text{bend}}(\theta_{ijk}) |
| 1561 |
> |
+ \sum_{\phi_{ijkl} \in I} V_{\text{torsion}}(\phi_{ijkl}) |
| 1562 |
> |
+ \sum_{\omega_{ijkl} \in I} V_{\text{inversion}}(\omega_{ijkl}) \\ |
| 1563 |
> |
& + \sum_{i \in I} \sum_{(j>i+4) \in I} |
| 1564 |
> |
\biggl[ V_{\text{dispersion}}(r_{ij}) + V_{\text{electrostatic}} |
| 1565 |
> |
(\mathbf{r}_{ij},\boldsymbol{\Omega}_{i},\boldsymbol{\Omega}_{j}) |
| 1566 |
> |
\biggr]. |
| 1567 |
> |
\label{eq:internalPotential} |
| 1568 |
> |
\end{align*} |
| 1569 |
> |
Here $V_{\text{bond}}, V_{\text{bend}}, |
| 1570 |
> |
V_{\text{torsion}},\mathrm{~and~} V_{\text{inversion}}$ represent the |
| 1571 |
> |
bond, bend, torsion, and inversion potentials for all |
| 1572 |
> |
topologically-connected sets of atoms within the molecule. Bonds are |
| 1573 |
> |
the primary way of specifying how the atoms are connected together to |
| 1574 |
> |
form the molecule (i.e. they define the molecular topology). The |
| 1575 |
> |
other short-range interactions may be specified explicitly in the |
| 1576 |
> |
molecule definition, or they may be deduced from bonding information. |
| 1577 |
> |
For example, bends can be implicitly deduced from two bonds which |
| 1578 |
> |
share an atom. Torsions can be deduced from two bends that share a |
| 1579 |
> |
bond. Inversion potentials are utilized primarily to enforce |
| 1580 |
> |
planarity around $sp^2$-hybridized sites, and these are specified with |
| 1581 |
> |
central atoms and satellites (or an atom with bonds to exactly three |
| 1582 |
> |
satellites). Non-bonded interactions are usually excluded for atom |
| 1583 |
> |
pairs that are involved in the same bond, bend, or torsion, but all |
| 1584 |
> |
other atom pairs are included in the standard non-bonded interactions. |
| 1585 |
> |
|
| 1586 |
> |
Bond lengths, angles, and torsions (dihedrals) are ``natural'' |
| 1587 |
> |
coordinates to treat molecular motion, as it is usually in these |
| 1588 |
> |
coordinates that most chemists understand the behavior of molecules. |
| 1589 |
> |
The bond lengths and angles are often considered ``hard'' degrees of |
| 1590 |
> |
freedom. That is, we can't deform them very much without a |
| 1591 |
> |
significant energetic penalty. On the other hand, dihedral angles or |
| 1592 |
> |
torsions are ``soft'' and typically undergo significant deformation |
| 1593 |
> |
under normal conditions. |
| 1594 |
> |
|
| 1595 |
> |
\subsection{\label{section:ffBond}The BondTypes block} |
| 1596 |
> |
|
| 1597 |
> |
Bonds are the primary way to specify how the atoms are connected |
| 1598 |
> |
together to form the molecule. In general, bonds exert strong |
| 1599 |
> |
restoring forces to keep the molecule compact. The bond energy |
| 1600 |
> |
functions are relatively simple functions of the distance between two |
| 1601 |
> |
atomic sites, |
| 1602 |
> |
\begin{equation} |
| 1603 |
> |
b = \left| \vec{r}_{ij} \right| = \sqrt{(x_j - x_i)^2 + (y_j - y_i)^2 |
| 1604 |
> |
+ (z_j - z_i)^2}. |
| 1605 |
> |
\end{equation} |
| 1606 |
> |
All BondTypes must specify two AtomType names ($i$ and $j$) that |
| 1607 |
> |
describe when that bond should be applied, as well as an equilibrium |
| 1608 |
> |
bond length, $b_{ij}^0$, in units of \AA. The most common forms for |
| 1609 |
> |
bonding potentials are {\tt Harmonic} bonds, |
| 1610 |
> |
\begin{equation} |
| 1611 |
> |
V_{\text{bond}}(b) = \frac{k_{ij}}{2} \left(b - |
| 1612 |
> |
b_{ij}^0 \right)^2 |
| 1613 |
> |
\end{equation} |
| 1614 |
> |
and {\tt Morse} bonds, |
| 1615 |
> |
\begin{equation} |
| 1616 |
> |
V_{\text{bond}}(b) = D_{ij} \left[ 1 - e^{-\beta_{ij} (b - b_{ij}^0)} \right]^2 |
| 1617 |
> |
\end{equation} |
| 1618 |
> |
|
| 1619 |
> |
\begin{figure}[h] |
| 1620 |
> |
\centering |
| 1621 |
> |
\includegraphics[width=2.5in]{bond.pdf} |
| 1622 |
> |
\caption[Bond coordinates]{The coordinate describing a |
| 1623 |
> |
a bond between atoms $i$ and $j$ is $|r_{ij}|$, the length of the |
| 1624 |
> |
$\vec{r}_{ij}$ vector. } |
| 1625 |
> |
\label{fig:bond} |
| 1626 |
> |
\end{figure} |
| 1627 |
> |
|
| 1628 |
> |
OpenMD can also simulate some less common types of bond potentials, |
| 1629 |
> |
including {\tt Fixed} bonds (which are constrained to be at a |
| 1630 |
> |
specified bond length), |
| 1631 |
> |
\begin{equation} |
| 1632 |
> |
V_{\text{bond}}(b) = 0. |
| 1633 |
> |
\end{equation} |
| 1634 |
> |
{\tt Cubic} bonds include terms to model anharmonicity, |
| 1635 |
> |
\begin{equation} |
| 1636 |
> |
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, |
| 1637 |
> |
\end{equation} |
| 1638 |
> |
and {\tt Quartic} bonds, include another term in the Taylor |
| 1639 |
> |
expansion around $b_{ij}^0$, |
| 1640 |
> |
\begin{equation} |
| 1641 |
> |
V_{\text{bond}}(b) = K_4 (b - b_{ij}^0)^4 + K_3 (b - b_{ij}^0)^3 + |
| 1642 |
> |
K_2 (b - b_{ij}^0)^2 + K_1 (b - b_{ij}^0) + K_0, |
| 1643 |
> |
\end{equation} |
| 1644 |
> |
can also be simulated. Note that the factor of $1/2$ that is included |
| 1645 |
> |
in the {\tt Harmonic} bond type force constant is {\it not} present in |
| 1646 |
> |
either the {\tt Cubic} or {\tt Quartic} bond types. |
| 1647 |
> |
|
| 1648 |
> |
{\tt Polynomial} bonds which can have any number of terms, |
| 1649 |
> |
\begin{equation} |
| 1650 |
> |
V_{\text{bond}}(b) = \sum_n K_n (b - b_{ij}^0)^n. |
| 1651 |
> |
\end{equation} |
| 1652 |
> |
can also be specified by giving a sequence of integer ($n$) and force |
| 1653 |
> |
constant ($K_n$) pairs. |
| 1654 |
> |
|
| 1655 |
> |
The order of terms in the BondTypes block is: |
| 1656 |
> |
\begin{itemize} |
| 1657 |
> |
\item {\tt AtomType} 1 |
| 1658 |
> |
\item {\tt AtomType} 2 |
| 1659 |
> |
\item {\tt BondType} (one of {\tt Harmonic}, {\tt Morse}, {\tt Fixed}, {\tt |
| 1660 |
> |
Cubic}, {\tt Quartic}, or {\tt Polynomial}) |
| 1661 |
> |
\item $b_{ij}^0$, the equilibrium bond length in \AA |
| 1662 |
> |
\item any other parameters required by the {\tt BondType} |
| 1663 |
> |
\end{itemize} |
| 1664 |
> |
|
| 1665 |
> |
\begin{code}[caption={[An example of a BondTypes block.] A |
| 1666 |
> |
simple example of a BondTypes block. Distances ($b_0$) |
| 1667 |
> |
are given in \AA\ and force constants are given in |
| 1668 |
> |
units so that when multiplied by the correct power of distance they |
| 1669 |
> |
return energies in kcal/mol. For example $k$ for a Harmonic bond is |
| 1670 |
> |
in units of kcal/mol/\AA$^2$.}, |
| 1671 |
> |
label={sch:BondTypes}] |
| 1672 |
> |
begin BondTypes |
| 1673 |
> |
//Atom1 Atom2 Harmonic b0 k (kcal/mol/A^2) |
| 1674 |
> |
CH2 CH2 Harmonic 1.526 260 |
| 1675 |
> |
//Atom1 Atom2 Morse b0 D beta (A^-1) |
| 1676 |
> |
CN NC Morse 1.157437 212.95 2.5802 |
| 1677 |
> |
//Atom1 Atom2 Fixed b0 |
| 1678 |
> |
CT HC Fixed 1.09 |
| 1679 |
> |
//Atom1 Atom2 Cubic b0 K3 K2 K1 K0 |
| 1680 |
> |
//Atom1 Atom2 Quartic b0 K4 K3 K2 K1 K0 |
| 1681 |
> |
//Atom1 Atom2 Polynomial b0 n Kn [m Km] |
| 1682 |
> |
end BondTypes |
| 1683 |
> |
\end{code} |
| 1684 |
> |
|
| 1685 |
> |
There are advantages and disadvantages of all of the different types |
| 1686 |
> |
of bonds, but specific simulation tasks may call for specific |
| 1687 |
> |
behaviors. |
| 1688 |
> |
|
| 1689 |
> |
\subsection{\label{section:ffBend}The BendTypes block} |
| 1690 |
> |
The equilibrium geometries and energy functions for bending motions in |
| 1691 |
> |
a molecule are strongly dependent on the bonding environment of the |
| 1692 |
> |
central atomic site. For example, different types of hybridized |
| 1693 |
> |
carbon centers require different bending angles and force constants to |
| 1694 |
> |
describe the local geometry. |
| 1695 |
> |
|
| 1696 |
> |
The bending potential energy functions used in most force fields are |
| 1697 |
> |
often simple functions of the angle between two bonds, |
| 1698 |
> |
\begin{equation} |
| 1699 |
> |
\theta_{ijk} = \cos^{-1} \left(\frac{\vec{r}_{ji} \cdot |
| 1700 |
> |
\vec{r}_{jk}}{\left| \vec{r}_{ji} \right| \left| \vec{r}_{jk} |
| 1701 |
> |
\right|} \right) |
| 1702 |
> |
\end{equation} |
| 1703 |
> |
Here atom $j$ is the central atom that is bonded to two partners $i$ |
| 1704 |
> |
and $k$. |
| 1705 |
> |
|
| 1706 |
> |
\begin{figure}[h] |
| 1707 |
> |
\centering |
| 1708 |
> |
\includegraphics[width=3.5in]{bend.pdf} |
| 1709 |
> |
\caption[Bend angle coordinates]{The coordinate describing a bend |
| 1710 |
> |
between atoms $i$, $j$, and $k$ is the angle $\theta_{ijk} = |
| 1711 |
> |
\cos^{-1} \left(\hat{r}_{ji} \cdot \hat{r}_{jk}\right)$ where $\hat{r}_{ji}$ is |
| 1712 |
> |
the unit vector between atoms $j$ and $i$. } |
| 1713 |
> |
\label{fig:bend} |
| 1714 |
> |
\end{figure} |
| 1715 |
> |
|
| 1716 |
> |
|
| 1717 |
> |
All BendTypes must specify three AtomType names ($i$, $j$ and $k$) |
| 1718 |
> |
that describe when that bend potential should be applied, as well as |
| 1719 |
> |
an equilibrium bending angle, $\theta_{ijk}^0$, in units of |
| 1720 |
> |
degrees. The most common forms for bending potentials are {\tt |
| 1721 |
> |
Harmonic} bends, |
| 1722 |
> |
\begin{equation} |
| 1723 |
> |
V_{\text{bend}}(\theta_{ijk}) = \frac{k_{ijk}}{2}( \theta_{ijk} - \theta_{ijk}^0 |
| 1724 |
> |
)^2, \label{eq:bendPot} |
| 1725 |
> |
\end{equation} |
| 1726 |
> |
where $k_{ijk}$ is the force constant which determines the strength of |
| 1727 |
> |
the harmonic bend. {\tt UreyBradley} bends utilize an additional 1-3 |
| 1728 |
> |
bond-type interaction in addition to the harmonic bending potential: |
| 1729 |
> |
\begin{equation} |
| 1730 |
> |
V_{\text{bend}}(\vec{r}_i , \vec{r}_j, \vec{r}_k) |
| 1731 |
> |
=\frac{k_{ijk}}{2}( \theta_{ijk} - \theta_{ijk}^0)^2 |
| 1732 |
> |
+ \frac{k_{ub}}{2}( r_{ik} - s_0 )^2. \label{eq:ubBend} |
| 1733 |
> |
\end{equation} |
| 1734 |
> |
|
| 1735 |
> |
A {\tt Cosine} bend is a variant on the harmonic bend which utilizes |
| 1736 |
> |
the cosine of the angle instead of the angle itself, |
| 1737 |
> |
\begin{equation} |
| 1738 |
> |
V_{\text{bend}}(\theta_{ijk}) = \frac{k_{ijk}}{2}( \cos\theta_{ijk} - |
| 1739 |
> |
\cos \theta_{ijk}^0 )^2. \label{eq:cosBend} |
| 1740 |
> |
\end{equation} |
| 1741 |
> |
|
| 1742 |
> |
OpenMD can also simulate some less common types of bend potentials, |
| 1743 |
> |
including {\tt Cubic} bends, which include terms to model anharmonicity, |
| 1744 |
> |
\begin{equation} |
| 1745 |
> |
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, |
| 1746 |
> |
\end{equation} |
| 1747 |
> |
and {\tt Quartic} bends, which include another term in the Taylor |
| 1748 |
> |
expansion around $\theta_{ijk}^0$, |
| 1749 |
> |
\begin{equation} |
| 1750 |
> |
V_{\text{bend}}(\theta_{ijk}) = K_4 (\theta_{ijk} - \theta_{ijk}^0)^4 + K_3 (\theta_{ijk} - \theta_{ijk}^0)^3 + |
| 1751 |
> |
K_2 (\theta_{ijk} - \theta_{ijk}^0)^2 + K_1 (\theta_{ijk} - |
| 1752 |
> |
\theta_{ijk}^0) + K_0, |
| 1753 |
> |
\end{equation} |
| 1754 |
> |
can also be simulated. Note that the factor of $1/2$ that is included |
| 1755 |
> |
in the {\tt Harmonic} bend type force constant is {\it not} present in |
| 1756 |
> |
either the {\tt Cubic} or {\tt Quartic} bend types. |
| 1757 |
> |
|
| 1758 |
> |
{\tt Polynomial} bends which can have any number of terms, |
| 1759 |
> |
\begin{equation} |
| 1760 |
> |
V_{\text{bend}}(\theta_{ijk}) = \sum_n K_n (\theta_{ijk} - \theta_{ijk}^0)^n. |
| 1761 |
> |
\end{equation} |
| 1762 |
> |
can also be specified by giving a sequence of integer ($n$) and force |
| 1763 |
> |
constant ($K_n$) pairs. |
| 1764 |
> |
|
| 1765 |
> |
The order of terms in the BendTypes block is: |
| 1766 |
> |
\begin{itemize} |
| 1767 |
> |
\item {\tt AtomType} 1 |
| 1768 |
> |
\item {\tt AtomType} 2 (this is the central atom) |
| 1769 |
> |
\item {\tt AtomType} 3 |
| 1770 |
> |
\item {\tt BendType} (one of {\tt Harmonic}, {\tt UreyBradley}, {\tt |
| 1771 |
> |
Cosine}, {\tt Cubic}, {\tt Quartic}, or {\tt Polynomial}) |
| 1772 |
> |
\item $\theta_{ijk}^0$, the equilibrium bending angle in degrees. |
| 1773 |
> |
\item any other parameters required by the {\tt BendType} |
| 1774 |
> |
\end{itemize} |
| 1775 |
> |
|
| 1776 |
> |
\begin{code}[caption={[An example of a BendTypes block.] A |
| 1777 |
> |
simple example of a BendTypes block. By convention, equilibrium angles |
| 1778 |
> |
($\theta_0$) are given in degrees but force constants are given in |
| 1779 |
> |
units so that when multiplied by the correct power of angle (in |
| 1780 |
> |
radians) they return energies in kcal/mol. For example $k$ for a |
| 1781 |
> |
Harmonic bend is in units of kcal/mol/radians$^2$.}, |
| 1782 |
> |
label={sch:BendTypes}] |
| 1783 |
> |
begin BendTypes |
| 1784 |
> |
//Atom1 Atom2 Atom3 Harmonic theta0(deg) Ktheta(kcal/mol/radians^2) |
| 1785 |
> |
CT CT CT Harmonic 109.5 80.000000 |
| 1786 |
> |
CH2 CH CH2 Harmonic 112.0 117.68 |
| 1787 |
> |
CH3 CH2 SH Harmonic 96.0 67.220 |
| 1788 |
> |
//UreyBradley |
| 1789 |
> |
//Atom1 Atom2 Atom3 UreyBradley theta0 Ktheta s0 Kub |
| 1790 |
> |
//Cosine |
| 1791 |
> |
//Atom1 Atom2 Atom3 Cosine theta0 Ktheta(kcal/mol) |
| 1792 |
> |
//Cubic |
| 1793 |
> |
//Atom1 Atom2 Atom3 Cubic theta0 K3 K2 K1 K0 |
| 1794 |
> |
//Quartic |
| 1795 |
> |
//Atom1 Atom2 Atom3 Quartic theta0 K4 K3 K2 K1 K0 |
| 1796 |
> |
//Polynomial |
| 1797 |
> |
//Atom1 Atom2 Atom3 Polynomial theta0 n Kn [m Km] |
| 1798 |
> |
end BendTypes |
| 1799 |
> |
\end{code} |
| 1800 |
> |
|
| 1801 |
> |
Note that the parameters for a particular bend type are the same for |
| 1802 |
> |
any bending triplet of the same atomic types (in the same or reversed |
| 1803 |
> |
order). Changing the AtomType in the Atom2 position will change the |
| 1804 |
> |
matched bend types in the force field. |
| 1805 |
> |
|
| 1806 |
> |
\subsection{\label{section:ffTorsion}The TorsionTypes block} |
| 1807 |
> |
The torsion potential for rotation of groups around a central bond can |
| 1808 |
> |
often be represented with various cosine functions. For two |
| 1809 |
> |
tetrahedral ($sp^3$) carbons connected by a single bond, the torsion |
| 1810 |
> |
potential might be |
| 1811 |
> |
\begin{equation*} |
| 1812 |
> |
V_{\text{torsion}} = \frac{v}{2} \left[ 1 + \cos( 3 \phi ) \right] |
| 1813 |
> |
\end{equation*} |
| 1814 |
> |
where $v$ is the barrier for going from {\em staggered} $\rightarrow$ |
| 1815 |
> |
{\em eclipsed} conformations, while for $sp^2$ carbons connected by a |
| 1816 |
> |
double bond, the potential might be |
| 1817 |
> |
\begin{equation*} |
| 1818 |
> |
V_{\text{torsion}} = \frac{w}{2} \left[ 1 - \cos( 2 \phi ) \right] |
| 1819 |
> |
\end{equation*} |
| 1820 |
> |
where $w$ is the barrier for going from {\em cis} $\rightarrow$ {\em |
| 1821 |
> |
trans} conformations. |
| 1822 |
> |
|
| 1823 |
> |
A general torsion potentials can be represented as a cosine series of |
| 1824 |
> |
the form: |
| 1825 |
> |
\begin{equation} |
| 1826 |
> |
V_{\text{torsion}}(\phi_{ijkl}) = c_1[1 + \cos \phi_{ijkl}] |
| 1827 |
> |
+ c_2[1 - \cos(2\phi_{ijkl})] |
| 1828 |
> |
+ c_3[1 + \cos(3\phi_{ijkl})], |
| 1829 |
> |
\label{eq:origTorsionPot} |
| 1830 |
> |
\end{equation} |
| 1831 |
> |
where the angle $\phi_{ijkl}$ is defined |
| 1832 |
> |
\begin{equation} |
| 1833 |
> |
\cos\phi_{ijkl} = (\hat{\mathbf{r}}_{ij} \times \hat{\mathbf{r}}_{jk}) \cdot |
| 1834 |
> |
(\hat{\mathbf{r}}_{jk} \times \hat{\mathbf{r}}_{kl}). |
| 1835 |
> |
\label{eq:torsPhi} |
| 1836 |
> |
\end{equation} |
| 1837 |
> |
Here, $\hat{\mathbf{r}}_{\alpha\beta}$ are the set of unit bond |
| 1838 |
> |
vectors between atoms $i$, $j$, $k$, and $l$. Note that some force |
| 1839 |
> |
fields define the zero of the $\phi_{ijkl}$ angle when atoms $i$ and |
| 1840 |
> |
$l$ are in the {\em trans} configuration, while most define the zero |
| 1841 |
> |
angle for when $i$ and $l$ are in the fully eclipsed orientation. The |
| 1842 |
> |
behavior of the torsion parser can be altered with the {\tt |
| 1843 |
> |
TorsionAngleConvention} keyword in the Options block. The default |
| 1844 |
> |
behavior is {\tt "180\_is\_trans"} while the opposite behavior can be |
| 1845 |
> |
invoked by setting this keyword to {\tt "0\_is\_trans"}. |
| 1846 |
> |
|
| 1847 |
> |
\begin{figure}[h] |
| 1848 |
> |
\centering |
| 1849 |
> |
\includegraphics[width=4.5in]{torsion.pdf} |
| 1850 |
> |
\caption[Torsion or dihedral angle coordinates]{The coordinate |
| 1851 |
> |
describing a torsion between atoms $i$, $j$, $k$, and $l$ is the |
| 1852 |
> |
dihedral angle $\phi_{ijkl}$ which measures the relative rotation of |
| 1853 |
> |
the two terminal atoms around the axis defined by the middle bond. } |
| 1854 |
> |
\label{fig:torsion} |
| 1855 |
> |
\end{figure} |
| 1856 |
> |
|
| 1857 |
> |
For computational efficiency, OpenMD recasts torsion potential in the |
| 1858 |
> |
method of {\sc charmm},\cite{Brooks83} in which the angle series is |
| 1859 |
> |
converted to a power series of the form: |
| 1860 |
> |
\begin{equation} |
| 1861 |
> |
V_{\text{torsion}}(\phi_{ijkl}) = |
| 1862 |
> |
k_3 \cos^3 \phi_{ijkl} + k_2 \cos^2 \phi_{ijkl} + k_1 \cos \phi_{ijkl} + k_0, |
| 1863 |
> |
\label{eq:torsionPot} |
| 1864 |
> |
\end{equation} |
| 1865 |
> |
where: |
| 1866 |
> |
\begin{align*} |
| 1867 |
> |
k_0 &= c_1 + 2 c_2 + c_3, \\ |
| 1868 |
> |
k_1 &= c_1 - 3c_3, \\ |
| 1869 |
> |
k_2 &= - 2 c_2, \\ |
| 1870 |
> |
k_3 &= 4 c_3. |
| 1871 |
> |
\end{align*} |
| 1872 |
> |
By recasting the potential as a power series, repeated trigonometric |
| 1873 |
> |
evaluations are avoided during the calculation of the potential |
| 1874 |
> |
energy. |
| 1875 |
> |
|
| 1876 |
> |
Using this framework, OpenMD implements a variety of different |
| 1877 |
> |
potential energy functions for torsions: |
| 1878 |
> |
\begin{itemize} |
| 1879 |
> |
\item {\tt Cubic}: |
| 1880 |
> |
\begin{equation*} |
| 1881 |
> |
V_{\text{torsion}}(\phi) = |
| 1882 |
> |
k_3 \cos^3 \phi + k_2 \cos^2 \phi + k_1 \cos \phi + k_0, |
| 1883 |
> |
\end{equation*} |
| 1884 |
> |
\item {\tt Quartic}: |
| 1885 |
> |
\begin{equation*} |
| 1886 |
> |
V_{\text{torsion}}(\phi) = k_4 \cos^4 \phi + |
| 1887 |
> |
k_3 \cos^3 \phi + k_2 \cos^2 \phi + k_1 \cos \phi + k_0, |
| 1888 |
> |
\end{equation*} |
| 1889 |
> |
\item {\tt Polynomial}: |
| 1890 |
> |
\begin{equation*} |
| 1891 |
> |
V_{\text{torsion}}(\phi) = \sum_n k_n \cos^n \phi , |
| 1892 |
> |
\end{equation*} |
| 1893 |
> |
\item {\tt Charmm}: |
| 1894 |
> |
\begin{equation*} |
| 1895 |
> |
V_{\text{torsion}}(\phi) = \sum_n K_n \left( 1 + cos(n |
| 1896 |
> |
\phi - \delta_n) \right), |
| 1897 |
> |
\end{equation*} |
| 1898 |
> |
\item {\tt Opls}: |
| 1899 |
> |
\begin{equation*} |
| 1900 |
> |
V_{\text{torsion}}(\phi) = \frac{1}{2} \left(v_1 (1 + \cos \phi) \right) |
| 1901 |
> |
+ v_2 (1 - \cos 2 \phi) + v_3 (1 + \cos 3 \phi), |
| 1902 |
> |
\end{equation*} |
| 1903 |
> |
\item {\tt Trappe}:\cite{Siepmann1998} |
| 1904 |
> |
\begin{equation*} |
| 1905 |
> |
V_{\text{torsion}}(\phi) = c_0 + c_1 (1 + \cos \phi) + c_2 (1 - \cos 2 \phi) + |
| 1906 |
> |
c_3 (1 + \cos 3 \phi), |
| 1907 |
> |
\end{equation*} |
| 1908 |
> |
\item {\tt Harmonic}: |
| 1909 |
> |
\begin{equation*} |
| 1910 |
> |
V_{\text{torsion}}(\phi) = \frac{d_0}{2} \left(\phi - \phi^0\right). |
| 1911 |
> |
\end{equation*} |
| 1912 |
> |
\end{itemize} |
| 1913 |
> |
|
| 1914 |
> |
Most torsion types don't require specific angle information in the |
| 1915 |
> |
parameters since they are typically expressed in cosine polynomials. |
| 1916 |
> |
{\tt Charmm} and {\tt Harmonic} torsions are a bit different. {\tt |
| 1917 |
> |
Charmm} torsion types require a set of phase angles, $\delta_n$ that |
| 1918 |
> |
are expressed in degrees, and associated periodicity numbers, $n$. |
| 1919 |
> |
{\tt Harmonic} torsions have an equilibrium torsion angle, $\phi_0$ |
| 1920 |
> |
that is measured in degrees, while $d_0$ has units of |
| 1921 |
> |
kcal/mol/degrees$^2$. All other torsion parameters are measured in |
| 1922 |
> |
units of kcal/mol. |
| 1923 |
> |
|
| 1924 |
> |
\begin{code}[caption={[An example of a TorsionTypes block.] A |
| 1925 |
> |
simple example of a TorsionTypes block. Energy constants are given in |
| 1926 |
> |
kcal / mol, and when required by the form, $\delta$ angles are given |
| 1927 |
> |
in degrees.}, |
| 1928 |
> |
label={sch:TorsionTypes}] |
| 1929 |
> |
begin TorsionTypes |
| 1930 |
> |
//Cubic |
| 1931 |
> |
//Atom1 Atom2 Atom3 Atom4 Cubic k3 k2 k1 k0 |
| 1932 |
> |
CH2 CH2 CH2 CH2 Cubic 5.9602 -0.2568 -3.802 2.1586 |
| 1933 |
> |
CH2 CH CH CH2 Cubic 3.3254 -0.4215 -1.686 1.1661 |
| 1934 |
> |
//Trappe |
| 1935 |
> |
//Atom1 Atom2 Atom3 Atom4 Trappe c0 c1 c2 c3 |
| 1936 |
> |
CH3 CH2 CH2 SH Trappe 0.10507 -0.10342 0.03668 0.60874 |
| 1937 |
> |
//Charmm |
| 1938 |
> |
//Atom1 Atom2 Atom3 Atom4 Charmm Kchi n delta [Kchi n delta] |
| 1939 |
> |
CT CT CT C Charmm 0.156 3 0.00 |
| 1940 |
> |
OH CT CT OH Charmm 0.144 3 0.00 1.175 2 0 |
| 1941 |
> |
HC CT CM CM Charmm 1.150 1 0.00 0.38 3 180 |
| 1942 |
> |
//Quartic |
| 1943 |
> |
//Atom1 Atom2 Atom3 Atom4 Quartic k4 k3 k2 k1 k0 |
| 1944 |
> |
//Polynomial |
| 1945 |
> |
//Atom1 Atom2 Atom3 Atom4 Polynomial n Kn [m Km] |
| 1946 |
> |
S CH2 CH2 C Polynomial 0 2.218 1 2.905 2 -3.136 3 -0.7313 4 6.272 5 -7.528 |
| 1947 |
> |
end TorsionTypes |
| 1948 |
> |
\end{code} |
| 1949 |
> |
|
| 1950 |
> |
Note that the parameters for a particular torsion type are the same |
| 1951 |
> |
for any torsional quartet of the same atomic types (in the same or |
| 1952 |
> |
reversed order). |
| 1953 |
> |
|
| 1954 |
> |
\subsection{\label{section:ffInversion}The InversionTypes block} |
| 1955 |
> |
|
| 1956 |
> |
Inversion potentials are often added to force fields to enforce |
| 1957 |
> |
planarity around $sp^2$-hybridized carbons or to correct vibrational |
| 1958 |
> |
frequencies for umbrella-like vibrational modes for central atoms |
| 1959 |
> |
bonded to exactly three satellite groups. |
| 1960 |
> |
|
| 1961 |
> |
In OpenMD's version of an inversion, the central atom is entered {\it |
| 1962 |
> |
first} in each line in the {\tt InversionTypes} block. Note that |
| 1963 |
> |
this is quite different than how other codes treat Improper torsional |
| 1964 |
> |
potentials to mimic inversion behavior. In some other widely-used |
| 1965 |
> |
simulation packages, the central atom is treated as atom 3 in a |
| 1966 |
> |
standard torsion form: |
| 1967 |
> |
\begin{itemize} |
| 1968 |
> |
\item OpenMD: I - (J - K - L) (e.g. I is $sp^2$ hybridized carbon) |
| 1969 |
> |
\item AMBER: I - J - K - L (e.g. K is $sp^2$ hybridized carbon) |
| 1970 |
> |
\end{itemize} |
| 1971 |
> |
|
| 1972 |
> |
The inversion angle itself is defined as: |
| 1973 |
|
\begin{equation} |
| 1974 |
< |
\label{eq:SCP1} |
| 1975 |
< |
U_{tot}=\sum _{i}\left[ \frac{1}{2}\sum _{j\neq |
| 1976 |
< |
i}D_{ij}V^{pair}_{ij}(r_{ij})-c_{i}D_{ii}\sqrt{\rho_{i}}\right] , |
| 1974 |
> |
\cos\omega_{i-jkl} = \left(\hat{\mathbf{r}}_{il} \times |
| 1975 |
> |
\hat{\mathbf{r}}_{ij}\right)\cdot\left( \hat{\mathbf{r}}_{il} \times |
| 1976 |
> |
\hat{\mathbf{r}}_{ik}\right) |
| 1977 |
|
\end{equation} |
| 1978 |
< |
where $V^{pair}_{ij}$ and $\rho_{i}$ are given by |
| 1979 |
< |
\begin{equation} |
| 1980 |
< |
\label{eq:SCP2} |
| 1981 |
< |
V^{pair}_{ij}(r)=\left( |
| 1982 |
< |
\frac{\alpha_{ij}}{r_{ij}}\right)^{n_{ij}}, \rho_{i}=\sum_{j\neq i}\left( |
| 1308 |
< |
\frac{\alpha_{ij}}{r_{ij}}\right) ^{m_{ij}} |
| 1309 |
< |
\end{equation} |
| 1978 |
> |
Here, $\hat{\mathbf{r}}_{\alpha\beta}$ are the set of unit bond |
| 1979 |
> |
vectors between the central atoms $i$, and the satellite atoms $j$, |
| 1980 |
> |
$k$, and $l$. Note that other definitions of inversion angles are |
| 1981 |
> |
possible, so users are encouraged to be particularly careful when |
| 1982 |
> |
converting other force field files for use with OpenMD. |
| 1983 |
|
|
| 1984 |
< |
$V^{pair}_{ij}$ is a repulsive pairwise potential that accounts for |
| 1985 |
< |
interactions of the pseudo-atom cores. The $\sqrt{\rho_i}$ term in |
| 1986 |
< |
Eq. (\ref{eq:SCP1}) is an attractive many-body potential that models |
| 1987 |
< |
the interactions between the valence electrons and the cores of the |
| 1988 |
< |
pseudo-atoms. $D_{ij}$, $D_{ii}$, $c_i$ and $\alpha_{ij}$ are |
| 1989 |
< |
parameters used to tune the potential for different transition |
| 1990 |
< |
metals. |
| 1984 |
> |
There are many common ways to create planarity or umbrella behavior in |
| 1985 |
> |
a potential energy function, and OpenMD implements a number of the |
| 1986 |
> |
more common functions: |
| 1987 |
> |
\begin{itemize} |
| 1988 |
> |
\item {\tt ImproperCosine}: |
| 1989 |
> |
\begin{equation*} |
| 1990 |
> |
V_{\text{torsion}}(\omega) = \sum_n \frac{K_n}{2} \left( 1 + cos(n |
| 1991 |
> |
\omega - \delta_n) \right), |
| 1992 |
> |
\end{equation*} |
| 1993 |
> |
\item {\tt AmberImproper}: |
| 1994 |
> |
\begin{equation*} |
| 1995 |
> |
V_{\text{torsion}}(\omega) = \frac{v}{2} (1 - \cos\left(2 \left(\omega - \omega_0\right)\right), |
| 1996 |
> |
\end{equation*} |
| 1997 |
> |
\item {\tt Harmonic}: |
| 1998 |
> |
\begin{equation*} |
| 1999 |
> |
V_{\text{torsion}}(\omega) = \frac{d}{2} \left(\omega - \omega_0\right). |
| 2000 |
> |
\end{equation*} |
| 2001 |
> |
\end{itemize} |
| 2002 |
> |
\begin{code}[caption={[An example of an InversionTypes block.] A |
| 2003 |
> |
simple example of a InversionTypes block. Angles ($\delta_n$ and |
| 2004 |
> |
$\omega_0$) angles are given in degrees, while energy parameters ($v, |
| 2005 |
> |
K_n$) are given in kcal / mol. The Harmonic Inversion type has a |
| 2006 |
> |
force constant that must be given in kcal/mol/degrees$^2$.}, |
| 2007 |
> |
label={sch:InversionTypes}] |
| 2008 |
> |
begin InversionTypes |
| 2009 |
> |
//Harmonic |
| 2010 |
> |
//Atom1 Atom2 Atom3 Atom4 Harmonic d(kcal/mol/deg^2) omega0 |
| 2011 |
> |
RCHar3 X X X Harmonic 1.21876e-2 180.0 |
| 2012 |
> |
//AmberImproper |
| 2013 |
> |
//Atom1 Atom2 Atom3 Atom4 AmberImproper v(kcal/mol) |
| 2014 |
> |
C CT N O AmberImproper 10.500000 |
| 2015 |
> |
CA CA CA CT AmberImproper 1.100000 |
| 2016 |
> |
//ImproperCosine |
| 2017 |
> |
//Atom1 Atom2 Atom3 Atom4 ImproperCosine Kn n delta_n [Kn n delta_n] |
| 2018 |
> |
end InversionTypes |
| 2019 |
> |
\end{code} |
| 2020 |
|
|
| 2021 |
< |
The {\sc sc} potential form has also been parameterized by Qi {\it et |
| 1320 |
< |
al.}\cite{Qi99} These parameters were obtained via empirical and {\it |
| 1321 |
< |
ab initio} calculations to match structural features of the FCC |
| 1322 |
< |
crystal. To specify the original Sutton-Chen variant of the {\sc sc} |
| 1323 |
< |
force field, the user would add the {\tt forceFieldVariant = "SC";} |
| 1324 |
< |
line to the meta-data file, while specification of the Qi {\it et al.} |
| 1325 |
< |
quantum-adapted variant of the {\sc sc} potential, the user would add |
| 1326 |
< |
the {\tt forceFieldVariant = "QSC";} line to the meta-data file. |
| 2021 |
> |
\section{\label{section::ffLongRange}Long Range Interactions} |
| 2022 |
|
|
| 2023 |
< |
\section{\label{section:clay}The CLAY force field} |
| 2023 |
> |
Calculating the long-range (non-bonded) potential involves a sum over |
| 2024 |
> |
all pairs of atoms (except for those atoms which are involved in a |
| 2025 |
> |
bond, bend, or torsion with each other). Many of these interactions |
| 2026 |
> |
can be inferred from the AtomTypes, |
| 2027 |
|
|
| 2028 |
< |
The {\sc clay} force field is based on an ionic (nonbonded) |
| 2029 |
< |
description of the metal-oxygen interactions associated with hydrated |
| 2030 |
< |
phases. All atoms are represented as point charges and are allowed |
| 2031 |
< |
complete translational freedom. Metal-oxygen interactions are based on |
| 2032 |
< |
a simple Lennard-Jones potential combined with electrostatics. The |
| 2033 |
< |
empirical parameters were optimized by Cygan {\it et |
| 2034 |
< |
al.}\cite{Cygan04} on the basis of known mineral structures, and |
| 2035 |
< |
partial atomic charges were derived from periodic DFT quantum chemical |
| 2036 |
< |
calculations of simple oxide, hydroxide, and oxyhydroxide model |
| 2037 |
< |
compounds with well-defined structures. |
| 2028 |
> |
\subsection{\label{section:ffNBinteraction}The NonBondedInteractions |
| 2029 |
> |
block} |
| 2030 |
> |
|
| 2031 |
> |
The user might want like to specify explicit or special interactions |
| 2032 |
> |
that override the default non-bodned interactions that are inferred |
| 2033 |
> |
from the AtomTypes. To do this, OpenMD implements a |
| 2034 |
> |
NonBondedInteractions block to allow the user to specify the following |
| 2035 |
> |
(pair-wise) non-bonded interactions: |
| 2036 |
> |
|
| 2037 |
> |
\begin{itemize} |
| 2038 |
> |
\item {\tt LennardJones}: |
| 2039 |
> |
\begin{equation*} |
| 2040 |
> |
V_{\text{NB}}(r) = 4 \epsilon_{ij} \left( |
| 2041 |
> |
\left(\frac{\sigma_{ij}}{r} \right)^{12} - |
| 2042 |
> |
\left(\frac{\sigma_{ij}}{r} \right)^{6} \right), |
| 2043 |
> |
\end{equation*} |
| 2044 |
> |
\item {\tt ShiftedMorse}: |
| 2045 |
> |
\begin{equation*} |
| 2046 |
> |
V_{\text{NB}}(r) = D_{ij} \left( e^{-2 \beta_{ij} (r - |
| 2047 |
> |
r^0)} - 2 e^{- \beta_{ij} (r - |
| 2048 |
> |
r^0)} \right), |
| 2049 |
> |
\end{equation*} |
| 2050 |
> |
\item {\tt RepulsiveMorse}: |
| 2051 |
> |
\begin{equation*} |
| 2052 |
> |
V_{\text{NB}}(r) = D_{ij} \left( e^{-2 \beta_{ij} (r - |
| 2053 |
> |
r^0)} \right), |
| 2054 |
> |
\end{equation*} |
| 2055 |
> |
\item {\tt RepulsivePower}: |
| 2056 |
> |
\begin{equation*} |
| 2057 |
> |
V_{\text{NB}}(r) = \epsilon_{ij} |
| 2058 |
> |
\left(\frac{\sigma_{ij}}{r} \right)^{n_{ij}}. |
| 2059 |
> |
\end{equation*} |
| 2060 |
> |
\end{itemize} |
| 2061 |
> |
|
| 2062 |
> |
\begin{code}[caption={[An example of a NonBondedInteractions block.] A |
| 2063 |
> |
simple example of a NonBondedInteractions block. Distances ($\sigma, |
| 2064 |
> |
r_0$) are given in \AA, while energies ($\epsilon, D0$) are in |
| 2065 |
> |
kcal/mol. The Morse potentials have an additional parameter $\beta_0$ |
| 2066 |
> |
which is in units of \AA$^{-1}$.}, |
| 2067 |
> |
label={sch:NonBondedInteractionTypes}] |
| 2068 |
> |
begin NonBondedInteractions |
| 2069 |
> |
|
| 2070 |
> |
//Lennard-Jones |
| 2071 |
> |
//Atom1 Atom2 LennardJones sigma epsilon |
| 2072 |
> |
Au CH3 LennardJones 3.54 0.2146 |
| 2073 |
> |
Au CH2 LennardJones 3.54 0.1749 |
| 2074 |
> |
Au CH LennardJones 3.54 0.1749 |
| 2075 |
> |
Au S LennardJones 2.40 8.465 |
| 2076 |
|
|
| 2077 |
+ |
//Shifted Morse |
| 2078 |
+ |
//Atom1 Atom2 ShiftedMorse r0 D0 beta0 |
| 2079 |
+ |
Au O_SPCE ShiftedMorse 3.70 0.0424 0.769 |
| 2080 |
|
|
| 2081 |
+ |
//Repulsive Morse |
| 2082 |
+ |
//Atom1 Atom2 RepulsiveMorse r0 D0 beta0 |
| 2083 |
+ |
Au H_SPCE RepulsiveMorse -1.00 0.00850 0.769 |
| 2084 |
+ |
|
| 2085 |
+ |
//Repulsive Power |
| 2086 |
+ |
//Atom1 Atom2 RepulsivePower sigma epsilon n |
| 2087 |
+ |
Au ON RepulsivePower 3.47005 0.186208 11 |
| 2088 |
+ |
Au NO RepulsivePower 3.53955 0.168629 11 |
| 2089 |
+ |
end NonBondedInteractions |
| 2090 |
+ |
\end{code} |
| 2091 |
+ |
|
| 2092 |
|
\section{\label{section:electrostatics}Electrostatics} |
| 2093 |
|
|
| 2094 |
< |
To aid in performing simulations in more traditional force fields, we |
| 2095 |
< |
have added routines to carry out electrostatic interactions using a |
| 2096 |
< |
number of different electrostatic summation methods. These methods |
| 2097 |
< |
are extended from the damped and cutoff-neutralized Coulombic sum |
| 2098 |
< |
originally proposed by Wolf, {\it et al.}\cite{Wolf99} One of these, |
| 2099 |
< |
the damped shifted force method, shows a remarkable ability to |
| 2100 |
< |
reproduce the energetic and dynamic characteristics exhibited by |
| 2101 |
< |
simulations employing lattice summation techniques. The basic idea is |
| 2102 |
< |
to construct well-behaved real-space summation methods using two tricks: |
| 2094 |
> |
Because nearly all force fields involve electrostatic interactions in |
| 2095 |
> |
one form or another, OpenMD implements a number of different |
| 2096 |
> |
electrostatic summation methods. These methods are extended from the |
| 2097 |
> |
damped and cutoff-neutralized Coulombic sum originally proposed by |
| 2098 |
> |
Wolf, {\it et al.}\cite{Wolf99} One of these, the damped shifted force |
| 2099 |
> |
method, shows a remarkable ability to reproduce the energetic and |
| 2100 |
> |
dynamic characteristics exhibited by simulations employing lattice |
| 2101 |
> |
summation techniques. The basic idea is to construct well-behaved |
| 2102 |
> |
real-space summation methods using two tricks: |
| 2103 |
|
\begin{enumerate} |
| 2104 |
|
\item shifting through the use of image charges, and |
| 2105 |
|
\item damping the electrostatic interaction. |
| 2218 |
|
\end{equation} |
| 2219 |
|
the shifted potential (eq. (\ref{eq:SPPot})) becomes |
| 2220 |
|
\begin{equation} |
| 2221 |
< |
V_{\textrm{DSP}}(r) = q_iq_j\left(\frac{\textrm{erfc}\left(\alpha r\right)}{r}-\ |
| 1472 |
< |
frac{\textrm{erfc}\left(\alpha R_\textrm{c}\right)}{R_\textrm{c}}\right) \quad r |
| 2221 |
> |
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 |
| 2222 |
|
\leqslant R_\textrm{c}, |
| 2223 |
|
\label{eq:DSPPot} |
| 2224 |
|
\end{equation} |
| 2262 |
|
this reason, the default electrostatic summation method utilized by |
| 2263 |
|
{\sc OpenMD} is the DSF (Eq. \ref{eq:DSFPot}) with a damping parameter |
| 2264 |
|
($\alpha$) that is set algorithmically from the cutoff radius. |
| 2265 |
+ |
|
| 2266 |
+ |
|
| 2267 |
+ |
\section{\label{section:cutoffGroups}Switching Functions} |
| 2268 |
+ |
|
| 2269 |
+ |
Calculating the the long-range interactions for $N$ atoms involves |
| 2270 |
+ |
$N(N-1)/2$ evaluations of atomic distances if it is done in a brute |
| 2271 |
+ |
force manner. To reduce the number of distance evaluations between |
| 2272 |
+ |
pairs of atoms, {\sc OpenMD} allows the use of hard or switched |
| 2273 |
+ |
cutoffs with Verlet neighbor lists.\cite{Allen87} Neutral groups which |
| 2274 |
+ |
contain charges can exhibit pathological forces unless the cutoff is |
| 2275 |
+ |
applied to the neutral groups evenly instead of to the individual |
| 2276 |
+ |
atoms.\cite{leach01:mm} {\sc OpenMD} allows users to specify cutoff |
| 2277 |
+ |
groups which may contain an arbitrary number of atoms in the molecule. |
| 2278 |
+ |
Atoms in a cutoff group are treated as a single unit for the |
| 2279 |
+ |
evaluation of the switching function: |
| 2280 |
+ |
\begin{equation} |
| 2281 |
+ |
V_{\mathrm{long-range}} = \sum_{a} \sum_{b>a} s(r_{ab}) \sum_{i \in a} \sum_{j \in b} V_{ij}(r_{ij}), |
| 2282 |
+ |
\end{equation} |
| 2283 |
+ |
where $r_{ab}$ is the distance between the centers of mass of the two |
| 2284 |
+ |
cutoff groups ($a$ and $b$). |
| 2285 |
+ |
|
| 2286 |
+ |
The sums over $a$ and $b$ are over the cutoff groups that are present |
| 2287 |
+ |
in the simulation. Atoms which are not explicitly defined as members |
| 2288 |
+ |
of a {\tt cutoffGroup} are treated as a group consisting of only one |
| 2289 |
+ |
atom. The switching function, $s(r)$ is the standard cubic switching |
| 2290 |
+ |
function, |
| 2291 |
+ |
\begin{equation} |
| 2292 |
+ |
S(r) = |
| 2293 |
+ |
\begin{cases} |
| 2294 |
+ |
1 & \text{if $r \le r_{\text{sw}}$},\\ |
| 2295 |
+ |
\frac{(r_{\text{cut}} + 2r - 3r_{\text{sw}})(r_{\text{cut}} - r)^2} |
| 2296 |
+ |
{(r_{\text{cut}} - r_{\text{sw}})^3} |
| 2297 |
+ |
& \text{if $r_{\text{sw}} < r \le r_{\text{cut}}$}, \\ |
| 2298 |
+ |
0 & \text{if $r > r_{\text{cut}}$.} |
| 2299 |
+ |
\end{cases} |
| 2300 |
+ |
\label{eq:dipoleSwitching} |
| 2301 |
+ |
\end{equation} |
| 2302 |
+ |
Here, $r_{\text{sw}}$ is the {\tt switchingRadius}, or the distance |
| 2303 |
+ |
beyond which interactions are reduced, and $r_{\text{cut}}$ is the |
| 2304 |
+ |
{\tt cutoffRadius}, or the distance at which interactions are |
| 2305 |
+ |
truncated. |
| 2306 |
+ |
|
| 2307 |
+ |
Users of {\sc OpenMD} do not need to specify the {\tt cutoffRadius} or |
| 2308 |
+ |
{\tt switchingRadius}. |
| 2309 |
+ |
If the {\tt cutoffRadius} was not explicitly set, OpenMD will attempt |
| 2310 |
+ |
to guess an appropriate choice. If the system contains electrostatic |
| 2311 |
+ |
atoms, the default cutoff radius is 12 \AA. In systems without |
| 2312 |
+ |
electrostatic (charge or multipolar) atoms, the atom types present in the simulation will be |
| 2313 |
+ |
polled for suggested cutoff values (e.g. $2.5 max(\left\{ \sigma |
| 2314 |
+ |
\right\})$ for Lennard-Jones atoms. The largest suggested value |
| 2315 |
+ |
that was found will be used. |
| 2316 |
+ |
|
| 2317 |
+ |
By default, OpenMD uses shifted force potentials to force the |
| 2318 |
+ |
potential energy and forces to smoothly approach zero at the cutoff |
| 2319 |
+ |
radius. If the user would like to use another cutoff method |
| 2320 |
+ |
they may do so by setting the {\tt cutoffMethod} parameter to: |
| 2321 |
+ |
\begin{itemize} |
| 2322 |
+ |
\item {\tt HARD} |
| 2323 |
+ |
\item {\tt SWITCHED} |
| 2324 |
+ |
\item {\tt SHIFTED\_FORCE} (default): |
| 2325 |
+ |
\item {\tt TAYLOR\_SHIFTED} |
| 2326 |
+ |
\item {\tt SHIFTED\_POTENTIAL} |
| 2327 |
+ |
\end{itemize} |
| 2328 |
+ |
|
| 2329 |
+ |
The {\tt switchingRadius} is set to a default value of 95\% of the |
| 2330 |
+ |
{\tt cutoffRadius}. In the special case of a simulation containing |
| 2331 |
+ |
{\it only} Lennard-Jones atoms, the default switching radius takes the |
| 2332 |
+ |
same value as the cutoff radius, and {\sc OpenMD} will use a shifted |
| 2333 |
+ |
potential to remove discontinuities in the potential at the cutoff. |
| 2334 |
+ |
Both radii may be specified in the meta-data file. |
| 2335 |
|
|
| 2336 |
|
\section{\label{section:pbc}Periodic Boundary Conditions} |
| 2337 |
|
|
| 2712 |
|
& (with changes to box shape) & \\ |
| 2713 |
|
NPTxyz & approximate isobaric-isothermal & {\tt ensemble = NPTxyz;} \\ |
| 2714 |
|
& (with separate barostats on each box dimension) & \\ |
| 2715 |
+ |
N$\gamma$T & constant lateral surface tension & {\tt ensemble = NgammaT;} \\ |
| 2716 |
+ |
& (must specify a surfaceTension) & \\ |
| 2717 |
+ |
NP$\gamma$T & constant normal pressure and lateral surface tension & {\tt ensemble = NPrT;} \\ |
| 2718 |
+ |
& (must specify a targetPressure and surfaceTension) & \\ |
| 2719 |
|
LD & Langevin Dynamics & {\tt ensemble = LD;} \\ |
| 2720 |
|
& (approximates the effects of an implicit solvent) & \\ |
| 2721 |
|
LangevinHull & Non-periodic Langevin Dynamics & {\tt ensemble = LangevinHull;} \\ |
| 2756 |
|
default value} \\ |
| 2757 |
|
$T_{\mathrm{target}}$ & {\tt targetTemperature = 300;} & K & none \\ |
| 2758 |
|
$P_{\mathrm{target}}$ & {\tt targetPressure = 1;} & atm & none \\ |
| 2759 |
+ |
$\gamma$ & {\tt surfaceTension = 0.015;} & Newtons / meter & none \\ |
| 2760 |
+ |
$\eta$ & {\tt viscosity = 0.0089;} & Poise & none \\ |
| 2761 |
|
$\tau_T$ & {\tt tauThermostat = 1e3;} & fs & none \\ |
| 2762 |
|
$\tau_B$ & {\tt tauBarostat = 5e3;} & fs & none \\ |
| 2763 |
|
& {\tt resetTime = 200;} & fs & none \\ |
| 3446 |
|
|
| 3447 |
|
\section{Constant Pressure without periodic boundary conditions (The LangevinHull)} |
| 3448 |
|
|
| 3449 |
< |
The Langevin Hull uses an external bath at a fixed constant pressure |
| 3449 |
> |
The Langevin Hull\cite{Vardeman2011} uses an external bath at a fixed constant pressure |
| 3450 |
|
($P$) and temperature ($T$) with an effective solvent viscosity |
| 3451 |
|
($\eta$). This bath interacts only with the objects on the exterior |
| 3452 |
|
hull of the system. Defining the hull of the atoms in a simulation is |
| 3758 |
|
\label{table:zconParams} |
| 3759 |
|
\end{longtable} |
| 3760 |
|
|
| 3761 |
< |
\chapter{\label{section:restraints}Restraints} |
| 3762 |
< |
Restraints are external potentials that are added to a system to keep |
| 3763 |
< |
particular molecules or collections of particles close to some |
| 3764 |
< |
reference structure. A restraint can be a collective |
| 3761 |
> |
\chapter{\label{chapter:restraints}Restraints} |
| 3762 |
> |
Restraints are external potential energy functions that are added to a |
| 3763 |
> |
system to keep particular molecules or collections of particles close |
| 3764 |
> |
to a reference structure. A \textbf{Molecular} restraint is a |
| 3765 |
> |
collective force applied to all atoms in a molecule, while an |
| 3766 |
> |
\textbf{Object} restraint is a simple harmonic spring connecting the |
| 3767 |
> |
position of a StuntDouble (or its orientation) close to a fixed |
| 3768 |
> |
reference geometry. |
| 3769 |
|
|
| 3770 |
+ |
Restraints require the specification of a reference geometry in the |
| 3771 |
+ |
{\tt Restraint\_file} parameter. These files are standard OpenMD {\tt |
| 3772 |
+ |
md} or {\tt eor} files which must have the same component |
| 3773 |
+ |
specification and StuntDouble indices as the simulation itself. |
| 3774 |
+ |
|
| 3775 |
+ |
Restraint potentials in OpenMD are harmonic, |
| 3776 |
+ |
\begin{equation} |
| 3777 |
+ |
V_\mathrm{trans} = \frac{k_\mathrm{trans}}{2} \left( \mathbf{i} - \mathbf{r} |
| 3778 |
+ |
\right)^2 |
| 3779 |
+ |
\end{equation} |
| 3780 |
+ |
where $k_\mathrm{trans}$ is a spring constant for translational |
| 3781 |
+ |
motion, $\mathbf{i}$ is the instantaneous position of the object, and |
| 3782 |
+ |
$\mathbf{r}$ is the position of the same object in the reference |
| 3783 |
+ |
structure. Alternatively, one might restrain the orientations of an object, |
| 3784 |
+ |
\begin{equation} |
| 3785 |
+ |
V_\mathrm{swing} = \frac{k_\mathrm{swing}}{2} \left( \theta - \theta_0 \right)^2 |
| 3786 |
+ |
\end{equation} |
| 3787 |
+ |
where $k_\mathrm{swing}$ is the force constant for swing motion of the |
| 3788 |
+ |
long axis of a molecule, $\theta$ is the instantaneous swing angle |
| 3789 |
+ |
relative to the reference structure, and $\theta_0$ is an optional |
| 3790 |
+ |
angle that the user can specify that is also measured relative to the |
| 3791 |
+ |
orientation of the reference structure. |
| 3792 |
+ |
|
| 3793 |
+ |
\begin{code}[caption={[Specifying restraints for all objects matching a |
| 3794 |
+ |
selection] Sample keywords defining object restraints (here the |
| 3795 |
+ |
object is the first rigid body associated with SPCE molecules},label={sch:restObj}] |
| 3796 |
+ |
restraint{ |
| 3797 |
+ |
restraintType = "object"; |
| 3798 |
+ |
objectSelection = "select SPCE_RB_0"; |
| 3799 |
+ |
displacementSpringConstant = 4.3; |
| 3800 |
+ |
twistSpringConstant = 750; |
| 3801 |
+ |
swingXSpringConstant = 700; |
| 3802 |
+ |
swingYSpringConstant = 700; |
| 3803 |
+ |
print = "false"; |
| 3804 |
+ |
} |
| 3805 |
+ |
|
| 3806 |
+ |
useRestraints = "true"; |
| 3807 |
+ |
Restraint_file = "idealStructure.in"; |
| 3808 |
+ |
\end{code} |
| 3809 |
+ |
|
| 3810 |
+ |
When the restraint is added to an entire molecule, the forces and |
| 3811 |
+ |
torques must be applied atom-by-atom. Translational restraints are |
| 3812 |
+ |
simple to apply, but torques for angular restraints require some |
| 3813 |
+ |
care. |
| 3814 |
+ |
|
| 3815 |
+ |
Defining the rotation angles for an instantaneous geometry of a |
| 3816 |
+ |
molecule relative to a reference structure is a difficult problem |
| 3817 |
+ |
because there are multiple combinations of three-angle rotations can |
| 3818 |
+ |
lead to the same structure. To tackle this problem, OpenMD combines |
| 3819 |
+ |
two methods, singular value decomposition (SVD) and twist-swing |
| 3820 |
+ |
decomposition. |
| 3821 |
+ |
|
| 3822 |
+ |
The core of the difficulty is in identifying the rotation matrix |
| 3823 |
+ |
($\mathsf{A}$) that relates the instantaneous geometry to the |
| 3824 |
+ |
reference structure. Because molecules generally are not rigid, the |
| 3825 |
+ |
internal dynamics makes this computationally demanding. Fortunately a |
| 3826 |
+ |
method that is widely used for protein alignment provides a reasonable |
| 3827 |
+ |
characterization of this rotation matrix. |
| 3828 |
+ |
|
| 3829 |
+ |
The molecule is represented as a $N \times 3$ matrix of coordinates. |
| 3830 |
+ |
The difference between the instantaneous and reference conformations |
| 3831 |
+ |
is encoded in the transformation between two of these matrices. The |
| 3832 |
+ |
transformation consists of a ``best fit'' translation vector and |
| 3833 |
+ |
rotation matrix such that after translation and rotation, the two |
| 3834 |
+ |
configurations will have the highest degree of overlap with each |
| 3835 |
+ |
other. I.e., the translation vector and rotation matrix minimize the |
| 3836 |
+ |
root mean-square distance (RMSD) between the two structures, |
| 3837 |
+ |
\begin{equation} |
| 3838 |
+ |
\mathrm{RMSD} = \sqrt{\frac{1}{N} \sum_n \left(i_n - r_n\right)^2 }, |
| 3839 |
+ |
\end{equation} |
| 3840 |
+ |
where $\left\{i_n\right\}$ and $\left\{r_n\right\}$ are the sets of |
| 3841 |
+ |
coordinates for the instantaneous and reference structures, and $N$ is |
| 3842 |
+ |
the total number of atoms in the molecule. A singular value |
| 3843 |
+ |
decomposition (SVD) is carried out in order to minimize the RMSD. |
| 3844 |
+ |
This operation provides the best-fitting translation vector $\vec{v}$ |
| 3845 |
+ |
and rotation matrix $\mathsf{A}$ that align the instantaneous and |
| 3846 |
+ |
reference structures. |
| 3847 |
+ |
|
| 3848 |
+ |
Twist-swing decomposition, a technique that has been widely used in |
| 3849 |
+ |
computer animation of articulated figures, is then employed to |
| 3850 |
+ |
calculate the relative rotational angles between the instantaneous and |
| 3851 |
+ |
reference structures. This decomposition regards the motion as a |
| 3852 |
+ |
``twist'' about one axis followed by a ``swing'' about another axis, |
| 3853 |
+ |
where the second axis is perpendicular to the |
| 3854 |
+ |
first.\cite{Kallmann:2008fk,Shoemake:1994uq} This model of rotation |
| 3855 |
+ |
provides a convenient way to define a unique relative rotational |
| 3856 |
+ |
angle. A simple example helps clarify: Suppose a cylinder moves from |
| 3857 |
+ |
original position $\mathbf{O}$ to end position $\mathbf{E}$ (Figure |
| 3858 |
+ |
\ref{fig:twistSwing}). Although there are many paths that can |
| 3859 |
+ |
accomplish this movement, the simplest path is to rotate the reference |
| 3860 |
+ |
configuration by $\theta$ (i.e. the swing angle) in the plane |
| 3861 |
+ |
$\mathbf{O} \times \mathbf{E}$ (the cross product of the two |
| 3862 |
+ |
orientation vectors). Then the reference structure is rotated by |
| 3863 |
+ |
$\phi$ (the twist angle) around the central axis. |
| 3864 |
+ |
|
| 3865 |
+ |
\begin{figure} |
| 3866 |
+ |
\includegraphics[width=3.5in]{twistSwing} |
| 3867 |
+ |
\caption[The twist-swing decomposition used in orientational |
| 3868 |
+ |
restraints] {The twist-swing decomposition defines relative rotational |
| 3869 |
+ |
angles using the simplest path to rotate the reference |
| 3870 |
+ |
configuration. The reference is rotated by a ``swing'' angle |
| 3871 |
+ |
$\theta$ in the plane of $\mathbf{O} \times \mathbf{E}$, then |
| 3872 |
+ |
rotated by a ``twist'' angle $\phi$ around the swing axis.} |
| 3873 |
+ |
\label{fig:twistSwing} |
| 3874 |
+ |
\end{figure} |
| 3875 |
+ |
|
| 3876 |
+ |
This illustrates the basic idea of the twist-swing decomposition, |
| 3877 |
+ |
which is a general operation that can be performed on the best-fitting |
| 3878 |
+ |
relative rotation matrix ($\mathsf{A}$) between the instantaneous and |
| 3879 |
+ |
reference structures. For objects that are not cylindrically |
| 3880 |
+ |
symmetric, there are two swing angles (one in X and one in Y) in |
| 3881 |
+ |
addition to the twist angle. |
| 3882 |
+ |
|
| 3883 |
+ |
\begin{code}[caption={[Specifying Restraints for a Molecule] Sample |
| 3884 |
+ |
keywords defining a molecular restraint and the associated force constants},label={sch:restMol}] |
| 3885 |
+ |
restraint{ |
| 3886 |
+ |
restraintType = "Molecular"; |
| 3887 |
+ |
molIndex = 0; |
| 3888 |
+ |
twistSpringConstant = 0.25; |
| 3889 |
+ |
swingXSpringConstant = 0.02; |
| 3890 |
+ |
restrainedSwingXAngle = -90.0; |
| 3891 |
+ |
swingYSpringConstant = 0.02; |
| 3892 |
+ |
print = "true"; |
| 3893 |
+ |
} |
| 3894 |
+ |
|
| 3895 |
+ |
useRestraints = "true"; |
| 3896 |
+ |
Restraint_file = "prism_ref.inc"; |
| 3897 |
+ |
\end{code} |
| 3898 |
+ |
|
| 3899 |
+ |
To specify a molecular restraint, it is necessary to give an exact |
| 3900 |
+ |
index of this molecule in the {\tt molIndex} parameter. In the |
| 3901 |
+ |
example in schem \ref{sch:restMol}, the restrained SwingX angle is -90 |
| 3902 |
+ |
degrees offset from the reference structure, so the molecule will be |
| 3903 |
+ |
restrained in a different orientation from the reference geometry. |
| 3904 |
+ |
|
| 3905 |
+ |
\begin{longtable}[c]{ABCD} |
| 3906 |
+ |
\caption{Meta-data Keywords: Restraint Parameters} |
| 3907 |
+ |
\\ |
| 3908 |
+ |
{\bf keyword} & {\bf units} & {\bf use} & {\bf remarks} \\ \hline |
| 3909 |
+ |
\endhead |
| 3910 |
+ |
\hline |
| 3911 |
+ |
\endfoot |
| 3912 |
+ |
{\tt restraintType} & string & What kind of restraint is this? & |
| 3913 |
+ |
choose either ``Object'' or ``Molecular'' \\ |
| 3914 |
+ |
{\tt molIndex} & integer & Which molecule to restrain & \\ |
| 3915 |
+ |
{\tt objectSelection} & string & Selection script for Object |
| 3916 |
+ |
restraints & \\ |
| 3917 |
+ |
{\tt displacementSpringConstant} & kcal/mol/\AA$^2$ & $k_\mathrm{trans}$ & \\ |
| 3918 |
+ |
{\tt twistSpringConstant} & kcal/mol/radian$^2$ & $k_\mathrm{twist}$ & \\ |
| 3919 |
+ |
{\tt swingXSpringConstant} & kcal/mol/radian$^2$ & $k_\mathrm{swingX}$ & \\ |
| 3920 |
+ |
{\tt swingYSpringConstant} & kcal/mol/radian$^2$ & $k_\mathrm{swingY}$ & \\ |
| 3921 |
+ |
{\tt restrainedTwistAngle} & degrees & $\omega_0$ (optional) & |
| 3922 |
+ |
defaults to 0\\ |
| 3923 |
+ |
{\tt restrainedSwingXAngle} & degrees & $\theta_0^x$ (optional) & defaults to 0 \\ |
| 3924 |
+ |
{\tt restrainedSwingYAngle} & degrees & $\theta_0^y$ (optional) & defaults to 0\\ |
| 3925 |
+ |
{\tt print} & logical & whether or not to print restraint value and energy & |
| 3926 |
+ |
defaults to true |
| 3927 |
+ |
\label{table:restraintParams} |
| 3928 |
+ |
\end{longtable} |
| 3929 |
+ |
|
| 3930 |
+ |
The various parameters for a {\tt restraint} are shown in table |
| 3931 |
+ |
\ref{table:restraintParams}. Because restraints have a large number of |
| 3932 |
+ |
parameters, these must be enclosed in a {\it separate} {\tt |
| 3933 |
+ |
restraint\{...\}} block. |
| 3934 |
+ |
|
| 3935 |
+ |
\chapter{\label{chapter:perturbations}Perturbations} |
| 3936 |
+ |
OpenMD allows the user to specify two external perturbations that |
| 3937 |
+ |
interact with the electrostatic properties of the atoms. |
| 3938 |
+ |
|
| 3939 |
+ |
\section{\label{sec:UniformField}Uniform Fields} |
| 3940 |
+ |
To apply a uniform |
| 3941 |
+ |
(vector) electric field to the system, the user adds the {\tt |
| 3942 |
+ |
uniformField} parameter to the MetaData section of the .md file. |
| 3943 |
+ |
|
| 3944 |
+ |
\begin{code}[caption={Specifying a uniform electric field.},label={sch:uniformField}] |
| 3945 |
+ |
uniformField = (a, b, c); |
| 3946 |
+ |
\end{code} |
| 3947 |
+ |
|
| 3948 |
+ |
The values of $a$, $b$, and $c$ are in units of V~/~\AA. The |
| 3949 |
+ |
electrostatic potential corresponding to this uniform field is |
| 3950 |
+ |
\begin{equation} |
| 3951 |
+ |
\phi(\mathbf{r}) = - a x - b y - c z |
| 3952 |
+ |
\end{equation} |
| 3953 |
+ |
which grows unbounded and is not periodic. For these reasons, care |
| 3954 |
+ |
should be taken in using a uniform field with point charges. |
| 3955 |
+ |
The field itself is |
| 3956 |
+ |
\begin{equation} |
| 3957 |
+ |
\mathbf{E} = \left( \begin{array}{c} a \\ b \\ c \end{array} \right). |
| 3958 |
+ |
\end{equation} |
| 3959 |
+ |
The uniform field applies a force on charged atoms, $ \mathbf{F} = C |
| 3960 |
+ |
\mathbf{E}$. For dipolar atoms, the uniform field applies both a |
| 3961 |
+ |
potential, $ U = - \mathbf{D} \cdot \mathbf{E}$ and a torque, $ |
| 3962 |
+ |
\mathbf{\tau} = \mathbf{D} \times \mathbf{E}$ that depends on the |
| 3963 |
+ |
instantaneous dipole ($\mathbf{D}$) of that atom. |
| 3964 |
+ |
|
| 3965 |
+ |
\section{\label{sec:UniformGradient}Uniform Field Gradients} |
| 3966 |
+ |
To apply a uniform electric field gradient to the system, the user |
| 3967 |
+ |
adds three parameters to the MetaData section |
| 3968 |
+ |
|
| 3969 |
+ |
\begin{code}[caption={Specifying a uniform electric field gradient.},label={sch:uniformFieldGradient}] |
| 3970 |
+ |
uniformGradientStrength = g; |
| 3971 |
+ |
uniformGradientDirection1 = (a1, a2, a3) |
| 3972 |
+ |
uniformGradientDirection2 = (b1, b2, b3); |
| 3973 |
+ |
\end{code} |
| 3974 |
+ |
|
| 3975 |
+ |
The two direction vectors, $\mathbf{a}$ and $\mathbf{b}$ are unit |
| 3976 |
+ |
vectors, and the value of $g$ is in units of |
| 3977 |
+ |
V~/~\AA\textsuperscript{2}. The electrostatic potential corresponding |
| 3978 |
+ |
to this uniform gradient is |
| 3979 |
+ |
\begin{align} |
| 3980 |
+ |
\phi(\mathbf{r}) = - \frac{g}{2} \Big[& |
| 3981 |
+ |
\left(a_1 b_1 - \frac{\cos\psi}{3}\right) x^2 |
| 3982 |
+ |
+ (a_1 b_2 + a_2 b_1) x y + (a_1 b_3 + a_3 b_1) x z \\ |
| 3983 |
+ |
& + (a_2 b_1 + a_1 b_2) y x |
| 3984 |
+ |
+ \left(a_2 b_2 - \frac{\cos\psi}{3}\right) y^2 |
| 3985 |
+ |
+ (a_2 b_3 + a_3 b_2) y z \\ |
| 3986 |
+ |
& \left. + (a_3 b_1 + a_1 b_3) z x |
| 3987 |
+ |
+ (a_3 b_2 + a_2 b_3) z y |
| 3988 |
+ |
+ \left(a_3 b_3 - \frac{\cos\psi}{3}\right) z^2 \right]. \\ |
| 3989 |
+ |
\end{align} |
| 3990 |
+ |
where $\cos \psi = \mathbf{a} \cdot \mathbf{b}$. Note that this |
| 3991 |
+ |
potential grows unbounded and is not periodic. For these reasons, |
| 3992 |
+ |
care should be taken in using a Uniform Gradient with point |
| 3993 |
+ |
charges. The corresponding field for this potential is: |
| 3994 |
+ |
\begin{equation} |
| 3995 |
+ |
\mathbf{E} = \frac{g}{2} \left( \begin{array}{c} |
| 3996 |
+ |
2\left(a_1 b_1 - \frac{\cos\psi}{3}\right) x + (a_1 b_2 + a_2 b_1) y |
| 3997 |
+ |
+ (a_1 b_3 + a_3 b_1) z \\ |
| 3998 |
+ |
(a_2 b_1 + a_1 b_2) x + 2 \left(a_2 b_2 - \frac{\cos\psi}{3}\right) y |
| 3999 |
+ |
+ (a_2 b_3 + a_3 b_2) z \\ |
| 4000 |
+ |
(a_3 b_1 + a_1 b_3) x + (a_3 b_2 + a_2 b_3) y |
| 4001 |
+ |
+ 2 \left(a_3 b_3 - \frac{\cos\psi}{3}\right) z \end{array} |
| 4002 |
+ |
\right). |
| 4003 |
+ |
\end{equation} |
| 4004 |
+ |
The field also grows unbounded and is not periodic. For these reasons, |
| 4005 |
+ |
care should be taken in using a Uniform Gradient with point dipoles. |
| 4006 |
+ |
|
| 4007 |
+ |
The corresponding field gradient, |
| 4008 |
+ |
\begin{equation} |
| 4009 |
+ |
\nabla \mathbf{E} = \frac{g}{2} \left( \begin{array}{ccc} |
| 4010 |
+ |
2\left(a_1 b_1 - \frac{\cos\psi}{3}\right) & |
| 4011 |
+ |
(a_1 b_2 + a_2 b_1) & (a_1 b_3 + a_3 b_1) \\ |
| 4012 |
+ |
(a_2 b_1 + a_1 b_2) & 2 \left(a_2 b_2 - \frac{\cos\psi}{3}\right) & |
| 4013 |
+ |
(a_2 b_3 + a_3 b_2) \\ |
| 4014 |
+ |
(a_3 b_1 + a_1 b_3) & (a_3 b_2 + a_2 b_3) & |
| 4015 |
+ |
2 \left(a_3 b_3 - \frac{\cos\psi}{3}\right) \end{array} \right) |
| 4016 |
+ |
\end{equation} |
| 4017 |
+ |
is uniform everywhere. The uniform field gradient applies a force on |
| 4018 |
+ |
charged atoms, $\mathbf{F} = C~\mathbf{E}(\mathbf{r})$. For dipolar |
| 4019 |
+ |
atoms, the gradient applies a potential, $U = -\mathbf{D} \cdot |
| 4020 |
+ |
\mathbf{E}(\mathbf{r})$, force, $\mathbf{F} = \mathbf{D} \cdot |
| 4021 |
+ |
\nabla \mathbf{E}$, and torque, $ \mathbf{\tau} = \mathbf{D} \times |
| 4022 |
+ |
\mathbf{E}(\mathbf{r})$. For quadrupolar atoms, the uniform field |
| 4023 |
+ |
gradient exerts a potential, $U = - \mathsf{Q}:\nabla \mathbf{E}$, and |
| 4024 |
+ |
a torque $ \mathbf{F} = 2 \mathsf{Q} \times \nabla \mathbf{E}$. |
| 4025 |
+ |
|
| 4026 |
+ |
Here, the $:$ indicates a tensor contraction (double dot product) of |
| 4027 |
+ |
two matrices, and the $\times$ for the quadrupole indicates a vector |
| 4028 |
+ |
(cross) product of two matrices, defined as |
| 4029 |
+ |
\begin{equation} |
| 4030 |
+ |
\left[\mathsf{A} \times \mathsf{B}\right]_\alpha = \sum_\beta |
| 4031 |
+ |
\left[\mathsf{A}_{\alpha+1,\beta} \mathsf{B}_{\alpha+2,\beta} |
| 4032 |
+ |
-\mathsf{A}_{\alpha+2,\beta} \mathsf{B}_{\alpha+1,\beta} |
| 4033 |
+ |
\right] |
| 4034 |
+ |
\label{eq:matrixCross} |
| 4035 |
+ |
\end{equation} |
| 4036 |
+ |
where $\alpha+1$ and $\alpha+2$ are regarded as cyclic |
| 4037 |
+ |
permuations of the matrix indices. |
| 4038 |
+ |
|
| 4039 |
|
\chapter{\label{section:thermInt}Thermodynamic Integration} |
| 4040 |
|
|
| 4041 |
|
Thermodynamic integration is an established technique that has been |
| 4070 |
|
\end{equation} |
| 4071 |
|
where $K_\mathrm{v}$, $K_\mathrm{\theta}$, and $K_\mathrm{\omega}$ are |
| 4072 |
|
the spring constants restraining translational motion and deflection |
| 4073 |
< |
of and rotation around the principle axis of the molecule |
| 4074 |
< |
respectively. The values of $\theta$ range from $0$ to $\pi$, while |
| 4075 |
< |
$\omega$ ranges from $-\pi$ to $\pi$. |
| 4073 |
> |
of (swing) and rotation around (twist) the principle axis of the |
| 4074 |
> |
molecule respectively. The values of $\theta$ range from $0$ to |
| 4075 |
> |
$\pi$, while $\omega$ ranges from $-\pi$ to $\pi$. |
| 4076 |
|
|
| 4077 |
|
The partition function for a molecular crystal restrained in this |
| 4078 |
|
fashion can be evaluated analytically, and the Helmholtz Free Energy |
| 4108 |
|
molecular centers of mass, and two for displacements from the ideal |
| 4109 |
|
orientations of the molecules. |
| 4110 |
|
|
| 4111 |
+ |
\begin{code}[caption={[Specifying Restraints for Thermodynamic |
| 4112 |
+ |
Integration to an Einstein Crystal]Sample keywords defining restraints |
| 4113 |
+ |
and their force constants for use in Thermodynamic |
| 4114 |
+ |
Integration to an Einstein Crystal},label={sch:tiSolid}] |
| 4115 |
+ |
useThermodynamicIntegration = "true"; |
| 4116 |
+ |
thermodynamicIntegrationLambda = 0.0; |
| 4117 |
+ |
thermodynamicIntegrationK = 1.0; |
| 4118 |
+ |
|
| 4119 |
+ |
restraint{ |
| 4120 |
+ |
restraintType = "object"; |
| 4121 |
+ |
objectSelection = "select SSD_E"; |
| 4122 |
+ |
displacementSpringConstant = 4.3; |
| 4123 |
+ |
twistSpringConstant = 750; |
| 4124 |
+ |
swingXSpringConstant = 700; |
| 4125 |
+ |
swingYSpringConstant = 700; |
| 4126 |
+ |
print = "false"; |
| 4127 |
+ |
} |
| 4128 |
+ |
|
| 4129 |
+ |
useRestraints = "true"; |
| 4130 |
+ |
Restraint_file = "idealCrystal.in"; |
| 4131 |
+ |
\end{code} |
| 4132 |
+ |
|
| 4133 |
|
To construct a thermodynamic integration path, the user would run a |
| 4134 |
< |
sequence of $N$ simulations, each with a different value of lambda |
| 4135 |
< |
between $0$ and $1$. When {\tt useSolidThermInt} is set to {\tt true} |
| 4136 |
< |
in the meta-data file, two additional energy columns are reported in |
| 4137 |
< |
the {\tt .stat} file for the simulation. The first, {\tt vRaw}, is |
| 4138 |
< |
the unperturbed energy for the configuration, and the second, {\tt |
| 4139 |
< |
vHarm}, is the energy of the harmonic (Einstein) system in an |
| 4140 |
< |
identical configuration. The total potential energy of the |
| 4141 |
< |
configuration is a linear combination of {\tt vRaw} and {\tt vHarm} |
| 4142 |
< |
weighted by the value of $\lambda$. |
| 4134 |
> |
sequence of $N$ simulations, each with a different value of $\lambda$ |
| 4135 |
> |
between $0$ and $1$. When {\tt useThermodynamicIntegration} is set to |
| 4136 |
> |
{\tt true} in the meta-data file and restraints are present, two |
| 4137 |
> |
additional energy columns are reported in the {\tt .stat} file for the |
| 4138 |
> |
simulation. The first, {\tt vRaw}, is the unperturbed energy for the |
| 4139 |
> |
configuration, and the second, {\tt vHarm}, is the energy of the |
| 4140 |
> |
harmonic (Einstein) system in an identical configuration. The total |
| 4141 |
> |
potential energy of the configuration is a linear combination of {\tt |
| 4142 |
> |
vRaw} and {\tt vHarm} weighted by the value of $\lambda$. |
| 4143 |
|
|
| 4144 |
|
From a running average of the difference between {\tt vRaw} and {\tt |
| 4145 |
|
vHarm}, the user can obtain the integrand in Eq. (\ref{eq:thermInt}) |
| 4146 |
|
for fixed value of $\lambda$. |
| 4147 |
|
|
| 3028 |
– |
There are two additional files with the suffixes {\tt .zang0} and {\tt |
| 3029 |
– |
.zang} generated by {\sc OpenMD} during the first run of a solid |
| 3030 |
– |
thermodynamic integration. These files contain the initial ({\tt |
| 3031 |
– |
.zang0}) and final ({\tt .zang}) values of the angular displacement |
| 3032 |
– |
coordinates for each of the molecules. These are particularly useful |
| 3033 |
– |
when chaining a number of simulations (with successive values of |
| 3034 |
– |
$\lambda$) together. |
| 3035 |
– |
|
| 4148 |
|
For {\it liquid} thermodynamic integrations, the reference system is |
| 4149 |
|
the ideal gas (with a potential exactly equal to 0), so the {\tt |
| 4150 |
|
.stat} file contains only the standard columns. The potential energy |
| 4153 |
|
potential energy directly as the $\Delta V$ in the integrand of |
| 4154 |
|
Eq. (\ref{eq:thermInt}). |
| 4155 |
|
|
| 4156 |
+ |
\begin{code}[caption={[Specifying Restraints for Thermodynamic |
| 4157 |
+ |
Integration to an Ideal Gas]Sample keywords for use in Thermodynamic |
| 4158 |
+ |
Integration to an Ideal Gas},label={sch:tiLiquid}] |
| 4159 |
+ |
useThermodynamicIntegration = "true"; |
| 4160 |
+ |
thermodynamicIntegrationLambda = 1.0; |
| 4161 |
+ |
thermodynamicIntegrationK = 1.0; |
| 4162 |
+ |
\end{code} |
| 4163 |
+ |
|
| 4164 |
|
Meta-data parameters concerning thermodynamic integrations are given in |
| 4165 |
|
Table~\ref{table:thermIntParams} |
| 4166 |
|
|
| 4171 |
|
\endhead |
| 4172 |
|
\hline |
| 4173 |
|
\endfoot |
| 4174 |
< |
{\tt useSolidThermInt} & logical & perform thermodynamic integration |
| 3055 |
< |
to an Einstein crystal? & default is ``false'' \\ |
| 3056 |
< |
{\tt useLiquidThermInt} & logical & perform thermodynamic integration |
| 3057 |
< |
to an ideal gas? & default is ``false'' \\ |
| 4174 |
> |
{\tt useThermodynamicIntegration} & logical & perform thermodynamic integration? & default is ``false'' \\ |
| 4175 |
|
{\tt thermodynamicIntegrationLambda} & & & \\ |
| 4176 |
|
& double & transformation |
| 4177 |
|
parameter & Sets how far along the thermodynamic integration path the |
| 4179 |
|
{\tt thermodynamicIntegrationK} & & & \\ |
| 4180 |
|
& double & & power of $\lambda$ |
| 4181 |
|
governing shape of integration pathway \\ |
| 3065 |
– |
{\tt thermIntDistSpringConst} & & & \\ |
| 3066 |
– |
& $\mbox{kcal~mol}^{-1} \mbox{\AA}^{-2}$ |
| 3067 |
– |
& & spring constant for translations in Einstein crystal \\ |
| 3068 |
– |
{\tt thermIntThetaSpringConst} & & & \\ |
| 3069 |
– |
& $\mbox{kcal~mol}^{-1} |
| 3070 |
– |
\mbox{rad}^{-2}$ & & spring constant for deflection away from z-axis |
| 3071 |
– |
in Einstein crystal \\ |
| 3072 |
– |
{\tt thermIntOmegaSpringConst} & & & \\ |
| 3073 |
– |
& $\mbox{kcal~mol}^{-1} |
| 3074 |
– |
\mbox{rad}^{-2}$ & & spring constant for rotation around z-axis in |
| 3075 |
– |
Einstein crystal |
| 4182 |
|
\label{table:thermIntParams} |
| 4183 |
|
\end{longtable} |
| 4184 |
+ |
|
| 4185 |
+ |
\chapter{\label{section:rnemd}Reverse Non-Equilibrium Molecular Dynamics (RNEMD)} |
| 4186 |
+ |
|
| 4187 |
+ |
There are many ways to compute transport properties from molecular |
| 4188 |
+ |
dynamics simulations. Equilibrium Molecular Dynamics (EMD) |
| 4189 |
+ |
simulations can be used by computing relevant time correlation |
| 4190 |
+ |
functions and assuming linear response theory holds. For some transport properties (notably thermal conductivity), EMD approaches |
| 4191 |
+ |
are subject to noise and poor convergence of the relevant |
| 4192 |
+ |
correlation functions. Traditional Non-equilibrium Molecular Dynamics |
| 4193 |
+ |
(NEMD) methods impose a gradient (e.g. thermal or momentum) on a |
| 4194 |
+ |
simulation. However, the resulting flux is often difficult to |
| 4195 |
+ |
measure. Furthermore, problems arise for NEMD simulations of |
| 4196 |
+ |
heterogeneous systems, such as phase-phase boundaries or interfaces, |
| 4197 |
+ |
where the type of gradient to enforce at the boundary between |
| 4198 |
+ |
materials is unclear. |
| 4199 |
+ |
|
| 4200 |
+ |
{\it Reverse} Non-Equilibrium Molecular Dynamics (RNEMD) methods adopt |
| 4201 |
+ |
a different approach in that an unphysical {\it flux} is imposed |
| 4202 |
+ |
between different regions or ``slabs'' of the simulation box. The |
| 4203 |
+ |
response of the system is to develop a temperature or momentum {\it |
| 4204 |
+ |
gradient} between the two regions. Since the amount of the applied |
| 4205 |
+ |
flux is known exactly, and the measurement of gradient is generally |
| 4206 |
+ |
less complicated, imposed-flux methods typically take shorter |
| 4207 |
+ |
simulation times to obtain converged results for transport properties. |
| 4208 |
+ |
|
| 4209 |
+ |
\begin{figure} |
| 4210 |
+ |
\includegraphics[width=\linewidth]{rnemdDemo} |
| 4211 |
+ |
\caption[Illustration of energy exchange in the VSS RNEMD method]{The (VSS) RNEMD approach imposes unphysical transfer of both |
| 4212 |
+ |
linear momentum and kinetic energy between a ``hot'' slab and a |
| 4213 |
+ |
``cold'' slab in the simulation box. The system responds to this |
| 4214 |
+ |
imposed flux by generating both momentum and temperature gradients. |
| 4215 |
+ |
The slope of the gradients can then be used to compute transport |
| 4216 |
+ |
properties (e.g. shear viscosity and thermal conductivity).} |
| 4217 |
+ |
\label{rnemdDemo} |
| 4218 |
+ |
\end{figure} |
| 4219 |
+ |
|
| 4220 |
+ |
\section{\label{section:algo}Three algorithms for carrying out RNEMD simulations} |
| 4221 |
+ |
\subsection{\label{subsection:swapping}The swapping algorithm} |
| 4222 |
+ |
The original ``swapping'' approaches by M\"{u}ller-Plathe {\it et |
| 4223 |
+ |
al.}\cite{ISI:000080382700030,MullerPlathe:1997xw} can be understood |
| 4224 |
+ |
as a sequence of imaginary elastic collisions between particles in |
| 4225 |
+ |
opposite slabs. In each collision, the entire momentum vectors of |
| 4226 |
+ |
both particles may be exchanged to generate a thermal |
| 4227 |
+ |
flux. Alternatively, a single component of the momentum vectors may be |
| 4228 |
+ |
exchanged to generate a shear flux. This algorithm turns out to be |
| 4229 |
+ |
quite useful in many simulations. However, the M\"{u}ller-Plathe |
| 4230 |
+ |
swapping approach perturbs the system away from ideal |
| 4231 |
+ |
Maxwell-Boltzmann distributions, and this may leads to undesirable |
| 4232 |
+ |
side-effects when the applied flux becomes large.\cite{Maginn:2010} |
| 4233 |
+ |
This limits the applicability of the swapping algorithm, so in OpenMD, |
| 4234 |
+ |
we have implemented two additional algorithms for RNEMD in addition to the |
| 4235 |
+ |
original swapping approach. |
| 4236 |
+ |
|
| 4237 |
+ |
\subsection{\label{subsection:nivs}Non-Isotropic Velocity Scaling (NIVS)} |
| 4238 |
+ |
Instead of having momentum exchange between {\it individual particles} |
| 4239 |
+ |
in each slab, the NIVS algorithm applies velocity scaling to all of |
| 4240 |
+ |
the selected particles in both slabs.\cite{kuang:164101} A combination of linear |
| 4241 |
+ |
momentum, kinetic energy, and flux constraint equations governs the |
| 4242 |
+ |
amount of velocity scaling performed at each step. Interested readers |
| 4243 |
+ |
should consult ref. \citealp{kuang:164101} for further details on the |
| 4244 |
+ |
methodology. |
| 4245 |
+ |
|
| 4246 |
+ |
NIVS has been shown to be very effective at producing thermal |
| 4247 |
+ |
gradients and for computing thermal conductivities, particularly for |
| 4248 |
+ |
heterogeneous interfaces. Although the NIVS algorithm can also be |
| 4249 |
+ |
applied to impose a directional momentum flux, thermal anisotropy was |
| 4250 |
+ |
observed in relatively high flux simulations, and the method is not |
| 4251 |
+ |
suitable for imposing a shear flux or for computing shear viscosities. |
| 4252 |
+ |
|
| 4253 |
+ |
\subsection{\label{subsection:vss}Velocity Shearing and Scaling (VSS)} |
| 4254 |
+ |
The third RNEMD algorithm implemented in OpenMD utilizes a series of |
| 4255 |
+ |
simultaneous velocity shearing and scaling exchanges between the two |
| 4256 |
+ |
slabs.\cite{2012MolPh.110..691K} This method results in a set of simpler equations to satisfy |
| 4257 |
+ |
the conservation constraints while creating a desired flux between the |
| 4258 |
+ |
two slabs. |
| 4259 |
+ |
|
| 4260 |
+ |
The VSS approach is versatile in that it may be used to implement both |
| 4261 |
+ |
thermal and shear transport either separately or simultaneously. |
| 4262 |
+ |
Perturbations of velocities away from the ideal Maxwell-Boltzmann |
| 4263 |
+ |
distributions are minimal, and thermal anisotropy is kept to a |
| 4264 |
+ |
minimum. This ability to generate simultaneous thermal and shear |
| 4265 |
+ |
fluxes has been utilized to map out the shear viscosity of SPC/E water |
| 4266 |
+ |
over a wide range of temperatures (90~K) just with a single simulation. |
| 4267 |
+ |
VSS-RNEMD also allows the directional momentum flux to have |
| 4268 |
+ |
arbitary directions, which could aid in the study of anisotropic solid |
| 4269 |
+ |
surfaces in contact with liquid environments. |
| 4270 |
+ |
|
| 4271 |
+ |
\section{\label{section:usingRNEMD}Using OpenMD to perform a RNEMD simulation} |
| 4272 |
+ |
\subsection{\label{section:rnemdParams} What the user needs to specify} |
| 4273 |
+ |
To carry out a RNEMD simulation, |
| 4274 |
+ |
a user must specify a number of parameters in the MetaData (.md) file. |
| 4275 |
+ |
Because the RNEMD methods have a large number of parameters, these |
| 4276 |
+ |
must be enclosed in a {\it separate} {\tt RNEMD\{...\}} block. The most important |
| 4277 |
+ |
parameters to specify are the {\tt useRNEMD}, {\tt fluxType} and flux |
| 4278 |
+ |
parameters. Most other parameters (summarized in table |
| 4279 |
+ |
\ref{table:rnemd}) have reasonable default values. {\tt fluxType} |
| 4280 |
+ |
sets up the kind of exchange that will be carried out between the two |
| 4281 |
+ |
slabs (either Kinetic Energy ({\tt KE}) or momentum ({\tt Px, Py, Pz, |
| 4282 |
+ |
Pvector}), or some combination of these). The flux is specified |
| 4283 |
+ |
with the use of three possible parameters: {\tt kineticFlux} for |
| 4284 |
+ |
kinetic energy exchange, as well as {\tt momentumFlux} or {\tt |
| 4285 |
+ |
momentumFluxVector} for simulations with directional exchange. |
| 4286 |
+ |
|
| 4287 |
+ |
\subsection{\label{section:rnemdResults} Processing the results} |
| 4288 |
+ |
OpenMD will generate a {\tt .rnemd} |
| 4289 |
+ |
file with the same prefix as the original {\tt .md} file. This file |
| 4290 |
+ |
contains a running average of properties of interest computed within a |
| 4291 |
+ |
set of bins that divide the simulation cell along the $z$-axis. The |
| 4292 |
+ |
first column of the {\tt .rnemd} file is the $z$ coordinate of the |
| 4293 |
+ |
center of each bin, while following columns may contain the average |
| 4294 |
+ |
temperature, velocity, or density within each bin. The output format |
| 4295 |
+ |
in the {\tt .rnemd} file can be altered with the {\tt outputFields}, |
| 4296 |
+ |
{\tt outputBins}, and {\tt outputFileName} parameters. A report at |
| 4297 |
+ |
the top of the {\tt .rnemd} file contains the current exchange totals |
| 4298 |
+ |
as well as the average flux applied during the simulation. Using the |
| 4299 |
+ |
slope of the temperature or velocity gradient obtaine from the {\tt |
| 4300 |
+ |
.rnemd} file along with the applied flux, the user can very simply |
| 4301 |
+ |
arrive at estimates of thermal conductivities ($\lambda$), |
| 4302 |
+ |
\begin{equation} |
| 4303 |
+ |
J_z = -\lambda \frac{\partial T}{\partial z}, |
| 4304 |
+ |
\end{equation} |
| 4305 |
+ |
and shear viscosities ($\eta$), |
| 4306 |
+ |
\begin{equation} |
| 4307 |
+ |
j_z(p_x) = -\eta \frac{\partial \langle v_x \rangle}{\partial z}. |
| 4308 |
+ |
\end{equation} |
| 4309 |
+ |
Here, the quantities on the left hand side are the actual flux values |
| 4310 |
+ |
(in the header of the {\tt .rnemd} file), while the slopes are |
| 4311 |
+ |
obtained from linear fits to the gradients observed in the {\tt |
| 4312 |
+ |
.rnemd} file. |
| 4313 |
+ |
|
| 4314 |
+ |
More complicated simulations (including interfaces) require a bit more |
| 4315 |
+ |
care. Here the second derivative may be required to compute the |
| 4316 |
+ |
interfacial thermal conductance, |
| 4317 |
+ |
\begin{align} |
| 4318 |
+ |
G^\prime &= \left(\nabla\lambda \cdot \mathbf{\hat{n}}\right)_{z_0} \\ |
| 4319 |
+ |
&= \frac{\partial}{\partial z}\left(-\frac{J_z}{ |
| 4320 |
+ |
\left(\frac{\partial T}{\partial z}\right)}\right)_{z_0} \\ |
| 4321 |
+ |
&= J_z\left(\frac{\partial^2 T}{\partial z^2}\right)_{z_0} \Big/ |
| 4322 |
+ |
\left(\frac{\partial T}{\partial z}\right)_{z_0}^2. |
| 4323 |
+ |
\label{derivativeG} |
| 4324 |
+ |
\end{align} |
| 4325 |
+ |
where $z_0$ is the location of the interface between two materials and |
| 4326 |
+ |
$\mathbf{\hat{n}}$ is a unit vector normal to the interface. We |
| 4327 |
+ |
suggest that users interested in interfacial conductance consult |
| 4328 |
+ |
reference \citealp{kuang:AuThl} for other approaches to computing $G$. |
| 4329 |
+ |
Users interested in {\it friction coefficients} at heterogeneous |
| 4330 |
+ |
interfaces may also find reference \citealp{2012MolPh.110..691K} |
| 4331 |
+ |
useful. |
| 4332 |
+ |
|
| 4333 |
+ |
\newpage |
| 4334 |
+ |
|
| 4335 |
+ |
\begin{longtable}[c]{JKLM} |
| 4336 |
+ |
\caption{Meta-data Keywords: Parameters for RNEMD simulations}\\ |
| 4337 |
+ |
\multicolumn{4}{c}{The following keywords must be enclosed inside a {\tt RNEMD\{...\}} block.} |
| 4338 |
+ |
\\ \hline |
| 4339 |
+ |
{\bf keyword} & {\bf units} & {\bf use} & {\bf remarks} \\ \hline |
| 4340 |
+ |
\endhead |
| 4341 |
+ |
\hline |
| 4342 |
+ |
\endfoot |
| 4343 |
+ |
{\tt useRNEMD} & logical & perform RNEMD? & default is ``false'' \\ |
| 4344 |
+ |
{\tt objectSelection} & string & see section \ref{section:syntax} |
| 4345 |
+ |
for selection syntax & default is ``select all'' \\ |
| 4346 |
+ |
{\tt method} & string & exchange method & one of the following: |
| 4347 |
+ |
{\tt Swap, NIVS,} or {\tt VSS} (default is {\tt VSS}) \\ |
| 4348 |
+ |
{\tt fluxType} & string & what is being exchanged between slabs? & one |
| 4349 |
+ |
of the following: {\tt KE, Px, Py, Pz, Pvector, KE+Px, KE+Py, KE+Pvector} \\ |
| 4350 |
+ |
{\tt kineticFlux} & kcal mol$^{-1}$ \AA$^{-2}$ fs$^{-1}$ & specify the kinetic energy flux & \\ |
| 4351 |
+ |
{\tt momentumFlux} & amu \AA$^{-1}$ fs$^{-2}$ & specify the momentum flux & \\ |
| 4352 |
+ |
{\tt momentumFluxVector} & amu \AA$^{-1}$ fs$^{-2}$ & specify the momentum flux when |
| 4353 |
+ |
{\tt Pvector} is part of the exchange & Vector3d input\\ |
| 4354 |
+ |
{\tt exchangeTime} & fs & how often to perform the exchange & default is 100 fs\\ |
| 4355 |
|
|
| 4356 |
+ |
{\tt slabWidth} & $\mbox{\AA}$ & width of the two exchange slabs & default is $\mathsf{H}_{zz} / 10.0$ \\ |
| 4357 |
+ |
{\tt slabAcenter} & $\mbox{\AA}$ & center of the end slab & default is 0 \\ |
| 4358 |
+ |
{\tt slabBcenter} & $\mbox{\AA}$ & center of the middle slab & default is $\mathsf{H}_{zz} / 2$ \\ |
| 4359 |
+ |
{\tt outputFileName} & string & file name for output histograms & default is the same prefix as the |
| 4360 |
+ |
.md file, but with the {\tt .rnemd} extension \\ |
| 4361 |
+ |
{\tt outputBins} & int & number of $z$-bins in the output histogram & |
| 4362 |
+ |
default is 20 \\ |
| 4363 |
+ |
{\tt outputFields} & string & columns to print in the {\tt .rnemd} |
| 4364 |
+ |
file where each column is separated by a pipe ($\mid$) symbol. & Allowed column names are: {\sc z, temperature, velocity, density} \\ |
| 4365 |
+ |
\label{table:rnemd} |
| 4366 |
+ |
\end{longtable} |
| 4367 |
|
|
| 4368 |
|
\chapter{\label{section:minimizer}Energy Minimization} |
| 4369 |
|
|
| 4370 |
< |
As one of the basic procedures of molecular modeling, energy |
| 3083 |
< |
minimization is used to identify local configurations that are stable |
| 4370 |
> |
Energy minimization is used to identify local configurations that are stable |
| 4371 |
|
points on the potential energy surface. There is a vast literature on |
| 4372 |
|
energy minimization algorithms have been developed to search for the |
| 4373 |
|
global energy minimum as well as to find local structures which are |
| 4417 |
|
the steepest descent method can be superior to the conjugate gradient |
| 4418 |
|
method. Hence, the steepest descent method is often used for the first |
| 4419 |
|
10-100 steps of minimization. Another useful feature of minimization |
| 4420 |
< |
methods in {\sc OpenMD} is that a modified {\sc shake} algorithm can be |
| 4421 |
< |
applied during the minimization to constraint the bond lengths if this |
| 4422 |
< |
is required by the force field. Meta-data parameters concerning the |
| 4423 |
< |
minimizer are given in Table~\ref{table:minimizeParams} |
| 4420 |
> |
methods in {\sc OpenMD} is that a modified {\sc shake} algorithm can |
| 4421 |
> |
be applied during the minimization to constraint the bond lengths if |
| 4422 |
> |
this is required by the force field. Meta-data parameters concerning |
| 4423 |
> |
the minimizer are given in Table~\ref{table:minimizeParams} Because |
| 4424 |
> |
the minimizer methods have a large number of parameters, these must be |
| 4425 |
> |
enclosed in a {\it separate} {\tt minimizer\{...\}} block. |
| 4426 |
|
|
| 4427 |
|
\begin{longtable}[c]{ABCD} |
| 4428 |
|
\caption{Meta-data Keywords: Energy Minimizer Parameters} |
| 4460 |
|
|
| 4461 |
|
\begin{description} |
| 4462 |
|
\item[{\bf Dump2XYZ}] Converts an {\sc OpenMD} dump file into a file |
| 4463 |
< |
suitable for viewing in a molecular dynamics viewer like Jmol |
| 4463 |
> |
suitable for viewing in a molecular dynamics viewer like Jmol or VMD |
| 4464 |
|
\item[{\bf StaticProps}] Computes static properties like the pair |
| 4465 |
|
distribution function, $g(r)$. |
| 4466 |
+ |
\item[{\bf SequentialProps}] Computes a time history of static |
| 4467 |
+ |
properties from a dump file. |
| 4468 |
|
\item[{\bf DynamicProps}] Computes time correlation functions like the |
| 4469 |
|
velocity autocorrelation function, $\langle v(t) \cdot v(0)\rangle$, |
| 4470 |
|
or the mean square displacement $\langle |r(t) - r(0)|^{2} \rangle$. |
| 4498 |
|
\begin{figure} |
| 4499 |
|
\centering |
| 4500 |
|
\includegraphics[width=3in]{heirarchy.pdf} |
| 4501 |
< |
\caption[Class heirarchy for StuntDoubles in {\sc OpenMD}-4]{ \\ The |
| 4502 |
< |
class heirarchy of StuntDoubles in {\sc OpenMD}-4. The selection |
| 4501 |
> |
\caption[Class heirarchy for StuntDoubles in {\sc OpenMD}]{ \\ The |
| 4502 |
> |
class heirarchy of StuntDoubles in {\sc OpenMD}. The selection |
| 4503 |
|
syntax allows the user to select any of the objects that are descended |
| 4504 |
|
from a StuntDouble.} |
| 4505 |
|
\label{fig:heirarchy} |
| 4514 |
|
DirectionalAtom}s which behaves as a single unit. |
| 4515 |
|
\end{itemize} |
| 4516 |
|
|
| 4517 |
< |
Every Molecule, Atom and DirectionalAtom in {\sc OpenMD} have their own names |
| 4518 |
< |
which are specified in the {\tt .md} file. In contrast, RigidBodies are |
| 4519 |
< |
denoted by their membership and index inside a particular molecule: |
| 4520 |
< |
[MoleculeName]\_RB\_[index] (the contents inside the brackets |
| 4521 |
< |
depend on the specifics of the simulation). The names of rigid bodies are |
| 4517 |
> |
Every Molecule, Atom and DirectionalAtom in {\sc OpenMD} have their |
| 4518 |
> |
own names which are specified in the {\tt .md} file. In contrast, |
| 4519 |
> |
other groupings of atoms are denoted by their membership and index |
| 4520 |
> |
inside a particular molecule. For example, RigidBodies are denoted |
| 4521 |
> |
[MoleculeName]\_RB\_[index] (the contents inside the brackets depend |
| 4522 |
> |
on the specifics of the simulation). The names of rigid bodies are |
| 4523 |
|
generated automatically. For example, the name of the first rigid body |
| 4524 |
< |
in a DMPC molecule is DMPC\_RB\_0. |
| 4524 |
> |
in a DMPC molecule is DMPC\_RB\_0. Similarly, bonds can be denoted |
| 4525 |
> |
as: [MoleculeName]\_Bond\_[index], bends as |
| 4526 |
> |
[MoleculeName]\_Bend\_[index], torsions as |
| 4527 |
> |
[MoleculeName]\_Torsion\_[index], and inversions as |
| 4528 |
> |
[MoleculeName]\_Inversion\_[index]. These selection names will select |
| 4529 |
> |
all of the atoms involved in that group, as well as the grouping |
| 4530 |
> |
itself. |
| 4531 |
|
|
| 4532 |
|
\section{\label{section:syntax}Syntax of the Select Command} |
| 4533 |
|
|
| 4638 |
|
\begin{center} |
| 4639 |
|
\begin{tabular}{|l|l|} |
| 4640 |
|
\hline |
| 4641 |
< |
{\bf property} & mass, charge \\ |
| 4641 |
> |
{\bf property} & mass, charge, x, y, z, r, wrappedx, wrappedy, wrappedz \\ |
| 4642 |
|
{\bf comparison operator} & ``$>$'', ``$<$'', ``$=$'', ``$>=$'', |
| 4643 |
|
``$<=$'', ``$!=$'' \\ |
| 4644 |
|
\hline |
| 4686 |
|
-z & {\tt -{}-zconstraint} & replace the atom types of zconstraint molecules (default=off) \\ |
| 4687 |
|
-r & {\tt -{}-rigidbody} & add a pseudo COM atom to rigidbody (default=off) \\ |
| 4688 |
|
-t & {\tt -{}-watertype} & replace the atom type of water model (default=on) \\ |
| 4689 |
< |
-b & {\tt -{}-basetype} & using base atom type (default=off) \\ |
| 4689 |
> |
-b & {\tt -{}-basetype} & using base atom type |
| 4690 |
> |
(default=off) \\ |
| 4691 |
> |
-v& {\tt -{}-velocities} & Print velocities in xyz file (default=off)\\ |
| 4692 |
> |
-f& {\tt -{}-forces} & Print forces xyz file (default=off)\\ |
| 4693 |
> |
-u& {\tt -{}-vectors} & Print vectors (dipoles, etc) in xyz file |
| 4694 |
> |
(default=off)\\ |
| 4695 |
> |
-c& {\tt -{}-charges} & Print charges in xyz file (default=off)\\ |
| 4696 |
> |
-e& {\tt -{}-efield} & Print electric field vector in xyz file |
| 4697 |
> |
(default=off)\\ |
| 4698 |
|
& {\tt -{}-repeatX=INT} & The number of images to repeat in the x direction (default=`0') \\ |
| 4699 |
|
& {\tt -{}-repeatY=INT} & The number of images to repeat in the y direction (default=`0') \\ |
| 4700 |
|
& {\tt -{}-repeatZ=INT} & The number of images to repeat in the z direction (default=`0') \\ |
| 4786 |
|
& {\tt -{}-sele1=selection script} & select the first StuntDouble set \\ |
| 4787 |
|
& {\tt -{}-sele2=selection script} & select the second StuntDouble set \\ |
| 4788 |
|
& {\tt -{}-sele3=selection script} & select the third StuntDouble set \\ |
| 4789 |
< |
& {\tt -{}-refsele=selection script} & select reference (can only be used with {\tt -{}-gxyz}) \\ |
| 4789 |
> |
& {\tt -{}-refsele=selection script} & select reference (can only |
| 4790 |
> |
be used with {\tt -{}-gxyz}) \\ |
| 4791 |
> |
& {\tt -{}-comsele=selection script} |
| 4792 |
> |
& select stunt doubles for center-of-mass |
| 4793 |
> |
reference point\\ |
| 4794 |
> |
& {\tt -{}-seleoffset=INT} & global index offset for a second object (used |
| 4795 |
> |
to define a vector between sites in molecule)\\ |
| 4796 |
> |
|
| 4797 |
|
& {\tt -{}-molname=STRING} & molecule name \\ |
| 4798 |
|
& {\tt -{}-begin=INT} & begin internal index \\ |
| 4799 |
|
& {\tt -{}-end=INT} & end internal index \\ |
| 4800 |
+ |
& {\tt -{}-radius=DOUBLE} & nanoparticle radius\\ |
| 4801 |
|
\hline |
| 4802 |
|
\multicolumn{3}{|l|}{One option from the following group of options is required:} \\ |
| 4803 |
|
\hline |
| 4804 |
< |
& {\tt -{}-gofr} & $g(r)$ \\ |
| 4805 |
< |
& {\tt -{}-r\_theta} & $g(r, \cos(\theta))$ \\ |
| 4806 |
< |
& {\tt -{}-r\_omega} & $g(r, \cos(\omega))$ \\ |
| 4807 |
< |
& {\tt -{}-theta\_omega} & $g(\cos(\theta), \cos(\omega))$ \\ |
| 4804 |
> |
& {\tt -{}-bo} & bond order parameter ({\tt -{}-rcut} must be specified) \\ |
| 4805 |
> |
& {\tt -{}-bor} & bond order parameter as a function of |
| 4806 |
> |
radius ({\tt -{}-rcut} must be specified) \\ |
| 4807 |
> |
& {\tt -{}-bad} & $N(\theta)$ bond angle density within ({\tt -{}-rcut} must be specified) \\ |
| 4808 |
> |
& {\tt -{}-count} & count of molecules matching selection |
| 4809 |
> |
criteria (and associated statistics) \\ |
| 4810 |
> |
-g& {\tt -{}-gofr} & $g(r)$ \\ |
| 4811 |
> |
& {\tt -{}-gofz} & $g(z)$ \\ |
| 4812 |
> |
& {\tt -{}-r\_theta} & $g(r, \cos(\theta))$ \\ |
| 4813 |
> |
& {\tt -{}-r\_omega} & $g(r, \cos(\omega))$ \\ |
| 4814 |
> |
& {\tt -{}-r\_z} & $g(r, z)$ \\ |
| 4815 |
> |
& {\tt -{}-theta\_omega} & $g(\cos(\theta), \cos(\omega))$ \\ |
| 4816 |
|
& {\tt -{}-gxyz} & $g(x, y, z)$ \\ |
| 4817 |
< |
& {\tt -{}-p2} & $P_2$ order parameter ({\tt -{}-sele1} and {\tt -{}-sele2} must be specified) \\ |
| 4817 |
> |
& {\tt -{}-twodgofr} & 2D $g(r)$ (Slab width {\tt -{}-dz} must be specified)\\ |
| 4818 |
> |
-p& {\tt -{}-p2} & $P_2$ order parameter ({\tt -{}-sele1} must be specified, {\tt -{}-sele2} is optional) \\ |
| 4819 |
> |
& {\tt -{}-rp2} & Ripple order parameter ({\tt -{}-sele1} and {\tt -{}-sele2} must be specified) \\ |
| 4820 |
|
& {\tt -{}-scd} & $S_{CD}$ order parameter(either {\tt -{}-sele1}, {\tt -{}-sele2}, {\tt -{}-sele3} are specified or {\tt -{}-molname}, {\tt -{}-begin}, {\tt -{}-end} are specified) \\ |
| 4821 |
< |
& {\tt -{}-density} & density plot ({\tt -{}-sele1} must be specified) \\ |
| 4822 |
< |
& {\tt -{}-slab\_density} & slab density ({\tt -{}-sele1} must be specified) |
| 4821 |
> |
-d& {\tt -{}-density} & density plot \\ |
| 4822 |
> |
& {\tt -{}-slab\_density} & slab density \\ |
| 4823 |
> |
& {\tt -{}-p\_angle} & $p(\cos(\theta))$ ($\theta$ |
| 4824 |
> |
is the angle between molecular axis and radial vector from origin\\ |
| 4825 |
> |
& {\tt -{}-hxy} & Calculates the undulation spectrum, $h(x,y)$, of an interface \\ |
| 4826 |
> |
& {\tt -{}-rho\_r} & $\rho(r)$\\ |
| 4827 |
> |
& {\tt -{}-angle\_r} & $\theta(r)$ (spatially resolves the |
| 4828 |
> |
angle between the molecular axis and the radial vector from the |
| 4829 |
> |
origin)\\ |
| 4830 |
> |
& {\tt -{}-hullvol} & hull volume of nanoparticle\\ |
| 4831 |
> |
& {\tt -{}-rodlength} & length of nanorod\\ |
| 4832 |
> |
-Q& {\tt -{}-tet\_param} & tetrahedrality order parameter ($Q$)\\ |
| 4833 |
> |
& {\tt -{}-tet\_param\_z} & spatially-resolved tetrahedrality order |
| 4834 |
> |
parameter $Q(z)$\\ |
| 4835 |
> |
& {\tt -{}-rnemdz} & slab-resolved RNEMD statistics (temperature, |
| 4836 |
> |
density, velocity)\\ |
| 4837 |
> |
& {\tt -{}-rnemdr} & shell-resolved RNEMD statistics (temperature, |
| 4838 |
> |
density, angular velocity) |
| 4839 |
|
\end{longtable} |
| 4840 |
|
|
| 4841 |
|
\subsection{\label{section:DynamicProps}DynamicProps} |
| 4876 |
|
-o& {\tt -{}-output=filename} & output file name \\ |
| 4877 |
|
& {\tt -{}-sele1=selection script} & select first StuntDouble set \\ |
| 4878 |
|
& {\tt -{}-sele2=selection script} & select second StuntDouble set (if sele2 is not set, use script from sele1) \\ |
| 4879 |
+ |
& {\tt -{}-order=INT} & Lengendre Polynomial Order\\ |
| 4880 |
+ |
-z& {\tt -{}-nzbins=INT} & Number of $z$ bins (default=`100`)\\ |
| 4881 |
+ |
-m& {\tt -{}-memory=memory specification} |
| 4882 |
+ |
&Available memory |
| 4883 |
+ |
(default=`2G`)\\ |
| 4884 |
|
\hline |
| 4885 |
|
\multicolumn{3}{|l|}{One option from the following group of options is required:} \\ |
| 4886 |
|
\hline |
| 4887 |
< |
-r& {\tt -{}-rcorr} & compute mean square displacement \\ |
| 4888 |
< |
-v& {\tt -{}-vcorr} & compute velocity correlation function \\ |
| 4889 |
< |
-d& {\tt -{}-dcorr} & compute dipole correlation function |
| 4887 |
> |
-s& {\tt -{}-selecorr} & selection correlation function \\ |
| 4888 |
> |
-r& {\tt -{}-rcorr} & compute mean squared displacement \\ |
| 4889 |
> |
-v& {\tt -{}-vcorr} & velocity autocorrelation function \\ |
| 4890 |
> |
-d& {\tt -{}-dcorr} & dipole correlation function \\ |
| 4891 |
> |
-l& {\tt -{}-lcorr} & Lengendre correlation function \\ |
| 4892 |
> |
& {\tt -{}-lcorrZ} & Lengendre correlation function binned by $z$ \\ |
| 4893 |
> |
& {\tt -{}-cohZ} & Lengendre correlation function for OH bond vectors binned by $z$\\ |
| 4894 |
> |
-M& {\tt -{}-sdcorr} & System dipole correlation function\\ |
| 4895 |
> |
& {\tt -{}-r\_rcorr} & Radial mean squared displacement\\ |
| 4896 |
> |
& {\tt -{}-thetacorr} & Angular mean squared displacement\\ |
| 4897 |
> |
& {\tt -{}-drcorr} & Directional mean squared displacement for particles with unit vectors\\ |
| 4898 |
> |
& {\tt -{}-helfandEcorr} & Helfand moment for thermal conductvity\\ |
| 4899 |
> |
-p& {\tt -{}-momentum} & Helfand momentum for viscosity\\ |
| 4900 |
> |
& {\tt -{}-stresscorr} & Stress tensor correlation function |
| 4901 |
|
\end{longtable} |
| 4902 |
|
|
| 4903 |
|
\chapter{\label{section:PreparingInput} Preparing Input Configurations} |
| 4904 |
|
|
| 4905 |
< |
{\sc OpenMD} version 4 comes with a few utility programs to aid in |
| 4906 |
< |
setting up initial configuration and meta-data files. Usually, a user |
| 4907 |
< |
is interested in either importing a structure from some other format |
| 4905 |
> |
{\sc OpenMD} comes with a few utility programs to aid in setting up |
| 4906 |
> |
initial configuration and meta-data files. Usually, a user is |
| 4907 |
> |
interested in either importing a structure from some other format |
| 4908 |
|
(usually XYZ or PDB), or in building an initial configuration in some |
| 4909 |
< |
perfect crystalline lattice. The programs bundled with {\sc OpenMD} |
| 4910 |
< |
which import coordinate files are {\tt atom2md}, {\tt xyz2md}, and |
| 4911 |
< |
{\tt pdb2md}. The programs which generate perfect crystals are called |
| 4912 |
< |
{\tt SimpleBuilder} and {\tt RandomBuilder} |
| 4909 |
> |
perfect crystalline lattice or nanoparticle geometry. The program |
| 4910 |
> |
bundled with {\sc OpenMD} that imports coordinate files is {\tt |
| 4911 |
> |
atom2md}, which is built if the initial CMake configuration can find |
| 4912 |
> |
the openbabel libraries. The programs which generate perfect crystals |
| 4913 |
> |
are called {\tt SimpleBuilder} and {\tt RandomBuilder}. There are |
| 4914 |
> |
programs to construct nanoparticles of various sizes and geometries |
| 4915 |
> |
also. These are {\tt nanoparticleBuilder}, {\tt icosahedralBuilder}, |
| 4916 |
> |
and {\tt nanorodBuilder}. |
| 4917 |
|
|
| 4918 |
< |
\section{\label{section:atom2md}atom2md, xyz2md, and pdb2md} |
| 4918 |
> |
\section{\label{section:atom2md}atom2md} |
| 4919 |
|
|
| 4920 |
< |
{\tt atom2md}, {\tt xyz2md}, and {\tt pdb2md} attempt to construct |
| 4921 |
< |
{\tt .md} files from a single file containing only atomic coordinate |
| 4922 |
< |
information. To do this task, they make reasonable guesses about |
| 4923 |
< |
bonding from the distance between atoms in the coordinate, and attempt |
| 4924 |
< |
to identify other terms in the potential energy from the topology of |
| 4925 |
< |
the graph of discovered bonds. This procedure is not perfect, and the |
| 4926 |
< |
user should check the discovered bonding topology that is contained in |
| 4927 |
< |
the {\tt $<$MetaData$>$} block in the file that is generated. |
| 4920 |
> |
{\tt atom2md} attempts to construct {\tt .md} files from a single file |
| 4921 |
> |
containing only atomic coordinate information. To do this task, they |
| 4922 |
> |
make reasonable guesses about bonding from the distance between atoms |
| 4923 |
> |
in the coordinate, and attempt to identify other terms in the |
| 4924 |
> |
potential energy from the topology of the graph of discovered bonds. |
| 4925 |
> |
This procedure is not perfect, and the user should check the |
| 4926 |
> |
discovered bonding topology that is contained in the {\tt |
| 4927 |
> |
$<$MetaData$>$} block in the file that is generated. |
| 4928 |
|
|
| 4929 |
|
Typically, the user would run: |
| 4930 |
|
|
| 4964 |
|
using -H$<$format-type$>$, e.g. -Hpdb} |
| 4965 |
|
\end{longtable} |
| 4966 |
|
|
| 3607 |
– |
The specific programs {\tt xyz2md} and {\tt pdb2md} are identical |
| 3608 |
– |
to {\tt atom2md}, but they use a specific input format and do not |
| 3609 |
– |
expect the the input specifier on the command line. |
| 3610 |
– |
|
| 4967 |
|
\section{\label{section:SimpleBuilder}SimpleBuilder} |
| 4968 |
|
|
| 4969 |
|
{\tt SimpleBuilder} creates simple lattice structures. It requires an |
| 4990 |
|
& {\tt -{}-nz=INT} & number of unit cells in z |
| 4991 |
|
\end{longtable} |
| 4992 |
|
|
| 4993 |
+ |
\section{\label{section:icosahedralBuilder}icosahedralBuilder} |
| 4994 |
+ |
|
| 4995 |
+ |
{\tt icosahedralBuilder} creates single-component geometric solids |
| 4996 |
+ |
that can be useful in simulating nanostructures. Like the other |
| 4997 |
+ |
builders, it requires an initial, but skeletal {\sc OpenMD} file to |
| 4998 |
+ |
specify the component that is to be placed on the lattice. The total |
| 4999 |
+ |
number of placed molecules will be shown at the top of the |
| 5000 |
+ |
configuration file that is generated, and that number may not match |
| 5001 |
+ |
the original meta-data file, so a new meta-data file is also generated |
| 5002 |
+ |
which matches the lattice structure. |
| 5003 |
+ |
|
| 5004 |
+ |
The options available for icosahedralBuilder are as follows: |
| 5005 |
+ |
\begin{longtable}[c]{|EFG|} |
| 5006 |
+ |
\caption{icosahedralBuilder Command-line Options} |
| 5007 |
+ |
\\ \hline |
| 5008 |
+ |
{\bf option} & {\bf verbose option} & {\bf behavior} \\ \hline |
| 5009 |
+ |
\endhead |
| 5010 |
+ |
\hline |
| 5011 |
+ |
\endfoot |
| 5012 |
+ |
-h& {\tt -{}-help} & Print help and exit\\ |
| 5013 |
+ |
-V& {\tt -{}-version} & Print version and exit\\ |
| 5014 |
+ |
-o& {\tt -{}-output=STRING} & Output file name\\ |
| 5015 |
+ |
-n& {\tt -{}-shells=INT} & Nanoparticle shells\\ |
| 5016 |
+ |
-d& {\tt -{}-latticeConstant=DOUBLE} & Lattice spacing in Angstroms for cubic lattice.\\ |
| 5017 |
+ |
-c& {\tt -{}-columnAtoms=INT} & Number of atoms along central |
| 5018 |
+ |
column (Decahedron only)\\ |
| 5019 |
+ |
-t& {\tt -{}-twinAtoms=INT} & Number of atoms along twin |
| 5020 |
+ |
boundary (Decahedron only) \\ |
| 5021 |
+ |
-p& {\tt -{}-truncatedPlanes=INT} & Number of truncated planes (Curling-stone Decahedron only)\\ |
| 5022 |
+ |
\hline |
| 5023 |
+ |
\multicolumn{3}{|l|}{One option from the following group of options is required:} \\ |
| 5024 |
+ |
\hline |
| 5025 |
+ |
& {\tt -{}-ico} & Create an Icosahedral cluster \\ |
| 5026 |
+ |
& {\tt -{}-deca} & Create a regualar Decahedral cluster\\ |
| 5027 |
+ |
& {\tt -{}-ino} & Create an Ino Decahedral cluster\\ |
| 5028 |
+ |
& {\tt -{}-marks} & Create a Marks Decahedral cluster\\ |
| 5029 |
+ |
& {\tt -{}-stone} & Create a Curling-stone Decahedral cluster |
| 5030 |
+ |
\end{longtable} |
| 5031 |
+ |
|
| 5032 |
+ |
|
| 5033 |
|
\section{\label{section:Hydro}Hydro} |
| 5034 |
|
{\tt Hydro} generates resistance tensor ({\tt .diff}) files which are |
| 5035 |
|
required when using the Langevin integrator using complex rigid |
| 5135 |
|
encourage other researchers to contribute their own code and make it a |
| 5136 |
|
more powerful package for everyone in the molecular dynamics community |
| 5137 |
|
to use. All source code for {\sc OpenMD} is available for download at |
| 5138 |
< |
{\tt http://openmd.net}. |
| 5138 |
> |
{\tt http://openmd.org}. |
| 5139 |
|
|
| 5140 |
|
\chapter{Acknowledgments} |
| 5141 |
|
|
| 5142 |
|
Development of {\sc OpenMD} was funded by a New Faculty Award from the |
| 5143 |
|
Camille and Henry Dreyfus Foundation and by the National Science |
| 5144 |
< |
Foundation under grant CHE-0134881. Computation time was provided by |
| 5145 |
< |
the Notre Dame Bunch-of-Boxes (B.o.B) computer cluster under NSF grant |
| 5146 |
< |
DMR-0079647. |
| 5144 |
> |
Foundation under grants CHE-0134881, CHE-0848243, and |
| 5145 |
> |
CHE-1362211. Computation time was provided by the Notre Dame |
| 5146 |
> |
Bunch-of-Boxes (B.o.B) computer cluster under NSF grant DMR-0079647. |
| 5147 |
|
|
| 5148 |
< |
|
| 3753 |
< |
\bibliographystyle{jcc} |
| 5148 |
> |
\bibliographystyle{aip} |
| 5149 |
|
\bibliography{openmdDoc} |
| 5150 |
|
|
| 5151 |
|
\end{document} |
