group/stokes/stokes.tex

\documentclass[11pt]{article}
\usepackage{amsmath}
\usepackage{amssymb}
\usepackage{setspace}
\usepackage{endfloat}
\usepackage{caption}
%\usepackage{tabularx}
\usepackage{graphicx}
\usepackage{multirow}
%\usepackage{booktabs}
%\usepackage{bibentry}
%\usepackage{mathrsfs}
%\usepackage[ref]{overcite}
\usepackage[square, comma, sort&compress]{natbib}
\usepackage{url}
\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 reference note
\bibpunct{[}{]}{,}{n}{}{;}
\bibliographystyle{achemso}

\begin{document}

\title{Velocity Shearing and Scaling RNEMD: a minimally perturbing
  method for producing temperature and momentum gradients}

\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 nonequilibrium
  velocity and temperature gradients in molecular dynamics simulations
  of heterogeneous systems. This method conserves both the linear
  momentum and total energy of the system and improves previous
  reverse non-equilibrium molecular dynamics (RNEMD) methods while
  retaining equilibrium thermal velocity distributions in each region
  of the system.  The new method avoids the thermal anisotropy
  produced by previous methods by using isotropic velocity scaling and
  shearing on all of the molecules in specific regions. To test the
  method, we have computed the thermal conductivity and shear
  viscosity of model liquid systems as well as the interfacial
  friction coeefficients of a series of solid / liquid interfaces. The
  method's ability to combine the thermal and momentum gradients
  allows us to obtain shear viscosity data for a range of temperatures
  from a single trajectory.
\end{abstract}

\newpage

%\narrowtext

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                          BODY OF TEXT
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\section{Introduction}
One of the standard ways to compute transport coefficients such as the
viscosities and thermal conductivities of liquids is to use 
imposed-flux non-equilibrium molecular dynamics
methods.\cite{MullerPlathe:1997xw,ISI:000080382700030,kuang:164101}
These methods establish stable non-equilibrium conditions in a
simulation box using an applied momentum or thermal flux between
different regions of the box.  The corresponding temperature or
velocity gradients which develop in response to the applied flux is
related (via linear response theory) to the transport coefficient of
interest.  These methods are quite efficient, in that they do not need
many trajectories to provide information about transport properties. To
date, they have been utilized in computing thermal and mechanical
transfer of both homogeneous liquids as well as heterogeneous systems
such as solid-liquid
interfaces.\cite{garde:nl2005,garde:PhysRevLett2009,kuang:AuThl}

The reverse non-equilibrium molecular dynamics (RNEMD) methods utilize
additional constraints that ensure conservation of linear momentum and
total energy of the system while imposing the desired flux.  The RNEMD
methods are therefore capable of sampling various thermodynamically
relevent ensembles, including the micro-canonical (NVE) ensemble,
without resorting to an external thermostat. The original approaches
proposed by M\"{u}ller-Plathe {\it et
  al.}\cite{MullerPlathe:1997xw,ISI:000080382700030} utilize simple
momentum swapping moves for generating energy/momentum fluxes.  The
swapping moves can also be made compatible with particles of different
identities. Although the swapping moves are simple to implement in
molecular simulations, Tenney and Maginn have shown that they create
nonthermal velocity distributions.\cite{Maginn:2010} Furthermore, this
approach is not particularly efficient for kinetic energy transfer
between particles of different identities, particularly when the mass
difference between the particles becomes significant.  This problem
makes applying swapping-move RNEMD methods on heterogeneous
interfacial systems somewhat difficult.

Recently, we developed a somewhat different approach to applying
thermal fluxes in RNEMD simulation using a Non-Isotropic Velocity
Scaling (NIVS) algorithm.\cite{kuang:164101} This algorithm scales all
atomic velocity vectors in two separate regions of a simulated system
using two diagonal scaling matrices.  The scaling matrices are
determined by solving single quartic equation which includes linear
momentum and kinetic energy conservation constraints as well as the
target thermal flux between the regions.  The NIVS method is able to
effectively impose a wide range of kinetic energy fluxes without
significant perturbation to the velocity distributions away from ideal
Maxwell-Boltzmann distributions, even in the presence of heterogeneous
interfaces. We successfully applied this approach in studying the
interfacial thermal conductance at chemically-capped metal-solvent
interfaces.\cite{kuang:AuThl}

The NIVS approach works very well for preparing stable {\it thermal}
gradients.  However, as we pointed out in the original paper, it has
limited application in imposing {\it linear} momentum fluxes (which
are required for measuring shear viscosities). The reason for this is
that linear momentum flux was being imposed by scaling random
fluctuations of the center of the velocity distributions.  Repeated
application of the original NIVS approach resulted in temperature
anisotropy, i.e. the width of the velocity distributions depended on
coordinates perpendicular to the desired gradient direction.  For this
reason, combining thermal and momentum fluxes was difficult with the
original NIVS algorithm.  However, combinations of thermal and
velocity gradients would provide a means to simulate thermal-linear
coupled processes such as Soret effect in liquid flows. Therefore,
developing improved approaches to the scaling imposed-flux methods
would be useful.

