| 63 |
|
\newcolumntype{M}{p{1.55in}} |
| 64 |
|
|
| 65 |
|
|
| 66 |
< |
\title{{\sc OpenMD-2.2}: Molecular Dynamics in the Open} |
| 66 |
> |
\title{{\sc OpenMD-2.3}: Molecular Dynamics in the Open} |
| 67 |
|
|
| 68 |
< |
\author{Joseph Michalka, James Marr, Kelsey Stocker, Madan Lamichhane, |
| 69 |
< |
Patrick Louden, \\ |
| 70 |
< |
Teng Lin, Charles F. Vardeman II, Christopher J. Fennell, Shenyu |
| 71 |
< |
Kuang, Xiuquan Sun, \\ |
| 72 |
< |
Chunlei Li, Kyle Daily, Yang Zheng, Matthew A. Meineke, and \\ |
| 73 |
< |
J. Daniel Gezelter \\ |
| 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 |
| 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{code}[caption={[An example of a complete OpenMD |
| 255 |
|
file] An example showing a complete OpenMD file.}, |
| 292 |
|
</OpenMD> |
| 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 |
| 299 |
< |
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{code}[caption={An example molecule definition in an |
| 307 |
|
include file.},label={sch:mdIncludeExample}] |
| 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 |
| 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. |
| 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} \\ |
| 560 |
– |
{\tt cutoffPolicy} & string & one of mix, max, or |
| 561 |
– |
traditional & the traditional cutoff policy is to set the cutoff |
| 562 |
– |
radius for all atoms in the system to the same value (governed by the |
| 563 |
– |
largest atom). mix and max are pair-dependent cutoff |
| 564 |
– |
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 |
| 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} |
| 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} |
| 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} |
| 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:InversionTypes}] |
| 2067 |
> |
label={sch:NonBondedInteractionTypes}] |
| 2068 |
|
begin NonBondedInteractions |
| 2069 |
|
|
| 2070 |
|
//Lennard-Jones |
| 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. |
| 2333 |
– |
|
| 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 \\ |
| 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 |
| 3763 |
< |
% keep 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 |
|
|
| 3846 |
– |
There are two additional files with the suffixes {\tt .zang0} and {\tt |
| 3847 |
– |
.zang} generated by {\sc OpenMD} during the first run of a solid |
| 3848 |
– |
thermodynamic integration. These files contain the initial ({\tt |
| 3849 |
– |
.zang0}) and final ({\tt .zang}) values of the angular displacement |
| 3850 |
– |
coordinates for each of the molecules. These are particularly useful |
| 3851 |
– |
when chaining a number of simulations (with successive values of |
| 3852 |
– |
$\lambda$) together. |
| 3853 |
– |
|
| 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 |
| 3873 |
< |
to an Einstein crystal? & default is ``false'' \\ |
| 3874 |
< |
{\tt useLiquidThermInt} & logical & perform thermodynamic integration |
| 3875 |
< |
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 \\ |
| 3883 |
– |
{\tt thermIntDistSpringConst} & & & \\ |
| 3884 |
– |
& $\mbox{kcal~mol}^{-1} \mbox{\AA}^{-2}$ |
| 3885 |
– |
& & spring constant for translations in Einstein crystal \\ |
| 3886 |
– |
{\tt thermIntThetaSpringConst} & & & \\ |
| 3887 |
– |
& $\mbox{kcal~mol}^{-1} |
| 3888 |
– |
\mbox{rad}^{-2}$ & & spring constant for deflection away from z-axis |
| 3889 |
– |
in Einstein crystal \\ |
| 3890 |
– |
{\tt thermIntOmegaSpringConst} & & & \\ |
| 3891 |
– |
& $\mbox{kcal~mol}^{-1} |
| 3892 |
– |
\mbox{rad}^{-2}$ & & spring constant for rotation around z-axis in |
| 3893 |
– |
Einstein crystal |
| 4182 |
|
\label{table:thermIntParams} |
| 4183 |
|
\end{longtable} |
| 4184 |
|
|
| 4208 |
|
|
| 4209 |
|
\begin{figure} |
| 4210 |
|
\includegraphics[width=\linewidth]{rnemdDemo} |
| 4211 |
< |
\caption{The (VSS) RNEMD approach imposes unphysical transfer of both |
| 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. |
| 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$. |
| 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 |
| 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 |
|
|
| 4664 |
– |
The specific programs {\tt xyz2md} and {\tt pdb2md} are identical |
| 4665 |
– |
to {\tt atom2md}, but they use a specific input format and do not |
| 4666 |
– |
expect the the input specifier on the command line. |
| 4667 |
– |
|
| 4668 |
– |
|
| 4967 |
|
\section{\label{section:SimpleBuilder}SimpleBuilder} |
| 4968 |
|
|
| 4969 |
|
{\tt SimpleBuilder} creates simple lattice structures. It requires an |
| 5063 |
|
\end{longtable} |
| 5064 |
|
|
| 5065 |
|
|
| 4768 |
– |
|
| 4769 |
– |
|
| 4770 |
– |
|
| 5066 |
|
\chapter{\label{section:parallelization} Parallel Simulation Implementation} |
| 5067 |
|
|
| 5068 |
|
Although processor power is continually improving, it is still |
| 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 |
|
|
| 4853 |
– |
|
| 5148 |
|
\bibliographystyle{aip} |
| 5149 |
|
\bibliography{openmdDoc} |
| 5150 |
|
|
