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}
\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

\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 distributions in molecular dynamics
  simulations of heterogeneous systems.  This method extends earlier
  Reverse Non-Equilibrium Molecular Dynamics (RNEMD) methods which use
  momentum exchange swapping moves that 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 and stable
  thermal and momentum gradients can then be established.  The scaling
  method we have developed conserves the total linear momentum and
  total energy of the system.  To test the methods, 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 difficult to measure quantities.

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.

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 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 \left[\vec{v}_i\right]_\alpha \\
P_h^\alpha & = & \sum_{j = 1}^{N_h} m_j \left[\vec{v}_j\right]_\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 momentum 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} K_h^\alpha + K_c^\alpha = \sum_{\alpha = x,y,z}
\left(\alpha^\prime\right)^2 K_h^\alpha + \alpha^2 K_c^\alpha
\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, [CITATIONS NEEDED] but numerically finding the roots
of high-degree polynomials is generally an ill-conditioned
problem.[CITATION NEEDED] One simplification is to maintain velocity
scalings that are {\it as isotropic as possible}.  To do this, we
impose $x=y$, and to 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 as gentle as possible, i.e. ${x,y,z
  \rightarrow 1}$, in order to avoid large perturbation 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
after each MD step.  We have tested the method in a variety of
different systems, including homogeneous fluids (Lennard-Jones and
SPC/E water), crystalline solids ({\sc eam}~\cite{PhysRevB.33.7983} and
quantum Sutton-Chen ({\sc q-sc})~\cite{PhysRevB.59.3527}
models for Gold), and heterogeneous interfaces (QSC gold - SPC/E
water). The last of these systems would have been difficult to study
using previous RNEMD methods, but using velocity scaling moves, we can
even obtain estimates of the interfacial thermal conductivities ($G$).

\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 ``cold'' slab, while the
central slab ($n= N/2 + 1$) was designated the ``hot'' 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.

\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 a the same momentum or energy flux exhibited in the
swapping simulation.

To test the temperature homogeneity (and to compute transport
coefficients), directional momentum and temperature distributions were
accumulated for molecules in each of the slabs.

\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 occurs in two directions due to our periodic
boundary conditions.  The velocity gradient ${\langle \partial v_x
  /\partial z \rangle}$ was obtained using linear regression of the
velocity profiles in the bins.  For Lennard-Jones simulations, shear
viscosities are reporte in reduced units (${\eta^* = \eta \sigma^2
  (\varepsilon m)^{-1/2}}$).

\subsection{Thermal Conductivities}

The energy flux was calculated similarly 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 materials with a relatively low interfacial conductance, and in
cases where the flux between the materials is small, 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.
In calculating the interfacial thermal conductivity $G$, this
assumption was made, and the conductance can be approximated as:

\begin{equation}
G = \frac{E_{total}}{2 t L_x L_y \left( \langle T_{gold}\rangle -
    \langle T_{water}\rangle \right)}
\label{interfaceCalc}
\end{equation}
where ${E_{total}}$ is the imposed non-physical kinetic energy
transfer and ${\langle T_{gold}\rangle}$ and ${\langle
  T_{water}\rangle}$ are the average observed temperature of gold and
water phases respectively.

