| 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, % |
| 41 |
|
\newcolumntype{H}{p{0.75in}} |
| 42 |
|
\newcolumntype{I}{p{5in}} |
| 43 |
|
|
| 44 |
+ |
\newcolumntype{J}{p{1.5in}} |
| 45 |
+ |
\newcolumntype{K}{p{1.2in}} |
| 46 |
+ |
\newcolumntype{L}{p{1.5in}} |
| 47 |
+ |
\newcolumntype{M}{p{1.55in}} |
| 48 |
|
|
| 42 |
– |
\title{{\sc OpenMD}: Molecular Dynamics in the Open} |
| 49 |
|
|
| 50 |
< |
\author{Shenyu Kuang, Chunlei Li, Charles F. Vardeman II, \\ |
| 45 |
< |
Teng Lin, Christopher J. Fennell, Xiuquan Sun, \\ |
| 46 |
< |
Kyle Daily, Yang Zheng, Matthew A. Meineke, and J. Daniel Gezelter\\ |
| 47 |
< |
Department of Chemistry and Biochemistry\\ |
| 48 |
< |
University of Notre Dame\\ |
| 49 |
< |
Notre Dame, Indiana 46556} |
| 50 |
> |
\title{{\sc OpenMD-2}: Molecular Dynamics in the Open} |
| 51 |
|
|
| 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\\ |
| 57 |
+ |
University of Notre Dame\\ |
| 58 |
+ |
Notre Dame, Indiana 46556} |
| 59 |
+ |
|
| 60 |
|
\maketitle |
| 61 |
|
|
| 62 |
|
\section*{Preface} |
| 73 |
|
that is easy to learn. |
| 74 |
|
|
| 75 |
|
\tableofcontents |
| 76 |
< |
%\listoffigures |
| 77 |
< |
%\listoftables |
| 76 |
> |
\listoffigures |
| 77 |
> |
\listoftables |
| 78 |
|
|
| 79 |
|
\mainmatter |
| 80 |
|
|
| 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 |
| 508 |
|
are SD and CG. Either {\tt ensemble} or {\tt minimizer} must be specified. \\ |
| 509 |
|
{\tt ensemble} & string & Sets the ensemble. & Possible ensembles are |
| 510 |
< |
NVE, NVT, NPTi, NPAT, NPTf, NPTxyz, and LD. Either {\tt ensemble} |
| 510 |
> |
NVE, NVT, NPTi, NPAT, NPTf, NPTxyz, LD and LangevinHull. Either {\tt ensemble} |
| 511 |
|
or {\tt minimizer} must be specified. \\ |
| 512 |
|
{\tt dt} & fs & Sets the time step. & Selection of {\tt dt} should be |
| 513 |
|
small enough to sample the fastest motion of the simulation. ({\tt |
| 717 |
|
<MetaData> |
| 718 |
|
molecule{ |
| 719 |
|
name = "I2"; |
| 720 |
< |
atom[0]{ |
| 721 |
< |
type = "I"; |
| 722 |
< |
} |
| 713 |
< |
atom[1]{ |
| 714 |
< |
type = "I"; |
| 715 |
< |
} |
| 716 |
< |
bond{ |
| 717 |
< |
members( 0, 1); |
| 718 |
< |
} |
| 720 |
> |
atom[0]{ type = "I"; } |
| 721 |
> |
atom[1]{ type = "I"; } |
| 722 |
> |
bond{ members( 0, 1); } |
| 723 |
|
} |
| 724 |
|
molecule{ |
| 725 |
|
name = "HCl" |
| 726 |
< |
atom[0]{ |
| 727 |
< |
type = "H"; |
| 728 |
< |
} |
| 725 |
< |
atom[1]{ |
| 726 |
< |
type = "Cl"; |
| 727 |
< |
} |
| 728 |
< |
bond{ |
| 729 |
< |
members( 0, 1); |
| 730 |
< |
} |
| 726 |
> |
atom[0]{ type = "H";} |
| 727 |
> |
atom[1]{ type = "Cl";} |
| 728 |
> |
bond{ members( 0, 1); } |
| 729 |
|
} |
| 730 |
|
component{ |
| 731 |
|
type = "HCl"; |
| 1904 |
|
& (with separate barostats on each box dimension) & \\ |
| 1905 |
|
LD & Langevin Dynamics & {\tt ensemble = LD;} \\ |
| 1906 |
|
& (approximates the effects of an implicit solvent) & \\ |
| 1907 |
+ |
LangevinHull & Non-periodic Langevin Dynamics & {\tt ensemble = LangevinHull;} \\ |
| 1908 |
+ |
& (Langevin Dynamics for molecules on convex hull;\\ |
| 1909 |
+ |
& Newtonian for interior molecules) & \\ |
| 1910 |
|
\end{tabular} |
| 1911 |
|
\end{center} |
| 1912 |
|
|
| 2392 |
|
in the body-fixed frame) and ${\bf V}$ is a generalized velocity, |
| 2393 |
|
${\bf V} = |
| 2394 |
|
\left\{{\bf v},{\bf \omega}\right\}$. The right side of |
| 2395 |
< |
Eq.~\ref{LDGeneralizedForm} consists of three generalized forces: a |
| 2395 |
> |
Eq. \ref{LDGeneralizedForm} consists of three generalized forces: a |
| 2396 |
|
system force (${\bf F}_{s}$), a frictional or dissipative force (${\bf |
| 2397 |
|
F}_{f}$) and a stochastic force (${\bf F}_{r}$). While the evolution |
| 2398 |
|
of the system in Newtonian mechanics is typically done in the lab |
| 2614 |
|
\endhead |
| 2615 |
|
\hline |
| 2616 |
|
\endfoot |
| 2617 |
< |
{\tt viscosity} & centipoise & Sets the value of viscosity of the implicit |
| 2617 |
> |
{\tt viscosity} & poise & Sets the value of viscosity of the implicit |
| 2618 |
|
solvent \\ |
| 2619 |
|
{\tt targetTemp} & K & Sets the target temperature of the system. |
| 2620 |
|
This parameter must be specified to use Langevin dynamics. \\ |
| 2621 |
|
{\tt HydroPropFile} & string & Specifies the name of the resistance |
| 2622 |
|
tensor (usually a {\tt .diff} file) which is precalculated by {\tt |
| 2623 |
< |
Hydro}. This keyworkd is not necessary if the simulation contains only |
| 2623 |
> |
Hydro}. This keyword is not necessary if the simulation contains only |
| 2624 |
|
simple bodies (spheres and ellipsoids). \\ |
| 2625 |
|
{\tt beadSize} & $\mbox{\AA}$ & Sets the diameter of the beads to use |
| 2626 |
|
when the {\tt RoughShell} model is used to approximate the resistance |
| 2627 |
|
tensor. |
| 2628 |
|
\label{table:ldParameters} |
| 2629 |
+ |
\end{longtable} |
| 2630 |
+ |
|
| 2631 |
+ |
\section{Constant Pressure without periodic boundary conditions (The LangevinHull)} |
| 2632 |
+ |
|
| 2633 |
+ |
The Langevin Hull\cite{Vardeman2011} uses an external bath at a fixed constant pressure |
| 2634 |
+ |
($P$) and temperature ($T$) with an effective solvent viscosity |
| 2635 |
+ |
($\eta$). This bath interacts only with the objects on the exterior |
| 2636 |
+ |
hull of the system. Defining the hull of the atoms in a simulation is |
| 2637 |
+ |
done in a manner similar to the approach of Kohanoff, Caro and |
| 2638 |
+ |
Finnis.\cite{Kohanoff:2005qm} That is, any instantaneous configuration |
| 2639 |
+ |
of the atoms in the system is considered as a point cloud in three |
| 2640 |
+ |
dimensional space. Delaunay triangulation is used to find all facets |
| 2641 |
+ |
between coplanar |
| 2642 |
+ |
neighbors.\cite{delaunay,springerlink:10.1007/BF00977785} In highly |
| 2643 |
+ |
symmetric point clouds, facets can contain many atoms, but in all but |
| 2644 |
+ |
the most symmetric of cases, the facets are simple triangles in |
| 2645 |
+ |
3-space which contain exactly three atoms. |
| 2646 |
+ |
|
| 2647 |
+ |
The convex hull is the set of facets that have {\it no concave |
| 2648 |
+ |
corners} at an atomic site.\cite{Barber96,EDELSBRUNNER:1994oq} This |
| 2649 |
+ |
eliminates all facets on the interior of the point cloud, leaving only |
| 2650 |
+ |
those exposed to the bath. Sites on the convex hull are dynamic; as |
| 2651 |
+ |
molecules re-enter the cluster, all interactions between atoms on that |
| 2652 |
+ |
molecule and the external bath are removed. Since the edge is |
| 2653 |
+ |
determined dynamically as the simulation progresses, no {\it a priori} |
| 2654 |
+ |
geometry is defined. The pressure and temperature bath interacts only |
| 2655 |
+ |
with the atoms on the edge and not with atoms interior to the |
| 2656 |
+ |
simulation. |
| 2657 |
+ |
|
| 2658 |
+ |
Atomic sites in the interior of the simulation move under standard |
| 2659 |
+ |
Newtonian dynamics, |
| 2660 |
+ |
\begin{equation} |
| 2661 |
+ |
m_i \dot{\mathbf v}_i(t)=-{\mathbf \nabla}_i U, |
| 2662 |
+ |
\label{eq:Newton} |
| 2663 |
+ |
\end{equation} |
| 2664 |
+ |
where $m_i$ is the mass of site $i$, ${\mathbf v}_i(t)$ is the |
| 2665 |
+ |
instantaneous velocity of site $i$ at time $t$, and $U$ is the total |
| 2666 |
+ |
potential energy. For atoms on the exterior of the cluster |
| 2667 |
+ |
(i.e. those that occupy one of the vertices of the convex hull), the |
| 2668 |
+ |
equation of motion is modified with an external force, ${\mathbf |
| 2669 |
+ |
F}_i^{\mathrm ext}$: |
| 2670 |
+ |
\begin{equation} |
| 2671 |
+ |
m_i \dot{\mathbf v}_i(t)=-{\mathbf \nabla}_i U + {\mathbf F}_i^{\mathrm ext}. |
| 2672 |
+ |
\end{equation} |
| 2673 |
+ |
|
| 2674 |
+ |
The external bath interacts indirectly with the atomic sites through |
| 2675 |
+ |
the intermediary of the hull facets. Since each vertex (or atom) |
| 2676 |
+ |
provides one corner of a triangular facet, the force on the facets are |
| 2677 |
+ |
divided equally to each vertex. However, each vertex can participate |
| 2678 |
+ |
in multiple facets, so the resultant force is a sum over all facets |
| 2679 |
+ |
$f$ containing vertex $i$: |
| 2680 |
+ |
\begin{equation} |
| 2681 |
+ |
{\mathbf F}_{i}^{\mathrm ext} = \sum_{\begin{array}{c}\mathrm{facets\ |
| 2682 |
+ |
} f \\ \mathrm{containing\ } i\end{array}} \frac{1}{3}\ {\mathbf |
| 2683 |
+ |
F}_f^{\mathrm ext} |
| 2684 |
+ |
\end{equation} |
| 2685 |
+ |
|
| 2686 |
+ |
The external pressure bath applies a force to the facets of the convex |
| 2687 |
+ |
hull in direct proportion to the area of the facet, while the thermal |
| 2688 |
+ |
coupling depends on the solvent temperature, viscosity and the size |
| 2689 |
+ |
and shape of each facet. The thermal interactions are expressed as a |
| 2690 |
+ |
standard Langevin description of the forces, |
| 2691 |
+ |
\begin{equation} |
| 2692 |
+ |
\begin{array}{rclclcl} |
| 2693 |
+ |
{\mathbf F}_f^{\text{ext}} & = & \text{external pressure} & + & \text{drag force} & + & \text{random force} \\ |
| 2694 |
+ |
& = & -\hat{n}_f P A_f & - & \Xi_f(t) {\mathbf v}_f(t) & + & {\mathbf R}_f(t) |
| 2695 |
+ |
\end{array} |
| 2696 |
+ |
\end{equation} |
| 2697 |
+ |
Here, $A_f$ and $\hat{n}_f$ are the area and (outward-facing) normal |
| 2698 |
+ |
vectors for facet $f$, respectively. ${\mathbf v}_f(t)$ is the |
| 2699 |
+ |
velocity of the facet centroid, |
| 2700 |
+ |
\begin{equation} |
| 2701 |
+ |
{\mathbf v}_f(t) = \frac{1}{3} \sum_{i=1}^{3} {\mathbf v}_i, |
| 2702 |
+ |
\end{equation} |
| 2703 |
+ |
and $\Xi_f(t)$ is an approximate ($3 \times 3$) resistance tensor that |
| 2704 |
+ |
depends on the geometry and surface area of facet $f$ and the |
| 2705 |
+ |
viscosity of the bath. The resistance tensor is related to the |
| 2706 |
+ |
fluctuations of the random force, $\mathbf{R}(t)$, by the |
| 2707 |
+ |
fluctuation-dissipation theorem (see Eq. \ref{eq:randomForce}). |
| 2708 |
+ |
|
| 2709 |
+ |
Once the resistance tensor is known for a given facet, a stochastic |
| 2710 |
+ |
vector that has the properties in Eq. (\ref{eq:randomForce}) can be |
| 2711 |
+ |
calculated efficiently by carrying out a Cholesky decomposition to |
| 2712 |
+ |
obtain the square root matrix of the resistance tensor (see |
| 2713 |
+ |
Eq. \ref{eq:Cholesky}). |
| 2714 |
+ |
|
| 2715 |
+ |
Our treatment of the resistance tensor for the Langevin Hull facets is |
| 2716 |
+ |
approximate. $\Xi_f$ for a rigid triangular plate would normally be |
| 2717 |
+ |
treated as a $6 \times 6$ tensor that includes translational and |
| 2718 |
+ |
rotational drag as well as translational-rotational coupling. The |
| 2719 |
+ |
computation of resistance tensors for rigid bodies has been detailed |
| 2720 |
+ |
elsewhere,\cite{JoseGarciadelaTorre02012000,Garcia-de-la-Torre:2001wd,GarciadelaTorreJ2002,Sun:2008fk} |
| 2721 |
+ |
but the standard approach involving bead approximations would be |
| 2722 |
+ |
prohibitively expensive if it were recomputed at each step in a |
| 2723 |
+ |
molecular dynamics simulation. |
| 2724 |
+ |
|
| 2725 |
+ |
Instead, we are utilizing an approximate resistance tensor obtained by |
| 2726 |
+ |
first constructing the Oseen tensor for the interaction of the |
| 2727 |
+ |
centroid of the facet ($f$) with each of the subfacets $\ell=1,2,3$, |
| 2728 |
+ |
\begin{equation} |
| 2729 |
+ |
T_{\ell f}=\frac{A_\ell}{8\pi\eta R_{\ell f}}\left(I + |
| 2730 |
+ |
\frac{\mathbf{R}_{\ell f}\mathbf{R}_{\ell f}^T}{R_{\ell f}^2}\right) |
| 2731 |
+ |
\end{equation} |
| 2732 |
+ |
Here, $A_\ell$ is the area of subfacet $\ell$ which is a triangle |
| 2733 |
+ |
containing two of the vertices of the facet along with the centroid. |
| 2734 |
+ |
$\mathbf{R}_{\ell f}$ is the vector between the centroid of facet $f$ |
| 2735 |
+ |
and the centroid of sub-facet $\ell$, and $I$ is the ($3 \times 3$) |
| 2736 |
+ |
identity matrix. $\eta$ is the viscosity of the external bath. |
| 2737 |
+ |
|
| 2738 |
+ |
The tensors for each of the sub-facets are added together, and the |
| 2739 |
+ |
resulting matrix is inverted to give a $3 \times 3$ resistance tensor |
| 2740 |
+ |
for translations of the triangular facet, |
| 2741 |
+ |
\begin{equation} |
| 2742 |
+ |
\Xi_f(t) =\left[\sum_{i=1}^3 T_{if}\right]^{-1}. |
| 2743 |
+ |
\end{equation} |
| 2744 |
+ |
Note that this treatment ignores rotations (and |
| 2745 |
+ |
translational-rotational coupling) of the facet. In compact systems, |
| 2746 |
+ |
the facets stay relatively fixed in orientation between |
| 2747 |
+ |
configurations, so this appears to be a reasonably good approximation. |
| 2748 |
+ |
|
| 2749 |
+ |
At each |
| 2750 |
+ |
molecular dynamics time step, the following process is carried out: |
| 2751 |
+ |
\begin{enumerate} |
| 2752 |
+ |
\item The standard inter-atomic forces ($\nabla_iU$) are computed. |
| 2753 |
+ |
\item Delaunay triangulation is carried out using the current atomic |
| 2754 |
+ |
configuration. |
| 2755 |
+ |
\item The convex hull is computed and facets are identified. |
| 2756 |
+ |
\item For each facet: |
| 2757 |
+ |
\begin{itemize} |
| 2758 |
+ |
\item[a.] The force from the pressure bath ($-\hat{n}_fPA_f$) is |
| 2759 |
+ |
computed. |
| 2760 |
+ |
\item[b.] The resistance tensor ($\Xi_f(t)$) is computed using the |
| 2761 |
+ |
viscosity ($\eta$) of the bath. |
| 2762 |
+ |
\item[c.] Facet drag ($-\Xi_f(t) \mathbf{v}_f(t)$) forces are |
| 2763 |
+ |
computed. |
| 2764 |
+ |
\item[d.] Random forces ($\mathbf{R}_f(t)$) are computed using the |
| 2765 |
+ |
resistance tensor and the temperature ($T$) of the bath. |
| 2766 |
+ |
\end{itemize} |
| 2767 |
+ |
\item The facet forces are divided equally among the vertex atoms. |
| 2768 |
+ |
\item Atomic positions and velocities are propagated. |
| 2769 |
+ |
\end{enumerate} |
| 2770 |
+ |
The Delaunay triangulation and computation of the convex hull are done |
| 2771 |
+ |
using calls to the qhull library,\cite{Qhull} and for this reason, if |
| 2772 |
+ |
qhull is not detected during the build, the Langevin Hull integrator |
| 2773 |
+ |
will not be available. There is a minimal penalty for computing the |
| 2774 |
+ |
convex hull and resistance tensors at each step in the molecular |
| 2775 |
+ |
dynamics simulation (roughly 0.02 $\times$ cost of a single force |
| 2776 |
+ |
evaluation). |
| 2777 |
+ |
|
| 2778 |
+ |
\begin{longtable}[c]{GBF} |
| 2779 |
+ |
\caption{Meta-data Keywords: Required parameters for the Langevin Hull integrator} |
| 2780 |
+ |
\\ |
| 2781 |
+ |
{\bf keyword} & {\bf units} & {\bf use} \\ \hline |
| 2782 |
+ |
\endhead |
| 2783 |
+ |
\hline |
| 2784 |
+ |
\endfoot |
| 2785 |
+ |
{\tt viscosity} & poise & Sets the value of viscosity of the implicit |
| 2786 |
+ |
solven . \\ |
| 2787 |
+ |
{\tt targetTemp} & K & Sets the target temperature of the system. |
| 2788 |
+ |
This parameter must be specified to use Langevin Hull dynamics. \\ |
| 2789 |
+ |
{\tt targetPressure} & atm & Sets the target pressure of the system. |
| 2790 |
+ |
This parameter must be specified to use Langevin Hull dynamics. \\ |
| 2791 |
+ |
{\tt usePeriodicBoundaryConditions} & logical & Turns off periodic boundary conditions. |
| 2792 |
+ |
This parameter must be set to \tt false \\ |
| 2793 |
+ |
\label{table:lhullParameters} |
| 2794 |
|
\end{longtable} |
| 2795 |
|
|
| 2796 |
+ |
|
| 2797 |
|
\section{\label{sec:constraints}Constraint Methods} |
| 2798 |
|
|
| 2799 |
|
\subsection{\label{section:rattle}The {\sc rattle} Method for Bond |
| 2942 |
|
\label{table:zconParams} |
| 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 |
| 2948 |
< |
reference structure. A restraint can be a collective |
| 2945 |
> |
% \chapter{\label{section:restraints}Restraints} |
| 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} |
| 2951 |
|
|
| 3083 |
|
\mbox{rad}^{-2}$ & & spring constant for rotation around z-axis in |
| 3084 |
|
Einstein crystal |
| 3085 |
|
\label{table:thermIntParams} |
| 3086 |
+ |
\end{longtable} |
| 3087 |
+ |
|
| 3088 |
+ |
\chapter{\label{section:rnemd}RNEMD} |
| 3089 |
+ |
|
| 3090 |
+ |
There are many ways to compute transport properties from molecular |
| 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. These approaches |
| 3094 |
+ |
are generally 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 |
| 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 |
+ |
\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 |
+ |
The original ``swapping'' approaches by M\"{u}ller-Plathe {\it et |
| 3124 |
+ |
al.}\cite{ISI:000080382700030,MullerPlathe:1997xw} can be understood |
| 3125 |
+ |
as a sequence of imaginary elastic collisions between particles in |
| 3126 |
+ |
opposite slabs. In each collision, the entire momentum vectors of |
| 3127 |
+ |
both particles may be exchanged to generate a thermal |
| 3128 |
+ |
flux. Alternatively, a single component of the momentum vectors may be |
| 3129 |
+ |
exchanged to generate a shear flux. This algorithm turns out to be |
| 3130 |
+ |
quite useful in many simulations. However, the M\"{u}ller-Plathe |
| 3131 |
+ |
swapping approach perturbs the system away from ideal |
| 3132 |
+ |
Maxwell-Boltzmann distributions, and this may leads to undesirable |
| 3133 |
+ |
side-effects when the applied flux becomes large.\cite{Maginn:2010} |
| 3134 |
+ |
This limits the application of the swapping algorithm, so in OpenMD, |
| 3135 |
+ |
we implement two additional algorithms for RNEMD in addition to the |
| 3136 |
+ |
original swapping approach. |
| 3137 |
+ |
|
| 3138 |
+ |
{\bf Non-Isotropic Velocity Scaling (NIVS):}\cite{kuang:164101} |
| 3139 |
+ |
Instead of having momentum exchange between {\it individual particles} |
| 3140 |
+ |
in each slab, the NIVS algorithm applies velocity scaling to all of |
| 3141 |
+ |
the selected particles in both slabs. A combination of linear |
| 3142 |
+ |
momentum, kinetic energy, and flux constraint equations governs the |
| 3143 |
+ |
amount of velocity scaling performed at each step. Interested readers |
| 3144 |
+ |
should consult ref. \citealp{kuang:164101} for further details on the |
| 3145 |
+ |
methodology. |
| 3146 |
+ |
|
| 3147 |
+ |
NIVS has been shown to be very effective at producing thermal |
| 3148 |
+ |
gradients and for computing thermal conductivities, particularly for |
| 3149 |
+ |
heterogeneous interfaces. Although the NIVS algorithm can also be |
| 3150 |
+ |
applied to impose a directional momentum flux, thermal anisotropy was |
| 3151 |
+ |
observed in relatively high flux simulations, and the method is not |
| 3152 |
+ |
suitable for imposing a shear flux. |
| 3153 |
+ |
|
| 3154 |
+ |
{\bf Velocity Shearing and Scaling (VSS)}:\cite{2012MolPh.110..691K} |
| 3155 |
+ |
The third RNEMD algorithm implemented in OpenMD utilizes a series of |
| 3156 |
+ |
simultaneous velocity shearing and scaling exchanges between the two |
| 3157 |
+ |
slabs. This method results in a set of simpler equations to satisfy |
| 3158 |
+ |
the conservation constraints while creating a desired flux between the |
| 3159 |
+ |
two slabs. |
| 3160 |
+ |
|
| 3161 |
+ |
The VSS approach is versatile in that it may be used to implement both |
| 3162 |
+ |
thermal and shear transport either separately or simultaneously. |
| 3163 |
+ |
Perturbations of velocities away from the ideal Maxwell-Boltzmann |
| 3164 |
+ |
distributions are minimal, and thermal anisotropy is kept to a |
| 3165 |
+ |
minimum. This ability to generate simultaneous thermal and shear |
| 3166 |
+ |
fluxes has been utilized to map out the shear viscosity of SPC/E water |
| 3167 |
+ |
in a wide range of temperature (90~K) just with a single simulation. |
| 3168 |
+ |
VSS-RNEMD also allows the directional momentum flux to have quite |
| 3169 |
+ |
arbitary directions, which could aid in the study of anisotropic solid |
| 3170 |
+ |
surfaces in contact with liquid environments. |
| 3171 |
+ |
|
| 3172 |
+ |
{\bf What the user needs to specify:} To carry out a RNEMD simulation, |
| 3173 |
+ |
a user must specify a number of parameters in the MetaData (.md) file. |
| 3174 |
+ |
Because the RNEMD methods have a large number of parameters, these |
| 3175 |
+ |
must be enclosed in a {\tt RNEMD\{...\}} block. The most important |
| 3176 |
+ |
parameters to specify are the {\tt useRNEMD}, {\tt fluxType} and flux |
| 3177 |
+ |
parameters. Most other parameters (summarized in table |
| 3178 |
+ |
\ref{table:rnemd}) have reasonable default values. {\tt fluxType} |
| 3179 |
+ |
sets up the kind of exchange that will be carried out between the two |
| 3180 |
+ |
slabs (either Kinetic Energy ({\tt KE}) or momentum ({\tt Px, Py, Pz, |
| 3181 |
+ |
Pvector}), or some combination of these). The flux is specified |
| 3182 |
+ |
with the use of three possible parameters: {\tt kineticFlux} for |
| 3183 |
+ |
kinetic energy exchange, as well as {\tt momentumFlux} or {\tt |
| 3184 |
+ |
momentumFluxVector} for simulations with directional exchange. |
| 3185 |
+ |
|
| 3186 |
+ |
{\bf How to process the results:} OpenMD will generate a {\tt .rnemd} |
| 3187 |
+ |
file with the same prefix as the original {\tt .md} file. This file |
| 3188 |
+ |
contains a running average of properties of interest computed within a |
| 3189 |
+ |
set of bins that divide the simulation cell along the $z$-axis. The |
| 3190 |
+ |
first column of the {\tt .rnemd} file is the $z$ coordinate of the |
| 3191 |
+ |
center of each bin, while following columns may contain the average |
| 3192 |
+ |
temperature, velocity, or density within each bin. The output format |
| 3193 |
+ |
in the {\tt .rnemd} file can be altered with the {\tt outputFields}, |
| 3194 |
+ |
{\tt outputBins}, and {\tt outputFileName} parameters. A report at |
| 3195 |
+ |
the top of the {\tt .rnemd} file contains the current exchange totals |
| 3196 |
+ |
as well as the average flux applied during the simulation. Using the |
| 3197 |
+ |
slope of the temperature or velocity gradient obtaine from the {\tt |
| 3198 |
+ |
.rnemd} file along with the applied flux, the user can very simply |
| 3199 |
+ |
arrive at estimates of thermal conductivities ($\lambda$), |
| 3200 |
+ |
\begin{equation} |
| 3201 |
+ |
J_z = -\lambda \frac{\partial T}{\partial z}, |
| 3202 |
+ |
\end{equation} |
| 3203 |
+ |
and shear viscosities ($\eta$), |
| 3204 |
+ |
\begin{equation} |
| 3205 |
+ |
j_z(p_x) = -\eta \frac{\partial \langle v_x \rangle}{\partial z}. |
| 3206 |
+ |
\end{equation} |
| 3207 |
+ |
Here, the quantities on the left hand side are the actual flux values |
| 3208 |
+ |
(in the header of the {\tt .rnemd} file), while the slopes are |
| 3209 |
+ |
obtained from linear fits to the gradients observed in the {\tt |
| 3210 |
+ |
.rnemd} file. |
| 3211 |
+ |
|
| 3212 |
+ |
More complicated simulations (including interfaces) require a bit more |
| 3213 |
+ |
care. Here the second derivative may be required to compute the |
| 3214 |
+ |
interfacial thermal conductance, |
| 3215 |
+ |
\begin{align} |
| 3216 |
+ |
G^\prime &= \left(\nabla\lambda \cdot \mathbf{\hat{n}}\right)_{z_0} \\ |
| 3217 |
+ |
&= \frac{\partial}{\partial z}\left(-\frac{J_z}{ |
| 3218 |
+ |
\left(\frac{\partial T}{\partial z}\right)}\right)_{z_0} \\ |
| 3219 |
+ |
&= J_z\left(\frac{\partial^2 T}{\partial z^2}\right)_{z_0} \Big/ |
| 3220 |
+ |
\left(\frac{\partial T}{\partial z}\right)_{z_0}^2. |
| 3221 |
+ |
\label{derivativeG} |
| 3222 |
+ |
\end{align} |
| 3223 |
+ |
where $z_0$ is the location of the interface between two materials and |
| 3224 |
+ |
$\mathbf{\hat{n}}$ is a unit vector normal to the interface. We |
| 3225 |
+ |
suggest that users interested in interfacial conductance consult |
| 3226 |
+ |
reference \citealp{kuang:AuThl} for other approaches to computing $G$. |
| 3227 |
+ |
Users interested in {\it friction coefficients} at heterogeneous |
| 3228 |
+ |
interfaces may also find reference \citealp{2012MolPh.110..691K} |
| 3229 |
+ |
useful. |
| 3230 |
+ |
|
| 3231 |
+ |
\newpage |
| 3232 |
+ |
|
| 3233 |
+ |
\begin{longtable}[c]{JKLM} |
| 3234 |
+ |
\caption{The following keywords must be enclosed inside a {\tt RNEMD\{\}} block} |
| 3235 |
+ |
\\ |
| 3236 |
+ |
{\bf keyword} & {\bf units} & {\bf use} & {\bf remarks} \\ \hline |
| 3237 |
+ |
\endhead |
| 3238 |
+ |
\hline |
| 3239 |
+ |
\endfoot |
| 3240 |
+ |
{\tt useRNEMD} & logical & perform RNEMD? & default is ``false'' \\ |
| 3241 |
+ |
{\tt objectSelection} & string & see section \ref{section:syntax} |
| 3242 |
+ |
for selection syntax & default is ``select all'' \\ |
| 3243 |
+ |
{\tt method} & string & exchange method & one of the following: |
| 3244 |
+ |
{\tt Swap, NIVS,} or {\tt VSS} (default is {\tt VSS}) \\ |
| 3245 |
+ |
{\tt fluxType} & string & what is being exchanged between slabs? & one |
| 3246 |
+ |
of the following: {\tt KE, Px, Py, Pz, Pvector, KE+Px, KE+Py, KE+Pvector} \\ |
| 3247 |
+ |
{\tt kineticFlux} & kcal mol$^{-1}$ \AA$^{-2}$ fs$^{-1}$ & specify the kinetic energy flux & \\ |
| 3248 |
+ |
{\tt momentumFlux} & amu \AA$^{-1}$ fs$^{-2}$ & specify the momentum flux & \\ |
| 3249 |
+ |
{\tt momentumFluxVector} & amu \AA$^{-1}$ fs$^{-2}$ & specify the momentum flux when |
| 3250 |
+ |
{\tt Pvector} is part of the exchange & Vector3d input\\ |
| 3251 |
+ |
{\tt exchangeTime} & fs & how often to perform the exchange & default is 100 fs\\ |
| 3252 |
+ |
|
| 3253 |
+ |
{\tt slabWidth} & $\mbox{\AA}$ & width of the two exchange slabs & default is $\mathsf{H}_{zz} / 10.0$ \\ |
| 3254 |
+ |
{\tt slabAcenter} & $\mbox{\AA}$ & center of the end slab & default is 0 \\ |
| 3255 |
+ |
{\tt slabBcenter} & $\mbox{\AA}$ & center of the middle slab & default is $\mathsf{H}_{zz} / 2$ \\ |
| 3256 |
+ |
{\tt outputFileName} & string & file name for output histograms & default is the same prefix as the |
| 3257 |
+ |
.md file, but with the {\tt .rnemd} extension \\ |
| 3258 |
+ |
{\tt outputBins} & int & number of $z$-bins in the output histogram & |
| 3259 |
+ |
default is 20 \\ |
| 3260 |
+ |
{\tt outputFields} & string & columns to print in the {\tt .rnemd} |
| 3261 |
+ |
file where each column is separated by a pipe ($\mid$) symbol. & Allowed column names are: {\sc z, temperature, velocity, density} \\ |
| 3262 |
+ |
\label{table:rnemd} |
| 3263 |
|
\end{longtable} |
| 3264 |
|
|
| 3265 |
|
|
| 3534 |
|
\end{center} |
| 3535 |
|
|
| 3536 |
|
For example, the phrase {\tt select mass > 16.0 and charge < -2} |
| 3537 |
< |
wouldselect StuntDoubles which have mass greater than 16.0 and charges |
| 3537 |
> |
would select StuntDoubles which have mass greater than 16.0 and charges |
| 3538 |
|
less than -2. |
| 3539 |
|
|
| 3540 |
|
\subsection{\label{section:within}Within expressions} |
| 3936 |
|
DMR-0079647. |
| 3937 |
|
|
| 3938 |
|
|
| 3939 |
< |
\bibliographystyle{jcc} |
| 3939 |
> |
\bibliographystyle{aip} |
| 3940 |
|
\bibliography{openmdDoc} |
| 3941 |
|
|
| 3942 |
|
\end{document} |
