trunk/nivsRnemd/nivsRnemd.tex

\documentclass[11pt]{article}
\usepackage{amsmath}
\usepackage{amssymb}
\usepackage{setspace}
\usepackage{endfloat}
\usepackage{caption}
%\usepackage{tabularx}
\usepackage{graphicx}
\usepackage{multirow}
%\usepackage{booktabs}
%\usepackage{bibentry}
%\usepackage{mathrsfs}
%\usepackage[ref]{overcite}
\usepackage[square, comma, sort&compress]{natbib}
\usepackage{url}
\pagestyle{plain} \pagenumbering{arabic} \oddsidemargin 0.0cm
\evensidemargin 0.0cm \topmargin -21pt \headsep 10pt \textheight
9.0in \textwidth 6.5in \brokenpenalty=10000

% double space list of tables and figures
\AtBeginDelayedFloats{\renewcommand{\baselinestretch}{1.66}}
\setlength{\abovecaptionskip}{20 pt}
\setlength{\belowcaptionskip}{30 pt}

%\renewcommand\citemid{\ } % no comma in optional referenc note
\bibpunct{[}{]}{,}{s}{}{;}
\bibliographystyle{aip}

\begin{document}

\title{A gentler approach to RNEMD: Non-isotropic Velocity Scaling for computing Thermal Conductivity and Shear Viscosity}

\author{Shenyu Kuang and J. Daniel
Gezelter\footnote{Corresponding author. \ Electronic mail: gezelter@nd.edu} \\
Department of Chemistry and Biochemistry,\\
University of Notre Dame\\
Notre Dame, Indiana 46556}

\date{\today}

\maketitle 

\begin{doublespace}

\begin{abstract}
  We present a new method for introducing stable non-equilibrium
  velocity and temperature gradients in molecular dynamics simulations
  of heterogeneous systems.  This method extends earlier Reverse
  Non-Equilibrium Molecular Dynamics (RNEMD) methods which use
  momentum exchange swapping moves. The standard swapping moves can
  create non-thermal velocity distributions and are difficult to use
  for interfacial calculations.  By using non-isotropic velocity
  scaling (NIVS) on the molecules in specific regions of a system, it
  is possible to impose momentum or thermal flux between regions of a
  simulation while conserving the linear momentum and total energy of
  the system.  To test the method, we have computed the thermal
  conductivity of model liquid and solid systems as well as the
  interfacial thermal conductivity of a metal-water interface.  We
  find that the NIVS-RNEMD improves the problematic velocity
  distributions that develop in other RNEMD methods.
\end{abstract}

\newpage

%\narrowtext

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                          BODY OF TEXT
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\section{Introduction}
The original formulation of Reverse Non-equilibrium Molecular Dynamics
(RNEMD) obtains transport coefficients (thermal conductivity and shear
viscosity) in a fluid by imposing an artificial momentum flux between
two thin parallel slabs of material that are spatially separated in
the simulation cell.\cite{MullerPlathe:1997xw,ISI:000080382700030} The
artificial flux is typically created by periodically ``swapping''
either the entire momentum vector $\vec{p}$ or single components of
this vector ($p_x$) between molecules in each of the two slabs.  If
the two slabs are separated along the $z$ coordinate, the imposed flux
is either directional ($j_z(p_x)$) or isotropic ($J_z$), and the
response of a simulated system to the imposed momentum flux will
typically be a velocity or thermal gradient (Fig. \ref{thermalDemo}).
The shear viscosity ($\eta$) and thermal conductivity ($\lambda$) are
easily obtained by assuming linear response of the system,
\begin{eqnarray}
j_z(p_x) & = & -\eta \frac{\partial v_x}{\partial z}\\
J_z & = & \lambda \frac{\partial T}{\partial z}
\end{eqnarray}
RNEMD has been widely used to provide computational estimates of
thermal conductivities and shear viscosities in a wide range of
materials, from liquid copper to both monatomic and molecular fluids
(e.g.  ionic
liquids).\cite{Bedrov:2000-1,Bedrov:2000,Muller-Plathe:2002,ISI:000184808400018,ISI:000231042800044,Maginn:2007,Muller-Plathe:2008,ISI:000258460400020,ISI:000258840700015,ISI:000261835100054}

\begin{figure}
\includegraphics[width=\linewidth]{thermalDemo}
\caption{RNEMD methods impose an unphysical transfer of momentum or
  kinetic energy between a ``hot'' slab and a ``cold'' slab in the
  simulation box.  The molecular system responds to this imposed flux
  by generating a momentum or temperature gradient.  The slope of the
  gradient can then be used to compute transport properties (e.g.
  shear viscosity and thermal conductivity).}
\label{thermalDemo}
\end{figure}

RNEMD is preferable in many ways to the forward NEMD
methods\cite{ISI:A1988Q205300014,hess:209,Vasquez:2004fk,backer:154503,ISI:000266247600008}
because it imposes what is typically difficult to measure (a flux or
stress) and it is typically much easier to compute the response
(momentum gradients or strains). For similar reasons, RNEMD is also
preferable to slowly-converging equilibrium methods for measuring
thermal conductivity and shear viscosity (using Green-Kubo
relations\cite{daivis:541,mondello:9327} or the Helfand moment
approach of Viscardy {\it et
  al.}\cite{Viscardy:2007bh,Viscardy:2007lq}) because these rely on
computing and integrating long-time correlation functions that are
subject to noise issues.