\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 first run under reduced
temperature ${\langle T^*\rangle = 0.72}$ in the microcanonical
ensemble, but other temperatures ([XXX, YYY, and ZZZ]) were also
sampled.  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. Four different swap intervals were
tested (Table \ref{thermalLJRes}). With a fixed 10 fs [WHY NOT REDUCED
UNITS???]  scaling interval, the target exchange kinetic energy
produced equivalent kinetic energy flux as in the swapping method.
Similar thermal gradients were observed with similar thermal flux
under the two different methods (Figure \ref{thermalGrad}).

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

      \caption{Thermal conductivity (in reduced units) of a
        Lennard-Jones fluid at ${\langle T^* \rangle = 0.72}$ and
        ${\rho^* = 0.85}$ for the swapping and scaling methods at
        various kinetic energy exchange rates. Uncertainties are
        indicated in parentheses.}
      
      \begin{tabular}{|cc|cc|}
        \hline
        \multicolumn{2}{|c|}{Swapping RNEMD} &
        \multicolumn{2}{|c|}{NIVS-RNEMD} \\
        \hline
        Swap Interval (fs) & $\lambda^*_{swap}$ & Equilvalent $J_z^*$ & $\lambda^*_{scale}$\\
        \hline 
        250  & 7.03(0.34) & 0.16  & 7.30(0.10)\\
        500  & 7.03(0.14) & 0.09  & 6.95(0.09)\\
        1000 & 6.91(0.42) & 0.047 & 7.19(0.07)\\
        2000 & 7.52(0.15) & 0.024 & 7.19(0.28)\\
        \hline
      \end{tabular}
      \label{thermalLJRes}
    \end{center}
  \end{minipage}
\end{table*}

\begin{figure}
\includegraphics[width=\linewidth]{thermalGrad}
\caption{NIVS-RNEMD method creates similar temperature gradients
  compared with the swapping method under a variety of imposed kinetic
  energy flux values.}
\label{thermalGrad}
\end{figure}

During these simulations, velocities were recorded every 1000 steps
and was used to produce distributions of both velocity and speed in
each of the slabs. From these distributions, we observed that under
relatively high non-physical kinetic energy flux, the spee of
molecules in slabs where 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 MB distribution.  In
the cold slab, the opposite notching and shouldering occurs.  This
phenomenon is more obvious at higher swapping rates.  

In the velocity distributions, the ideal Gaussian peak 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 diffusions of low (high) speed particles from neighboring slabs,
unless the swapping rate is sufficiently small. Simutaneously, surplus
low speed particles in the low temperature 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 (Figure \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 an unphysical thermal flux
as the previous RNEMD methods but without large perturbations to the
velocity distributions in the two slabs.

\begin{figure}
\includegraphics[width=\linewidth]{thermalHist}
\caption{Speed distribution for thermal conductivity using a)
  ``swapping'' and b) NIVS- RNEMD methods. Shown is from the
  simulations with an exchange or equilvalent exchange interval of 250
  fs. In circled areas, distributions from ``swapping'' RNEMD
  simulation have deviation from ideal Maxwell-Boltzmann distribution
  (curves fit for each distribution).}
\label{thermalHist}
\end{figure}


\subsubsection*{Shear Viscosity}
Our calculations (Table \ref{shearRate}) show that velocity-scaling
RNEMD predicted comparable shear viscosities to swap RNEMD method.  All
the scale method results were from simulations that had a scaling
interval of 10 time steps. The average molecular momentum gradients of
these samples are shown in Figure \ref{shear} (a) and (b).

\begin{table*}
  \begin{minipage}{\linewidth}
    \begin{center}
      
      \caption{Shear viscosities of Lennard-Jones fluid at ${T^* =
          0.72}$ and ${\rho^* = 0.85}$ using swapping and NIVS methods
        at various momentum exchange rates.  Uncertainties are
        indicated in parentheses. }
      
      \begin{tabular}{ccccc}
        Swapping method & & & NIVS-RNEMD & \\
        \hline
        Swap Interval (fs) & $\eta^*_{swap}$ & & Equilvalent $j_p^*(v_x)$ &
        $\eta^*_{scale}$\\
        \hline
        500  & 3.64(0.05) & & 0.09  & 3.76(0.09)\\
        1000 & 3.52(0.16) & & 0.046 & 3.66(0.06)\\
        2000 & 3.72(0.05) & & 0.024 & 3.32(0.18)\\
        2500 & 3.42(0.06) & & 0.019 & 3.43(0.08)\\
        \hline
      \end{tabular}
      \label{shearRate}
    \end{center}
  \end{minipage}
\end{table*}

