| 38 |
|
\newcolumntype{H}{p{0.75in}} |
| 39 |
|
\newcolumntype{I}{p{5in}} |
| 40 |
|
|
| 41 |
+ |
\newcolumntype{J}{p{1.5in}} |
| 42 |
+ |
\newcolumntype{K}{p{1.2in}} |
| 43 |
+ |
\newcolumntype{L}{p{1.5in}} |
| 44 |
+ |
\newcolumntype{M}{p{1.55in}} |
| 45 |
|
|
| 46 |
+ |
|
| 47 |
|
\title{{\sc OpenMD}: Molecular Dynamics in the Open} |
| 48 |
|
|
| 49 |
< |
\author{Shenyu Kuang, Chunlei Li, Charles F. Vardeman II, \\ |
| 50 |
< |
Teng Lin, Christopher J. Fennell, Xiuquan Sun, \\ |
| 51 |
< |
Kyle Daily, Yang Zheng, Matthew A. Meineke, and J. Daniel Gezelter\\ |
| 52 |
< |
Department of Chemistry and Biochemistry\\ |
| 53 |
< |
University of Notre Dame\\ |
| 54 |
< |
Notre Dame, Indiana 46556} |
| 49 |
> |
\author{Shenyu Kuang, Charles F. Vardeman II, \\ |
| 50 |
> |
Teng Lin, Christopher J. Fennell, Xiuquan Sun, \\ |
| 51 |
> |
Chunlei Li, Kyle Daily, Yang Zheng, Matthew A. Meineke, and \\ |
| 52 |
> |
J. Daniel Gezelter \\ |
| 53 |
> |
Department of Chemistry and Biochemistry\\ |
| 54 |
> |
University of Notre Dame\\ |
| 55 |
> |
Notre Dame, Indiana 46556} |
| 56 |
|
|
| 57 |
|
\maketitle |
| 58 |
|
|
| 503 |
|
{\tt minimizer} & string & Chooses a minimizer & Possible minimizers |
| 504 |
|
are SD and CG. Either {\tt ensemble} or {\tt minimizer} must be specified. \\ |
| 505 |
|
{\tt ensemble} & string & Sets the ensemble. & Possible ensembles are |
| 506 |
< |
NVE, NVT, NPTi, NPAT, NPTf, NPTxyz, LD and LHull. Either {\tt ensemble} |
| 506 |
> |
NVE, NVT, NPTi, NPAT, NPTf, NPTxyz, LD and LangevinHull. Either {\tt ensemble} |
| 507 |
|
or {\tt minimizer} must be specified. \\ |
| 508 |
|
{\tt dt} & fs & Sets the time step. & Selection of {\tt dt} should be |
| 509 |
|
small enough to sample the fastest motion of the simulation. ({\tt |
| 713 |
|
<MetaData> |
| 714 |
|
molecule{ |
| 715 |
|
name = "I2"; |
| 716 |
< |
atom[0]{ |
| 717 |
< |
type = "I"; |
| 718 |
< |
} |
| 713 |
< |
atom[1]{ |
| 714 |
< |
type = "I"; |
| 715 |
< |
} |
| 716 |
< |
bond{ |
| 717 |
< |
members( 0, 1); |
| 718 |
< |
} |
| 716 |
> |
atom[0]{ type = "I"; } |
| 717 |
> |
atom[1]{ type = "I"; } |
| 718 |
> |
bond{ members( 0, 1); } |
| 719 |
|
} |
| 720 |
|
molecule{ |
| 721 |
|
name = "HCl" |
| 722 |
< |
atom[0]{ |
| 723 |
< |
type = "H"; |
| 724 |
< |
} |
| 725 |
< |
atom[1]{ |
| 726 |
< |
type = "Cl"; |
| 727 |
< |
} |
| 728 |
< |
bond{ |
| 729 |
< |
members( 0, 1); |
| 730 |
< |
} |
| 722 |
> |
atom[0]{ type = "H";} |
| 723 |
> |
atom[1]{ type = "Cl";} |
| 724 |
> |
bond{ members( 0, 1); } |
| 725 |
|
} |
| 726 |
|
component{ |
| 727 |
|
type = "HCl"; |
| 1900 |
|
& (with separate barostats on each box dimension) & \\ |
| 1901 |
|
LD & Langevin Dynamics & {\tt ensemble = LD;} \\ |
| 1902 |
|
& (approximates the effects of an implicit solvent) & \\ |
| 1903 |
< |
LangevinHull & Non-periodic Langevin Dynamics & {\tt ensemble = LHull;} \\ |
| 1903 |
> |
LangevinHull & Non-periodic Langevin Dynamics & {\tt ensemble = LangevinHull;} \\ |
| 1904 |
|
& (Langevin Dynamics for molecules on convex hull;\\ |
| 1905 |
|
& Newtonian for interior molecules) & \\ |
| 1906 |
|
\end{tabular} |
| 2624 |
|
\label{table:ldParameters} |
| 2625 |
|
\end{longtable} |
| 2626 |
|
|
| 2627 |
< |
\section{Langevin Hull Dynamics (LHull)} |
| 2627 |
> |
\section{Constant Pressure without periodic boundary conditions (The LangevinHull)} |
| 2628 |
|
|
| 2629 |
< |
The Langevin Hull uses an external bath at a fixed constant pressure |
| 2629 |
> |
The Langevin Hull\cite{Vardeman2011} uses an external bath at a fixed constant pressure |
| 2630 |
|
($P$) and temperature ($T$) with an effective solvent viscosity |
| 2631 |
|
($\eta$). This bath interacts only with the objects on the exterior |
| 2632 |
|
hull of the system. Defining the hull of the atoms in a simulation is |
| 2700 |
|
depends on the geometry and surface area of facet $f$ and the |
| 2701 |
|
viscosity of the bath. The resistance tensor is related to the |
| 2702 |
|
fluctuations of the random force, $\mathbf{R}(t)$, by the |
| 2703 |
< |
fluctuation-dissipation theorem, |
| 2710 |
< |
\begin{eqnarray} |
| 2711 |
< |
\left< {\mathbf R}_f(t) \right> & = & 0 \\ |
| 2712 |
< |
\left<{\mathbf R}_f(t) {\mathbf R}_f^T(t^\prime)\right> & = & 2 k_B T\ |
| 2713 |
< |
\Xi_f(t)\delta(t-t^\prime). |
| 2714 |
< |
\label{eq:randomForce} |
| 2715 |
< |
\end{eqnarray} |
| 2703 |
> |
fluctuation-dissipation theorem (see Eq. \ref{eq:randomForce}). |
| 2704 |
|
|
| 2705 |
|
Once the resistance tensor is known for a given facet, a stochastic |
| 2706 |
|
vector that has the properties in Eq. (\ref{eq:randomForce}) can be |
| 2707 |
|
calculated efficiently by carrying out a Cholesky decomposition to |
| 2708 |
< |
obtain the square root matrix of the resistance tensor, |
| 2709 |
< |
\begin{equation} |
| 2722 |
< |
\Xi_f = {\bf S} {\bf S}^{T}, |
| 2723 |
< |
\label{eq:Cholesky} |
| 2724 |
< |
\end{equation} |
| 2725 |
< |
where ${\bf S}$ is a lower triangular matrix.\cite{Schlick2002} A |
| 2726 |
< |
vector with the statistics required for the random force can then be |
| 2727 |
< |
obtained by multiplying ${\bf S}$ onto a random 3-vector ${\bf Z}$ which |
| 2728 |
< |
has elements chosen from a Gaussian distribution, such that: |
| 2729 |
< |
\begin{equation} |
| 2730 |
< |
\langle {\bf Z}_i \rangle = 0, \hspace{1in} \langle {\bf Z}_i \cdot |
| 2731 |
< |
{\bf Z}_j \rangle = \frac{2 k_B T}{\delta t} \delta_{ij}, |
| 2732 |
< |
\end{equation} |
| 2733 |
< |
where $\delta t$ is the timestep in use during the simulation. The |
| 2734 |
< |
random force, ${\bf R}_{f} = {\bf S} {\bf Z}$, can be shown to |
| 2735 |
< |
have the correct properties required by Eq. (\ref{eq:randomForce}). |
| 2708 |
> |
obtain the square root matrix of the resistance tensor (see |
| 2709 |
> |
Eq. \ref{eq:Cholesky}). |
| 2710 |
|
|
| 2711 |
< |
Our treatment of the resistance tensor is approximate. $\Xi_f$ for a |
| 2712 |
< |
rigid triangular plate would normally be treated as a $6 \times 6$ |
| 2713 |
< |
tensor that includes translational and rotational drag as well as |
| 2714 |
< |
translational-rotational coupling. The computation of resistance |
| 2715 |
< |
tensors for rigid bodies has been detailed |
| 2711 |
> |
Our treatment of the resistance tensor for the Langevin Hull facets is |
| 2712 |
> |
approximate. $\Xi_f$ for a rigid triangular plate would normally be |
| 2713 |
> |
treated as a $6 \times 6$ tensor that includes translational and |
| 2714 |
> |
rotational drag as well as translational-rotational coupling. The |
| 2715 |
> |
computation of resistance tensors for rigid bodies has been detailed |
| 2716 |
|
elsewhere,\cite{JoseGarciadelaTorre02012000,Garcia-de-la-Torre:2001wd,GarciadelaTorreJ2002,Sun:2008fk} |
| 2717 |
|
but the standard approach involving bead approximations would be |
| 2718 |
|
prohibitively expensive if it were recomputed at each step in a |
| 2764 |
|
\item Atomic positions and velocities are propagated. |
| 2765 |
|
\end{enumerate} |
| 2766 |
|
The Delaunay triangulation and computation of the convex hull are done |
| 2767 |
< |
using calls to the qhull library.\cite{Qhull} There is a minimal |
| 2768 |
< |
penalty for computing the convex hull and resistance tensors at each |
| 2769 |
< |
step in the molecular dynamics simulation (roughly 0.02 $\times$ cost |
| 2770 |
< |
of a single force evaluation), and the convex hull is remarkably easy |
| 2771 |
< |
to parallelize on distributed memory machines. |
| 2772 |
< |
|
| 2767 |
> |
using calls to the qhull library,\cite{Qhull} and for this reason, if |
| 2768 |
> |
qhull is not detected during the build, the Langevin Hull integrator |
| 2769 |
> |
will not be available. There is a minimal penalty for computing the |
| 2770 |
> |
convex hull and resistance tensors at each step in the molecular |
| 2771 |
> |
dynamics simulation (roughly 0.02 $\times$ cost of a single force |
| 2772 |
> |
evaluation). |
| 2773 |
|
|
| 2774 |
|
\begin{longtable}[c]{GBF} |
| 2775 |
|
\caption{Meta-data Keywords: Required parameters for the Langevin Hull integrator} |
| 2784 |
|
This parameter must be specified to use Langevin Hull dynamics. \\ |
| 2785 |
|
{\tt targetPressure} & atm & Sets the target pressure of the system. |
| 2786 |
|
This parameter must be specified to use Langevin Hull dynamics. \\ |
| 2787 |
< |
{\tt usePeriodicBoundaryConditions = false} & logical & Turns off periodic boundary conditions. |
| 2787 |
> |
{\tt usePeriodicBoundaryConditions} & logical & Turns off periodic boundary conditions. |
| 2788 |
|
This parameter must be set to \tt false \\ |
| 2789 |
|
\label{table:lhullParameters} |
| 2790 |
|
\end{longtable} |
| 2938 |
|
\label{table:zconParams} |
| 2939 |
|
\end{longtable} |
| 2940 |
|
|
| 2941 |
< |
\chapter{\label{section:restraints}Restraints} |
| 2942 |
< |
Restraints are external potentials that are added to a system to keep |
| 2943 |
< |
particular molecules or collections of particles close to some |
| 2944 |
< |
reference structure. A restraint can be a collective |
| 2941 |
> |
% \chapter{\label{section:restraints}Restraints} |
| 2942 |
> |
% Restraints are external potentials that are added to a system to keep |
| 2943 |
> |
% particular molecules or collections of particles close to some |
| 2944 |
> |
% reference structure. A restraint can be a collective |
| 2945 |
|
|
| 2946 |
|
\chapter{\label{section:thermInt}Thermodynamic Integration} |
| 2947 |
|
|
| 3079 |
|
\mbox{rad}^{-2}$ & & spring constant for rotation around z-axis in |
| 3080 |
|
Einstein crystal |
| 3081 |
|
\label{table:thermIntParams} |
| 3082 |
+ |
\end{longtable} |
| 3083 |
+ |
|
| 3084 |
+ |
\chapter{\label{section:rnemd}RNEMD} |
| 3085 |
+ |
|
| 3086 |
+ |
There are many ways to compute transport properties from molecular |
| 3087 |
+ |
dynamic simulations. Equilibrium Molecular Dynamics (EMD) simulations |
| 3088 |
+ |
can be used by computing relevant time correlation functions and |
| 3089 |
+ |
assuming linear response theory holds. These approaches are generally |
| 3090 |
+ |
subject to noise and poor convergence of the relevant correlation |
| 3091 |
+ |
functions. Traditional Non-equilibrium Molecular Dynamics (NEMD) |
| 3092 |
+ |
methods impose a gradient (e.g. thermal or momentum) on a simulation. |
| 3093 |
+ |
However, the resulting flux is often difficult to |
| 3094 |
+ |
measure. Furthermore, problems arise for NEMD simulations of |
| 3095 |
+ |
heterogeneous systems, such as phase-phase boundaries or interfaces, |
| 3096 |
+ |
where the type of gradient to enforce at the boundary between |
| 3097 |
+ |
materials is unclear. |
| 3098 |
+ |
|
| 3099 |
+ |
{\it Reverse} Non-Equilibrium Molecular Dynamics (RNEMD) methods adopt a |
| 3100 |
+ |
different approach in that an unphysical {\it flux} is imposed between |
| 3101 |
+ |
different regions or ``slabs'' of the simulation box. The response of |
| 3102 |
+ |
the system is to develop a temperature or momentum {\it gradient} |
| 3103 |
+ |
between the two regions. Since the amount of the applied flux is known |
| 3104 |
+ |
exactly, and the measurement of gradient is generally less |
| 3105 |
+ |
complicated, imposed-flux methods typically take shorter simulation |
| 3106 |
+ |
times to obtain converged results for transport properties. |
| 3107 |
+ |
|
| 3108 |
+ |
%RNEMD figure |
| 3109 |
+ |
|
| 3110 |
+ |
|
| 3111 |
+ |
RNEMD methods further its advantages by utilizing momentum- and |
| 3112 |
+ |
energy-conserving approaches to apply fluxes. The original |
| 3113 |
+ |
``swapping'' approach by Muller-Plathe {\it et al.} %CITATIONS |
| 3114 |
+ |
can be seen as an imaginary elastic collision between selected |
| 3115 |
+ |
particles for each momentum exchange. This simple to implement |
| 3116 |
+ |
algorithm turned out to be quite useful in many simulations. However, |
| 3117 |
+ |
the approach inherently perturbs the ideal Maxwell-Boltzmann |
| 3118 |
+ |
distributions, which leads to undesirable side-effects when the |
| 3119 |
+ |
applied exchanged flux becomes quite large. %CITATION |
| 3120 |
+ |
This limits the range of flux available to the method, and also its |
| 3121 |
+ |
applications. |
| 3122 |
+ |
|
| 3123 |
+ |
In OpenMD, we improve the above method by introducing two alternative |
| 3124 |
+ |
approaches: |
| 3125 |
+ |
|
| 3126 |
+ |
Non-Isotropic Velocity Scaling (NIVS): %CITATION |
| 3127 |
+ |
Instead of have two individual particles involved in momentum |
| 3128 |
+ |
exchange, this algorithm applies scaling to all the particles in |
| 3129 |
+ |
particular regions: |
| 3130 |
+ |
|
| 3131 |
+ |
%NIVS equations |
| 3132 |
+ |
|
| 3133 |
+ |
Although the above matrices can be diagonal as shown, these |
| 3134 |
+ |
coefficients cannot be always the same, in order to satisfy the linear |
| 3135 |
+ |
momentum and kinetic energy conservation constraints: |
| 3136 |
+ |
|
| 3137 |
+ |
%Conservation equations |
| 3138 |
+ |
|
| 3139 |
+ |
And to apply a kinetic energy exchange between the two regions, the |
| 3140 |
+ |
following should be satisfied as well: |
| 3141 |
+ |
|
| 3142 |
+ |
%Flux equations |
| 3143 |
+ |
|
| 3144 |
+ |
Mathematically, any points in the 3-dimensional space of the solution |
| 3145 |
+ |
set would satisfy the equations. However, to avoid solving an |
| 3146 |
+ |
ill-conditioned high-order polynomial in actual practice, another |
| 3147 |
+ |
constraint, ${x_c=y_c}$, is applied, taking into consideration of its |
| 3148 |
+ |
physical relevance. Therefore, a quartic equation is solved in actual |
| 3149 |
+ |
practice to determine the sets of possible coefficients. To determine |
| 3150 |
+ |
which set is actually used to perform the scaling, two criteria are |
| 3151 |
+ |
mainly considered: 1. ${x,y,z\rightarrow 1}$ so that the perturbation |
| 3152 |
+ |
could be as gentle as possible. 2. ${K^x, K^y, K^z}$ have minimal |
| 3153 |
+ |
difference among each other, so that the anisotropy introduced by the |
| 3154 |
+ |
algorithm can be offset to some extend. One set of scaling |
| 3155 |
+ |
coefficients is chosen against these criteria, and the best one is |
| 3156 |
+ |
used to perform the scaling for that particular step. However, if no |
| 3157 |
+ |
solution found, the NIVS move is not performed in that step. |
| 3158 |
+ |
|
| 3159 |
+ |
Although the NIVS algorithm can also be applied to impose a |
| 3160 |
+ |
directional momentum flux, thermal anisotropy was observed in |
| 3161 |
+ |
relatively high flux simulations. %This is because... |
| 3162 |
+ |
However, the gentleness and ability to apply a wide range of kinetic |
| 3163 |
+ |
energy flux makes the method useful in thermal transport simulations, |
| 3164 |
+ |
particularly for complex and heterogeneous systems including |
| 3165 |
+ |
interfaces. %CITATION |
| 3166 |
+ |
|
| 3167 |
+ |
Velocity Shearing and Scaling (VSS): %CITATION |
| 3168 |
+ |
Learning from NIVS that imposing directional momentum flux by velocity |
| 3169 |
+ |
scaling could cause problem, we shift the approach to combine the move |
| 3170 |
+ |
of velocity shearing and scaling: |
| 3171 |
+ |
|
| 3172 |
+ |
%VSS equations |
| 3173 |
+ |
|
| 3174 |
+ |
It turned out that this approach results in a set of simpler-to-solve |
| 3175 |
+ |
equations for conservation and to satisfy momentum exchange: |
| 3176 |
+ |
|
| 3177 |
+ |
%conservation equations |
| 3178 |
+ |
|
| 3179 |
+ |
Furthermore, isotropic scaling is now possible, with the presence of |
| 3180 |
+ |
velocity shearing quantities. Only a set of simple quadratic equations |
| 3181 |
+ |
need to be solved, and the positive set of coefficients are chosen, in |
| 3182 |
+ |
order to reach minimal perturbations. Similar to the NIVS method, no |
| 3183 |
+ |
VSS is performed in a step given that no solution can be found. |
| 3184 |
+ |
|
| 3185 |
+ |
The VSS approach turned out to be versatile in both thermal and |
| 3186 |
+ |
directional momentum transport simulations. It is found that the |
| 3187 |
+ |
perturbation is minimal and undesired side-effects like thermal |
| 3188 |
+ |
anisotropy can be avoided. Another nice feature of VSS is its ability |
| 3189 |
+ |
to combine a thermal and a directional momentum flux. This feature has |
| 3190 |
+ |
been utilized to map out the shear viscosity of SPC/E water in a wide |
| 3191 |
+ |
range of temperature (90~K) just with one single simulation. Possible |
| 3192 |
+ |
applications may also include the studies of thermal-momentum coupled |
| 3193 |
+ |
transport phenomena. VSS also allows the directional momentum flux to |
| 3194 |
+ |
have quite arbitary directions, which could benefit researches of |
| 3195 |
+ |
anisotropic systems. |
| 3196 |
+ |
|
| 3197 |
+ |
Table \ref{table:rnemd} summarizes the parameters used in RNEMD |
| 3198 |
+ |
simulations. |
| 3199 |
+ |
|
| 3200 |
+ |
\begin{longtable}[c]{JKLM} |
| 3201 |
+ |
\caption{The following keywords must be enclosed inside a {\tt RNEMD\{\}} block} |
| 3202 |
+ |
\\ |
| 3203 |
+ |
{\bf keyword} & {\bf units} & {\bf use} & {\bf remarks} \\ \hline |
| 3204 |
+ |
\endhead |
| 3205 |
+ |
\hline |
| 3206 |
+ |
\endfoot |
| 3207 |
+ |
{\tt useRNEMD} & logical & perform RNEMD? & default is ``false'' \\ |
| 3208 |
+ |
{\tt objectSelection} & string & see section \ref{section:syntax} |
| 3209 |
+ |
for selection syntax & default is ``select all'' \\ |
| 3210 |
+ |
{\tt method} & string & exchange method & one of the following: |
| 3211 |
+ |
{\tt Swap, NIVS,} or {\tt VSS} (default is {\tt VSS}) \\ |
| 3212 |
+ |
{\tt fluxType} & string & what is being exchanged between slabs? & one |
| 3213 |
+ |
of the following: {\tt KE, Px, Py, Pz, Pvector, KE+Px, KE+Py, KE+Pvector} \\ |
| 3214 |
+ |
{\tt kineticFlux} & kcal mol$^{-1}$ \AA$^{-2}$ fs$^{-1}$ & specify the kinetic energy flux & \\ |
| 3215 |
+ |
{\tt momentumFlux} & amu \AA$^{-1}$ fs$^{-2}$ & specify the momentum flux & \\ |
| 3216 |
+ |
{\tt momentumFluxVector} & amu \AA$^{-1}$ fs$^{-2}$ & specify the momentum flux when |
| 3217 |
+ |
{\tt Pvector} is part of the exchange & Vector3d input\\ |
| 3218 |
+ |
{\tt exchangeTime} & fs & how often to perform the exchange & default is 100 fs\\ |
| 3219 |
+ |
|
| 3220 |
+ |
{\tt slabWidth} & $\mbox{\AA}$ & width of the two exchange slabs & default is $\mathsf{H}_{zz} / 10.0$ \\ |
| 3221 |
+ |
{\tt slabAcenter} & $\mbox{\AA}$ & center of the end slab & default is 0 \\ |
| 3222 |
+ |
{\tt slabBcenter} & $\mbox{\AA}$ & center of the middle slab & default is $\mathsf{H}_{zz} / 2$ \\ |
| 3223 |
+ |
{\tt outputFileName} & string & file name for output histograms & default is the same prefix as the |
| 3224 |
+ |
.md file, but with the {\tt .rnemd} extension \\ |
| 3225 |
+ |
{\tt outputBins} & int & number of $z$-bins in the output histogram & |
| 3226 |
+ |
default is 20 \\ |
| 3227 |
+ |
{\tt outputFields} & string & columns to print in the {\tt .rnemd} |
| 3228 |
+ |
file where each column is separated by a pipe ($\mid$) symbol. & Allowed column names are: {\sc z, temperature, velocity, density}} \\ |
| 3229 |
+ |
\label{table:rnemd} |
| 3230 |
|
\end{longtable} |
| 3231 |
|
|
| 3232 |
|
|
