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}

\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 transport
coefficients (shear viscosity and thermal conductivity) 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{ISI:000246190100032}

\begin{figure}
\includegraphics[width=\linewidth]{thermalDemo}
\caption{Demostration of thermal gradient estalished by RNEMD method.}
\label{thermalDemo}
\end{figure}

RNEMD is preferable in many ways to the forward NEMD methods 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 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{ISI:000273472300004} 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
develop 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 Muller-Plathe's RNEMD method; the periodic
system is partitioned into a series of thin slabs along a particular
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 bin 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 equation}:
\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 equation defines the set of cold slab scaling
parameters which can be applied while preserving both linear momentum
in all three directions as well as 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 equation} 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,
it is sufficient to 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{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 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 equation} is essentially a plane which 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.

Similarly, to satisfy a momentum flux as well as the conservation
constraints, it is sufficient to 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.

To summarize, by solving respective equation sets, one can determine
possible sets of scaling variables for cold slab. And corresponding
sets of scaling variables for hot slab can be determine as well.

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}
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 methods.

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 shear
rate. Furthermore, various scaling frequencies can be tested for one
single swapping rate. To compare the performance between swapping and
scaling method, 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 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}
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 calculation, simulations were first run under
reduced temperature ${T^* = 0.72}$ in NVE ensemble. Muller-Plathe's
algorithm was adopted in the swapping method. Under identical
simulation box parameters, in each swap, the top slab exchange the
molecule with least kinetic energy with 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}}$.

Another series of our simulation is to calculate the interfacial
thermal conductivity of a Au/H$_2$O system. Respective calculations of
water (SPC/E) and gold (QSC) 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 (NPT), followed by equilibration
in fixed volume (NVT). The system was then equilibrated in
microcanonical ensemble (NVE). Also in NVE ensemble, establishing
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. And the face centered cubic crystal simulation box consists
of 2880 Au atoms.

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. A thermal gradient was
built. It was found out that compared to homogeneous systems, the two
phases could have large temperature difference under a relatively low
thermal flux. Therefore, under our low flux condition, it is assumed
that the metal and water phases have respectively homogeneous
temperature. 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 Discussion}
\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{shearGrad}.

\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]{shearGrad}
\caption{Average momentum gradients of shear viscosity simulations}
\label{shearGrad}
\end{figure}

\begin{figure}
\includegraphics[width=\linewidth]{shearTempScale}
\caption{Temperature profile for scaling RNEMD simulation.}
\label{shearTempScale}
\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{shearTempScale} 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{Thermal Conductivity}
\subsubsection{Lennard-Jones Fluid}
Our thermal conductivity calculations also show that scaling method results
agree with swapping method. Table \ref{thermal} lists our simulation
results with similar manner we used in shear viscosity
calculation. All the data reported from scaling method were obtained
by simulations of 10-step exchange frequency, and the target exchange
kinetic energy were set to produce equivalent kinetic energy flux as
in swapping method. Figure \ref{thermalGrad} exhibits similar thermal
gradients of respective similar kinetic energy flux.

\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
Series & $\lambda^*_{swap}$ & $\lambda^*_{scale}$\\
\hline 
20-250  & 7.03(0.34) & 7.30(0.10)\\
20-500  & 7.03(0.14) & 6.95(0.09)\\
20-1000 & 6.91(0.42) & 7.19(0.07)\\
20-2000 & 7.52(0.15) & 7.19(0.28)\\
\hline
\end{tabular}
\label{thermal}
\end{center}
\end{minipage}
\end{table*}

\begin{figure}
\includegraphics[width=\linewidth]{thermalGrad}
\caption{Temperature gradients of thermal conductivity simulations}
\label{thermalGrad}
\end{figure}

During these simulations, molecule velocities were recorded in 1000 of
all the snapshots. These velocity data were used to produce histograms
of velocity and speed distribution in different slabs. From these
histograms, it is observed that with increasing unphysical kinetic
energy flux, speed and velocity distribution of molecules in slabs
where swapping occured could deviate from Maxwell-Boltzmann
distribution. Figure \ref{histSwap} indicates how these distributions
deviate from ideal condition. In high temperature slabs, probability
density in low speed is confidently smaller than ideal distribution;
in low temperature slabs, probability density in high speed is smaller
than ideal. This phenomenon is observable even in our relatively low
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.

\begin{figure}
\includegraphics[width=\linewidth]{histSwap}
\caption{Speed distribution for thermal conductivity using swapping RNEMD.}
\label{histSwap}
\end{figure}

\begin{figure}
\includegraphics[width=\linewidth]{histScale}
\caption{Speed distribution for thermal conductivity using scaling RNEMD.}
\label{histScale}
\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. Our simulations
were able to produce a similar temperature gradient 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 two series of results. 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 for SPC/E water thermal conductivity.}
\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 used QSC force field are
shown in Table \ref{AuThermal}. Although our calculation is smaller
than experimental value by an order of more than 100, this difference
is mainly attributed to the lack of electron interaction
representation in our force field parameters. Richardson {\it et
  al.}\cite{ISI:A1992HX37800010} used similar force field parameters
in their metal thermal conductivity calculations. The EMD method they
employed in their simulations produced comparable results to
ours. 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 crystal gold thermal conductivity.}
\label{AuGrad}
\end{figure}

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

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

\begin{tabular}{ccc}
\hline
$\langle dT/dz\rangle$(K/\AA) & $\lambda$(W/m/K)\\
& This work & Previous simulations\cite{ISI:A1992HX37800010} \\
\hline
1.44 & 1.10(0.01) & \\
2.86 & 1.08(0.02) & \\
5.14 & 1.15(0.01) & \\
\hline
\end{tabular}
\label{AuThermal}
\end{center}
\end{minipage}
\end{table*}

\subsection{Interfaciel Thermal Conductivity}
After valid simulations of homogeneous water and gold systems using
NIVS-RNEMD method, calculation of gold/water interfacial thermal
conductivity was followed. It is found out that the interfacial
conductance is low due to a hydrophobic surface in our system. Figure
\ref{interfaceDensity} demonstrates this observance. Consequently, our
reported results (Table \ref{interfaceRes}) are of two orders of
magnitude smaller than our calculations on homogeneous systems.

\begin{figure}
\includegraphics[width=\linewidth]{interfaceDensity}
\caption{Density profile for interfacial thermal conductivity
  simulation box.}
\label{interfaceDensity}
\end{figure}

\begin{figure}
\includegraphics[width=\linewidth]{interfaceGrad}
\caption{Temperature profiles for interfacial thermal conductivity
  simulation box.}
\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*}

\section{Conclusions}
NIVS-RNEMD simulation method is developed and tested on various
systems. Simulation results demonstrate its validity of thermal
conductivity calculations. NIVS-RNEMD improves non-Boltzmann-Maxwell
distributions existing in previous RNEMD methods, and extends its
applicability to interfacial systems. NIVS-RNEMD has also limited
application on shear viscosity calculations, but under high momentum
flux, it  could cause temperature difference among different
dimensions. 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{jcp2}
\bibliography{nivsRnemd}
\end{doublespace}
\end{document}

Revision:	3575
Committed:	Sat Mar 13 02:08:11 2010 UTC (15 years, 6 months ago) by skuang
Content type:	application/x-tex
File size:	29920 byte(s)
Log Message:	a few changes in intro and method part.