--- stokes/stokes.tex	2011/12/09 05:08:56	3776
+++ stokes/stokes.tex	2011/12/11 03:22:47	3780
@@ -43,25 +43,17 @@ Notre Dame, Indiana 46556}
 \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.
+  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}
 
@@ -226,7 +218,6 @@ section will illustrate this effect.
 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.
-%NEW METHOD (AND NIVS) HAVE LESS PERTURBATION THAN MOMENTUM SWAPPING
 
 \section{Computational Details}
 The algorithm has been implemented in our MD simulation code,
@@ -392,8 +383,9 @@ solid-liquid interface. [MAY INCLUDE A SNAPSHOT IN FIG
 \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.}
+  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}
 
@@ -405,8 +397,6 @@ data.
 interface can be measured respectively. Further calculations and
 characterizations of the interface can be carried out using these
 data.
-[MENTION IN RESULTS THAT ETA OBTAINED HERE DOES NOT NECESSARILY EQUAL
-TO BULK VALUES]
 
 \section{Results and Discussions}
 \subsection{Lennard-Jones fluid}
@@ -437,7 +427,7 @@ calculations with various fluxes in reduced units.
         \multicolumn{2}{c}{$\lambda^*$} &
         \multicolumn{2}{c}{$\eta^*$} \\
         \hline
-        Swap Interval & Equivalent $J_z^*$ or $j_z^*(p_x)$ &
+        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)\\
@@ -469,7 +459,7 @@ in swapping ($\vec{p}_i$ in slab ``c'' with $\vec{p}_j
 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 ``j'')
+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
@@ -498,7 +488,12 @@ computations.
 
 \begin{figure}
 \includegraphics[width=\linewidth]{tempXyz}
-\caption{.}
+\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}
 
@@ -525,136 +520,232 @@ computations.
 
 \begin{figure}
 \includegraphics[width=\linewidth]{velDist}
-\caption{.}
+\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}
-[WATER COMPARED TO RNEMD NIVS AND EMD]
+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}
-[VSIS DOES AS WELL AS NIVS]
+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.
 
-\subsubsection{Shear viscosity}
-[COMPARE W EMD]
+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.
 
-[MAY HAVE A FIRURE FOR DATA]
-
-[MAY COMBINE JZPX AND JZKE TO RUN ONE JOB BUT GET ETA=F(T)]
-[PUT RESULT AND FIGURE HERE IF IT WORKS]
-\subsection{Interfacial frictions}
-[SLIP BOUNDARY VS STICK BOUNDARY; ICE-WATER INTERFACES]
+\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*}
 
-qualitative agreement w interfacial thermal conductance
+\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.
 
-[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.}
+      \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}{llccc}
+      \begin{tabular}{cccc}
         \hline\hline
-        Butanethiol model & Solvent & $G$ & $G^\prime$ \\
-        (or bare surface) & model & \multicolumn{2}{c}{(MW/m$^2$/K)} \\
+        $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
-        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) \\
+        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{modelTest}
+      \label{spceShear}
     \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.
+[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.
 
-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.
+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]
 
-\subsubsection{Effects due to average temperature}
+\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*}
 
-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.
+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.
 
-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.
+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.
 
-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!]
+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.
 
-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.
+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.
 
-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.
+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.
 
-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