Another attractive feature of RNEMD is that it conserves both total
linear momentum and total energy during the swaps (as long as the two
molecules have the same identity), so the swapped configurations are
typically samples from the same manifold of states in the
microcanonical ensemble. Furthermore, the method is applicable with
different ensembles, unlike the heat-exchange (HeX) algorithm proposed
by Hafskjold {\it et al.} \cite{HeX:1994,HeX:1993}, which is
incompatible with non-microcanonical integrators.

Recently, Tenney and Maginn\cite{Maginn:2010} have discovered some
problems with the original RNEMD swap technique.  Notably, large
momentum fluxes (equivalent to frequent momentum swaps between the
slabs) can result in ``notched'', ``peaked'' and generally non-thermal
momentum distributions in the two slabs, as well as non-linear thermal
and velocity distributions along the direction of the imposed flux
($z$). Tenney and Maginn obtained reasonable limits on imposed flux
and proposed self-adjusting metrics for retaining the usability of the
method.

In this paper, we develop and test a method for non-isotropic velocity
scaling (NIVS) which retains the desirable features of RNEMD
(conservation of linear momentum and total energy, compatibility with
periodic boundary conditions) while establishing true thermal
distributions in each of the two slabs. In the next section, we
present the method for determining the scaling constraints.  We then
test the method on both liquids and solids as well as a non-isotropic
liquid-solid interface and show that it is capable of providing
reasonable estimates of the thermal conductivity and shear viscosity
in all of these cases.

\section{Methodology}
We retain the basic idea of M\"{u}ller-Plathe's RNEMD method; the
periodic system is partitioned into a series of thin slabs along one
axis ($z$).  One of the slabs at the end of the periodic box is
designated the ``hot'' slab, while the slab in the center of the box
is designated the ``cold'' slab.  The artificial momentum flux will be
established by transferring momentum from the cold slab and into the
hot slab.

Rather than using momentum swaps, we use a series of velocity scaling
moves.  For molecules $\{i\}$  located within the cold slab,
\begin{equation}
\vec{v}_i \leftarrow \left( \begin{array}{ccc}
x & 0 & 0 \\
0 & y & 0 \\
0 & 0 & z \\
\end{array} \right) \cdot \vec{v}_i
\end{equation}
where ${x, y, z}$ are a set of 3 velocity-scaling variables for each
of the three directions in the system.  Likewise, the molecules
$\{j\}$ located in the hot slab will see a concomitant scaling of
velocities,
\begin{equation}
\vec{v}_j \leftarrow \left( \begin{array}{ccc}
x^\prime & 0 & 0 \\
0 & y^\prime & 0 \\
0 & 0 & z^\prime \\
\end{array} \right) \cdot \vec{v}_j
\end{equation}

Conservation of linear momentum in each of the three directions
($\alpha = x,y,z$) ties the values of the hot and cold scaling
parameters together:
\begin{equation}
P_h^\alpha + P_c^\alpha = \alpha^\prime P_h^\alpha + \alpha P_c^\alpha
\end{equation}
where
\begin{eqnarray}
P_c^\alpha & = & \sum_{i = 1}^{N_c} m_i v_{i\alpha} \\
P_h^\alpha & = & \sum_{j = 1}^{N_h} m_j v_{j\alpha}
\label{eq:momentumdef}
\end{eqnarray}
Therefore, for each of the three directions, the hot scaling
parameters are a simple function of the cold scaling parameters and
the instantaneous linear momenta in each of the two slabs.
\begin{equation}
\alpha^\prime = 1 + (1 - \alpha) p_\alpha
\label{eq:hotcoldscaling}
\end{equation}
where
\begin{equation}
p_\alpha = \frac{P_c^\alpha}{P_h^\alpha}
\end{equation}
for convenience.

Conservation of total energy also places constraints on the scaling:
\begin{equation}
\sum_{\alpha = x,y,z} \left\{ K_h^\alpha + K_c^\alpha \right\} = \sum_{\alpha = x,y,z}
\left\{ \left(\alpha^\prime\right)^2 K_h^\alpha + \alpha^2 K_c^\alpha \right\}
\end{equation}
where the translational kinetic energies, $K_h^\alpha$ and
$K_c^\alpha$, are computed for each of the three directions in a
similar manner to the linear momenta (Eq. \ref{eq:momentumdef}).
Substituting in the expressions for the hot scaling parameters
($\alpha^\prime$) from Eq. (\ref{eq:hotcoldscaling}), we obtain the
{\it constraint ellipsoid}:
\begin{equation}
\sum_{\alpha = x,y,z} \left( a_\alpha \alpha^2 + b_\alpha \alpha +
  c_\alpha \right) = 0
\label{eq:constraintEllipsoid}
\end{equation}
where the constants are obtained from the instantaneous values of the
linear momenta and kinetic energies for the hot and cold slabs,
\begin{eqnarray}
a_\alpha & = &\left(K_c^\alpha + K_h^\alpha
  \left(p_\alpha\right)^2\right) \\
b_\alpha & = & -2 K_h^\alpha p_\alpha \left( 1 + p_\alpha\right) \\
c_\alpha & = & K_h^\alpha p_\alpha^2 + 2 K_h^\alpha p_\alpha - K_c^\alpha
\label{eq:constraintEllipsoidConsts}
\end{eqnarray}
This ellipsoid (Eq. \ref{eq:constraintEllipsoid}) defines the set of
cold slab scaling parameters which, when applied, preserve the linear
momentum of the system in all three directions as well as total
kinetic energy.