In this paper, we improve the RNEMD methods by introducing a novel
approach to creating imposed fluxes. This approach separates the
shearing and scaling of the velocity distributions in different
spatial regions and can apply both transformations within a single
time step.  The approach retains desirable features of previous RNEMD
approaches and is simpler to implement compared to the earlier NIVS
method. In what follows, we first present the Shearing-and-Scaling
(SS) RNEMD method and its implementation in a simulation. Then we
compare the SS-RNEMD method in bulk fluids to previous methods.  We
also compute interfacial frictions are computed for a series of
heterogeneous interfaces.

\section{Methodology}
In an approach similar to the earlier NIVS method,\cite{kuang:164101}
we consider a periodic system which has been divided into a series of
slabs along a single axis (e.g. $z$). The unphysical thermal and
momentum fluxes are applied from one of the end slabs to the center
slab, and thus the thermal flux produces a higher temperature in the
center slab than in the end slab, and the momentum flux results in a
center slab moving along the positive $y$ axis (see Fig. \ref{cartoon}).
The applied fluxes can be set to negative values if the reverse
gradients are desired.  For convenience the center slab is denoted as
the {\it hot} or {\it ``H''} slab, while the end slab is denoted {\it
  ``C''} (or {\it cold}).

\begin{figure}
\includegraphics[width=\linewidth]{cartoon}
\caption{The SS-RNEMD approach we are introducing imposes unphysical
  transfer of both momentum and kinetic energy between a ``hot'' slab
  and a ``cold'' slab in the simulation box.  The molecular system
  responds to this imposed flux by generating both momentum and
  temperature gradients.  The slope of the gradients can then be used
  to compute transport properties (e.g. shear viscosity and thermal
  conductivity).}
\label{cartoon}
\end{figure}

To impose these fluxes, we periodically apply a set of operations on
the velocities of particles {$i$} within the cold slab and a separate
operation on particles {$j$} within the hot slab.
\begin{eqnarray}
\vec{v}_i & \leftarrow & c\cdot\left(\vec{v}_i - \langle\vec{v}_c
  \rangle\right) + \left(\langle\vec{v}_c\rangle + \vec{a}_c\right) \\
\vec{v}_j & \leftarrow & h\cdot\left(\vec{v}_j - \langle\vec{v}_h
  \rangle\right) + \left(\langle\vec{v}_h\rangle + \vec{a}_h\right)
\end{eqnarray}
where $\langle\vec{v}_c\rangle$ and $\langle\vec{v}_h\rangle$ denote
the instantaneous average velocity of the molecules within slabs $C$
and $H$ respectively.  When a momentum flux $\vec{j}_z(\vec{p})$ is
present, these slab-averaged velocities also get corresponding
incremental changes ($\vec{a}_c$ and $\vec{a}_h$ respectively) that
are applied to all particles within each slab.  The incremental
changes are obtained using Newton's second law:
\begin{eqnarray}
M_c \vec{a}_c & = & -\vec{j}_z(\vec{p}) \Delta t \label{eq:newton1} \\
M_h \vec{a}_h & = &  \vec{j}_z(\vec{p}) \Delta t
\label{eq:newton2}
\end{eqnarray}
where $M$ denotes total mass of particles within a slab:
\begin{eqnarray}
M_c & = & \sum_{i = 1}^{N_c} m_i \\
M_h & = & \sum_{j = 1}^{N_h} m_j
\end{eqnarray}
and $\Delta t$ is the interval between two separate operations.

The operations in Eqs. \ref{eq:newton1} and \ref{eq:newton2} already
conserve the linear momentum of a periodic system.  To further satisfy
total energy conservation as well as to impose the thermal flux $J_z$,
the following constraint equations must be solved for the two scaling
variables $c$ and $h$:
\begin{eqnarray}
K_c - J_z\Delta t & = & c^2 (K_c - \frac{1}{2}M_c \langle\vec{v}_c
\rangle^2) + \frac{1}{2}M_c (\langle \vec{v}_c \rangle + \vec{a}_c)^2 \\
K_h + J_z\Delta t & = & h^2 (K_h - \frac{1}{2}M_h \langle\vec{v}_h
\rangle^2) + \frac{1}{2}M_h (\langle \vec{v}_h \rangle + \vec{a}_h)^2.
\label{constraint}
\end{eqnarray}
Here $K_c$ and $K_h$ denote the translational kinetic energy of slabs
$C$ and $H$ respectively. These conservation equations are sufficient
to ensure total energy conservation, as the operations applied in our
method do not change the kinetic energy related to orientational
degrees of freedom or the potential energy of a system (as long as the
potential energy is independent of particle velocity).

Equations \ref{eq:newton1}-\ref{constraint} are sufficient to
determine the velocity scaling coefficients ($c$ and $h$) as well as
$\vec{a}_c$ and $\vec{a}_h$. Note that there are usually two roots
respectively for $c$ and $h$. However, the positive roots (which are
closer to 1) are chosen so that the perturbations to the system are
minimal.  Figure \ref{method} illustrates the implementation of this
algorithm in an individual step.

\begin{figure}
\centering
\includegraphics[width=5in]{method}
\caption{Illustration of a single step implementation of the
  algorithm.  Starting with the velocity distributions for the two
  slabs in a shearing fluid, the transformation is used to apply the
  effect of both a thermal and a momentum flux from the ``c'' slab to
  the ``h'' slab. As the figure shows, gaussian distributions are
  preserved by both the scaling and shearing operations.}
