group/stokes/stokes.tex

\documentclass{tMPH2e}
\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{tMPH}
\markboth{Kuang \& Gezelter}{Molecular Physics}
\begin{document}

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

\author{Shenyu Kuang and J. Daniel Gezelter$^{\ast}$\thanks{
    $^{\ast}$ Corresponding author. Email: gezelter@nd.edu} \\ \vspace{6pt}
{\em Department of Chemistry and Biochemistry,\\
University of Notre Dame\\
Notre Dame, Indiana 46556} \\ \vspace{6pt} 
\received{15 December 2011}
}

\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 and interfacial 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 (VSS) 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 velocity shearing and
scaling (VSS) RNEMD method and its implementation in a
simulation. Then we compare the VSS-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 VSS-RNEMD approach imposes unphysical transfer of both
  linear momentum and kinetic energy between a ``hot'' slab and a
  ``cold'' slab in the simulation box.  The 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) \label{eq:xformc} \\
\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) \label{eq:xformh}
\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{A single step implementation of the VSS algorithm starts with
  the velocity distributions for the two slabs in a shearing fluid
  (solid lines).  Equations \ref{eq:xformc} and \ref{eq:xformh} are
  used to scale and shear velocities in both slabs to mimic the effect
  of both a thermal and a momentum flux. 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 VSS-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 VSS-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 VSS 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 VSS 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 solvated 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 that produces 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 VSS-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-molecule organic
solvents included in our simulations are hexane and toluene. The
United-Atom models\cite{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, VSS-RNEMD imposes kinetic energy
transfer between the slabs, which can be used for computation of
thermal conductivities. As with other 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
then 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
{\it C} to slab {\it 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 the
simulation cell that is 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
be convenient 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 VSS-RNEMD method allows the user the ability to {\it
 simultaneously} impose both a thermal and a momentum flux during a
single simulation. This creates configurations 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 VSS-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 velocity of the fluid 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 ``stick'' boundary conditions, $\delta\rightarrow 0$, so $\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 along an axis perpendicular to the interface. The data shown
  is for an Au / hexane interface -- the crystalline region (Au) is
  moving as a block (lower dots), while the measured velocity gradient
  in the hexane phase is discontinuous at the interface.}
\label{slipLength}
\end{figure}

Since the VSS-RNEMD 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 previous work to facilitate
comparison.\cite{kuang:164101} Thermal conductivities and shear
viscosities were computed with the new VSS algorithm, and the results
were compared with the previous NIVS algorithm.  However, since the
NIVS algorithm produces temperature anisotropy under shear stress,
these results are also compared to the previous momentum swapping
approaches. Table \ref{LJ} lists these values with various fluxes in
reduced units.

\begin{table}
  \tbl{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 VSS method yields similar
    results to previous RNEMD methods. All results are reported in
    reduced unit. Uncertainties are indicated in parentheses.}
      
  {\begin{tabular}{cc|ccc|cc}\toprule
      Swap & Equivalent &
      \multicolumn{3}{|c}{$\lambda^*$} &
      \multicolumn{2}{|c}{$\eta^*$} \\
      Interval & $J_z^*$ or $j_z^*(p_x)$ & Swapping &
        NIVS\cite{kuang:164101} & This work & Swapping & This work \\
        \colrule
        0.116 & 0.16  & 7.03(0.34) & 7.30(0.10) & 6.25(0.23) & 3.57(0.06) & 3.53(0.16)\\
        0.232 & 0.09  & 7.03(0.14) & 6.95(0.09) & 8.02(0.56) & 3.64(0.05) & 3.43(0.17)\\
        0.463 & 0.047 & 6.91(0.42) & 7.19(0.07) & 6.48(0.15) & 3.52(0.16) & 3.51(0.08)\\
        0.926 & 0.024 & 7.52(0.15) & 7.19(0.28) & 7.02(0.13) & 3.72(0.05) & 3.79(0.11)\\
        1.158 & 0.019 & 7.41(0.29) & 7.98(0.33) & 7.37(0.23) & 3.42(0.06) & 3.53(0.06)\\
        \botrule
      \end{tabular}}
    \label{LJ}
\end{table}

\subsubsection{Thermal conductivity}
Our thermal conductivity calculations yield comparable results to the
previous NIVS algorithm. This indicates that the thermal gradients
produced using this method are also quite close to previous RNEMD
methods. Simulations with moderately large thermal fluxes tend to
yield more reliable thermal gradients than under NIVS and thus avoid
large errors.  Extreme values for the applied thermal flux can
introduce side effects such as non-linear temperature gradients and
inadvertent phase transitions, so these are avoided.

Since the scaling operation is isotropic in VSS, the temperature
anisotropy that was observed under the earlier NIVS approach should be
absent.  Furthermore, the VSS method avoids unintended momentum flux
when only thermal flux is being imposed.  This was not always possible
with swapping or NIVS approaches. The thermal energy exchange in
swapping ($\vec{p}_i$ in slab {\it C} swapped with $\vec{p}_j$ in slab
{\it H}) or NIVS (total slab momentum components $P^\alpha$ scaled to
$\alpha P^\alpha$) do not achieve this unless thermal flux vanishes
(i.e. $\vec{p}_i=\vec{p}_j$ or $\alpha=1$). In this sense, the VSS
method achieves minimal perturbation to a simulation while imposing
thermal flux.

\subsubsection{Shear viscosity}
Table \ref{LJ} also compares our shear viscosity results with the
momentum swapping approach. Our calculations show that the VSS method
predicted similar values for shear viscosities to the momentum
swapping approach, as well as providing similar velocity gradient
profiles. Moderately large momentum fluxes are helpful in reducing the
errors in measured velocity gradients and thus the shear viscosity
values. However, it should be noted that the momentum swapping
approach tends to produce non-thermal velocity distributions in the
two slabs.\cite{Maginn:2010}

To show that the temperature is isotropic within each slab under VSS,
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 (shearing) velocities of the slabs. The one-dimensional
temperature profiles show no observable differences between the three
box dimensions. This ensures that VSS automatically preserves
temperature isotropy during shear viscosity calculations.

\begin{figure}
\includegraphics[width=\linewidth]{tempXyz}
\caption{Unlike the previous NIVS algorithm, the VSS method does not
  produce thermal anisotropy under shearing stress.  Temperature
  differences along the different box axes were significantly smaller
  than the error bars. Note that the two regions of elevated
  temperature are caused by the shear stress.  This is the same
  frictional heating reported previously by Tenney and
  Maginn.\cite{Maginn:2010}}
\label{tempXyz}
\end{figure}

The velocity distribution profiles were further tested by imposing a
large shear stress. Figure \ref{vDist} demonstrates how the VSS method
is able to maintain nearly ideal Maxwell-Boltzmann velocity
distributions even under large imposed fluxes. Previous swapping
methods tend to deplete particles with positive velocities in the
negative velocity slab ({\it C}) and deplete particles with negative
velocities in the {\it H} slab, leaving significant `notches' in the
velocity distributions.  The problematic velocity distributions become
significant when the imposed-flux becomes so large that diffusion from
neighboring slabs cannot offset the depletion.  Simultaneously,
abnormal peaks appear in the other regions of the velocity
distributions as high (or low) momentum particles are introduced in to
the corresponding slab.  These nonthermal distributions limit
application of the swapping approach in shear stress simulations. The
VSS method avoids the above problematic distributions by altering the
means of applying momentum flux. Comparatively, velocity distributions
recorded from simulations with the VSS method is so close to the ideal
Maxwell-Boltzmann distributions that no obvious difference is visible
in Figure \ref{vDist}. 

\begin{figure}
\centering
\includegraphics[width=5.5in]{velDist}
\caption{Velocity distributions that develop under VSS more closely
  resemble ideal Maxwell-Boltzmann distributions than earlier
  momentum-swapping RNEMD approaches. This data was 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
  which demonstrates some of the non-thermal velocity distributions
  that can develop under swapping RNEMD.}
\label{vDist}
\end{figure}

\subsection{Bulk SPC/E water}
We also computed the thermal conductivity and shear viscosity of SPC/E
water. A simulation cell with 1000 molecules was set up in a similar
manner to Ref. \cite{kuang:164101}. For thermal conductivity
calculations, 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 equilibrium molecular dynamics
calculations.\cite{10.1063/1.3330544,Medina2011} Additionally, a
simulation with {\it both} applied thermal and momentum gradients was
carried out to map out shear viscosity as a function of temperature.

\subsubsection{Thermal conductivity}
Table \ref{spceThermal} summarizes the thermal conductivity of SPC/E
under different temperatures and thermal gradients compared with 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 the VSS method
was applied, which indicates that our algorithm conserves total energy
well for systems involving electrostatic interactions.

Measurements using the VSS method established similar temperature
gradients to the previous NIVS method. Our simulation results are in
fairly good agreement with those from previous simulations. Both
methods yield values in reasonable agreement with experimental
values. Simulations using larger thermal gradients and those with
longer gradient axes ($z$) have less measurable noise.

\begin{table}
  
  \tbl{Thermal conductivity of SPC/E water under various
    imposed thermal gradients. Uncertainties are indicated in
    parentheses.}
  
  {\begin{tabular}{cc|ccc}\toprule
        $\langle T\rangle$ & $\langle dT/dz\rangle$ & \multicolumn{3}{c}
        {$\lambda (\mathrm{W m}^{-1} \mathrm{K}^{-1})$} \\
        (K) & (K/\AA) & This work &  NIVS\cite{kuang:164101} &
        Experiment\cite{WagnerKruse} \\
        \colrule
        300 & 0.8 & 0.815(0.027)     & 0.770(0.008) & 0.61 \\ \colrule
        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)$^{\rm a}$ &              &      \\
        \botrule
      \end{tabular}}
\tabnote{$^{\rm a}$Simulation with 2000 SPC/E molecules}
\label{spceThermal}
\end{table}

\subsubsection{Shear viscosity}
Shear viscosity was computed for SPC/E water at a range of
temperatures that span the liquidus range for water under atmospheric
pressure.  VSS-RNEMD simulations predict a similar trend of $\eta$
vs. $T$ to equilibrium molecular dynamics (EMD) results and are
presented in table \ref{spceShear}.  Our results show no significant
differences from the earlier EMD calculations.  Since each value
reported using our method takes a single 1.5 ns trajectory, instead of
average from many trajectories when using EMD, the VSS method provides
an efficient means for computing the shear viscosity.

\begin{table}
      \tbl{Computed shear viscosity of SPC/E water under different
        temperatures. Results are compared to those obtained with
        equilibrium molecular dynamics.\cite{Medina2011} Uncertainties
        are indicated in parentheses.}
      
      {\begin{tabular}{cc|cc}\toprule
        $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} \\
        \colrule
        273 & 1.12 & 1.218(0.004) & 1.282(0.048) \\
            & 1.79 & 1.140(0.012) &              \\
        \colrule
        303 & 2.09 & 0.646(0.008) & 0.643(0.019) \\
        \colrule
        318 & 2.50 & 0.536(0.007) &              \\
            & 5.25 & 0.510(0.007) &              \\
            & 2.82 & 0.474(0.003)$^{\rm a}$ &              \\
        \colrule
        333 & 3.10 & 0.428(0.002) & 0.421(0.008) \\
        \colrule
        363 & 2.34 & 0.279(0.014) & 0.291(0.005) \\
            & 4.26 & 0.306(0.001) &              \\
        \botrule
      \end{tabular}}
\tabnote{$^{\rm a}$Simulation with 2000 SPC/E molecules}\label{spceShear}
\end{table}

A much more efficient way to map out $\eta$ vs $T$ is to combine the
momentum flux with a thermal flux. Figure \ref{Tvxdvdz} shows the
thermal and velocity gradient in one such simulation. The local
temperature varies along the $z$-axis, and the velocity gradient can
be computed locally via finite difference approximations. With the
data provided in Figure \ref{Tvxdvdz}, the entire temperature
dependent viscosity $\eta(T)$ can be calculated from a {\it single} 2
ns trajectory.  This data is shown in figure \ref{etaT}.  A linear fit
can be used to approximate $\partial v_x/\partial z$ vs. $z$ and the
fit velocity gradient produced the solid curve in Fig. \ref{etaT}.

\begin{figure}
\centering
\includegraphics[width=5.5in]{tvxdvdz}
\caption{With a combination of a thermal and a momentum flux, a
  simulation can exhibit both a temperature (top) and a velocity
  (middle) gradient. $\partial v_x/\partial z$ depends on the local
  temperature, so it varies along the gradient axis.  This derivative
  can be computed using finite difference approximations (lower) and
  can be used to calculate $\eta$ vs $T$ (Figure \ref{etaT}).}
\label{Tvxdvdz}
\end{figure}

In Figure \ref{etaT}, the fit gradient yields viscosity values that
are in agreement with VSS-RNEMD simulations carried out at fixed
temperatures, as well as results reported using EMD methods through
most of the simulated temperature
range.\cite{10.1063/1.3330544,Medina2011} However, there is larger
error in lower temperature regions and the single-trajectory approach
has difficulty predicting $\eta$ near 273K.  Longer trajectories and
larger box geometries would be required to converge these values with
single-temperature calculations. Since previous work already pointed
out that SPC/E water tends to yield lower viscosities than
experimental data,\cite{Medina2011} experimental comparisons are not
given here.

\begin{figure}
\includegraphics[width=\linewidth]{etaT}
\caption{The temperature dependence of the shear viscosity generated
  by a {\it single} VSS-RNEMD simulation.  Applying both thermal and
  momentum gradient predicts reasonable values in much of the
  temperature range being tested.}
\label{etaT}
\end{figure}

\subsection{Interfacial friction and slip length}
Another attractive aspect of the VSS-RNEMD method is the ability to
apply shear stress to interfacial and heterogeneous systems, where
molecules of different identities (or phases) are segregated in
different regions of space. We have previously studied the interfacial
thermal conductivity of a series of metal-liquid
interfaces,\cite{kuang:164101,kuang:AuThl} and would like to further
investigate the relationship between this phenomenon and the
interfacial friction.  Knowing the friction coefficients and slip
lengths of interfaces could also help predict the mass- and
heat-transport properties of colloidal suspensions of nanocrystalline
materials.

Table \ref{etaKappaDelta} includes results of a few model interfaces.
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. This 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}
  \tbl{Computed interfacial friction coefficients and slip
    lengths for various solid -- liquid interfaces.  All of the
    solvated bare metal surfaces appear to exist in the ``slip''
    boundary conditions under shear stress, while the ice/water
    interface is in the opposite or ``stick'' boundary
    limit. Error estimates are indicated in parentheses.}
  
  {\begin{tabular}{ccccccc}
      \toprule
        Solid & Liquid & $T$ & $j_z(p_x)$ & $\eta_{liquid}$ & $\kappa$
        & $\delta$  \\
        surface & phase & K & MPa & mPa$\cdot$s &
        10$^4$Pa$\cdot$s/m & nm  \\
        \colrule
        Au(111) & hexane & 200 & 1.08 & 0.197(0.009) & 5.30(0.36) &
        3.72(0.25)  \\
                &        &     & 2.15 & 0.141(0.002) & 5.31(0.26) &
        2.76(0.14)     \\
        \colrule
        Au(111) & toluene & 200 & 1.08 & 0.722(0.035) & 15.7(0.7) &
        4.60(0.22) \\
                &         &     & 2.16 & 0.544(0.030) & 11.2(0.5) &
        4.86(0.27)      \\
        \colrule
        Au(111) & water & 300 & 1.08 & 0.399(0.050) & 1.928(0.022) &
        20.7(2.6)   \\
                &       &     & 2.16 & 0.794(0.255) & 1.895(0.003) &
        41.9(13.5)       \\
        \colrule
        ice Ih (basal) & water & 225 & 19.4 & 15.8(0.2) & $\infty$ & 0  \\
        \botrule
      \end{tabular}}
    \label{etaKappaDelta}
\end{table}

In addition to these previously studied interfaces, we also
constructed ice Ih (basal face) -- water interfaces. In contrast with
the Au(111) -- water interface, where the friction coefficient is
small and large slip effect is present, the ice -- water interface
demonstrates strong solid-liquid interactions and appears to be in the
stick boundary limit. The supercooled liquid phase is an order of
magnitude more viscous than measurements of liquid water done in the
previous section. It would be of interst to investigate the effects of
different ice lattice planes (e.g. the prism, secondary prism, and
pyramidal faces of ice) on interfacial friction and the corresponding
liquid viscosity.

\section{Conclusions}
We have developed a novel method for simultaneously applying both
shear stress and thermal flux during molecular dynamics
simulations. This method, the VSS-RNEMD approach, enables efficient
computation of both thermal conductivity and shear viscosity in bulk
liquids. The method maintains thermal velocity distributions and
avoids thermal anisotropy in previous NIVS shear stress simulations,
while retaining the attractive features of previous RNEMD
methods. There is no {\it a priori} restrictions on the method to
particular statistical ensembles, so use of this method within
extended-system simulations (to maintain constant temperature and
pressure) is also possible.

Most importantly, the VSS-RNEMD method is capable of simultaneously
imposing thermal flux and shear stress accross a heterogeneous
interface. This will facilitate the study of interfacial properties
that depend on the chemical details of the interface. Further
investigations will be carried out to characterize interfacial
friction behavior using the method.

Since VSS-RNEMD can simultaneously introduce thermal and momentum
gradients, it can also be used to efficiently map out the viscosity as
a function of temperature with a single simulation.  Complex systems
that involve thermal and momentum gradients might potentially benefit
from this feature. For example, the behavior of the Soret effect
(solute migration across a thermal gradient) in a solution flowing
through a medium would be of interest to those simulating purification
and separation processes.

\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:	3791
Committed:	Wed Mar 21 17:34:37 2012 UTC (12 years, 5 months ago) by gezelter
Content type:	application/x-tex
File size:	42647 byte(s)
Log Message:	submitted version