The goal of using these velocity scaling variables is to transfer
kinetic energy from the cold slab to the hot slab.  If the hot and
cold slabs are separated along the z-axis, the energy flux is given
simply by the decrease in kinetic energy of the cold bin:
\begin{equation}
(1-x^2) K_c^x  + (1-y^2) K_c^y + (1-z^2) K_c^z = J_z\Delta t
\end{equation}
The expression for the energy flux can be re-written as another
ellipsoid centered on $(x,y,z) = 0$:
\begin{equation}
\sum_{\alpha = x,y,z} K_c^\alpha \alpha^2 = \sum_{\alpha = x,y,z}
K_c^\alpha -J_z \Delta t
\label{eq:fluxEllipsoid}
\end{equation}
The spatial extent of the {\it thermal flux ellipsoid} is governed
both by the target flux, $J_z$ as well as the instantaneous values of
the kinetic energy components in the cold bin.

To satisfy an energetic flux as well as the conservation constraints,
we must determine the points ${x,y,z}$ that lie on both the constraint
ellipsoid (Eq. \ref{eq:constraintEllipsoid}) and the flux ellipsoid
(Eq. \ref{eq:fluxEllipsoid}), i.e. the intersection of the two
ellipsoids in 3-dimensional space.

\begin{figure}
\includegraphics[width=\linewidth]{ellipsoids}
\caption{Velocity scaling coefficients which maintain both constant
  energy and constant linear momentum of the system lie on the surface
  of the {\it constraint ellipsoid} while points which generate the
  target momentum flux lie on the surface of the {\it flux ellipsoid}.
  The velocity distributions in the cold bin are scaled by only those
  points which lie on both ellipsoids.}
\label{ellipsoids}
\end{figure}

Since ellipsoids can be expressed as polynomials up to second order in
each of the three coordinates, finding the the intersection points of
two ellipsoids is isomorphic to finding the roots a polynomial of
degree 16.  There are a number of polynomial root-finding methods in
the literature,\cite{Sommese2005,384119} but numerically finding
the roots of high-degree polynomials is generally an ill-conditioned
problem.\cite{Sommese2005} One simplification is to maintain
velocity scalings that are {\it as isotropic as possible}.  To do
this, we impose $x=y$, and treat both the constraint and flux
ellipsoids as 2-dimensional ellipses.  In reduced dimensionality, the
intersecting-ellipse problem reduces to finding the roots of
polynomials of degree 4.

Depending on the target flux and current velocity distributions, the
ellipsoids can have between 0 and 4 intersection points.  If there are
no intersection points, it is not possible to satisfy the constraints
while performing a non-equilibrium scaling move, and no change is made
to the dynamics.  

With multiple intersection points, any of the scaling points will
conserve the linear momentum and kinetic energy of the system and will
generate the correct target flux.  Although this method is inherently
non-isotropic, the goal is still to maintain the system as close to an
isotropic fluid as possible.  With this in mind, we would like the
kinetic energies in the three different directions could become as
close as each other as possible after each scaling.  Simultaneously,
one would also like each scaling to be as gentle as possible, i.e.
${x,y,z \rightarrow 1}$, in order to avoid large perturbations to the
system.  To do this, we pick the intersection point which maintains
the three scaling variables ${x, y, z}$ as well as the ratio of
kinetic energies ${K_c^z/K_c^x, K_c^z/K_c^y}$ as close as possible to
1.

After the valid scaling parameters are arrived at by solving geometric
intersection problems in $x, y, z$ space in order to obtain cold slab
scaling parameters, Eq. (\ref{eq:hotcoldscaling}) is used to
determine the conjugate hot slab scaling variables.

\subsection{Introducing shear stress via velocity scaling}
It is also possible to use this method to magnify the random
fluctuations of the average momentum in each of the bins to induce a
momentum flux.  Doing this repeatedly will create a shear stress on
the system which will respond with an easily-measured strain.  The
momentum flux (say along the $x$-direction) may be defined as:
\begin{equation}
(1-x) P_c^x = j_z(p_x)\Delta t
\label{eq:fluxPlane}
\end{equation}
This {\it momentum flux plane} is perpendicular to the $x$-axis, with
its position governed both by a target value, $j_z(p_x)$ as well as
the instantaneous value of the momentum along the $x$-direction.

In order to satisfy a momentum flux as well as the conservation
constraints, we must determine the points ${x,y,z}$ which lie on both
the constraint ellipsoid (Eq. \ref{eq:constraintEllipsoid}) and the
flux plane (Eq. \ref{eq:fluxPlane}), i.e. the intersection of an
ellipsoid and a plane in 3-dimensional space.

In the case of momentum flux transfer, we also impose another
constraint to set the kinetic energy transfer as zero. In other
words, we apply Eq. \ref{eq:fluxEllipsoid} and let ${J_z = 0}$.  With
one variable fixed by Eq. \ref{eq:fluxPlane}, one now have a similar
set of quartic equations to the above kinetic energy transfer problem.

\section{Computational Details}

We have implemented this methodology in our molecular dynamics code,
OpenMD,\cite{Meineke:2005gd,openmd} performing the NIVS scaling moves
with a variable frequency after the molecular dynamics (MD) steps.  We
have tested the method in a variety of different systems, including:
homogeneous fluids (Lennard-Jones and SPC/E water), crystalline
solids, using both the embedded atom method
(EAM)~\cite{PhysRevB.33.7983} and quantum Sutton-Chen
(QSC)~\cite{PhysRevB.59.3527} models for Gold, and heterogeneous
interfaces (QSC gold - SPC/E water).  

\subsection{Simulation Cells}

In each of the systems studied, the dynamics was carried out in a
rectangular simulation cell using periodic boundary conditions in all
three dimensions.  The cells were longer along the $z$ axis and the
space was divided into $N$ slabs along this axis (typically $N=20$).
The top slab ($n=1$) was designated the ``hot'' slab, while the
central slab ($n= N/2 + 1$) was designated the ``cold'' slab. In all
cases, simulations were first thermalized in canonical ensemble (NVT)
using a Nos\'{e}-Hoover thermostat.\cite{Hoover85} then equilibrated in
microcanonical ensemble (NVE) before introducing any non-equilibrium
method.