\label{method}
\end{figure}

By implementing these operations at a fixed frequency, stable thermal
and momentum fluxes can both be applied and the corresponding
temperature and momentum gradients can be established.

Compared to the previous NIVS method, the SS-RNEMD approach is
computationally simpler in that only quadratic equations are involved
to determine a set of scaling coefficients, while the NIVS method
required solution of quartic equations. Furthermore, this method
implements {\it isotropic} scaling of velocities in respective slabs,
unlike NIVS, which required extra flexibility in the choice of scaling
coefficients to allow for the energy and momentum constraints.  Most
importantly, separating the linear momentum flux from velocity scaling
avoids the underlying cause of the thermal anisotropy in NIVS.  In
later sections, we demonstrate that this can improve the calculated
shear viscosities from RNEMD simulations.

The SS-RNEMD approach has advantages over the original momentum
swapping in many respects. In the swapping method, the velocity
vectors involved are usually very different (or the generated flux
would be quite small), thus the swapping tends to create strong
perturbations in the neighborhood of the particles involved. In
comparison, the SS approach distributes the flux widely to all
particles in a slab so that perturbations in the flux generating
region are minimized. Additionally, momentum swapping steps tend to
result in nonthermal velocity distributions when the imposed flux is
large and diffusion from the neighboring slabs cannot carry momentum
away quickly enough.  In comparison, the scaling and shearing moves
made in the SS approach preserve the shapes of the equilibrium
velocity disributions (e.g. Maxwell-Boltzmann).  The results presented
in later sections will illustrate that this is quite helpful in
retaining reasonable thermal distributions in a simulation.

\section{Computational Details}
The algorithm has been implemented in our MD simulation code,
OpenMD.\cite{Meineke:2005gd,openmd} We will first compare the method
with previous RNEMD methods and equilibrium MD in homogeneous
fluids (Lennard-Jones and SPC/E water).  We have also used the new
method to simulate the interfacial friction of different heterogeneous
interfaces (Au (111) with organic solvents, Au(111) with SPC/E water,
and the Ice Ih - liquid water interface).

\subsection{Simulation Protocols}
The systems we investigated were set up in orthorhombic simulation
cells with periodic boundary conditions in all three dimensions. The
$z$-axes of these cells were typically quite long and served as the
temperature and/or momentum gradient axes.  The cells were evenly
divided into $N$ slabs along this axis, with $N$ varying for the
individual system. The $x$ and $y$ axes were of similar lengths in all
simulations. In all cases, before introducing a nonequilibrium method
to establish steady thermal and/or momentum gradients, equilibration
simulations were run under the canonical ensemble with a Nos\'e-Hoover
thermostat\cite{hoover85} followed by further equilibration using
standard constant energy (NVE) conditions.  For SPC/E water
simulations, isobaric-isothermal equilibrations\cite{melchionna93}
were performed before equilibration to reach standard densities at
atmospheric pressure (1 bar); for interfacial systems, similar
equilibrations with anisotropic box relaxations are used to relax the
surface tension in the $xy$ plane.

While homogeneous fluid systems can be set up with random
configurations, interfacial systems are typically prepared with a
single crystal face presented perpendicular to the $z$-axis. This
crystal face is aligned in the x-y plane of the periodic cell, and the
solvent occupies the region directly above and below a crystalline
slab.  The preparation and equilibration of butanethiol covered
Au(111) surfaces, as well as the solvation and equilibration processes
used for these interfaces are described in detail in reference
\cite{kuang:AuThl}.

For the ice / liquid water interfaces, the basal surface of ice
lattice was first constructed. Hirsch and Ojam\"{a}e
\cite{doi:10.1021/jp048434u} explored the energetics of ice lattices
with all possible proton ordered configurations. We utilized Hirsch
and Ojam\"{a}e's structure 6 ($P2_12_12_1$) which is an orthorhombic
cell giving a proton-ordered version of Ice Ih.  The basal face of ice
in this unit cell geometry is the $\{0~0~1\}$ face.  Although
experimental solid/liquid coexistant temperature under normal pressure
are close to 273K, Bryk and Haymet's simulations of ice/liquid water
interfaces with different models suggest that for SPC/E, the most
stable interface is observed at 225$\pm$5K.\cite{bryk:10258}
Therefore, our ice/liquid water simulations were carried out at an
average temperature of 225K.  Molecular translation and orientational
restraints were applied in the early stages of equilibration to
prevent melting of the ice slab.  These restraints were removed during
NVT equilibration, well before data collection was carried out.

\subsection{Force Field Parameters}
For comparing the SS-RNEMD method with previous work, we retained
force field parameters consistent with previous simulations. Argon is
the Lennard-Jones fluid used here, and its results are reported in
reduced units for purposes of direct comparison with previous
calculations.

For our water simulations, we utilized the SPC/E model throughout this
work.  Previous work for transport properties of SPC/E water model is
available\cite{Bedrov:2000,10.1063/1.3330544,Medina2011} so direct
comparison with previous calculation methods is possible.

