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{A minimal perturbation approach to RNEMD able to simultaneously
impose thermal 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 the linear momentum
  and total energy of the system and improves previous Reverse
  Non-Equilibrium Molecular Dynamics (RNEMD) methods while maintaining
  thermal velocity distributions. It also avoid thermal anisotropy
  occured in previous NIVS simulations by using isotropic velocity
  scaling on the molecules in specific regions of a system. To test
  the method, we have computed the thermal conductivity and shear
  viscosity of model liquid systems as well as the interfacial
  frictions of a series of metal/liquid interfaces. Its ability to
  combine the thermal and momentum gradients allows us to obtain shear
  viscosity data for a range of temperatures in only one trajectory.

\end{abstract}

\newpage

%\narrowtext

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

\section{Introduction}
Imposed-flux methods in Molecular Dynamics (MD)
simulations\cite{MullerPlathe:1997xw,ISI:000080382700030,kuang:164101}
can establish steady state systems with an applied flux set vs a
corresponding gradient that can be measured. These methods does not
need many trajectories to provide information of transport properties
of a given system. Thus, they are utilized in computing thermal and
mechanical transfer of homogeneous bulk systems as well as
heterogeneous systems such as solid-liquid
interfaces.\cite{garde:nl2005,garde:PhysRevLett2009,kuang:AuThl}

The Reverse Non-Equilibrium MD (RNEMD) methods adopt constraints that
satisfy linear momentum and total energy conservation of a system when
imposing fluxes in a simulation. Thus they are compatible with various
ensembles, including the micro-canonical (NVE) ensemble, without the
need of an external thermostat. The original approaches proposed by
M\"{u}ller-Plathe {\it et
  al.}\cite{MullerPlathe:1997xw,ISI:000080382700030} utilize simple
momentum swapping for generating energy/momentum fluxes, which can
also be compatible with particles of different identities. Although
simple to implement in a simulation, this approach can create
nonthermal velocity distributions, as discovered by Tenney and
Maginn\cite{Maginn:2010}. Furthermore, this approach is less efficient
for kinetic energy transfer between particles of different identities,
especially when the mass difference between the particles becomes
significant. This also limits its applications on heterogeneous
interfacial systems.

Recently, we developed a different approach, using Non-Isotropic
Velocity Scaling (NIVS) \cite{kuang:164101} algorithm to impose
fluxes. Compared to the momentum swapping move, it scales the velocity
vectors in two separate regions of a simulated system with respective
diagonal scaling matrices. These matrices are determined by solving a
set of equations including linear momentum and kinetic energy
conservation constraints and target flux satisfaction. This method is
able to effectively impose a wide range of kinetic energy fluxes
without obvious perturbation to the velocity distributions of the
simulated systems, regardless of the presence of heterogeneous
interfaces. We have successfully applied this approach in studying the
interfacial thermal conductance at metal-solvent
interfaces.\cite{kuang:AuThl}

However, the NIVS approach has limited applications in imposing
momentum fluxes. Temperature anisotropy could happen under high
momentum fluxes due to the implementation of this algorithm. Thus,
combining thermal and momentum flux is also difficult to obtain with
this approach. However, such combination may provide a means to
simulate thermal/momentum gradient coupled processes such as Soret
effect in liquid flows. Therefore, developing improved approaches to
extend the applications of the imposed-flux method is desirable.

In this paper, we improve the RNEMD methods by proposing a novel
approach to impose fluxes. This approach separate the means of applying
momentum and thermal flux with operations in one time step and thus is
able to simutaneously impose thermal and momentum flux. Furthermore,
the approach retains desirable features of previous RNEMD approaches
and is simpler to implement compared to the NIVS method. In what
follows, we first present the method and its implementation in a
simulation. Then we compare the method on bulk fluids to previous
methods. Also, interfacial frictions are computed for a series of
interfaces.

\section{Methodology}
Similar to the NIVS method,\cite{kuang:164101} we consider a
periodic system divided into a series of slabs along a certain axis
(e.g. $z$). The unphysical thermal and/or momentum flux is designated
from the center slab to one of the end slabs, and thus the thermal
flux results in a lower temperature of the center slab than the end
slab, and the momentum flux results in negative center slab momentum
with positive end slab momentum (unless these fluxes are set
negative). Therefore, the center slab is denoted as ``$c$'', while the
end slab as ``$h$''.

To impose these fluxes, we periodically apply different set of
operations on velocities of particles {$i$} within the center slab and
those of particles {$j$} within the end 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$ denotes
the instantaneous bulk velocity of slabs $c$ and $h$ respectively
before an operation is applied. When a momentum flux $\vec{j}_z(\vec{p})$
presents, these bulk velocities would have a corresponding change
($\vec{a}_c$ and $\vec{a}_h$ respectively) according to Newton's
second law:
\begin{eqnarray}
M_c \vec{a}_c & = & -\vec{j}_z(\vec{p}) \Delta t \\
M_h \vec{a}_h & = &  \vec{j}_z(\vec{p}) \Delta t
\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 above operations 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 equations are
included as well:
[MAY PUT EXTRA MATH IN SUPPORT INFO OR APPENDIX]
\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
\end{eqnarray}
where $K_c$ and $K_h$ denotes translational kinetic energy of slabs
$c$ and $h$ respectively before an operation is applied. These
translational kinetic energy conservation equations are sufficient to
ensure total energy conservation, as the operations applied in our
method do not change the kinetic energy related to other degrees of
freedom or the potential energy of a system, given that its potential
energy does not depend on particle velocity.