\subsection{RNEMD with M\"{u}ller-Plathe swaps}

In order to compare our new methodology with the original
M\"{u}ller-Plathe swapping algorithm,\cite{ISI:000080382700030} we
first performed simulations using the original technique. At fixed
intervals, kinetic energy or momentum exchange moves were performed
between the hot and the cold slabs.  The interval between exchange
moves governs the effective momentum flux ($j_z(p_x)$) or energy flux
($J_z$) between the two slabs so to vary this quantity, we performed
simulations with a variety of delay intervals between the swapping moves.

For thermal conductivity measurements, the particle with smallest
speed in the hot slab and the one with largest speed in the cold slab
had their entire momentum vectors swapped.  In the test cases run
here, all particles had the same chemical identity and mass, so this
move preserves both total linear momentum and total energy.  It is
also possible to exchange energy by assuming an elastic collision
between the two particles which are exchanging energy.

For shear stress simulations, the particle with the most negative
$p_x$ in the hot slab and the one with the most positive $p_x$ in the
cold slab exchanged only this component of their momentum vectors.

\subsection{RNEMD with NIVS scaling}

For each simulation utilizing the swapping method, a corresponding
NIVS-RNEMD simulation was carried out using a target momentum flux set
to produce the same flux experienced in the swapping simulation.

To test the temperature homogeneity, momentum and temperature
distributions (for all three dimensions) were accumulated for
molecules in each of the slabs.  Transport coefficients were computed
using the temperature (and momentum) gradients across the $z$-axis as
well as the total momentum flux the system experienced during data
collection portion of the simulation.

\subsection{Shear viscosities}

The momentum flux was calculated using the total non-physical momentum
transferred (${P_x}$) and the data collection time ($t$):
\begin{equation}
j_z(p_x) = \frac{P_x}{2 t L_x L_y}
\end{equation}
where $L_x$ and $L_y$ denote the $x$ and $y$ lengths of the simulation
box.  The factor of two in the denominator is present because physical
momentum transfer between the slabs occurs in two directions ($+z$ and
$-z$).  The velocity gradient ${\langle \partial v_x /\partial z
  \rangle}$ was obtained using linear regression of the mean $x$
component of the velocity, $\langle v_x \rangle$, in each of the bins.
For Lennard-Jones simulations, shear viscosities are reported in
reduced units (${\eta^* = \eta \sigma^2 (\varepsilon m)^{-1/2}}$).

\subsection{Thermal Conductivities}

The energy flux was calculated in a similar manner to the momentum
flux, using the total non-physical energy transferred (${E_{total}}$)
and the data collection time $t$:
\begin{equation}
J_z = \frac{E_{total}}{2 t L_x L_y}
\end{equation}
The temperature gradient ${\langle\partial T/\partial z\rangle}$ was
obtained by a linear regression of the temperature profile. For
Lennard-Jones simulations, thermal conductivities are reported in
reduced units (${\lambda^*=\lambda \sigma^2 m^{1/2}
  k_B^{-1}\varepsilon^{-1/2}}$).

\subsection{Interfacial Thermal Conductivities}

For interfaces with a relatively low interfacial conductance, the bulk
regions on either side of an interface rapidly come to a state in
which the two phases have relatively homogeneous (but distinct)
temperatures.  The interfacial thermal conductivity $G$ can therefore
be approximated as:

\begin{equation}
G = \frac{E_{total}}{2 t L_x L_y \left( \langle T_\mathrm{hot}\rangle -
    \langle T_\mathrm{cold}\rangle \right)}
\label{interfaceCalc}
\end{equation}
where ${E_{total}}$ is the imposed non-physical kinetic energy
transfer and ${\langle T_\mathrm{hot}\rangle}$ and ${\langle
  T_\mathrm{cold}\rangle}$ are the average observed temperature of the
two separated phases.  If the interfacial conductance is {\it not}
small, it would also be possible to compute the interfacial thermal
conductivity using this method by computing the change in the slope of
the thermal gradient ($\partial^2 \langle T \rangle /
\partial z^2$) at the interface.

\section{Results}

\subsection{Lennard-Jones Fluid}
2592 Lennard-Jones atoms were placed in an orthorhombic cell
${10.06\sigma \times 10.06\sigma \times 30.18\sigma}$ on a side.  The
reduced density ${\rho^* = \rho\sigma^3}$ was thus 0.85, which enabled
direct comparison between our results and previous methods.  These
simulations were carried out with a reduced timestep ${\tau^* =
  4.6\times10^{-4}}$.  For the shear viscosity calculations, the mean
temperature was ${T^* = k_B T/\varepsilon = 0.72}$.  For thermal
conductivity calculations, simulations were run under reduced
temperature ${\langle T^*\rangle = 0.72}$ in the microcanonical
ensemble. The simulations included $10^5$ steps of equilibration
without any momentum flux, $10^5$ steps of stablization with an
imposed momentum transfer to create a gradient, and $10^6$ steps of
data collection under RNEMD.

\subsubsection*{Thermal Conductivity}

Our thermal conductivity calculations show that the NIVS method agrees
well with the swapping method. Five different swap intervals were
tested (Table \ref{LJ}). Similar thermal gradients were observed with
similar thermal flux under the two different methods (Fig.
\ref{thermalGrad}). Furthermore, the 1-d temperature profiles showed
no observable differences between the $x$, $y$ and $z$ axes (lower
panel of Fig. \ref{thermalGrad}), so even though we are using a
non-isotropic scaling method, none of the three directions are
experience disproportionate heating due to the velocity scaling.

\begin{table*}
  \begin{minipage}{\linewidth}
    \begin{center}

      \caption{Thermal conductivity ($\lambda^*$) and shear viscosity
        ($\eta^*$) (in reduced units) of a Lennard-Jones fluid at
        ${\langle T^* \rangle = 0.72}$ and ${\rho^* = 0.85}$ computed
        at various momentum fluxes.  The original swapping method and
        the velocity scaling method give similar results.
        Uncertainties are indicated in parentheses.}
      
      \begin{tabular}{|cc|cc|cc|}
        \hline
        \multicolumn{2}{|c}{Momentum Exchange} &
        \multicolumn{2}{|c}{Swapping RNEMD} &
        \multicolumn{2}{|c|}{NIVS-RNEMD} \\
        \hline
        Swap Interval & Equivalent $J_z^*$ or &
        \multirow{2}{*}{$\lambda^*_{swap}$} & 
        \multirow{2}{*}{$\eta^*_{swap}$}  & 
        \multirow{2}{*}{$\lambda^*_{scale}$} & 
        \multirow{2}{*}{$\eta^*_{scale}$} \\
        (reduced units) & $j_z^*(p_x)$ (reduced units) & & & & \\
        \hline 
        0.116  & 0.16  & 7.03(0.34) & 3.57(0.06) & 7.30(0.10) & 3.54(0.04)\\
        0.232  & 0.09  & 7.03(0.14) & 3.64(0.05) & 6.95(0.09) & 3.76(0.09)\\
        0.463 & 0.047 & 6.91(0.42) & 3.52(0.16) & 7.19(0.07) & 3.66(0.06)\\
        0.926 & 0.024 & 7.52(0.15) & 3.72(0.05) & 7.19(0.28) & 3.32(0.18)\\
        1.158 & 0.019 & 7.41(0.29) & 3.42(0.06) & 7.98(0.33) & 3.43(0.08)\\
        \hline
      \end{tabular}
      \label{LJ}
    \end{center}
  \end{minipage}
\end{table*}

\begin{figure}
  \includegraphics[width=\linewidth]{thermalGrad}
  \caption{The NIVS-RNEMD method creates similar temperature gradients
    compared with the swapping method under a variety of imposed
    kinetic energy flux values. Furthermore, the implementation of
    Non-Isotropic Velocity Scaling does not cause temperature
    anisotropy to develop in thermal conductivity calculations.}
  \label{thermalGrad}
\end{figure}

\subsubsection*{Velocity Distributions}

To test the effects on the velocity distributions, we accumulated
velocities every 100 steps and produced distributions of both velocity
and speed in each of the slabs. From these distributions, we observed
that under high non-physical kinetic energy flux, the speed of
molecules in slabs where {\it swapping} occured could deviate from the
Maxwell-Boltzmann distribution. This behavior was also noted by Tenney
and Maginn.\cite{Maginn:2010} Figure \ref{thermalHist} shows how these
distributions deviate from an ideal distribution. In the ``hot'' slab,
the probability density is notched at low speeds and has a substantial
shoulder at higher speeds relative to the ideal distribution. In the
cold slab, the opposite notching and shouldering occurs.  This
phenomenon is more obvious at high swapping rates.

The peak of the velocity distribution is substantially flattened in
the hot slab, and is overly sharp (with truncated wings) in the cold
slab. This problem is rooted in the mechanism of the swapping method.
Continually depleting low (high) speed particles in the high (low)
temperature slab is not complemented by diffusion of low (high) speed
particles from neighboring slabs unless the swapping rate is
sufficiently small. Simutaneously, surplus low speed particles in the
cold slab do not have sufficient time to diffuse to neighboring slabs.
Since the thermal exchange rate must reach a minimum level to produce
an observable thermal gradient, the swapping-method RNEMD has a
relatively narrow choice of exchange times that can be utilized.

For comparison, NIVS-RNEMD produces a speed distribution closer to the
Maxwell-Boltzmann curve (Fig. \ref{thermalHist}).  The reason for this
is simple; upon velocity scaling, a Gaussian distribution remains
Gaussian.  Although a single scaling move is non-isotropic in three
dimensions, our criteria for choosing a set of scaling coefficients
helps maintain the distributions as close to isotropic as possible.
Consequently, NIVS-RNEMD is able to impose a non-physical thermal flux
without large perturbations to the velocity distributions in the two
slabs.

\begin{figure}
\includegraphics[width=\linewidth]{thermalHist}
\caption{Velocity and speed distributions that develop under the
  swapping and NIVS-RNEMD methods at high flux.  The distributions for
  the hot bins (upper panels) and cold bins (lower panels) were
  obtained from Lennard-Jones simulations with $\langle T^* \rangle =
  4.5$ with a flux of $J_z^* \sim 5$ (equivalent to a swapping interval
  of 10 time steps).  This is a relatively large flux which shows the
  non-thermal distributions that develop under the swapping method.
  NIVS does a better job of producing near-thermal distributions in
  the bins.}
\label{thermalHist}
\end{figure}


\subsubsection*{Shear Viscosity}
Our calculations (Table \ref{LJ}) show that velocity-scaling RNEMD
predicted similar values for shear viscosities to the swapping RNEMD
method. The average molecular momentum gradients of these samples are
shown in the upper two panels of Fig. \ref{shear}.

\begin{figure}
  \includegraphics[width=\linewidth]{shear}
  \caption{Average momentum gradients in shear viscosity simulations,
    using ``swapping'' method (top panel) and NIVS-RNEMD method
    (middle panel). NIVS-RNEMD produces a thermal anisotropy artifact
    in the hot and cold bins when used for shear viscosity.  This
    artifact does not appear in thermal conductivity calculations.}
  \label{shear}
\end{figure}

Observations of the three one-dimensional temperatures in each of the
slabs shows that NIVS-RNEMD does produce some thermal anisotropy,
particularly in the hot and cold slabs.  (Note that these temperature
measurements include the kinetic energy contributed by the slab field
velocities.  This contribution is small relative to the total kinetic
energy in the slab, and is of a comparable magnitude to the size of
the error bars.)  The lower panel of Fig.  \ref{shear} indicates that
with a relatively large imposed momentum flux, $j_z(p_x)$, the $x$
direction approaches a different temperature from the $y$ and $z$
directions in both the hot and cold bins.  This is an artifact of the
scaling constraints.  After a momentum gradient has been established,
$P_c^x$ is always less than zero.  Consequently, the scaling factor
$x$ is always greater than one in order to satisfy the constraints.
This will continually increase the kinetic energy in $x$-dimension,
$K_c^x$.  If there is not enough time for the kinetic energy to
exchange among different directions and different slabs, the system
will exhibit the observed thermal anisotropy in the hot and cold bins.

Although results between scaling and swapping methods are comparable,
the inherent temperature anisotropy does make NIVS-RNEMD method less
attractive than swapping RNEMD for shear viscosity calculations.  We
note that this problem appears only when a large {\it linear} momentum
flux is applied, and does not appear in {\it thermal} flux
calculations.

\subsection{Bulk SPC/E water}

We compared the thermal conductivity of SPC/E water using NIVS-RNEMD
to the work of Bedrov {\it et al.}\cite{Bedrov:2000}, which employed
the original swapping RNEMD method.  Bedrov {\it et
  al.}\cite{Bedrov:2000} argued that exchange of the molecular
center-of-mass velocities instead of single atom velocities conserves
the total kinetic energy and linear momentum.  This principle is also
adopted in our simulations. Scaling was applied to the center-of-mass
velocities of SPC/E water molecules.

To construct the simulations, a simulation box consisting of 1000
molecules were first equilibrated under ambient pressure and
temperature conditions using the isobaric-isothermal (NPT)
ensemble.\cite{melchionna93} A fixed volume was chosen to match the
average volume observed in the NPT simulations, and this was followed
by equilibration, first in the canonical (NVT) ensemble, followed by a
100~ps period under constant-NVE conditions without any momentum flux.
Another 100~ps was allowed to stabilize the system with an imposed
momentum transfer to create a gradient, and 1~ns was allotted for data
collection under RNEMD.  In all simulations involving point charges,
we utilized damped shifted-force electrostatics\cite{Fennell06} which
is a variant of the Wolf summation\cite{wolf:8254} that has been shown
to provide good forces and torques on molecular models for water in a
computationally efficient manner.\cite{Fennell06} The damping
parameter ($\alpha$) was set to 0.18 \AA$^{-1}$, and the cutoff radius
was set to 12 \AA.  The total system energy and temperature were
recorded to ensure that there was no appreciable drift in either
quantity even without application of a thermostat.

In our simulations, the established temperature gradients were similar
to the previous work.  Our simulation results at 318K are in good
agreement with those from Bedrov {\it et al.} (Table
\ref{spceThermal}). And both methods yield values in reasonable
agreement with experimental values.  

\begin{table*}
  \begin{minipage}{\linewidth}
    \begin{center}
      
      \caption{Thermal conductivity of SPC/E water under various
        imposed thermal gradients. Uncertainties are indicated in
        parentheses.}
      
      \begin{tabular}{|c|c|ccc|}
        \hline
        \multirow{2}{*}{$\langle T\rangle$(K)} &
        \multirow{2}{*}{$\langle dT/dz\rangle$(K/\AA)} &
        \multicolumn{3}{|c|}{$\lambda (\mathrm{W m}^{-1}
          \mathrm{K}^{-1})$} \\
        & & This work & Previous simulations\cite{Bedrov:2000} &
        Experiment\cite{WagnerKruse}\\
        \hline
        \multirow{3}{*}{300} & 0.38 & 0.816(0.044) & &
        \multirow{3}{*}{0.61}\\
        & 0.81 & 0.770(0.008) & & \\
        & 1.54 & 0.813(0.007) & & \\
        \hline
        \multirow{2}{*}{318} & 0.75 & 0.750(0.032) & 0.784 &
        \multirow{2}{*}{0.64}\\
        & 1.59 & 0.778(0.019) & 0.730 & \\
        \hline
      \end{tabular}
      \label{spceThermal}
    \end{center}
  \end{minipage}
\end{table*}

\subsection{Crystalline Gold} 

To see how the method performed in a solid, we calculated thermal
conductivities using two atomistic models for gold.  Several different
potential models have been developed that reasonably describe
interactions in transition metals. In particular, the Embedded Atom
Model (EAM)~\cite{PhysRevB.33.7983} and Sutton-Chen (SC)~\cite{Chen90}
potential have been used to study a wide range of phenomena in both
bulk materials and
nanoparticles.\cite{ISI:000207079300006,Clancy:1992,Vardeman:2008fk,Vardeman-II:2001jn,ShibataT._ja026764r,Sankaranarayanan:2005lr,Chui:2003fk,Wang:2005qy,Medasani:2007uq}
Both potentials are based on a model of a metal which treats the
nuclei and core electrons as pseudo-atoms embedded in the electron
density due to the valence electrons on all of the other atoms in the
system. The SC potential has a simple form that closely resembles the
Lennard Jones potential,
\begin{equation}
\label{eq:SCP1} 
U_{tot}=\sum _{i}\left[ \frac{1}{2}\sum _{j\neq i}D_{ij}V^{pair}_{ij}(r_{ij})-c_{i}D_{ii}\sqrt{\rho_{i}}\right] , 
\end{equation}
where $V^{pair}_{ij}$ and $\rho_{i}$ are given by 
\begin{equation}
\label{eq:SCP2} 
V^{pair}_{ij}(r)=\left( \frac{\alpha_{ij}}{r_{ij}}\right)^{n_{ij}}, \rho_{i}=\sum_{j\neq i}\left( \frac{\alpha_{ij}}{r_{ij}}\right) ^{m_{ij}}. 
\end{equation}
$V^{pair}_{ij}$ is a repulsive pairwise potential that accounts for
interactions between the pseudoatom cores. The $\sqrt{\rho_i}$ term in
Eq. (\ref{eq:SCP1}) is an attractive many-body potential that models
the interactions between the valence electrons and the cores of the
pseudo-atoms. $D_{ij}$, $D_{ii}$ set the appropriate overall energy
scale, $c_i$ scales the attractive portion of the potential relative
to the repulsive interaction and $\alpha_{ij}$ is a length parameter
that assures a dimensionless form for $\rho$. These parameters are
tuned to various experimental properties such as the density, cohesive
energy, and elastic moduli for FCC transition metals. The quantum
Sutton-Chen (QSC) formulation matches these properties while including
zero-point quantum corrections for different transition
metals.\cite{PhysRevB.59.3527} The EAM functional forms differ
slightly from SC but the overall method is very similar.

In this work, we have utilized both the EAM and the QSC potentials to
test the behavior of scaling RNEMD.

A face-centered-cubic (FCC) lattice was prepared containing 2880 Au
atoms (i.e. ${6\times 6\times 20}$ unit cells). Simulations were run
both with and without isobaric-isothermal (NPT)~\cite{melchionna93}
pre-equilibration at a target pressure of 1 atm. When equilibrated
under NPT conditions, our simulation box expanded by approximately 1\%
in volume. Following adjustment of the box volume, equilibrations in
both the canonical and microcanonical ensembles were carried out. With
the simulation cell divided evenly into 10 slabs, different thermal
gradients were established by applying a set of target thermal fluxes.

The results for the thermal conductivity of gold are shown in Table
\ref{AuThermal}.  In these calculations, the end and middle slabs were
excluded from the thermal gradient linear regession. EAM predicts
slightly larger thermal conductivities than QSC.  However, both values
are smaller than experimental value by a factor of more than 200. This
behavior has been observed previously by Richardson and Clancy, and
has been attributed to the lack of electronic contribution in these
force fields.\cite{Clancy:1992} It should be noted that the density of
the metal being simulated has an effect on thermal conductance. With
an expanded lattice, lower thermal conductance is expected (and
observed). We also observed a decrease in thermal conductance at
higher temperatures, a trend that agrees with experimental
measurements.\cite{AshcroftMermin}

\begin{table*}
  \begin{minipage}{\linewidth}
    \begin{center}
      
      \caption{Calculated thermal conductivity of crystalline gold
        using two related force fields. Calculations were done at both
        experimental and equilibrated densities and at a range of
        temperatures and thermal flux rates. Uncertainties are
        indicated in parentheses. Richardson {\it et
          al.}\cite{Clancy:1992} give an estimate of 1.74 $\mathrm{W
          m}^{-1}\mathrm{K}^{-1}$ for EAM gold
         at a density of 19.263 g / cm$^3$.}
      
      \begin{tabular}{|c|c|c|cc|}
        \hline
        Force Field & $\rho$ (g/cm$^3$) & ${\langle T\rangle}$ (K) &
        $\langle dT/dz\rangle$ (K/\AA) & $\lambda (\mathrm{W m}^{-1} \mathrm{K}^{-1})$\\
        \hline
        \multirow{7}{*}{QSC} & \multirow{3}{*}{19.188} & \multirow{3}{*}{300} & 1.44 & 1.10(0.06)\\
        &        &     & 2.86 & 1.08(0.05)\\
        &        &     & 5.14 & 1.15(0.07)\\\cline{2-5}
        & \multirow{4}{*}{19.263} & \multirow{2}{*}{300} & 2.31 & 1.25(0.06)\\
        &        &     & 3.02 & 1.26(0.05)\\\cline{3-5}
        &        & \multirow{2}{*}{575} & 3.02 & 1.02(0.07)\\
        &        &     & 4.84 & 0.92(0.05)\\
        \hline
        \multirow{8}{*}{EAM} & \multirow{3}{*}{19.045} & \multirow{3}{*}{300} & 1.24 & 1.24(0.16)\\
        &        &     & 2.06 & 1.37(0.04)\\
        &        &     & 2.55 & 1.41(0.07)\\\cline{2-5}
        & \multirow{5}{*}{19.263} & \multirow{3}{*}{300} & 1.06 & 1.45(0.13)\\
        &        &     & 2.04 & 1.41(0.07)\\
        &        &     & 2.41 & 1.53(0.10)\\\cline{3-5}
        &        & \multirow{2}{*}{575} & 2.82 & 1.08(0.03)\\
        &        &     & 4.14 & 1.08(0.05)\\
        \hline
      \end{tabular}
      \label{AuThermal}
    \end{center}
  \end{minipage}
\end{table*}

\subsection{Thermal Conductance at the Au/H$_2$O interface}
The most attractive aspect of the scaling approach for RNEMD is the
ability to use the method in non-homogeneous systems, where molecules
of different identities are segregated in different slabs.  To test
this application, we simulated a Gold (111) / water interface.  To
construct the interface, a box containing a lattice of 1188 Au atoms
(with the 111 surface in the $+z$ and $-z$ directions) was allowed to
relax under ambient temperature and pressure.  A separate (but
identically sized) box of SPC/E water was also equilibrated at ambient
conditions.  The two boxes were combined by removing all water
molecules within 3 \AA~ radius of any gold atom.  The final
configuration contained 1862 SPC/E water molecules.

The Spohr potential was adopted in depicting the interaction between
metal atoms and water molecules.\cite{ISI:000167766600035} A similar
protocol of equilibration to our water simulations was followed.  We
observed that the two phases developed large temperature differences
even under a relatively low thermal flux.

The low interfacial conductance is probably due to an acoustic
impedance mismatch between the solid and the liquid
phase.\cite{Cahill:793,RevModPhys.61.605} Experiments on the thermal
conductivity of gold nanoparticles and nanorods in solvent and in
glass cages have predicted values for $G$ between 100 and 350
(MW/m$^2$/K), two orders of magnitude larger than the value reported
here.  The experiments typically have multiple surfaces that have been
protected by ionic surfactants, either
citrate\cite{Wilson:2002uq,plech:195423} or cetyltrimethylammonium
bromide (CTAB), or which are in direct contact with various glassy
solids.  In these cases, the acoustic impedance mismatch would be
substantially reduced, leading to much higher interfacial
conductances.  It is also possible, however, that the lack of
electronic effects that gives rise to the low thermal conductivity of
EAM gold is also causing a low reading for this particular interface.

Under this low thermal conductance, both gold and water phases have
sufficient time to eliminate local temperature differences (Fig.
\ref{interface}). With flat thermal profiles within each phase, it is
valid to assume that the temperature difference between gold and water
surfaces would be approximately the same as the difference between the
gold and water bulk regions.  This assumption enables convenient
calculation of $G$ using Eq. \ref{interfaceCalc}.

\begin{figure}
\includegraphics[width=\linewidth]{interface}
\caption{Temperature profiles of the Gold / Water interface at four
  different values for the thermal flux.  Temperatures for slabs
  either in the gold or in the water show no significant differences,
  although there is a large discontinuity between the materials
  because of the relatively low interfacial thermal conductivity.}
\label{interface}
\end{figure}

\begin{table*}
  \begin{minipage}{\linewidth}
    \begin{center}
      
      \caption{Computed interfacial thermal conductivity ($G$) values
        for the Au(111) / water interface at ${\langle T\rangle \sim}$
        300K using a range of energy fluxes. Uncertainties are
        indicated in parentheses. }
      
      \begin{tabular}{|cccc| }
        \hline
        $J_z$ (MW/m$^2$) & $\langle T_{gold} \rangle$ (K) & $\langle
        T_{water} \rangle$ (K) & $G$
        (MW/m$^2$/K)\\
        \hline
        98.0 & 355.2 & 295.8 & 1.65(0.21) \\
        78.8 & 343.8 & 298.0 & 1.72(0.32) \\
        73.6 & 344.3 & 298.0 & 1.59(0.24) \\
        49.2 & 330.1 & 300.4 & 1.65(0.35) \\
        \hline
      \end{tabular}
      \label{interfaceRes}
    \end{center}
  \end{minipage}
\end{table*}


\section{Conclusions}

Our simulations demonstrate that validity of non-isotropic velocity
scaling (NIVS) in RNEMD calculations of thermal conductivity in atomic
and molecular liquids and solids.  NIVS-RNEMD improves the problematic
velocity distributions which can develop in other RNEMD methods.
Furthermore, it provides a means for carrying out non-physical thermal
transfer between different species of molecules, and thus extends its
applicability to interfacial systems. 

We note that the swapping RNEMD approach can also be used for
heterogeneous interfaces by assuming an elastic collision between the
two particles which are exchanging energy.  However, energy transfer
from imaginary elastic collisions is less effective when the two
particles involved have a large difference in mass. In the particular
interface studied in this work, an imarginary elastic collision
between gold atom and water molecule has a kinetic energy transfer
efficiency of only about 30\%.  Under the NIVS-RNEMD approach, as long
as an intersection between the constraint and flux ellipsoids can be
found, the target amount of energy is transferred at every time step.

Although this method samples states from the microcanonical ensemble,
there is no {\it a priori} reason it could not be used in tandem with
extended-system molecular dynamics methods that sample constant
temperature or pressure ensembles.  This has been implemented with the
swapping RNEMD approach,\cite{muller:014102} and since the energy
drift from NIVS-RNEMD is minimal under NVE dynamics (a fairly
stringent test of integration methods), we expect similar good
behavior in simulations in other ensembles

NIVS-RNEMD also has limited applications for shear viscosity
calculations, but has the potential to cause temperature anisotropy
under high momentum fluxes. Further work will be necessary to
eliminate the one-dimensional heating if shear viscosities are
required.

\section{Acknowledgments}
The authors would like to thank Craig Tenney and Ed Maginn for many
helpful discussions.  Support for this project was provided by the
National Science Foundation under grant CHE-0848243. Computational
time was provided by the Center for Research Computing (CRC) at the
University of Notre Dame.  
\newpage

\bibliography{nivsRnemd}

\end{doublespace}
\end{document}

Revision:	3648
Committed:	Wed Sep 22 19:56:09 2010 UTC (13 years, 11 months ago) by gezelter
Content type:	application/x-tex
File size:	43942 byte(s)
Log Message:	final version!