\begin{figure}
  \includegraphics[width=\linewidth]{shear}
  \caption{Average momentum gradients in shear viscosity simulations,
    using (a) ``swapping'' method and (b) NIVS-RNEMD method
    respectively. (c) Temperature difference among x and y, z dimensions
    observed when using NIVS-RNEMD with equivalent exchange interval of
    500 fs.}
  \label{shear}
\end{figure}

However, observations of temperatures along three dimensions show that
inhomogeneity occurs in scaling RNEMD simulations, particularly in the
two slabs which were scaled. Figure \ref{shear} (c) indicate that with
relatively large imposed momentum flux, the temperature difference among $x$
and the other two dimensions was significant. This would result from the
algorithm of scaling method. From Eq. \ref{eq:fluxPlane}, after
momentum gradient is set up, $P_c^x$ would be roughly stable
($<0$). Consequently, scaling factor $x$ would most probably larger
than 1. Therefore, the kinetic energy in $x$-dimension $K_c^x$ would
keep increase after most scaling steps. And if there is not enough time
for the kinetic energy to exchange among different dimensions and
different slabs, the system would finally build up temperature
(kinetic energy) difference among the three dimensions. Also, between
$y$ and $z$ dimensions in the scaled slabs, temperatures of $z$-axis
are closer to neighbor slabs. This is due to momentum transfer along
$z$ dimension between slabs.

Although results between scaling and swapping methods are comparable,
the inherent temperature inhomogeneity even in relatively low imposed
exchange momentum flux simulations makes scaling RNEMD method less
attractive than swapping RNEMD in shear viscosity calculation. 


\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 molecule
center-of-mass velocities instead of single atom velocities in a
molecule 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 the rigid bodies of SPC/E model 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
[XXX ps] period under constant-NVE conditions without any momentum
flux.  [YYY ps] was allowed to stabilize the system with an imposed
momentum transfer to create a gradient, and [ZZZ ps] was alotted for
data collection under RNEMD.

As shown in Figure \ref{spceGrad}, temperature gradients were
established similar to the previous work.  However, the average
temperature of our system is 300K, while that in Bedrov {\it et al.}
is 318K, which would be attributed for part of the difference between
the final calculation results (Table \ref{spceThermal}). [WHY DIDN'T
WE DO 318 K?]  Both methods yield values in reasonable agreement with
experiment [DONE AT WHAT TEMPERATURE?]

\begin{figure}
  \includegraphics[width=\linewidth]{spceGrad}
  \caption{Temperature gradients in SPC/E water thermal conductivity
    simulations.}
  \label{spceGrad}
\end{figure}