The above sets of equations 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 two roots respectively for $c$ and
$h$. However, the positive roots (which are closer to 1) are chosen so
that the perturbations to a system can be reduced to a minimum. Figure
\ref{method} illustrates the implementation sketch of this algorithm
in an individual step.

\begin{figure}
\includegraphics[width=\linewidth]{method}
\caption{Illustration of the implementation of the algorithm in a
  single step. Starting from an ideal velocity distribution, 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, thermal distributions can preserve after this operation.}
\label{method}
\end{figure}

By implementing these operations at a certain frequency, a steady
thermal and/or momentum flux can be applied and the corresponding
temperature and/or momentum gradients can be established.

Compared to the previous NIVS method, this approach is computationally
more efficient in that only quadratic equations are involved to
determine a set of scaling coefficients, while the NIVS method needs
to solve quartic equations. Furthermore, this method implements
isotropic scaling of velocities in respective slabs, unlike the NIVS,
where an extra criteria function is necessary to choose a set of
coefficients that performs a scaling as isotropic as possible. More
importantly, separating the means of momentum flux imposing from
velocity scaling avoids the underlying cause to thermal anisotropy in
NIVS when applying a momentum flux. And later sections will
demonstrate that this can improve the performance in shear viscosity
simulations.

This approach is advantageous over the original momentum swapping in
many aspects. In one swapping, the velocity vectors involved are
usually very different (or the generated flux is trivial to obtain
gradients), thus the swapping tends to incur perturbations to the
neighbors of the particles involved. Comparatively, our approach
disperse the flux to every selected particle in a slab so that
perturbations in the flux generating region could be
minimized. Additionally, because the momentum swapping steps tend to
result in a nonthermal distribution, when an imposed flux is
relatively large and diffusions from the neighboring slabs could no
longer remedy this effect, problematic distributions would be
observed. In comparison, the operations of our approach has the nature
of preserving the equilibrium velocity distributions (commonly
Maxwell-Boltzmann), and results in later section will illustrate that
this is helpful to retain thermal distributions in a simulation.

\section{Computational Details}
The algorithm has been implemented in our MD simulation code,
OpenMD\cite{Meineke:2005gd,openmd}. We compare the method with
previous RNEMD methods or equilibrium MD (EMD) methods in homogeneous fluids
(Lennard-Jones and SPC/E water). And taking advantage of the method,
we simulate the interfacial friction of different heterogeneous
interfaces (gold-organic solvent and gold-SPC/E water and ice-liquid
water).

\subsection{Simulation Protocols}
The systems to be investigated are set up in orthorhombic simulation
cells with periodic boundary conditions in all three dimensions. The
$z$ axis of these cells were longer and set as the temperature and/or
momentum gradient axis. And the cells were evenly divided into $N$
slabs along this axis, with various $N$ depending on individual
system. The $x$ and $y$ axis were of the same length in homogeneous
systems or had length scale close to each other where heterogeneous
interfaces presents. In all cases, before introducing a nonequilibrium
method to establish steady thermal and/or momentum gradients for
further measurements and calculations, canonical ensemble with a
Nos\'e-Hoover thermostat\cite{hoover85} and microcanonical ensemble
equilibrations were used before data collections. For SPC/E water
simulations, isobaric-isothermal equilibrations\cite{melchionna93} are
performed before the above to reach normal pressure (1 bar); for
interfacial systems, similar equilibrations are used to relax the
surface tensions of the $xy$ plane.

While homogeneous fluid systems can be set up with rather random
configurations, our interfacial systems needs a series of steps to
ensure the interfaces be established properly for computations. The
preparation and equilibration of butanethiol covered gold (111)
surface and further solvation and equilibration process is described
in details as in reference \cite{kuang:AuThl}.

As for the ice/liquid water interfaces, the basal surface of ice
lattice was first constructed. Hirsch {\it et
  al.}\cite{doi:10.1021/jp048434u} explored the energetics of ice
lattices with all possible proton order configurations. We refer to
their results and choose the configuration of the lowest energy after
geometry optimization as the unit cell for our ice lattices. Although
experimental solid/liquid coexistant temperature under normal pressure
should be 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
225K. To have extra protection of the ice lattice during initial
equilibration (when the randomly generated liquid phase configuration
could release large amount of energy in relaxation), restraints were
applied to the ice lattice to avoid inadvertent melting by the heat
dissipated from the high enery configurations.
[MAY ADD A SNAPSHOT FOR BASAL PLANE]

