trunk/nivsRnemd/nivsRnemd.tex

\documentclass[11pt]{article}
\usepackage{amsmath}
\usepackage{amssymb}
\usepackage{setspace}
\usepackage{endfloat}
\usepackage{caption}
%\usepackage{tabularx}
\usepackage{graphicx}
%\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 some
  earlier Reverse Non-Equilibrium Molecular Dynamics (RNEMD) methods
  which use momentum exchange swapping moves that can create
  non-thermal velocity distributions (and which 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 monatomic liquids to 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 momentum gradients
or strains (the response).  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-RNEMD) 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 single component, multi-component, and
non-isotropic mixtures and show that it is capable of providing
reasonable estimates of the thermal conductivity and shear viscosity
in 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 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} a_\alpha \alpha^2 + b_\alpha \alpha + c_\alpha = 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 can be applied while preserving
both linear momentum in all three directions as well as total kinetic
energy.

The goal of using velocity scaling variables is to transfer linear
momentum or 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}
x^2 K_c^x + y^2 K_c^y + z^2 K_c^z = K_c^x + K_c^y + K_c^z - J_z\Delta t
\label{eq:fluxEllipsoid}
\end{equation}
The spatial extent of the {\it thermal flux ellipsoid} is governed
both by a targetted value, $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}$ which 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{Illustration from the perspective of a space having cold
  slab scaling coefficients as its coordinates. Scaling points 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}