The Au-Au interaction parameters in all simulations are described by
the quantum Sutton-Chen (QSC) formulation.\cite{PhysRevB.59.3527} The
QSC potentials include zero-point quantum corrections and are
reparametrized for accurate surface energies compared to the
Sutton-Chen potentials.\cite{Chen90} For gold/water interfaces, the
Spohr potential was adopted\cite{ISI:000167766600035} to depict
Au-H$_2$O interactions.

For our gold/organic liquid interfaces, the small organic molecules
included in our simulations are the Au surface capping agent
butanethiol as well as liquid hexane and liquid toluene. The
United-Atom
models\cite{TraPPE-UA.thiols,TraPPE-UA.alkanes,TraPPE-UA.alkylbenzenes}
for these components were used in this work for computational
efficiency, while maintaining good accuracy. We refer readers to our
previous work\cite{kuang:AuThl} for further details of these models,
as well as the interactions between Au and the above organic molecule
components.

\subsection{Thermal conductivities}
When the linear momentum flux $\vec{j}_z(\vec{p})$ is set to zero and
the target $J_z$ is non-zero, SS-RNEMD imposes kinetic energy transfer
between the slabs, which can be used for computation of thermal
conductivities. Similar to previous RNEMD methods, we assume that we
are in the linear response regime of the temperature gradient with
respect to the thermal flux.  The thermal conductivity ($\lambda$) can
be calculated using the imposed kinetic energy flux and the measured
thermal gradient:
\begin{equation}
J_z = -\lambda \frac{\partial T}{\partial z}
\end{equation}
Like other imposed-flux methods, the energy flux was calculated using
the total non-physical energy transferred (${E_{total}}$) from slab
``c'' to slab ``h'', which was recorded throughout the simulation, and
the total data collection time $t$:
\begin{equation}
J_z = \frac{E_{total}}{2 t L_x L_y}.
\end{equation}
Here, $L_x$ and $L_y$ denote the dimensions of the plane in a
simulation cell perpendicular to the thermal gradient, and a factor of
two in the denominator is necessary as the heat transport occurs in
both the $+z$ and $-z$ directions. The average temperature gradient
${\langle\partial T/\partial z\rangle}$ can be obtained by a linear
regression of the temperature profile, which is recorded during a
simulation for each slab in a cell. For Lennard-Jones simulations,
thermal conductivities are reported in reduced units
(${\lambda^*=\lambda \sigma^2 m^{1/2} k_B^{-1}\varepsilon^{-1/2}}$).

\subsection{Shear viscosities}
Alternatively, when the linear momentum flux $\vec{j}_z(\vec{p})$ is
non-zero in either the $x$ or $y$ directions, the method can be used
to compute the shear viscosity.  For isotropic systems, the direction
of $\vec{j}_z(\vec{p})$ on the $xy$ plane should not matter, but the
ability to arbitarily specify the vector direction in our method could
provide convenience when working with anisotropic interfaces.

In a manner similar to the thermal conductivity calculations, a linear
response of the momentum gradient with respect to the shear stress is
assumed, and the shear viscosity ($\eta$) can be obtained with the
imposed momentum flux (e.g. in $x$ direction) and the measured
velocity gradient:
\begin{equation}
j_z(p_x) = -\eta \frac{\partial v_x}{\partial z}
\end{equation}
where the flux is similarly defined:
\begin{equation}
j_z(p_x) = \frac{P_x}{2 t L_x L_y}
\end{equation}
with $P_x$ being the total non-physical momentum transferred within
the data collection time. Also, the averaged velocity gradient
${\langle\partial v_x/\partial z\rangle}$ can be obtained using linear
regression of the $x$ component of the mean velocity ($\langle
v_x\rangle$) in each of the bins. For Lennard-Jones simulations, shear
viscosities are also reported in reduced units
(${\eta^*=\eta\sigma^2(\varepsilon m)^{-1/2}}$).

Although $J_z$ may be switched off for shear viscosity simulations,
the SS-RNEMD method allows the user the ability to simultaneously
impose both a thermal and a momentum flux during a single
simulation. This can create system with coincident temperature and a
velocity gradients.  Since the viscosity is generally a function of
temperature, the local viscosity depends on the local temperature in
the fluid. Therefore, in a single simulation, viscosity at $z$
(corresponding to a certain $T$) can be computed with the applied
shear flux and the local velocity gradient (which can be obtained by
finite difference approximation). This means that the temperature
dependence of the viscosity can be mapped out in only one
trajectory. Results for shear viscosity computations of SPC/E water
will demonstrate SS-RNEMD's efficiency in this respect.

