| 17 |
|
\textwidth 6.5in |
| 18 |
|
\brokenpenalty=10000 |
| 19 |
|
\renewcommand{\baselinestretch}{1.2} |
| 20 |
+ |
\usepackage[square, comma, sort&compress]{natbib} |
| 21 |
+ |
\bibpunct{[}{]}{,}{n}{}{;} |
| 22 |
|
|
| 23 |
+ |
|
| 24 |
|
%\renewcommand\citemid{\ } % no comma in optional reference note |
| 25 |
|
\lstset{language=C,frame=TB,basicstyle=\tiny,basicstyle=\ttfamily, % |
| 26 |
|
xleftmargin=0.25in, xrightmargin=0.25in,captionpos=b, % |
| 47 |
|
\newcolumntype{M}{p{1.55in}} |
| 48 |
|
|
| 49 |
|
|
| 50 |
< |
\title{{\sc OpenMD}: Molecular Dynamics in the Open} |
| 50 |
> |
\title{{\sc OpenMD-2}: Molecular Dynamics in the Open} |
| 51 |
|
|
| 52 |
< |
\author{Shenyu Kuang, Charles F. Vardeman II, \\ |
| 53 |
< |
Teng Lin, Christopher J. Fennell, Xiuquan Sun, \\ |
| 52 |
> |
\author{Shenyu Kuang, Joseph Michalka, Kelsey Stocker, James Marr, \\ |
| 53 |
> |
Teng Lin, Charles F. Vardeman II, Christopher J. Fennell, Xiuquan Sun, \\ |
| 54 |
|
Chunlei Li, Kyle Daily, Yang Zheng, Matthew A. Meineke, and \\ |
| 55 |
|
J. Daniel Gezelter \\ |
| 56 |
|
Department of Chemistry and Biochemistry\\ |
| 499 |
|
\endhead |
| 500 |
|
\hline |
| 501 |
|
\endfoot |
| 502 |
< |
{\tt forceField} & string & Sets the force field. & Possible force |
| 503 |
< |
fields are DUFF, WATER, LJ, EAM, SC, and CLAY. \\ |
| 502 |
> |
{\tt forceField} & string & Sets the base name for the force field file & |
| 503 |
> |
OpenMD appends a {\tt .frc} to the end of this to look for a force |
| 504 |
> |
field file.\\ |
| 505 |
|
{\tt component} & & Defines the molecular components of the system & |
| 506 |
|
Every {\tt $<$MetaData$>$} block must have a component statement. \\ |
| 507 |
|
{\tt minimizer} & string & Chooses a minimizer & Possible minimizers |
| 2943 |
|
\end{longtable} |
| 2944 |
|
|
| 2945 |
|
% \chapter{\label{section:restraints}Restraints} |
| 2946 |
< |
% Restraints are external potentials that are added to a system to keep |
| 2947 |
< |
% particular molecules or collections of particles close to some |
| 2946 |
> |
% Restraints are external potentials that are added to a system to |
| 2947 |
> |
% keep particular molecules or collections of particles close to some |
| 2948 |
|
% reference structure. A restraint can be a collective |
| 2949 |
|
|
| 2950 |
|
\chapter{\label{section:thermInt}Thermodynamic Integration} |
| 3085 |
|
\label{table:thermIntParams} |
| 3086 |
|
\end{longtable} |
| 3087 |
|
|
| 3088 |
< |
\chapter{\label{section:rnemd}RNEMD} |
| 3088 |
> |
\chapter{\label{section:rnemd}Reverse Non-Equilibrium Molecular Dynamics (RNEMD)} |
| 3089 |
|
|
| 3090 |
|
There are many ways to compute transport properties from molecular |
| 3091 |
< |
dynamic simulations. Equilibrium Molecular Dynamics (EMD) simulations |
| 3092 |
< |
can be used by computing relevant time correlation functions and |
| 3093 |
< |
assuming linear response theory holds. These approaches are generally |
| 3094 |
< |
subject to noise and poor convergence of the relevant correlation |
| 3095 |
< |
functions. Traditional Non-equilibrium Molecular Dynamics (NEMD) |
| 3096 |
< |
methods impose a gradient (e.g. thermal or momentum) on a simulation. |
| 3097 |
< |
However, the resulting flux is often difficult to |
| 3091 |
> |
dynamics simulations. Equilibrium Molecular Dynamics (EMD) |
| 3092 |
> |
simulations can be used by computing relevant time correlation |
| 3093 |
> |
functions and assuming linear response theory holds. For some transport properties (notably thermal conductivity), EMD approaches |
| 3094 |
> |
are subject to noise and poor convergence of the relevant |
| 3095 |
> |
correlation functions. Traditional Non-equilibrium Molecular Dynamics |
| 3096 |
> |
(NEMD) methods impose a gradient (e.g. thermal or momentum) on a |
| 3097 |
> |
simulation. However, the resulting flux is often difficult to |
| 3098 |
|
measure. Furthermore, problems arise for NEMD simulations of |
| 3099 |
|
heterogeneous systems, such as phase-phase boundaries or interfaces, |
| 3100 |
|
where the type of gradient to enforce at the boundary between |
| 3101 |
|
materials is unclear. |
| 3102 |
|
|
| 3103 |
< |
{\it Reverse} Non-Equilibrium Molecular Dynamics (RNEMD) methods adopt a |
| 3104 |
< |
different approach in that an unphysical {\it flux} is imposed between |
| 3105 |
< |
different regions or ``slabs'' of the simulation box. The response of |
| 3106 |
< |
the system is to develop a temperature or momentum {\it gradient} |
| 3107 |
< |
between the two regions. Since the amount of the applied flux is known |
| 3108 |
< |
exactly, and the measurement of gradient is generally less |
| 3109 |
< |
complicated, imposed-flux methods typically take shorter simulation |
| 3110 |
< |
times to obtain converged results for transport properties. |
| 3103 |
> |
{\it Reverse} Non-Equilibrium Molecular Dynamics (RNEMD) methods adopt |
| 3104 |
> |
a different approach in that an unphysical {\it flux} is imposed |
| 3105 |
> |
between different regions or ``slabs'' of the simulation box. The |
| 3106 |
> |
response of the system is to develop a temperature or momentum {\it |
| 3107 |
> |
gradient} between the two regions. Since the amount of the applied |
| 3108 |
> |
flux is known exactly, and the measurement of gradient is generally |
| 3109 |
> |
less complicated, imposed-flux methods typically take shorter |
| 3110 |
> |
simulation times to obtain converged results for transport properties. |
| 3111 |
|
|
| 3112 |
< |
%RNEMD figure |
| 3112 |
> |
\begin{figure} |
| 3113 |
> |
\includegraphics[width=\linewidth]{rnemdDemo} |
| 3114 |
> |
\caption{The (VSS) RNEMD approach imposes unphysical transfer of both |
| 3115 |
> |
linear momentum and kinetic energy between a ``hot'' slab and a |
| 3116 |
> |
``cold'' slab in the simulation box. The system responds to this |
| 3117 |
> |
imposed flux by generating both momentum and temperature gradients. |
| 3118 |
> |
The slope of the gradients can then be used to compute transport |
| 3119 |
> |
properties (e.g. shear viscosity and thermal conductivity).} |
| 3120 |
> |
\label{rnemdDemo} |
| 3121 |
> |
\end{figure} |
| 3122 |
> |
|
| 3123 |
> |
\section{\label{section:algo}Three algorithms for carrying out RNEMD simulations} |
| 3124 |
> |
\subsection{\label{subsection:swapping}The swapping algorithm} |
| 3125 |
> |
The original ``swapping'' approaches by M\"{u}ller-Plathe {\it et |
| 3126 |
> |
al.}\cite{ISI:000080382700030,MullerPlathe:1997xw} can be understood |
| 3127 |
> |
as a sequence of imaginary elastic collisions between particles in |
| 3128 |
> |
opposite slabs. In each collision, the entire momentum vectors of |
| 3129 |
> |
both particles may be exchanged to generate a thermal |
| 3130 |
> |
flux. Alternatively, a single component of the momentum vectors may be |
| 3131 |
> |
exchanged to generate a shear flux. This algorithm turns out to be |
| 3132 |
> |
quite useful in many simulations. However, the M\"{u}ller-Plathe |
| 3133 |
> |
swapping approach perturbs the system away from ideal |
| 3134 |
> |
Maxwell-Boltzmann distributions, and this may leads to undesirable |
| 3135 |
> |
side-effects when the applied flux becomes large.\cite{Maginn:2010} |
| 3136 |
> |
This limits the applicability of the swapping algorithm, so in OpenMD, |
| 3137 |
> |
we have implemented two additional algorithms for RNEMD in addition to the |
| 3138 |
> |
original swapping approach. |
| 3139 |
|
|
| 3140 |
+ |
\subsection{\label{subsection:nivs}Non-Isotropic Velocity Scaling (NIVS)} |
| 3141 |
+ |
Instead of having momentum exchange between {\it individual particles} |
| 3142 |
+ |
in each slab, the NIVS algorithm applies velocity scaling to all of |
| 3143 |
+ |
the selected particles in both slabs.\cite{kuang:164101} A combination of linear |
| 3144 |
+ |
momentum, kinetic energy, and flux constraint equations governs the |
| 3145 |
+ |
amount of velocity scaling performed at each step. Interested readers |
| 3146 |
+ |
should consult ref. \citealp{kuang:164101} for further details on the |
| 3147 |
+ |
methodology. |
| 3148 |
|
|
| 3149 |
< |
RNEMD methods further its advantages by utilizing momentum- and |
| 3150 |
< |
energy-conserving approaches to apply fluxes. The original |
| 3151 |
< |
``swapping'' approach by Muller-Plathe {\it et al.} %CITATIONS |
| 3152 |
< |
can be seen as an imaginary elastic collision between selected |
| 3153 |
< |
particles for each momentum exchange. This simple to implement |
| 3154 |
< |
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. |
| 3149 |
> |
NIVS has been shown to be very effective at producing thermal |
| 3150 |
> |
gradients and for computing thermal conductivities, particularly for |
| 3151 |
> |
heterogeneous interfaces. Although the NIVS algorithm can also be |
| 3152 |
> |
applied to impose a directional momentum flux, thermal anisotropy was |
| 3153 |
> |
observed in relatively high flux simulations, and the method is not |
| 3154 |
> |
suitable for imposing a shear flux or for computing shear viscosities. |
| 3155 |
|
|
| 3156 |
< |
In OpenMD, we improve the above method by introducing two alternative |
| 3157 |
< |
approaches: |
| 3156 |
> |
\subsection{\label{subsection:vss}Velocity Shearing and Scaling (VSS)} |
| 3157 |
> |
The third RNEMD algorithm implemented in OpenMD utilizes a series of |
| 3158 |
> |
simultaneous velocity shearing and scaling exchanges between the two |
| 3159 |
> |
slabs.\cite{2012MolPh.110..691K} This method results in a set of simpler equations to satisfy |
| 3160 |
> |
the conservation constraints while creating a desired flux between the |
| 3161 |
> |
two slabs. |
| 3162 |
|
|
| 3163 |
< |
Non-Isotropic Velocity Scaling (NIVS): %CITATION |
| 3164 |
< |
Instead of have two individual particles involved in momentum |
| 3165 |
< |
exchange, this algorithm applies scaling to all the particles in |
| 3166 |
< |
particular regions: |
| 3163 |
> |
The VSS approach is versatile in that it may be used to implement both |
| 3164 |
> |
thermal and shear transport either separately or simultaneously. |
| 3165 |
> |
Perturbations of velocities away from the ideal Maxwell-Boltzmann |
| 3166 |
> |
distributions are minimal, and thermal anisotropy is kept to a |
| 3167 |
> |
minimum. This ability to generate simultaneous thermal and shear |
| 3168 |
> |
fluxes has been utilized to map out the shear viscosity of SPC/E water |
| 3169 |
> |
over a wide range of temperatures (90~K) just with a single simulation. |
| 3170 |
> |
VSS-RNEMD also allows the directional momentum flux to have |
| 3171 |
> |
arbitary directions, which could aid in the study of anisotropic solid |
| 3172 |
> |
surfaces in contact with liquid environments. |
| 3173 |
|
|
| 3174 |
< |
%NIVS equations |
| 3175 |
< |
|
| 3176 |
< |
Although the above matrices can be diagonal as shown, these |
| 3177 |
< |
coefficients cannot be always the same, in order to satisfy the linear |
| 3178 |
< |
momentum and kinetic energy conservation constraints: |
| 3179 |
< |
|
| 3180 |
< |
%Conservation equations |
| 3181 |
< |
|
| 3182 |
< |
And to apply a kinetic energy exchange between the two regions, the |
| 3183 |
< |
following should be satisfied as well: |
| 3174 |
> |
\section{\label{section:usingRNEMD}Using OpenMD to perform a RNEMD simulation} |
| 3175 |
> |
\subsection{\label{section:rnemdParams} What the user needs to specify} |
| 3176 |
> |
To carry out a RNEMD simulation, |
| 3177 |
> |
a user must specify a number of parameters in the MetaData (.md) file. |
| 3178 |
> |
Because the RNEMD methods have a large number of parameters, these |
| 3179 |
> |
must be enclosed in a {\it separate} {\tt RNEMD\{...\}} block. The most important |
| 3180 |
> |
parameters to specify are the {\tt useRNEMD}, {\tt fluxType} and flux |
| 3181 |
> |
parameters. Most other parameters (summarized in table |
| 3182 |
> |
\ref{table:rnemd}) have reasonable default values. {\tt fluxType} |
| 3183 |
> |
sets up the kind of exchange that will be carried out between the two |
| 3184 |
> |
slabs (either Kinetic Energy ({\tt KE}) or momentum ({\tt Px, Py, Pz, |
| 3185 |
> |
Pvector}), or some combination of these). The flux is specified |
| 3186 |
> |
with the use of three possible parameters: {\tt kineticFlux} for |
| 3187 |
> |
kinetic energy exchange, as well as {\tt momentumFlux} or {\tt |
| 3188 |
> |
momentumFluxVector} for simulations with directional exchange. |
| 3189 |
|
|
| 3190 |
< |
%Flux equations |
| 3190 |
> |
\subsection{\label{section:rnemdResults} Processing the results} |
| 3191 |
> |
OpenMD will generate a {\tt .rnemd} |
| 3192 |
> |
file with the same prefix as the original {\tt .md} file. This file |
| 3193 |
> |
contains a running average of properties of interest computed within a |
| 3194 |
> |
set of bins that divide the simulation cell along the $z$-axis. The |
| 3195 |
> |
first column of the {\tt .rnemd} file is the $z$ coordinate of the |
| 3196 |
> |
center of each bin, while following columns may contain the average |
| 3197 |
> |
temperature, velocity, or density within each bin. The output format |
| 3198 |
> |
in the {\tt .rnemd} file can be altered with the {\tt outputFields}, |
| 3199 |
> |
{\tt outputBins}, and {\tt outputFileName} parameters. A report at |
| 3200 |
> |
the top of the {\tt .rnemd} file contains the current exchange totals |
| 3201 |
> |
as well as the average flux applied during the simulation. Using the |
| 3202 |
> |
slope of the temperature or velocity gradient obtaine from the {\tt |
| 3203 |
> |
.rnemd} file along with the applied flux, the user can very simply |
| 3204 |
> |
arrive at estimates of thermal conductivities ($\lambda$), |
| 3205 |
> |
\begin{equation} |
| 3206 |
> |
J_z = -\lambda \frac{\partial T}{\partial z}, |
| 3207 |
> |
\end{equation} |
| 3208 |
> |
and shear viscosities ($\eta$), |
| 3209 |
> |
\begin{equation} |
| 3210 |
> |
j_z(p_x) = -\eta \frac{\partial \langle v_x \rangle}{\partial z}. |
| 3211 |
> |
\end{equation} |
| 3212 |
> |
Here, the quantities on the left hand side are the actual flux values |
| 3213 |
> |
(in the header of the {\tt .rnemd} file), while the slopes are |
| 3214 |
> |
obtained from linear fits to the gradients observed in the {\tt |
| 3215 |
> |
.rnemd} file. |
| 3216 |
|
|
| 3217 |
< |
Mathematically, any points in the 3-dimensional space of the solution |
| 3218 |
< |
set would satisfy the equations. However, to avoid solving an |
| 3219 |
< |
ill-conditioned high-order polynomial in actual practice, another |
| 3220 |
< |
constraint, ${x_c=y_c}$, is applied, taking into consideration of its |
| 3221 |
< |
physical relevance. Therefore, a quartic equation is solved in actual |
| 3222 |
< |
practice to determine the sets of possible coefficients. To determine |
| 3223 |
< |
which set is actually used to perform the scaling, two criteria are |
| 3224 |
< |
mainly considered: 1. ${x,y,z\rightarrow 1}$ so that the perturbation |
| 3225 |
< |
could be as gentle as possible. 2. ${K^x, K^y, K^z}$ have minimal |
| 3226 |
< |
difference among each other, so that the anisotropy introduced by the |
| 3227 |
< |
algorithm can be offset to some extend. One set of scaling |
| 3228 |
< |
coefficients is chosen against these criteria, and the best one is |
| 3229 |
< |
used to perform the scaling for that particular step. However, if no |
| 3230 |
< |
solution found, the NIVS move is not performed in that step. |
| 3231 |
< |
|
| 3232 |
< |
Although the NIVS algorithm can also be applied to impose a |
| 3233 |
< |
directional momentum flux, thermal anisotropy was observed in |
| 3234 |
< |
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. |
| 3217 |
> |
More complicated simulations (including interfaces) require a bit more |
| 3218 |
> |
care. Here the second derivative may be required to compute the |
| 3219 |
> |
interfacial thermal conductance, |
| 3220 |
> |
\begin{align} |
| 3221 |
> |
G^\prime &= \left(\nabla\lambda \cdot \mathbf{\hat{n}}\right)_{z_0} \\ |
| 3222 |
> |
&= \frac{\partial}{\partial z}\left(-\frac{J_z}{ |
| 3223 |
> |
\left(\frac{\partial T}{\partial z}\right)}\right)_{z_0} \\ |
| 3224 |
> |
&= J_z\left(\frac{\partial^2 T}{\partial z^2}\right)_{z_0} \Big/ |
| 3225 |
> |
\left(\frac{\partial T}{\partial z}\right)_{z_0}^2. |
| 3226 |
> |
\label{derivativeG} |
| 3227 |
> |
\end{align} |
| 3228 |
> |
where $z_0$ is the location of the interface between two materials and |
| 3229 |
> |
$\mathbf{\hat{n}}$ is a unit vector normal to the interface. We |
| 3230 |
> |
suggest that users interested in interfacial conductance consult |
| 3231 |
> |
reference \citealp{kuang:AuThl} for other approaches to computing $G$. |
| 3232 |
> |
Users interested in {\it friction coefficients} at heterogeneous |
| 3233 |
> |
interfaces may also find reference \citealp{2012MolPh.110..691K} |
| 3234 |
> |
useful. |
| 3235 |
|
|
| 3236 |
< |
Table \ref{table:rnemd} summarizes the parameters used in RNEMD |
| 3198 |
< |
simulations. |
| 3236 |
> |
\newpage |
| 3237 |
|
|
| 3238 |
|
\begin{longtable}[c]{JKLM} |
| 3239 |
< |
\caption{The following keywords must be enclosed inside a {\tt RNEMD\{\}} block} |
| 3240 |
< |
\\ |
| 3239 |
> |
\caption{Meta-data Keywords: Parameters for RNEMD simulations}\\ |
| 3240 |
> |
\multicolumn{4}{c}{The following keywords must be enclosed inside a {\tt RNEMD\{...\}} block.} |
| 3241 |
> |
\\ \hline |
| 3242 |
|
{\bf keyword} & {\bf units} & {\bf use} & {\bf remarks} \\ \hline |
| 3243 |
|
\endhead |
| 3244 |
|
\hline |
| 3264 |
|
{\tt outputBins} & int & number of $z$-bins in the output histogram & |
| 3265 |
|
default is 20 \\ |
| 3266 |
|
{\tt outputFields} & string & columns to print in the {\tt .rnemd} |
| 3267 |
< |
file where each column is separated by a pipe ($\mid$) symbol. & Allowed column names are: {\sc z, temperature, velocity, density}} \\ |
| 3267 |
> |
file where each column is separated by a pipe ($\mid$) symbol. & Allowed column names are: {\sc z, temperature, velocity, density} \\ |
| 3268 |
|
\label{table:rnemd} |
| 3269 |
|
\end{longtable} |
| 3270 |
|
|
| 3232 |
– |
|
| 3271 |
|
\chapter{\label{section:minimizer}Energy Minimization} |
| 3272 |
|
|
| 3273 |
< |
As one of the basic procedures of molecular modeling, energy |
| 3236 |
< |
minimization is used to identify local configurations that are stable |
| 3273 |
> |
Energy minimization is used to identify local configurations that are stable |
| 3274 |
|
points on the potential energy surface. There is a vast literature on |
| 3275 |
|
energy minimization algorithms have been developed to search for the |
| 3276 |
|
global energy minimum as well as to find local structures which are |
| 3940 |
|
DMR-0079647. |
| 3941 |
|
|
| 3942 |
|
|
| 3943 |
< |
\bibliographystyle{jcc} |
| 3943 |
> |
\bibliographystyle{aip} |
| 3944 |
|
\bibliography{openmdDoc} |
| 3945 |
|
|
| 3946 |
|
\end{document} |