\subsection{Force Field Parameters}
For comparison of our new method with previous work, we retain our
force field parameters consistent with previous simulations. Argon is
the Lennard-Jones fluid used here, and its results are reported in
reduced unit for direct comparison purpose.

As for our water simulations, SPC/E model is used throughout this work
for consistency. Previous work for transport properties of SPC/E water
model is available\cite{Bedrov:2000,10.1063/1.3330544,Medina2011} so
that unnecessary repetition of previous methods can be avoided.

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 and liquid hexane and toluene. The United-Atom
models\cite{TraPPE-UA.thiols,TraPPE-UA.alkanes,TraPPE-UA.alkylbenzenes}
for these components were used in this work for better 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 $\vec{j}_z(\vec{p})$ is set to zero and a target $J_z$ is set to
impose kinetic energy transfer, the method can be used for thermal
conductivity computations. Similar to previous RNEMD methods, we
assume linear response of the temperature gradient with respect to the
thermal flux in general case. And the thermal conductivity ($\lambda$)
can be obtained with 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 is recorded throughout a simulation, and
the time for data collection $t$:
\begin{equation}
J_z = \frac{E_{total}}{2 t L_x L_y}
\end{equation}
where $L_x$ and $L_y$ denotes the dimensions of the plane in a
simulation cell perpendicular to the thermal gradient, and a factor of
two in the denominator is necessary for the heat transport occurs in
both $+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, the method can carry out shear viscosity calculations
by specify a momentum flux. In our algorithm, one can specify the
three components of the flux vector $\vec{j}_z(\vec{p})$
respectively. For shear viscosity simulations, $j_z(p_z)$ is usually
set to zero. For isotropic systems, the direction of
$\vec{j}_z(\vec{p})$ on the $xy$ plane should not matter, but the
ability of arbitarily specifying the vector direction in our method
could provide convenience in anisotropic simulations.

Similar to thermal conductivity computations, for a homogeneous
system, 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 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 at
a certain temperature, our method's ability to impose both a thermal
and a momentum flux in one simulation allows the combination of a
temperature and a velocity gradient. In this case, since viscosity is
generally a function of temperature, the local viscosity also depends
on the local temperature. Therefore, in one such 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). As a whole, the
viscosity can be mapped out as the function of temperature in one
single trajectory of simulation. Results for shear viscosity
computations of SPC/E water will demonstrate its effectiveness in
detail.

\subsection{Interfacial friction and Slip length}
While the shear stress results in a velocity gradient within bulk
fluid phase, its effect at a solid-liquid interface could vary due to
the interaction strength between the two phases. The interfacial
friction coefficient $\kappa$ is defined to relate the shear stress
(e.g. along $x$-axis) with the relative fluid velocity tangent to the
interface:
\begin{equation}
j_z(p_x)|_{interface} = \kappa\Delta v_x|_{interface}
\end{equation}
Under ``stick'' boundary condition, $\Delta v_x|_{interface}
\rightarrow 0$, which leads to $\kappa\rightarrow\infty$. However, for
``slip'' boundary conditions at the solid-liquid interfaces, $\kappa$
becomes finite. To characterize the interfacial boundary conditions,
slip length ($\delta$) is defined using $\kappa$ and the shear
viscocity of liquid phase ($\eta$):
\begin{equation}
\delta = \frac{\eta}{\kappa}
\end{equation}
so that $\delta\rightarrow 0$ in the ``no-slip'' boundary condition,
and depicts how ``slippery'' an interface is. Figure \ref{slipLength}
illustrates how this quantity is defined and computed for a
solid-liquid interface. [MAY INCLUDE SNAPSHOT IN FIGURE]

\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. An example of Au/hexane interfaces is shown, and the
  calculation for the left side is illustrated. The calculation for
  the right side is similar to the left.}
\label{slipLength}
\end{figure}

In our method, a shear stress can be applied similar to shear
viscosity computations by applying an unphysical momentum flux
(e.g. $j_z(p_x)$). A corresponding velocity profile can be obtained as
shown in Figure \ref{slipLength}, in which the velocity gradients
within liquid phase and velocity difference at the liquid-solid
interface can be measured respectively. Further calculations and
characterizations of the interface can be carried out using these
data.

\section{Results and Discussions}
\subsection{Lennard-Jones fluid}
Our orthorhombic simulation cell of Lennard-Jones fluid has identical
parameters to our previous work\cite{kuang:164101} to facilitate
comparison. Thermal conductivitis and shear viscosities were computed
with the algorithm applied to the simulations. The results of thermal
conductivity are compared with our previous NIVS algorithm. However,
since the NIVS algorithm could produce temperature anisotropy for
shear viscocity computations, these results are instead compared to
the momentum swapping approaches. Table \ref{LJ} lists these
calculations 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 $t^*$ & 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 to achieve
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:	3785
Committed:	Tue Dec 13 23:31:35 2011 UTC (12 years, 6 months ago) by skuang
Content type:	application/x-tex
File size:	41347 byte(s)
Log Message:	a revised version ready