\subsection{Interfacial friction and slip length}
Shear stress creates a velocity gradient within bulk fluid phases, but
at solid-liquid interfaces, the effects of a shear stress depend on
the molecular details of the interface. The interfacial friction
coefficient, $\kappa$, relates the shear stress (e.g. along the
$x$-axis) with the relative fluid velocity tangent to the interface:
\begin{equation}
j_z(p_x) = \kappa \left[v_x(fluid) - v_x(solid)\right]
\end{equation}
where $v_x(fluid)$ and $v_x(solid)$ are the velocities measured
directly adjacent to the interface.  Under ``stick'' boundary
condition, $\Delta v_x|_\mathrm{interface} \rightarrow 0$, which leads
to $\kappa\rightarrow\infty$. However, for ``slip'' boundary
conditions at a solid-liquid interface, $\kappa$ becomes finite. To
characterize the interfacial boundary conditions, the slip length,
$\delta$, is defined by the ratio of the fluid-phase viscosity to the
friction coefficient of the interface:
\begin{equation}
\delta = \frac{\eta}{\kappa}
\end{equation}
In ``no-slip'' or ``stick'' boundary conditions, $\delta\rightarrow
0$, and $\delta$ is a measure of how ``slippery'' an interface is.
Figure \ref{slipLength} illustrates how this quantity is defined and
computed for a solid-liquid interface. 

\begin{figure}
\includegraphics[width=\linewidth]{defDelta}
\caption{The slip length $\delta$ can be obtained from a velocity
  profile of a solid-liquid interface simulation, when a momentum flux
  is applied. The data shown is for a simulated Au/hexane interface.
  The Au crystalline region is moving as a block (lower dots), while
  the measured velocity gradient in the hexane phase is discontinuous
  a the interface.}
\label{slipLength}
\end{figure}

Since the method can be applied for interfaces as well as for bulk
materials, the shear stress is applied in an identical manner to the
shear viscosity computations, e.g. by applying an unphysical momentum
flux, $j_z(\vec{p})$.  With the correct choice of $\vec{p}$ in the
$x-y$ plane, one can compute friction coefficients and slip lengths
for a number of different dragging vectors on a given slab.  The
corresponding velocity profiles can be obtained as shown in Figure
\ref{slipLength}, in which the velocity gradients within the liquid
phase and the velocity difference at the liquid-solid interface can be
easily measured from saved simulation data.

\section{Results and Discussions}
\subsection{Lennard-Jones fluid}
Our orthorhombic simulation cell for the Lennard-Jones fluid has
identical parameters to our previous work\cite{kuang:164101} to
facilitate comparison. Thermal conductivities and shear viscosities
were computed with the new algorithm applied to the simulations. The
results of thermal conductivity are compared with our previous NIVS
algorithm. However, since the NIVS algorithm produced temperature
anisotropy for shear viscocity computations, these results are instead
compared to the previous momentum swapping approaches. Table \ref{LJ}
lists these values with various fluxes in reduced units.

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

      \caption{Thermal conductivity ($\lambda^*$) and shear viscosity
        ($\eta^*$) (in reduced units) of a Lennard-Jones fluid at
        ${\langle T^* \rangle = 0.72}$ and ${\rho^* = 0.85}$ computed
        at various momentum fluxes.  The new method yields similar
        results to previous RNEMD methods. All results are reported in
        reduced unit. Uncertainties are indicated in parentheses.}
      
      \begin{tabular}{cccccc}
        \hline\hline
        \multicolumn{2}{c}{Momentum Exchange} &
        \multicolumn{2}{c}{$\lambda^*$} &
        \multicolumn{2}{c}{$\eta^*$} \\
        \hline
        Swap Interval & Equivalent $J_z^*$ or $j_z^*(p_x)$ &
        NIVS\cite{kuang:164101} & This work & Swapping & This work \\
        \hline
        0.116 & 0.16  & 7.30(0.10) & 6.25(0.23) & 3.57(0.06) & 3.53(0.16)\\
        0.232 & 0.09  & 6.95(0.09) & 8.02(0.56) & 3.64(0.05) & 3.43(0.17)\\
        0.463 & 0.047 & 7.19(0.07) & 6.48(0.15) & 3.52(0.16) & 3.51(0.08)\\
        0.926 & 0.024 & 7.19(0.28) & 7.02(0.13) & 3.72(0.05) & 3.79(0.11)\\
        1.158 & 0.019 & 7.98(0.33) & 7.37(0.23) & 3.42(0.06) & 3.53(0.06)\\
        \hline\hline
      \end{tabular}
      \label{LJ}
    \end{center}
  \end{minipage}
\end{table*}

\subsubsection{Thermal conductivity}
Our thermal conductivity calculations with this method yields
comparable results to the previous NIVS algorithm. This indicates that
the thermal gradients introduced using this method are also close to
previous RNEMD methods. Simulations with moderately higher thermal
fluxes tend to yield more reliable thermal gradients and thus avoid
large errors, while overly high thermal fluxes could introduce side
effects such as non-linear temperature gradient response or
inadvertent phase transitions.

Since the scaling operation is isotropic in this method, one does not
need extra care to ensure temperature isotropy between the $x$, $y$
and $z$ axes, while for NIVS, thermal anisotropy might happen if the
criteria function for choosing scaling coefficients does not perform
as expected. Furthermore, this method avoids inadvertent concomitant
momentum flux when only thermal flux is imposed, which could not be
achieved with swapping or NIVS approaches. The thermal energy exchange
in swapping ($\vec{p}_i$ in slab ``c'' with $\vec{p}_j$ in slab ``h'')
or NIVS (total slab momentum conponemts $P^\alpha$ scaled to $\alpha
P^\alpha$) would not achieve this effect unless thermal flux vanishes
(i.e. $\vec{p}_i=\vec{p}_j$ or $\alpha=1$ which do not contribute to
applying a thermal flux). In this sense, this method aids in achieving
minimal perturbation to a simulation while imposing a thermal flux.

