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 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
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}
\sum_{\alpha = x,y,z} K_c^\alpha \alpha^2 = -J_z \Delta t +
\sum_{\alpha = x,y,z} K_c^\alpha  
\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 way around this is to try 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 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}
Rather than using this method to induce a thermal flux, it is possible
to use 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 another
word, 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,\cite{Meineke:2005gd} by performing the NIVS scaling moves after
each MD step.  We have tested it for a variety of different
situations, including homogeneous fluids (Lennard-Jones and SPC/E
water), crystalline solids (EAM and Sutton-Chen models for Gold), and
heterogeneous interfaces (EAM gold - SPC/E water). The last of these
systems would have been very difficult to study using previous RNEMD
methods, but using velocity scaling moves, we can even obtain
estimates of the interfacial thermal conductivity.

\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 calculation, the mean
temperature was ${T^* = k_B T/\varepsilon = 0.72}$.  Simulations were
first thermalized in canonical ensemble (NVT), then equilibrated in
microcanonical ensemble (NVE) before introducing any non-equilibrium
method.

We have compared the momentum gradients established using our method
to those obtained using the original M\"{u}ller-Plathe swapping
algorithm.\cite{ISI:000080382700030} In both cases, the simulation box
was divided into ${N = 20}$ slabs.  In the swapping algorithm, the top
slab $(n = 1)$ exchanges the most negative $x$ momentum with the most
positive $x$ momentum in the center slab $(n = N/2 + 1)$.  The rate at
which the swapping moves are carried out defines the momentum or
thermal flux between the two slabs.  In their work, Tenney {\it et
  al.}\cite{Maginn:2010} found problematic behavior with large swap
frequencies.

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}{ccccc}
Swapping method & & & 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 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{Bedrov:2000}, which employed the
previous swapping RNEMD method for their calculation. 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 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{Bedrov:2000} &
Experiment[CITATION]\\
\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 in Table \ref{AuThermal}. 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 a factor of more than 200. 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. Another factor that affects the calculation results could be
the density of the metal. After equilibration under
isobaric-isothermal conditions, our crystall simulation cell expanded
by the order of 1\%. Under larger lattice constant than default,
lower thermal conductance were expected and observed. Furthermore,
from our simulations, larger thermal gradients tend to produce larger
results. Consequently, having a higher temperature gradient (3.6
K/\AA), Richardson {\it et al.} could obtain a larger
result. Therefore, it is still confident to apply NIVS-RNEMD to metal
force field systems.

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

\caption{Calculation results for thermal conductivity of crystal gold
  using different force fields at different temperatures and  various
  thermal exchange rates. Errors of calculations in parentheses.}

\begin{tabular}{ccccc}
\hline
Force Field Used & $\rho$(g/cm$^3$) & ${\langle T\rangle}$/K &
$\langle dT/dz\rangle$(K/\AA) & $\lambda$(W/m/K)\\
\hline
QSC & 19.188 & 300 & 1.44 & 1.10(0.06)\\
    &        &     & 2.86 & 1.08(0.05)\\
    &        &     & 5.14 & 1.15(0.07)\\
\\
    & 19.263 & 300 & 2.31 & 1.25(0.06)\\
    &        & 575 & 3.02 & 1.02(0.07)\\
\\
\hline
EAM & 19.045 & 300 & 1.24 & 1.24(0.16)\\
    &        &     & 2.06 & 1.37(0.04)\\
    &        &     & 2.55 & 1.41(0.07)\\
\\
    & 19.263 & 300 & 1.06 & 1.45(0.13)\\
    &        &     & 2.04 & 1.41(0.07)\\
    &        &     & 2.41 & 1.53(0.10)\\
    &        & 575 & 2.82 & 1.08(0.03)\\
\hline
\end{tabular}
\label{AuThermal}
\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{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}$ & $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}{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. 

\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:	3603
Committed:	Tue May 4 17:24:52 2010 UTC (14 years, 1 month ago) by skuang
Content type:	application/x-tex
File size:	38631 byte(s)
Log Message:	edited gold thermal part.