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{ENTER TITLE HERE}

\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 and maintains
  thermal velocity distributions. It also avoid thermal anisotropy
  occured in 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.

\end{abstract}

\newpage

%\narrowtext

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

\section{Introduction}
[REFINE LATER, ADD MORE REF.S]
Imposed-flux methods in Molecular Dynamics (MD)
simulations\cite{MullerPlathe:1997xw} can establish steady state
systems with a set applied flux 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
or bulk systems as well as heterogeneous systems such as liquid-solid
interfaces.\cite{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 by
M\"{u}ller-Plathe {\it et
  al.}\cite{MullerPlathe:1997xw,ISI:000080382700030} utilize simple
momentum swapping for generating energy/momentum fluxes, which is also
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 to kinetic energy
transfer between particles of different identities is less efficient
when the mass difference between the particles becomes significant,
which also limits its application 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 limits its application in imposing momentum
fluxes. Temperature anisotropy can happen under high momentum fluxes,
due to the nature of the algorithm. Thus, combining thermal and
momentum flux is also difficult to implement with this
approach. However, such combination may provide a means to simulate
thermal/momentum gradient coupled processes such as freeze
desalination. Therefore, developing novel approaches to extend the
application of imposed-flux method is desired.

In this paper, we improve the NIVS method and propose 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 to implement the method 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:
[SUPPORT INFO OR APPENDIX MIGHT BE NECESSARY TO PUT EXTRA MATH IN]
\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.

The advantages of this approach over the original momentum swapping
one lies in its nature of preserve a Maxwell-Boltzmann distribution
(mathimatically a Gaussian function). However, because the momentum
swapping steps tend to result in a nonthermal distribution, when an
imposed flux is relatively large, diffusions from the neighboring
slabs could no longer remedy this effect, problematic distributions
would be observed. Results in later section will illustrate this
effect in more detail.

\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 a orthorhombic simulation
cell with periodic boundary conditions in all three dimensions. The
$z$ axis of these cells were longer and was set as the gradient axis
of temperature and/or momentum. Thus the cells were divided into $N$
slabs along this axis, with various $N$ depending on individual
system. The $x$ and $y$ axis were usually of the same length in
homogeneous systems or close to each other where 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 to prepare systems ready for data
collections. Isobaric-isothermal equilibrations are performed before
this for SPC/E water systems to reach normal pressure (1 bar), while
similar equilibrations are used for interfacial systems to relax the
surface tensions.

While homogeneous fluid systems can be set up with random
configurations, our interfacial systems needs extra 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
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 different proton orders. We refer to their results and
choose the configuration of the lowest energy after geometry
optimization as the unit cells of our ice lattices. Although
experimental solid/liquid coexistant temperature near normal pressure
is 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. Therefore, all our ice/liquid
water simulations were carried out under 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 FIGURE HERE FOR BASAL PLANE]

\subsection{Force Field Parameters}
For comparison of our new method with previous work, we retain our
force field parameters consistent with the results we will compare
with. The Lennard-Jones fluid used here for argon , and reduced unit
results are reported 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.

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 present for the heat transport occurs in
both $+z$ and $-z$ directions. The 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 switching off $J_z$. One can specify the vector
$\vec{j}_z(\vec{p})$ by choosing the three components
respectively. For shear viscosity simulations, $j_z(p_z)$ is usually
set to zero. Although for isotropic systems, the direction of
$\vec{j}_z(\vec{p})$ on the $xy$ plane should not matter, the ability
of arbitarily specifying the vector direction in our method could
provide convenience in anisotropic simulations.

Similar to thermal conductivity computations, 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 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 reported in reduced units
(${\eta^*=\eta\sigma^2(\varepsilon m)^{-1/2}}$).

\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) and 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 condition at the solid-liquid interface, $\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 A 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. An example of
  Au/hexane interfaces is shown. Calculation for the left side is 
  illustrated. The right side is similar to the left side.}
\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 & 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 rendered 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 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 obtain this result unless thermal flux vanishes
(i.e. $\vec{p}_i=\vec{p}_j$ or $\alpha=1$ which are trivial to apply a
thermal flux). In this sense, this method contributes to having
minimal perturbation to a simulation while imposing thermal flux.

