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.  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 & = & \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}

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 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 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 hot 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 flux equation} is essentially a plane which is
perpendicular to the x-axis, with its position governed both by a
targetted 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, the mixture system was simulated.

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

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}

\subsection{Interfaciel Thermal Conductivity}

\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:	3569
Committed:	Tue Mar 9 23:14:21 2010 UTC (14 years, 3 months ago) by gezelter
Content type:	application/x-tex
File size:	23157 byte(s)
Log Message:	Added a figure