One may also define {\it momentum} flux (say along the $x$-direction) as:
\begin{equation}
(1-x) P_c^x = j_z(p_x)\Delta t
\label{eq:fluxPlane}
\end{equation}
The above {\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 both the momentum and energy flux scenarios, valid scaling
parameters are arrived at by solving geometric intersection problems
in $x, y, z$ space in order to obtain cold slab scaling parameters.
Once the scaling variables for the cold slab are known, the hot slab
scaling has also been determined.


The following problem will be choosing an optimal set of scaling
variables among the possible sets. Although this method is inherently
non-isotropic, the goal is still to maintain the system as isotropic
as possible. Under this consideration, one would like the kinetic
energies in different directions could become as close as each other
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. Therefore, one approach to obtain
the scaling variables would be constructing an criteria function, with
constraints as above equation sets, and solving the function's minimum
by method like Lagrange multipliers.

In order to save computation time, we have a different approach to a
relatively good set of scaling variables with much less calculation
than above. Here is the detail of our simplification of the problem.

In the case of kinetic energy transfer, we impose another constraint
${x = y}$, into the equation sets. Consequently, there are two
variables left. And now one only needs to solve a set of two {\it
  ellipses equations}. This problem would be transformed into solving
one quartic equation for one of the two variables. There are known
generic methods that solve real roots of quartic equations. Then one
can determine the other variable and obtain sets of scaling
variables. Among these sets, one can apply the above criteria to
choose the best set, while much faster with only a few sets to choose.

In the case of momentum flux transfer, we impose another constraint to
set the kinetic energy transfer as zero. In another word, we apply
Eq. \ref{eq:fluxEllipsoid} and let ${J_z = 0}$. After that, with one
variable fixed by Eq. \ref{eq:fluxPlane}, one now have a similar set
of equations on the above kinetic energy transfer problem. Therefore,
an approach similar to the above would be sufficient for this as well.

\section{Computational Details}
\subsection{Lennard-Jones Fluid}
Our simulation consists of a series of systems. All of these
simulations were run with the OpenMD simulation software
package\cite{Meineke:2005gd} integrated with RNEMD codes.

A Lennard-Jones fluid system was built and tested first. In order to
compare our method with swapping RNEMD, a series of simulations were
performed to calculate the shear viscosity and thermal conductivity of
argon. 2592 atoms were in a orthorhombic cell, which was ${10.06\sigma
  \times 10.06\sigma \times 30.18\sigma}$ by size. The reduced density
${\rho^* = \rho\sigma^3}$ was thus 0.85, which enabled direct
comparison between our results and others. These simulations used
velocity Verlet algorithm with reduced timestep ${\tau^* =
  4.6\times10^{-4}}$. 

For shear viscosity calculation, the reduced temperature was ${T^* =
  k_B T/\varepsilon = 0.72}$. Simulations were first equilibrated in canonical
ensemble (NVT), then equilibrated in microcanonical ensemble
(NVE). Establishing and stablizing momentum gradient were followed
also in NVE ensemble. For the swapping part, Muller-Plathe's algorithm was
adopted.\cite{ISI:000080382700030} The simulation box was under
periodic boundary condition, and devided into ${N = 20}$ slabs. In each swap,
the top slab ${(n = 1)}$ exchange the most negative $x$ momentum with the
most positive $x$ momentum in the center slab ${(n = N/2 + 1)}$. Referred
to Tenney {\it et al.}\cite{ISI:000273472300004}, a series of swapping
frequency were chosen. According to each result from swapping
RNEMD, scaling RNEMD simulations were run with the target momentum
flux set to produce a similar momentum flux, and consequently shear
rate. Furthermore, various scaling frequencies can be tested for one
single swapping rate. To test the temperature homogeneity in our
system of swapping and scaling methods, temperatures of different
dimensions in all the slabs were observed. Most of the simulations
include $10^5$ steps of equilibration without imposing momentum flux,
$10^5$ steps of stablization with imposing unphysical momentum
transfer, and $10^6$ steps of data collection under RNEMD. For
relatively high momentum flux simulations, ${5\times10^5}$ step data
collection is sufficient. For some low momentum flux simulations,
${2\times10^6}$ steps were necessary. 

After each simulation, the shear viscosity was calculated in reduced
unit. The momentum flux was calculated with total unphysical
transferred momentum ${P_x}$ and 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$ denotes $x$ and $y$ lengths of the simulation
box, and physical momentum transfer occurs in two ways due to our
periodic boundary condition settings. And the velocity gradient
${\langle \partial v_x /\partial z \rangle}$ can be obtained by a
linear regression of the velocity profile. From the shear viscosity
$\eta$ calculated with the above parameters, one can further convert
it into reduced unit ${\eta^* = \eta \sigma^2 (\varepsilon m)^{-1/2}}$.

For thermal conductivity calculations, simulations were first run under
reduced temperature ${\langle T^*\rangle = 0.72}$ in NVE
ensemble. Muller-Plathe's algorithm was adopted in the swapping
method. Under identical simulation box parameters with our shear
viscosity calculations, in each swap, the top slab exchanges all three
translational momentum components of the molecule with least kinetic
energy with the same components of the molecule in the center slab
with most kinetic energy, unless this ``coldest'' molecule in the
``hot'' slab is still ``hotter'' than the ``hottest'' molecule in the
``cold'' slab. According to swapping RNEMD results, target energy flux
for scaling RNEMD simulations can be set. Also, various scaling
frequencies can be tested for one target energy flux. To compare the
performance between swapping and scaling method, distributions of
velocity and speed in different slabs were observed.

For each swapping rate, thermal conductivity was calculated in reduced
unit. The energy flux was calculated similarly to the momentum flux,
with total unphysical transferred energy ${E_{total}}$ and data collection
time $t$:
\begin{equation}
J_z = \frac{E_{total}}{2 t L_x L_y}
\end{equation}
And the temperature gradient ${\langle\partial T/\partial z\rangle}$
can be obtained by a linear regression of the temperature
profile. From the thermal conductivity $\lambda$ calculated, one can
further convert it into reduced unit ${\lambda^*=\lambda \sigma^2
  m^{1/2} k_B^{-1}\varepsilon^{-1/2}}$.

\subsection{ Water / Metal Thermal Conductivity}
Another series of our simulation is the calculation of interfacial
thermal conductivity of a Au/H$_2$O system. Respective calculations of
liquid water (Extended Simple Point Charge model) and crystal gold
thermal conductivity were performed and compared with current results
to ensure the validity of NIVS-RNEMD. After that, a mixture system was
simulated.

For thermal conductivity calculation of bulk water, a simulation box
consisting of 1000 molecules were first equilibrated under ambient
pressure and temperature conditions using NPT ensemble, followed by
equilibration in fixed volume (NVT). The system was then equilibrated in
microcanonical ensemble (NVE). Also in NVE ensemble, establishing a
stable thermal gradient was followed. The simulation box was under
periodic boundary condition and devided into 10 slabs. Data collection
process was similar to Lennard-Jones fluid system.

Thermal conductivity calculation of bulk crystal gold used a similar
protocol. Two types of force field parameters, Embedded Atom Method
(EAM) and Quantum Sutten-Chen (QSC) force field were used
respectively. The face-centered cubic crystal simulation box consists of
2880 Au atoms. The lattice was first allowed volume change to relax
under ambient temperature and pressure. Equilibrations in canonical and
microcanonical ensemble were followed in order. With the simulation
lattice devided evenly into 10 slabs, different thermal gradients were
established by applying a set of target thermal transfer flux. Data of
the series of thermal gradients was collected for calculation.

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. Therefore, under our low flux conditions, it is assumed
that the metal and water phases have respectively homogeneous
temperature, excluding the surface regions. In calculating the
interfacial thermal conductivity $G$, this assumptioin was applied and
thus our formula becomes:

\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 unphysical 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 And Discussions}
\subsection{Thermal Conductivity}
\subsubsection{Lennard-Jones Fluid}
Our thermal conductivity calculations show that scaling method results
agree with swapping method. Four different exchange intervals were
tested (Table \ref{thermalLJRes}) using swapping method. With a fixed
10fs exchange interval, target exchange kinetic energy was set to
produce equivalent kinetic energy flux as in swapping method. And
similar thermal gradients were observed with similar thermal flux in
two simulation methods (Figure \ref{thermalGrad}).

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

