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}
  REPLACE ABSTRACT HERE
  With the Non-Isotropic Velocity Scaling (NIVS) approach to Reverse
  Non-Equilibrium Molecular Dynamics (RNEMD) it is possible to impose
  an unphysical thermal flux between different regions of
  inhomogeneous systems such as solid / liquid interfaces.  We have
  applied NIVS to compute the interfacial thermal conductance at a
  metal / organic solvent interface that has been chemically capped by
  butanethiol molecules.  Our calculations suggest that coupling
  between the metal and liquid phases is enhanced by the capping
  agents, leading to a greatly enhanced conductivity at the interface.
  Specifically, the chemical bond between the metal and the capping
  agent introduces a vibrational overlap that is not present without
  the capping agent, and the overlap between the vibrational spectra
  (metal to cap, cap to solvent) provides a mechanism for rapid
  thermal transport across the interface. Our calculations also
  suggest that this is a non-monotonic function of the fractional
  coverage of the surface, as moderate coverages allow diffusive heat
  transport of solvent molecules that have been in close contact with
  the capping agent.

\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 methodology,\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 center slab
would have a lower temperature than the end slab (unless the thermal
flux is negative). Therefore, the center slab is denoted as ``$c$''
while the end slab as ``$h$''.

To impose these fluxes, we periodically apply separate operations to
velocities of particles {$i$} within the center slab and 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 occurs. 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
\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 operations.

The above operations conserve the linear momentum of a periodic
system. To satisfy total energy conservation as well as to impose a
thermal flux $J_z$, one would have
[SUPPORT INFO 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 occurs. These
translational kinetic energy conservation equations are sufficient to
ensure total energy conservation, as the operations applied do not
change the potential energy of a system, given that the 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 two roots of $c$ and $h$ exist
respectively. However, to avoid dramatic perturbations to a system,
the positive roots (which are closer to 1) are chosen. Figure
\ref{method} illustrates the implementation 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 both thermal and momentum flux from
  the ``c'' slab to the ``h'' slab. As the figure shows, the thermal
  distributions 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.

This approach is more computationaly efficient compared to the
previous NIVS method, in that only quadratic equations are involved,
while the NIVS method needs to solve a quartic equations. Furthermore,
the 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 the most isotropic
scaling. More importantly, separating the momentum flux imposing from
velocity scaling avoids the underlying cause that NIVS produced
thermal anisotropy when applying a momentum flux.

The advantages of the approach over the original momentum swapping
approach lies in its nature to preserve a Gaussian
distribution. Because the momentum swapping tends to render a
nonthermal distribution, when the imposed flux is relatively large,
diffusion of the neighboring slabs could no longer remedy this effect,
and nonthermal distributions would be observed. Results in later
section will illustrate this effect.

\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 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), a constraint method (REF?) was
adopted until the high energy configuration was relaxed.
[MAY ADD A FIGURE HERE FOR BASAL PLANE, MAY INCLUDE PRISM IF POSSIBLE]

\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 provides
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. An example of Au/hexane
  interfaces is shown.}
\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 & 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)$^a$ &              &      \\
        \hline\hline
      \end{tabular}
      $^a$Simulation with $L_z$ twice as long.
      \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}). 

\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}
[SLIP BOUNDARY VS STICK BOUNDARY]

qualitative agreement w interfacial thermal conductance

[ETA OBTAINED HERE DOES NOT NECESSARILY EQUAL TO BULK VALUES]


[ATTEMPT TO CONSTRUCT BASAL PLANE ICE-WATER INTERFACE]

[FUTURE WORK HERE OR IN CONCLUSIONS]