\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|ccc|}
        \hline
        $\langle dT/dz\rangle$(K/\AA) & \multicolumn{3}{|c|}{$\lambda
          (\mathrm{W m}^{-1} \mathrm{K}^{-1})$} \\
        & This work (300K) & Previous simulations (318K)\cite{Bedrov:2000} &
        Experiment\cite{WagnerKruse}\\
        \hline
        0.38 & 0.816(0.044) & & 0.64\\
        0.81 & 0.770(0.008) & 0.784 & \\
        1.54 & 0.813(0.007) & 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 ({\sc eam})~\cite{PhysRevB.33.7983} and Sutton-Chen ({\sc
  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 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 ({\sc q-sc}) formulation matches these properties while
including zero-point quantum corrections for different transition
metals.\cite{PhysRevB.59.3527}  The {\sc eam} functional forms differ
slightly from {\sc sc} but the overall method is very similar.

In this work, we have utilized both the {\sc eam} and the {\sc q-sc}
potentials to test the behavior of scaling RNEMD. 

A face-centered-cubic (FCC) lattice was prepared containing 2880 Au
atoms.  [LxMxN 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 OR LINEAR DIMENSIONS ?].  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 transfer 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 in thermal gradient linear regession. {\sc eam} predicts
slightly larger thermal conductivities than {\sc q-sc}.  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 effects in
these force fields.\cite{Clancy:1992} The non-equilibrium MD method
employed in their simulations produced comparable results to ours.  It
should be noted that the density of the metal being simulated also
greatly affects the thermal conductivity.  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 [PAGE
NUMBERS?].\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. [CLANCY COMPARISON? SWAPPING
        COMPARISON?]}
      
      \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}{*}{\sc q-sc} & \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}{*}{\sc 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 withing 3 \AA radius of any gold atom.  The final
configuration contained 1862 SPC/E water molecules.

After simulations of bulk water and crystal gold, a mixture system was
constructed, consisting of 1188 Au atoms and 1862 H$_2$O
molecules. Spohr potential was adopted in depicting the interaction
between metal atom and water molecule.\cite{ISI:000167766600035} A
similar protocol of equilibration was followed. Several thermal
gradients was built under different target thermal flux. It was found
out that compared to our previous simulation systems, the two phases
could have large temperature difference even under a relatively low
thermal flux.


After simulations of homogeneous water and gold systems using
NIVS-RNEMD method were proved valid, calculation of gold/water
interfacial thermal conductivity was followed. It is found out that
the low interfacial conductance is probably due to the hydrophobic
surface in our system. Figure \ref{interface} (a) demonstrates mass
density change along $z$-axis, which is perpendicular to the
gold/water interface. It is observed that water density significantly
decreases when approaching the surface. Under this low thermal
conductance, both gold and water phase have sufficient time to
eliminate temperature difference inside respectively (Figure
\ref{interface} b). With indistinguishable temperature difference
within respective phase, it is valid to assume that the temperature
difference between gold and water on surface would be approximately
the same as the difference between the gold and water phase. This
assumption enables convenient calculation of $G$ using
Eq. \ref{interfaceCalc} instead of measuring temperatures of thin
layer of water and gold close enough to surface, which would have
greater fluctuation and lower accuracy. Reported results (Table
\ref{interfaceRes}) are of two orders of magnitude smaller than our
calculations on homogeneous systems, and thus have larger relative
errors than our calculation results on homogeneous systems.

\begin{figure}
\includegraphics[width=\linewidth]{interface}
\caption{Simulation results for Gold/Water interfacial thermal
  conductivity: (a) Significant water density decrease is observed on
  crystalline gold surface, which indicates low surface contact and
  leads to low thermal conductance. (b) Temperature profiles for a
  series of simulations. Temperatures of different slabs in the same
  phase show no significant differences.}
\label{interface}
\end{figure}

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

\caption{Calculation results for interfacial thermal conductivity
  at ${\langle T\rangle \sim}$ 300K at various thermal exchange
  rates. Errors of calculations in parentheses. }

\begin{tabular}{cccc}
\hline
$J_z$ (MW/m$^2$) & $T_{gold}$ (K) & $T_{water}$ (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}
NIVS-RNEMD simulation method is developed and tested on various
systems. Simulation results demonstrate its validity in thermal
conductivity calculations, from Lennard-Jones fluid to multi-atom
molecule like water and metal crystals. NIVS-RNEMD improves
non-Boltzmann-Maxwell distributions, which exist in previous RNEMD
methods. Furthermore, it develops a valid means for unphysical thermal
transfer between different species of molecules, and thus extends its
applicability to interfacial systems. Our calculation of gold/water
interfacial thermal conductivity demonstrates this advantage over
previous RNEMD methods. NIVS-RNEMD has also limited application on
shear viscosity calculations, but could cause temperature difference
among different dimensions under high momentum flux. Modification is
necessary to extend the applicability of NIVS-RNEMD in shear viscosity
calculations.

\section{Acknowledgments}
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

\bibliographystyle{aip}
\bibliography{nivsRnemd}

\end{doublespace}
\end{document}

Revision:	3611
Committed:	Wed Jul 14 15:20:31 2010 UTC (14 years, 1 month ago) by gezelter
Content type:	application/x-tex
File size:	40620 byte(s)
Log Message:	edits