\caption{Calculation results for thermal conductivity of Lennard-Jones
  fluid at ${\langle T^* \rangle = 0.72}$ and ${\rho^* = 0.85}$, with
  swap and scale methods at various kinetic energy exchange rates. Results
  in reduced unit. Errors of calculations in parentheses.} 

\begin{tabular}{ccc}
\hline
(Equilvalent) Exchange Interval (fs) & $\lambda^*_{swap}$ &
$\lambda^*_{scale}$\\
\hline 
250  & 7.03(0.34) & 7.30(0.10)\\
500  & 7.03(0.14) & 6.95(0.09)\\
1000 & 6.91(0.42) & 7.19(0.07)\\
2000 & 7.52(0.15) & 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 introduced similar temperature gradients
  compared to ``swapping'' method under various kinetic energy flux in
  thermal conductivity simulations.}
\label{thermalGrad}
\end{figure}

During these simulations, molecule velocities were recorded in 1000 of
all the snapshots of one single data collection process. These
velocity data were used to produce histograms of velocity and speed
distribution in different slabs. From these histograms, it is observed
that under relatively high unphysical kinetic energy flux, speed and
velocity distribution of molecules in slabs where swapping occured
could deviate from Maxwell-Boltzmann distribution. Figure
\ref{thermalHist} a) illustrates how these distributions deviate from an
ideal distribution. In high temperature slab, probability density in
low speed is confidently smaller than ideal curve fit; in low
temperature slab, probability density in high speed is smaller than
ideal, while larger than ideal in low speed. This phenomenon is more
obvious in our high swapping rate simulations. And this deviation
could also leads to deviation of distribution of velocity in various
dimensions. One feature of these deviated distribution is that in high
temperature slab, the ideal Gaussian peak was changed into a
relatively flat plateau; while in low temperature slab, that peak
appears sharper. This problem is rooted in the mechanism of the
swapping method. Continually depleting low (high) speed particles in
the high (low) temperature slab could not be complemented by
diffusions of low (high) speed particles from neighbor slabs, unless
in suffciently low swapping rate. Simutaneously, surplus low speed
particles in the low temperature slab do not have sufficient time to
diffuse to neighbor slabs. However, thermal exchange rate should reach
a minimum level to produce an observable thermal gradient under noise
interference. Consequently, swapping RNEMD has a relatively narrow
choice of swapping rate to satisfy these above restrictions.

Comparatively, NIVS-RNEMD has a speed distribution closer to the ideal
curve fit (Figure \ref{thermalHist} b). Essentially, after scaling, a
Gaussian distribution function would remain Gaussian. Although a
single scaling is non-isotropic in all three dimensions, our scaling
coefficient criteria could help maintian the scaling region as
isotropic as possible. On the other hand, scaling coefficients are
preferred to be as close to 1 as possible, which also helps minimize
the difference among different dimensions. This is possible if scaling
interval and one-time thermal transfer energy are well
chosen. Consequently, NIVS-RNEMD is able to impose an unphysical
thermal flux as the previous RNEMD method without large perturbation
to the distribution of velocity and speed in the exchange regions.

\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{SPC/E Water}
Our results of SPC/E water thermal conductivity are comparable to
Bedrov {\it et al.}\cite{ISI:000090151400044}, which employed the
previous swapping RNEMD method for their calculation. Bedrov {\it et
  al.}\cite{ISI:000090151400044} 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 adopted in our simulations. Scaling is applied to the
velocities of the rigid bodies of SPC/E model water molecules, instead
of each hydrogen and oxygen atoms in relevant water molecules. As
shown in Figure \ref{spceGrad}, temperature gradients were established
similar to their system. However, the average temperature of our
system is 300K, while theirs is 318K, which would be attributed for
part of the difference between the final calculation results (Table
\ref{spceThermal}). Both methods yields values in agreement with
experiment. And this shows the applicability of our method to
multi-atom molecular system.

\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{Calculation results for thermal conductivity of SPC/E water
  at ${\langle T\rangle}$ = 300K at various thermal exchange rates. Errors of
  calculations in parentheses. }