\subsubsection{Shear viscosity}
Table \ref{LJ} also compares our shear viscosity results with the
momentum swapping approach. Our calculations show that our method
predicted similar values of shear viscosities to the momentum swapping
approach, as well as the velocity gradient profiles. Moderately larger
momentum fluxes are helpful to reduce the errors of measured velocity
gradients and thus the final result. However, it is pointed out that
the momentum swapping approach tends to produce nonthermal velocity
distributions.\cite{Maginn:2010}

To examine that temperature isotropy holds in simulations using our
method, we measured the three one-dimensional temperatures in each of
the slabs (Figure \ref{tempXyz}). Note that the $x$-dimensional
temperatures were calculated after subtracting the contribution from
bulk velocities of the slabs. The one-dimensional temperature profiles
showed no observable difference between the three dimensions. This
ensures that isotropic scaling automatically preserves temperature
isotropy and that our method is useful in shear viscosity
computations.

\begin{figure}
\includegraphics[width=\linewidth]{tempXyz}
\caption{Unlike the previous NIVS algorithm, the new method does not
  produce a thermal anisotropy. No temperature difference between
  different dimensions were observed beyond the magnitude of the error
  bars. Note that the two ``hotter'' regions are caused by the shear
  stress, as reported by Tenney and Maginn\cite{Maginn:2010}, but not
  an effect that only observed in our methods.}
\label{tempXyz}
\end{figure}

Furthermore, the velocity distribution profiles are tested by imposing
a large shear stress into the simulations. Figure \ref{vDist}
demonstrates how our method is able to maintain thermal velocity
distributions against the momentum swapping approach even under large
imposed fluxes. Previous swapping methods tend to deplete particles of
positive velocities in the negative velocity slab (``c'') and vice
versa in slab ``h'', where the distributions leave notchs. This
problematic profiles become significant when the imposed-flux becomes
larger and diffusions from neighboring slabs could not offset the
depletions. Simutaneously, abnormal peaks appear corresponding to
excessive particles having velocity swapped from the other slab. These
nonthermal distributions limit applications of the swapping approach
in shear stress simulations. Our method avoids the above problematic
distributions by altering the means of applying momentum
fluxes. Comparatively, velocity distributions recorded from
simulations with our method is so close to the ideal thermal
prediction that no obvious difference is shown in Figure
\ref{vDist}. Conclusively, our method avoids problems that occurs in
previous RNEMD methods and provides a useful means for shear viscosity
computations.

\begin{figure}
\includegraphics[width=\linewidth]{velDist}
\caption{Velocity distributions that develop under the swapping and
  our methods at a large flux. These distributions were obtained from
  Lennard-Jones simulations with $j_z^*(p_x)\sim 0.4$ (equivalent to a
  swapping interval of 20 time steps). This is a relatively large flux
  to demonstrate the nonthermal distributions that develop under the
  swapping method. In comparison, distributions produced by our method
  are very close to the ideal thermal situations.}
\label{vDist}
\end{figure}

\subsection{Bulk SPC/E water}
We extend our applications of thermal conductivity and shear viscosity
computations to a complex fluid model of SPC/E water. A simulation
cell with 1000 molecules was set up in the similar manner as in
\cite{kuang:164101}. For thermal conductivity simulations,
measurements were taken to compare with previous RNEMD methods; for
shear viscosity computations, simulations were run under a series of
temperatures (with corresponding pressure relaxation using the
isobaric-isothermal ensemble\cite{melchionna93}), and results were
compared to available data from EMD
methods\cite{10.1063/1.3330544,Medina2011}. Besides, a simulation with
both thermal and momentum gradient was carried out to map out shear
viscosity as a function of temperature to see the effectiveness and
accuracy our method could reach.

\subsubsection{Thermal conductivity}
Table \ref{spceThermal} summarizes our thermal conductivity
computations under different temperatures and thermal gradients, in
comparison to the previous NIVS results\cite{kuang:164101} and
experimental measurements\cite{WagnerKruse}. Note that no appreciable
drift of total system energy or temperature was observed when our
method is applied, which indicates that our algorithm conserves total
energy well for systems involving electrostatic interactions.

Measurements using our method established similar temperature
gradients to the previous NIVS method. Our simulation results are in
good agreement with those from previous simulations. And both methods
yield values in reasonable agreement with experimental
values. Simulations using moderately higher thermal gradient or those
with longer gradient axis ($z$) for measurement seem to have better
accuracy, from our results.