\begin{table*}
  \begin{minipage}{\linewidth}
    \begin{center}
      
      \caption{Computed interfacial thermal conductance ($G$ and
        $G^\prime$) values for interfaces using various models for
        solvent and capping agent (or without capping agent) at
        $\langle T\rangle\sim$200K.  Here ``D'' stands for deuterated
        solvent or capping agent molecules. Error estimates are
        indicated in parentheses.}
      
      \begin{tabular}{llccc}
        \hline\hline
        Butanethiol model & Solvent & $G$ & $G^\prime$ \\
        (or bare surface) & model & \multicolumn{2}{c}{(MW/m$^2$/K)} \\
        \hline
        UA    & UA hexane    & 131(9)    & 87(10)    \\
              & UA hexane(D) & 153(5)    & 136(13)   \\
              & AA hexane    & 131(6)    & 122(10)   \\
              & UA toluene   & 187(16)   & 151(11)   \\
              & AA toluene   & 200(36)   & 149(53)   \\
        \hline
        bare  & UA hexane    & 46.5(3.2) & 49.4(4.5) \\
              & UA hexane(D) & 43.9(4.6) & 43.0(2.0) \\
              & AA hexane    & 31.0(1.4) & 29.4(1.3) \\
              & UA toluene   & 70.1(1.3) & 65.8(0.5) \\
        \hline\hline
      \end{tabular}
      \label{modelTest}
    \end{center}
  \end{minipage}
\end{table*}

On bare metal / solvent surfaces, different force field models for
hexane yield similar results for both $G$ and $G^\prime$, and these
two definitions agree with each other very well. This is primarily an
indicator of weak interactions between the metal and the solvent.

For the fully-covered surfaces, the choice of force field for the
capping agent and solvent has a large impact on the calculated values
of $G$ and $G^\prime$.  The AA thiol to AA solvent conductivities are
much larger than their UA to UA counterparts, and these values exceed
the experimental estimates by a large measure.  The AA force field
allows significant energy to go into C-H (or C-D) stretching modes,
and since these modes are high frequency, this non-quantum behavior is
likely responsible for the overestimate of the conductivity.  Compared
to the AA model, the UA model yields more reasonable conductivity
values with much higher computational efficiency.

\subsubsection{Effects due to average temperature}

We also studied the effect of average system temperature on the
interfacial conductance.  The simulations are first equilibrated in
the NPT ensemble at 1 atm.  The TraPPE-UA model for hexane tends to
predict a lower boiling point (and liquid state density) than
experiments.  This lower-density liquid phase leads to reduced contact
between the hexane and butanethiol, and this accounts for our
observation of lower conductance at higher temperatures.  In raising
the average temperature from 200K to 250K, the density drop of
$\sim$20\% in the solvent phase leads to a $\sim$40\% drop in the
conductance.

Similar behavior is observed in the TraPPE-UA model for toluene,
although this model has better agreement with the experimental
densities of toluene.  The expansion of the toluene liquid phase is
not as significant as that of the hexane (8.3\% over 100K), and this
limits the effect to $\sim$20\% drop in thermal conductivity.

Although we have not mapped out the behavior at a large number of
temperatures, is clear that there will be a strong temperature
dependence in the interfacial conductance when the physical properties
of one side of the interface (notably the density) change rapidly as a
function of temperature.

\section{Conclusions}
[VSIS WORKS! COMBINES NICE FEATURES OF PREVIOUS RNEMD METHODS AND
IMPROVEMENTS TO THEIR PROBLEMS! ROBUST AND VERSATILE!]

The NIVS algorithm has been applied to simulations of
butanethiol-capped Au(111) surfaces in the presence of organic
solvents. This algorithm allows the application of unphysical thermal
flux to transfer heat between the metal and the liquid phase. With the
flux applied, we were able to measure the corresponding thermal
gradients and to obtain interfacial thermal conductivities. Under
steady states, 2-3 ns trajectory simulations are sufficient for
computation of this quantity.

Our simulations have seen significant conductance enhancement in the
presence of capping agent, compared with the bare gold / liquid
interfaces. The vibrational coupling between the metal and the liquid
phase is enhanced by a chemically-bonded capping agent. Furthermore,
the coverage percentage of the capping agent plays an important role
in the interfacial thermal transport process. Moderately low coverages
allow higher contact between capping agent and solvent, and thus could
further enhance the heat transfer process, giving a non-monotonic
behavior of conductance with increasing coverage.

Our results, particularly using the UA models, agree well with
available experimental data.  The AA models tend to overestimate the
interfacial thermal conductance in that the classically treated C-H
vibrations become too easily populated. Compared to the AA models, the
UA models have higher computational efficiency with satisfactory
accuracy, and thus are preferable in modeling interfacial thermal
transport.

\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:	3778
Committed:	Sat Dec 10 02:42:55 2011 UTC (12 years, 9 months ago) by skuang
Content type:	application/x-tex
File size:	36115 byte(s)
Log Message:	finish much of SPC/E results.