\begin{tabular}{cccc}
\hline
$\langle dT/dz\rangle$(K/\AA) & & $\lambda$(W/m/K) & \\
& This work & Previous simulations\cite{ISI:000090151400044} &
Experiment$^a$\\
\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*}

\subsubsection{Crystal Gold}
Our results of gold thermal conductivity using two force fields are
shown separately in Table \ref{qscThermal} and \ref{eamThermal}. In
these calculations,the end and middle slabs were excluded in thermal
gradient regession and only used as heat source and drain in the
systems. Our yielded values using EAM force field are slightly larger
than those using QSC force field. However, both series are
significantly smaller than experimental value by an order of more than
100. It has been verified that this difference is mainly attributed to
the lack of electron interaction representation in these force field
parameters. Richardson {\it et al.}\cite{Clancy:1992} used EAM
force field parameters in their metal thermal conductivity
calculations. The Non-Equilibrium MD method they employed in their
simulations produced comparable results to ours. As Zhang {\it et
  al.}\cite{ISI:000231042800044} stated, thermal conductivity values
are influenced mainly by force field. Therefore, it is confident to
conclude that NIVS-RNEMD is applicable to metal force field system.

\begin{figure}
\includegraphics[width=\linewidth]{AuGrad}
\caption{Temperature gradients for thermal conductivity calculation of
  crystal gold using QSC force field.}
\label{AuGrad}
\end{figure}

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

\caption{Calculation results for thermal conductivity of crystal gold
  using QSC force field at ${\langle T\rangle}$ = 300K at various
  thermal exchange rates. Errors of calculations in parentheses. }

\begin{tabular}{cc}
\hline
$\langle dT/dz\rangle$(K/\AA) & $\lambda$(W/m/K)\\
\hline
1.44 & 1.10(0.01)\\
2.86 & 1.08(0.02)\\
5.14 & 1.15(0.01)\\
\hline
\end{tabular}
\label{qscThermal}
\end{center}
\end{minipage}
\end{table*}

\begin{figure}
\includegraphics[width=\linewidth]{eamGrad}
\caption{Temperature gradients for thermal conductivity calculation of
  crystal gold using EAM force field.}
\label{eamGrad}
\end{figure}

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

\caption{Calculation results for thermal conductivity of crystal gold
  using EAM force field at ${\langle T\rangle}$ = 300K at various
  thermal exchange rates. Errors of calculations in parentheses. }

\begin{tabular}{cc}
\hline
$\langle dT/dz\rangle$(K/\AA) & $\lambda$(W/m/K)\\
\hline
1.24 & 1.24(0.06)\\
2.06 & 1.37(0.04)\\
2.55 & 1.41(0.03)\\
\hline
\end{tabular}
\label{eamThermal}
\end{center}
\end{minipage}
\end{table*}


\subsection{Interfaciel Thermal Conductivity}
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{interfaceDensity} 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{interfaceGrad}). 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]{interfaceDensity}
\caption{Density profile for interfacial thermal conductivity
  simulation box. Significant water density decrease is observed on
  gold surface.}
\label{interfaceDensity}
\end{figure}

\begin{figure}
\includegraphics[width=\linewidth]{interfaceGrad}
\caption{Temperature profiles for interfacial thermal conductivity
  simulation box. Temperatures of different slabs in the same phase
  show no significant difference.}
\label{interfaceGrad}
\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}$ & $T_{water}$ & $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*}

\subsection{Shear Viscosity}
Our calculations (Table \ref{shearRate}) shows that scale RNEMD method
produced comparable shear viscosity to swap RNEMD method. In Table
\ref{shearRate}, the names of the calculated samples are devided into
two parts. The first number refers to total slabs in one simulation
box. The second number refers to the swapping interval in swap method, or
in scale method the equilvalent swapping interval that the same
momentum flux would theoretically result in swap 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{Calculation results for shear viscosity of Lennard-Jones
  fluid at ${T^* = 0.72}$ and ${\rho^* = 0.85}$, with swap and scale
  methods at various momentum exchange rates. Results in reduced
  unit. Errors of calculations in parentheses. }

\begin{tabular}{ccc}
\hline
Series & $\eta^*_{swap}$ & $\eta^*_{scale}$\\
\hline
20-500 & 3.64(0.05) & 3.76(0.09)\\
20-1000 & 3.52(0.16) & 3.66(0.06)\\
20-2000 & 3.72(0.05) & 3.32(0.18)\\
20-2500 & 3.42(0.06) & 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. 

\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:	3593
Committed:	Fri Apr 16 22:08:12 2010 UTC (14 years, 2 months ago) by skuang
Content type:	application/x-tex
File size:	37551 byte(s)
Log Message:	minor edit to some captions.