\begin{table*}
  \begin{minipage}{\linewidth}
    \begin{center}
      
      \caption{Thermal conductivity of SPC/E water under various
        imposed thermal gradients. Uncertainties are indicated in
        parentheses.}
      
      \begin{tabular}{ccccc}
        \hline\hline
        $\langle T\rangle$ & $\langle dT/dz\rangle$ & \multicolumn{3}{c}
        {$\lambda (\mathrm{W m}^{-1} \mathrm{K}^{-1})$} \\
        (K) & (K/\AA) & This work & Previous NIVS\cite{kuang:164101} &
        Experiment\cite{WagnerKruse} \\
        \hline
        300 & 0.8 & 0.815(0.027)     & 0.770(0.008) & 0.61 \\
        318 & 0.8 & 0.801(0.024)     & 0.750(0.032) & 0.64 \\
            & 1.6 & 0.766(0.007)     & 0.778(0.019) &      \\
            & 0.8 & 0.786(0.009)\footnote{Simulation with $L_z$
                     twice as long.} &              &      \\
        \hline\hline
      \end{tabular}
      \label{spceThermal}
    \end{center}
  \end{minipage}
\end{table*}

\subsubsection{Shear viscosity}
The improvement our method achieves for shear viscosity computations
enables us to apply it on SPC/E water models. The series of
temperatures under which our shear viscosity calculations were carried
out covers the liquid range under normal pressure. Our simulations
predict a similar trend of $\eta$ vs. $T$ to EMD results we refer to
(Table \ref{spceShear}). Considering subtlties such as temperature or
pressure/density errors in these two series of measurements, our
results show no significant difference from those with EMD
methods. Since each value reported using our method takes only one
single trajectory of simulation, instead of average from many
trajectories when using EMD, our method provides an effective means
for shear viscosity computations.

\begin{table*}
  \begin{minipage}{\linewidth}
    \begin{center}
      
      \caption{Computed shear viscosity of SPC/E water under different
        temperatures. Results are compared to those obtained with EMD
        method[CITATION]. Uncertainties are indicated in parentheses.}
      
      \begin{tabular}{cccc}
        \hline\hline
        $T$ & $\partial v_x/\partial z$ & \multicolumn{2}{c}
        {$\eta (\mathrm{mPa}\cdot\mathrm{s})$} \\
        (K) & (10$^{10}$s$^{-1}$) & This work & Previous
        simulations\cite{Medina2011} \\
        \hline
        273 & 1.12 & 1.218(0.004) & 1.282(0.048) \\
            & 1.79 & 1.140(0.012) &              \\
        303 & 2.09 & 0.646(0.008) & 0.643(0.019) \\
        318 & 2.50 & 0.536(0.007) &              \\
            & 5.25 & 0.510(0.007) &              \\
            & 2.82 & 0.474(0.003)\footnote{Simulation with $L_z$ twice
                        as long.} &              \\
        333 & 3.10 & 0.428(0.002) & 0.421(0.008) \\
        363 & 2.34 & 0.279(0.014) & 0.291(0.005) \\
            & 4.26 & 0.306(0.001) &              \\
        \hline\hline
      \end{tabular}
      \label{spceShear}
    \end{center}
  \end{minipage}
\end{table*}

A more effective way to map out $\eta$ vs $T$ is to combine a momentum
flux with a thermal flux. Figure \ref{Tvxdvdz} shows the thermal and
velocity gradient in one such simulation. At different positions with
different temperatures, the velocity gradient is not a constant but
can be computed locally. With the data provided in Figure
\ref{Tvxdvdz}, a series of $\eta$ is calculated as in Figure
\ref{etaT} and a linear fit was performed to $\partial v_x/\partial z$
vs. $z$ so that the resulted $\eta$ can be present as a curve as
well. For comparison, other results are also mapped in the figure.

\begin{figure}
\includegraphics[width=\linewidth]{tvxdvdz}
\caption{With a combination of a thermal and a momentum flux, a
  simulation can have both a temperature (top) and a velocity (middle)
  gradient. Due to the thermal gradient, $\partial v_x/\partial z$ is
  not constant but can be computed using finite difference
  approximations (lower). These data can be used further to calculate
  $\eta$ vs $T$ (Figure \ref{etaT}).}
\label{Tvxdvdz}
\end{figure}

From Figure \ref{etaT}, one can see that the generated curve agrees
well with the above RNEMD simulations at different temperatures, as
well as results reported using EMD
methods\cite{10.1063/1.3330544,Medina2011} in much of the temperature
range simulated. However, this curve has relatively large error in
lower temperature regions and has some difference in predicting $\eta$
near 273K. Provided that this curve only takes one trajectory to
generate, these results are of satisfactory efficiency and
accuracy. Since previous work already pointed out that the SPC/E model
tends to predict lower viscosity compared to experimental
data,\cite{Medina2011} experimental comparison are not given here.

\begin{figure}
\includegraphics[width=\linewidth]{etaT}
\caption{The curve generated by single simulation with thermal and
  momentum gradient predicts satisfatory values in much of the
  temperature range under test.}
\label{etaT}
\end{figure}

\subsection{Interfacial frictions and slip lengths}
Another attractive aspect of our method is the ability to apply
momentum and/or thermal flux in nonhomogeneous systems, where
molecules of different identities (or phases) are segregated in
different regions. We have previously studied the interfacial thermal
transport of a series of metal gold-liquid
surfaces\cite{kuang:164101,kuang:AuThl}, and would like to further
investigate the relationship between this phenomenon and the
interfacial frictions.