\subsubsection{Shear viscosity}
Table \ref{LJ} also compares our shear viscosity results with momentum
swapping approach. Our calculations show that our method predicted
similar values for 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 effects 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 a notch. This
problematic profiles become significant when the imposed-flux becomes
larger and diffusions from neighboring slabs could not offset the
depletion. Simutaneously, abnormal peaks appear corresponding to
excessive velocity swapped from the other slab. This 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 observable difference is shown in Figure
\ref{vDist}. Conclusively, our method avoids problems happened 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 high 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. Distributions produced by our method are very close
  to the ideal thermal situations.}
\label{vDist}
\end{figure}

\subsection{Bulk SPC/E water}
Since our method was in good performance of thermal conductivity and
shear viscosity computations for simple Lennard-Jones fluid, we extend
our applications of these simulations to complex fluid like SPC/E
water model. A simulation cell with 1000 molecules was set up in the
same 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 NIVS REF 32]), and results were
compared to available data from Equilibrium MD methods[CITATIONS].

\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 even 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[CITATION]\\
        \hline
        273 &  & 1.218(0.004) &  \\
            &  & 1.140(0.012) &  \\
        303 &  & 0.646(0.008) &  \\
        318 &  & 0.536(0.007) &  \\
            &  & 0.510(0.007) &  \\
            &  &  &  \\
        333 &  & 0.428(0.002) &  \\
        363 &  & 0.279(0.014) &  \\
            &  & 0.306(0.001) &  \\
        \hline\hline
      \end{tabular}
      \label{spceShear}
    \end{center}
  \end{minipage}
\end{table*}

[MAY COMBINE JZPX AND JZKE TO RUN ONE JOB BUT GET ETA=F(T)]
[PUT RESULTS AND FIGURE HERE IF IT WORKS]
\subsection{Interfacial frictions and slip lengths}
An 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 attemptions have been
made to 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. Instead, the butanethiol
covered Au(111) surface appeared to be sticky to the organic liquid
molecules 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 BARRAT, GARDE]

\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 & Pa$\cdot$s/m & nm
        & MW/m$^2$/K \\
        \hline
        Au(111) & hexane & 200 & 1.08 & 0.20() & 5.3$\times$10$^4$() &
        3.7 & 46.5 \\
                &        &     & 2.15 & 0.14() & 5.3$\times$10$^4$() &
        2.7 &      \\
        Au-SC$_4$H$_9$ & hexane & 200 & 2.16 & 0.29() & $\infty$ & 0 &
        131 \\
                       &        &     & 5.39 & 0.32() & $\infty$ & 0 &
            \\
        \hline
        Au(111) & toluene & 200 & 1.08 & 0.72() & 1.?$\times$10$^5$() &
        4.6 & 70.1 \\
                &         &     & 2.16 & 0.54() & 1.?$\times$10$^5$() &
        4.9 &      \\
        Au-SC$_4$H$_9$ & toluene & 200 & 5.39 & 0.98() & $\infty$ & 0
        & 187 \\
                       &         &     & 10.8 & 0.99() & $\infty$ & 0
        &     \\
        \hline
        Au(111) & water & 300 & 1.08 & 0.40() & 1.9$\times$10$^4$() &
        20.7 & 1.65 \\
                &       &     & 2.16 & 0.79() & 1.9$\times$10$^4$() &
        41.9 &      \\
        \hline
        ice(basal) & water & 225 & 19.4 & 15.8() & $\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. Note that $\eta$ measured near a ``slip'' surface tends
to be smaller than that near a ``stick'' surface. This suggests that
an interface could affect the dynamic properties on 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 SCHMIDT AND SKINNER] Our observations could provide
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
first studied. In contrast to the Au(111)/water interface, where the
friction coefficient is relatively small and large slip effect
presents, the ice/liquid water interface demonstrates strong
interactions and appears to be sticky. The supercooled liquid phase is
an order of magnitude viscous than measurements in previous
section. It would be of interst to investigate the effect of different
ice lattice planes (such as prism surface) on interfacial friction and
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.

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

Another attractive feature of our method is the ability of
simultaneously imposing thermal and momentum flux in a
system. potential researches that might be benefit include complex
systems that involve thermal and momentum gradients. For example, the
Soret effects under a velocity gradient would be 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:	3782
Committed:	Mon Dec 12 19:39:38 2011 UTC (12 years, 9 months ago) by skuang
Content type:	application/x-tex
File size:	36821 byte(s)
Log Message:	add a figure, plus a few edits.