Table \ref{etaKappaDelta} includes these computations and previous
calculations of corresponding interfacial thermal conductance. For
bare Au(111) surfaces, slip boundary conditions were observed for both
organic and aqueous liquid phases, corresponding to previously
computed low interfacial thermal conductance. In comparison, the
butanethiol covered Au(111) surface appeared to be sticky to the
organic liquid layers in our simulations. We have reported conductance
enhancement effect for this surface capping agent,\cite{kuang:AuThl}
and these observations have a qualitative agreement with the thermal
conductance results. This agreement also supports discussions on the
relationship between surface wetting and slip effect and thermal
conductance of the
interface.\cite{PhysRevLett.82.4671,doi:10.1080/0026897031000068578,garde:PhysRevLett2009}

\begin{table*}
  \begin{minipage}{\linewidth}
    \begin{center}
      
      \caption{Computed interfacial friction coefficient values for
        interfaces with various components for liquid and solid
        phase. Error estimates are indicated in parentheses.}
      
      \begin{tabular}{llcccccc}
        \hline\hline
        Solid & Liquid & $T$ & $j_z(p_x)$ & $\eta_{liquid}$ & $\kappa$
        & $\delta$ & $G$\footnote{References \cite{kuang:AuThl} and
          \cite{kuang:164101}.} \\
        surface & molecules & K & MPa & mPa$\cdot$s &
        10$^4$Pa$\cdot$s/m & nm & MW/m$^2$/K \\
        \hline
        Au(111) & hexane & 200 & 1.08 & 0.197(0.009) & 5.30(0.36) &
        3.72 & 46.5 \\
                &        &     & 2.15 & 0.141(0.002) & 5.31(0.26) &
        2.76 &      \\
        Au-SC$_4$H$_9$ & hexane & 200 & 2.16 & 0.286(0.019) & $\infty$
        & 0 & 131 \\
                       &        &     & 5.39 & 0.320(0.006) & $\infty$
        & 0 &     \\
        \hline
        Au(111) & toluene & 200 & 1.08 & 0.722(0.035) & 15.7(0.7) &
        4.60 & 70.1 \\
                &         &     & 2.16 & 0.544(0.030) & 11.2(0.5) &
        4.86 &      \\
        Au-SC$_4$H$_9$ & toluene & 200 & 5.39 & 0.980(0.057) &
        $\infty$ & 0 & 187 \\
                       &         &     & 10.8 & 0.995(0.005) &
        $\infty$ & 0 &     \\
        \hline
        Au(111) & water & 300 & 1.08 & 0.399(0.050) & 1.928(0.022) &
        20.7 & 1.65 \\
                &       &     & 2.16 & 0.794(0.255) & 1.895(0.003) &
        41.9 &      \\
        \hline
        ice(basal) & water & 225 & 19.4 & 15.8(0.2) & $\infty$ & 0 & \\
        \hline\hline
      \end{tabular}
      \label{etaKappaDelta}
    \end{center}
  \end{minipage}
\end{table*}

An interesting effect alongside the surface friction change is
observed on the shear viscosity of liquids in the regions close to the
solid surface. In our results, $\eta$ measured near a ``slip'' surface
tends to be smaller than that near a ``stick'' surface. This may
suggest the influence from an interface on the dynamic properties of
liquid within its neighbor regions. It is known that diffusions of
solid particles in liquid phase is affected by their surface
conditions (stick or slip boundary).\cite{10.1063/1.1610442} Our
observations could provide a support to this phenomenon.

In addition to these previously studied interfaces, we attempt to
construct ice-water interfaces and the basal plane of ice lattice was
studied here. In contrast to the Au(111)/water interface, where the
friction coefficient is substantially small and large slip effect
presents, the ice/liquid water interface demonstrates strong
solid-liquid interactions and appears to be sticky. The supercooled
liquid phase is an order of magnitude more viscous than measurements
in previous section. It would be of interst to investigate the effect
of different ice lattice planes (such as prism and other surfaces) on
interfacial friction and the corresponding liquid viscosity.

\section{Conclusions}
Our simulations demonstrate the validity of our method in RNEMD
computations of thermal conductivity and shear viscosity in atomic and
molecular liquids. Our method maintains thermal velocity distributions
and avoids thermal anisotropy in previous NIVS shear stress
simulations, as well as retains attractive features of previous RNEMD
methods. There is no {\it a priori} restrictions to the method to be
applied in various ensembles, so prospective applications to
extended-system methods are possible.

Our method is capable of effectively imposing thermal and/or momentum
flux accross an interface. This facilitates studies that relates
dynamic property measurements to the chemical details of an
interface. Therefore, investigations can be carried out to
characterize interfacial interactions using the method.

Another attractive feature of our method is the ability of
simultaneously introducing thermal and momentum gradients in a
system. This facilitates us to effectively map out the shear viscosity
with respect to a range of temperature in single trajectory of
simulation with satisafactory accuracy. Complex systems that involve
thermal and momentum gradients might potentially benefit from
this. For example, the Soret effects under a velocity gradient might
be models of interest to purification and separation researches.

\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

\bibliography{stokes}

\end{doublespace}
\end{document}
Revision:	3786
Committed:	Wed Dec 14 21:59:12 2011 UTC (12 years, 9 months ago) by gezelter
Content type:	application/x-tex
File size:	44054 byte(s)
Log Message: