--- stokes/stokes.tex	2011/12/11 03:22:47	3780
+++ stokes/stokes.tex	2011/12/13 23:31:35	3785
@@ -28,7 +28,8 @@
 
 \begin{document}
 
-\title{ENTER TITLE HERE}
+\title{A minimal perturbation approach to RNEMD able to simultaneously
+impose thermal and momentum gradients}
 
 \author{Shenyu Kuang and J. Daniel
 Gezelter\footnote{Corresponding author. \ Electronic mail: gezelter@nd.edu} \\
@@ -47,13 +48,15 @@ Notre Dame, Indiana 46556}
   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
+  Non-Equilibrium Molecular Dynamics (RNEMD) methods while maintaining
   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.
+  occured in previous NIVS simulations by using isotropic velocity
+  scaling on the molecules in specific regions of a system. To test
+  the method, we have computed the thermal conductivity and shear
+  viscosity of model liquid systems as well as the interfacial
+  frictions of a series of metal/liquid interfaces. Its ability to
+  combine the thermal and momentum gradients allows us to obtain shear
+  viscosity data for a range of temperatures in only one trajectory.
 
 \end{abstract}
 
@@ -66,32 +69,32 @@ Notre Dame, Indiana 46556}
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
 \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}
+simulations\cite{MullerPlathe:1997xw,ISI:000080382700030,kuang:164101}
+can establish steady state systems with an applied flux set vs a
+corresponding gradient that can be measured. These methods does not
+need many trajectories to provide information of transport properties
+of a given system. Thus, they are utilized in computing thermal and
+mechanical transfer of homogeneous bulk systems as well as
+heterogeneous systems such as solid-liquid
+interfaces.\cite{garde:nl2005,garde:PhysRevLett2009,kuang:AuThl}
 
 The Reverse Non-Equilibrium MD (RNEMD) methods adopt constraints that
 satisfy linear momentum and total energy conservation of a system when
 imposing fluxes in a simulation. Thus they are compatible with various
 ensembles, including the micro-canonical (NVE) ensemble, without the
-need of an external thermostat. The original approaches by
+need of an external thermostat. The original approaches proposed by
 M\"{u}ller-Plathe {\it et
   al.}\cite{MullerPlathe:1997xw,ISI:000080382700030} utilize simple
-momentum swapping for generating energy/momentum fluxes, which 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.
+momentum swapping for generating energy/momentum fluxes, which can
+also be compatible with particles of different identities. Although
+simple to implement in a simulation, this approach can create
+nonthermal velocity distributions, as discovered by Tenney and
+Maginn\cite{Maginn:2010}. Furthermore, this approach is less efficient
+for kinetic energy transfer between particles of different identities,
+especially when the mass difference between the particles becomes
+significant. This also limits its applications on heterogeneous
+interfacial systems.
 
 Recently, we developed a different approach, using Non-Isotropic
 Velocity Scaling (NIVS) \cite{kuang:164101} algorithm to impose
@@ -107,38 +110,40 @@ However, the NIVS approach limits its application in i
 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.
+However, the NIVS approach has limited applications in imposing
+momentum fluxes. Temperature anisotropy could happen under high
+momentum fluxes due to the implementation of this algorithm. Thus,
+combining thermal and momentum flux is also difficult to obtain with
+this approach. However, such combination may provide a means to
+simulate thermal/momentum gradient coupled processes such as Soret
+effect in liquid flows. Therefore, developing improved approaches to
+extend the applications of the imposed-flux method is desirable.
 
-In this paper, we improve the NIVS method and propose a novel approach
-to impose fluxes. This approach separate the means of applying
+In this paper, we improve the RNEMD methods by proposing a novel
+approach to impose fluxes. This approach separate the means of applying
 momentum and thermal flux with operations in one time step and thus is
 able to simutaneously impose thermal and momentum flux. Furthermore,
 the approach retains desirable features of previous RNEMD approaches
 and is simpler to implement compared to the NIVS method. In what
-follows, we first present the method to implement the method in a
+follows, we first present the method and its implementation in a
 simulation. Then we compare the method on bulk fluids to previous
 methods. Also, interfacial frictions are computed for a series of
 interfaces.
 
 \section{Methodology}
-Similar to the NIVS methodology,\cite{kuang:164101} we consider a
+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 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$''.
+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 separate operations to
-velocities of particles {$i$} within the center slab and of particles
-{$j$} within the end slab:
+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) \\
@@ -147,7 +152,7 @@ before an operation occurs. When a momentum flux $\vec
 \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})$
+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:
@@ -155,17 +160,18 @@ where
 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
+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 operations.
+and $\Delta t$ is the interval between two separate 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]
+The above operations already conserve the linear momentum of a
+periodic system. To further satisfy total energy conservation as well
+as to impose the thermal flux $J_z$, the following equations are
+included as well:
+[MAY PUT EXTRA MATH IN SUPPORT INFO OR APPENDIX]
 \begin{eqnarray}
 K_c - J_z\Delta t & = & c^2 (K_c - \frac{1}{2}M_c \langle\vec{v}_c
 \rangle^2) + \frac{1}{2}M_c (\langle \vec{v}_c \rangle + \vec{a}_c)^2 \\
@@ -173,27 +179,28 @@ $c$ and $h$ respectively before an operation occurs. T
 \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
+$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 do not
-change the potential energy of a system, given that the potential
+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 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.
+$\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 both thermal and momentum flux from
-  the ``c'' slab to the ``h'' slab. As the figure shows, the thermal
-  distributions preserve after this operation.}
+  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}
 
@@ -201,80 +208,92 @@ This approach is more computationaly efficient compare
 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.
+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 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.
+This approach is advantageous over the original momentum swapping in
+many aspects. In one swapping, the velocity vectors involved are
+usually very different (or the generated flux is trivial to obtain
+gradients), thus the swapping tends to incur perturbations to the
+neighbors of the particles involved. Comparatively, our approach
+disperse the flux to every selected particle in a slab so that
+perturbations in the flux generating region could be
+minimized. Additionally, because the momentum swapping steps tend to
+result in a nonthermal distribution, when an imposed flux is
+relatively large and diffusions from the neighboring slabs could no
+longer remedy this effect, problematic distributions would be
+observed. In comparison, the operations of our approach has the nature
+of preserving the equilibrium velocity distributions (commonly
+Maxwell-Boltzmann), and results in later section will illustrate that
+this is helpful to retain thermal distributions in a simulation.
 
 \section{Computational Details}
 The algorithm has been implemented in our MD simulation code,
 OpenMD\cite{Meineke:2005gd,openmd}. We compare the method with
-previous RNEMD methods or equilibrium MD methods in homogeneous fluids
+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$
+The systems to be investigated are set up in orthorhombic simulation
+cells with periodic boundary conditions in all three dimensions. The
+$z$ axis of these cells were longer and set as the temperature and/or
+momentum gradient axis. And the cells were evenly divided into $N$
 slabs along this axis, with various $N$ depending on individual
-system. The $x$ and $y$ axis were 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.
+system. The $x$ and $y$ axis were of the same length in homogeneous
+systems or had length scale close to each other where heterogeneous
+interfaces presents. In all cases, before introducing a nonequilibrium
+method to establish steady thermal and/or momentum gradients for
+further measurements and calculations, canonical ensemble with a
+Nos\'e-Hoover thermostat\cite{hoover85} and microcanonical ensemble
+equilibrations were used before data collections. For SPC/E water
+simulations, isobaric-isothermal equilibrations\cite{melchionna93} are
+performed before the above to reach normal pressure (1 bar); for
+interfacial systems, similar equilibrations are used to relax the
+surface tensions of the $xy$ plane.
 
-While homogeneous fluid systems can be set up with random
-configurations, our interfacial systems needs extra steps to ensure
-the interfaces be established properly for computations. The
+While homogeneous fluid systems can be set up with rather random
+configurations, our interfacial systems needs a series of steps to
+ensure the interfaces be established properly for computations. The
 preparation and equilibration of butanethiol covered gold (111)
 surface and further solvation and equilibration process is described
-as in reference \cite{kuang:AuThl}. 
+in details as in reference \cite{kuang:AuThl}.
 
 As for the ice/liquid water interfaces, the basal surface of ice
 lattice was first constructed. Hirsch {\it et
   al.}\cite{doi:10.1021/jp048434u} explored the energetics of ice
-lattices with 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]
+lattices with all possible proton order configurations. We refer to
+their results and choose the configuration of the lowest energy after
+geometry optimization as the unit cell for our ice lattices. Although
+experimental solid/liquid coexistant temperature under normal pressure
+should be close to 273K, Bryk and Haymet's simulations of ice/liquid
+water interfaces with different models suggest that for SPC/E, the
+most stable interface is observed at 225$\pm$5K.\cite{bryk:10258}
+Therefore, our ice/liquid water simulations were carried out at
+225K. To have extra protection of the ice lattice during initial
+equilibration (when the randomly generated liquid phase configuration
+could release large amount of energy in relaxation), restraints were
+applied to the ice lattice to avoid inadvertent melting by the heat
+dissipated from the high enery configurations.
+[MAY ADD A SNAPSHOT FOR BASAL PLANE]
 
 \subsection{Force Field Parameters}
 For comparison of our new method with previous work, we retain our
-force field parameters consistent with the results we will compare
-with. The Lennard-Jones fluid used here for argon , and reduced unit
-results are reported for direct comparison purpose.
+force field parameters consistent with previous simulations. Argon is
+the Lennard-Jones fluid used here, and its results are reported in
+reduced unit for direct comparison purpose.
 
 As for our water simulations, SPC/E model is used throughout this work
 for consistency. Previous work for transport properties of SPC/E water
@@ -289,9 +308,9 @@ The small organic molecules included in our simulation
 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
+For our gold/organic liquid interfaces, the small organic molecules
+included in our simulations are the Au surface capping agent
+butanethiol and liquid hexane and toluene. The United-Atom
 models\cite{TraPPE-UA.thiols,TraPPE-UA.alkanes,TraPPE-UA.alkylbenzenes}
 for these components were used in this work for better computational
 efficiency, while maintaining good accuracy. We refer readers to our
@@ -319,8 +338,8 @@ two in the denominator is present for the heat transpo
 \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
+two in the denominator is necessary for the heat transport occurs in
+both $+z$ and $-z$ directions. The average temperature gradient
 ${\langle\partial T/\partial z\rangle}$ can be obtained by a linear
 regression of the temperature profile, which is recorded during a
 simulation for each slab in a cell. For Lennard-Jones simulations,
@@ -329,18 +348,19 @@ by switching off $J_z$. One can specify the vector
 
 \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
+by specify a momentum flux. In our algorithm, one can specify the
+three components of the flux vector $\vec{j}_z(\vec{p})$
 respectively. For shear viscosity simulations, $j_z(p_z)$ is usually
-set to zero. 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.
+set to zero. For isotropic systems, the direction of
+$\vec{j}_z(\vec{p})$ on the $xy$ plane should not matter, but the
+ability of arbitarily specifying the vector direction in our method
+could provide convenience in anisotropic simulations.
 
-Similar to thermal conductivity computations, 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:
+Similar to thermal conductivity computations, for a homogeneous
+system, linear response of the momentum gradient with respect to the
+shear stress is assumed, and the shear viscosity ($\eta$) can be
+obtained with the imposed momentum flux (e.g. in $x$ direction) and
+the measured gradient:
 \begin{equation}
 j_z(p_x) = -\eta \frac{\partial v_x}{\partial z}
 \end{equation}
@@ -349,26 +369,40 @@ the data collection time. Also, the velocity gradient
 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
+the data collection time. Also, the averaged velocity gradient
 ${\langle\partial v_x/\partial z\rangle}$ can be obtained using linear
-regression of the $x$ component of the mean velocity, $\langle
-v_x\rangle$, in each of the bins. For Lennard-Jones simulations, shear
-viscosities are reported in reduced units
+regression of the $x$ component of the mean velocity ($\langle
+v_x\rangle$) in each of the bins. For Lennard-Jones simulations, shear
+viscosities are also reported in reduced units
 (${\eta^*=\eta\sigma^2(\varepsilon m)^{-1/2}}$).
+
+Although $J_z$ may be switched off for shear viscosity simulations at
+a certain temperature, our method's ability to impose both a thermal
+and a momentum flux in one simulation allows the combination of a
+temperature and a velocity gradient. In this case, since viscosity is
+generally a function of temperature, the local viscosity also depends
+on the local temperature. Therefore, in one such simulation, viscosity
+at $z$ (corresponding to a certain $T$) can be computed with the
+applied shear flux and the local velocity gradient (which can be
+obtained by finite difference approximation). As a whole, the
+viscosity can be mapped out as the function of temperature in one
+single trajectory of simulation. Results for shear viscosity
+computations of SPC/E water will demonstrate its effectiveness in
+detail.
 
 \subsection{Interfacial friction and Slip length}
 While the shear stress results in a velocity gradient within bulk
 fluid phase, its effect at a solid-liquid interface could vary due to
 the interaction strength between the two phases. The interfacial
 friction coefficient $\kappa$ is defined to relate the shear stress
-(e.g. along $x$-axis) and the relative fluid velocity tangent to the
+(e.g. along $x$-axis) with the relative fluid velocity tangent to the
 interface:
 \begin{equation}
 j_z(p_x)|_{interface} = \kappa\Delta v_x|_{interface}
 \end{equation}
 Under ``stick'' boundary condition, $\Delta v_x|_{interface}
 \rightarrow 0$, which leads to $\kappa\rightarrow\infty$. However, for
-``slip'' boundary condition at the solid-liquid interface, $\kappa$
+``slip'' boundary conditions at the solid-liquid interfaces, $\kappa$
 becomes finite. To characterize the interfacial boundary conditions,
 slip length ($\delta$) is defined using $\kappa$ and the shear
 viscocity of liquid phase ($\eta$):
@@ -378,14 +412,15 @@ solid-liquid interface. [MAY INCLUDE A SNAPSHOT IN FIG
 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]
+solid-liquid interface. [MAY INCLUDE SNAPSHOT IN FIGURE]
 
 \begin{figure}
 \includegraphics[width=\linewidth]{defDelta}
 \caption{The slip length $\delta$ can be obtained from a velocity
-  profile of a solid-liquid interface simulation. An example of
-  Au/hexane interfaces is shown. Calculation for the left side is 
-  illustrated. The right side is similar to the left side.}
+  profile of a solid-liquid interface simulation, when a momentum flux
+  is applied. An example of Au/hexane interfaces is shown, and the
+  calculation for the left side is illustrated. The calculation for
+  the right side is similar to the left.}
 \label{slipLength}
 \end{figure}
 
@@ -428,7 +463,7 @@ calculations with various fluxes in reduced units.
         \multicolumn{2}{c}{$\eta^*$} \\
         \hline
         Swap Interval $t^*$ & Equivalent $J_z^*$ or $j_z^*(p_x)$ &
-        NIVS & This work & Swapping & This work \\
+        NIVS\cite{kuang:164101} & This work & Swapping & This work \\
         \hline
         0.116 & 0.16  & 7.30(0.10) & 6.25(0.23) & 3.57(0.06) & 3.53(0.16)\\
         0.232 & 0.09  & 6.95(0.09) & 8.02(0.56) & 3.64(0.05) & 3.43(0.17)\\
@@ -445,7 +480,7 @@ the thermal gradients rendered using this method are a
 \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
+the thermal gradients introduced using this method are also close to
 previous RNEMD methods. Simulations with moderately higher thermal
 fluxes tend to yield more reliable thermal gradients and thus avoid
 large errors, while overly high thermal fluxes could introduce side
@@ -454,22 +489,22 @@ and $z$ axes, while thermal anisotropy might happen if
 
 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
+and $z$ axes, while for NIVS, thermal anisotropy might happen if the
+criteria function for choosing scaling coefficients does not perform
+as expected. Furthermore, this method avoids inadvertent concomitant
 momentum flux when only thermal flux is imposed, which could not be
 achieved with swapping or NIVS approaches. The thermal energy exchange
 in swapping ($\vec{p}_i$ in slab ``c'' with $\vec{p}_j$ in slab ``h'')
 or NIVS (total slab momentum conponemts $P^\alpha$ scaled to $\alpha
-P^\alpha$) would not 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.
+P^\alpha$) would not achieve this effect unless thermal flux vanishes
+(i.e. $\vec{p}_i=\vec{p}_j$ or $\alpha=1$ which do not contribute to
+applying a thermal flux). In this sense, this method aids to achieve
+minimal perturbation to a simulation while imposing a thermal flux.
 
 \subsubsection{Shear viscosity}
-Table \ref{LJ} also compares our shear viscosity results with momentum
-swapping approach. Our calculations show that our method predicted
-similar values for shear viscosities to the momentum swapping
+Table \ref{LJ} also compares our shear viscosity results with the
+momentum swapping approach. Our calculations show that our method
+predicted similar values of shear viscosities to the momentum swapping
 approach, as well as the velocity gradient profiles. Moderately larger
 momentum fluxes are helpful to reduce the errors of measured velocity
 gradients and thus the final result. However, it is pointed out that
@@ -479,8 +514,8 @@ temperatures were calculated after subtracting the eff
 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
+temperatures were calculated after subtracting the contribution from
+bulk velocities of the slabs. The one-dimensional temperature profiles
 showed no observable difference between the three dimensions. This
 ensures that isotropic scaling automatically preserves temperature
 isotropy and that our method is useful in shear viscosity
@@ -503,44 +538,47 @@ versa in slab ``h'', where the distributions leave a n
 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
+versa in slab ``h'', where the distributions leave notchs. This
 problematic profiles become significant when the imposed-flux becomes
 larger and diffusions from neighboring slabs could not offset the
-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
+depletions. Simutaneously, abnormal peaks appear corresponding to
+excessive particles having velocity swapped from the other slab. These
+nonthermal distributions limit applications of the swapping approach
+in shear stress simulations. Our method avoids the above problematic
 distributions by altering the means of applying momentum
 fluxes. Comparatively, velocity distributions recorded from
 simulations with our method is so close to the ideal thermal
-prediction that no observable difference is shown in Figure
-\ref{vDist}. Conclusively, our method avoids problems happened in
+prediction that no obvious difference is shown in Figure
+\ref{vDist}. Conclusively, our method avoids problems that occurs in
 previous RNEMD methods and provides a useful means for shear viscosity
 computations.
 
 \begin{figure}
 \includegraphics[width=\linewidth]{velDist}
 \caption{Velocity distributions that develop under the swapping and
-  our methods at high flux. These distributions were obtained from
+  our methods at a large flux. These distributions were obtained from
   Lennard-Jones simulations with $j_z^*(p_x)\sim 0.4$ (equivalent to a
   swapping interval of 20 time steps). This is a relatively large flux
   to demonstrate the nonthermal distributions that develop under the
-  swapping method. Distributions produced by our method are very close
-  to the ideal thermal situations.}
+  swapping method. In comparison, distributions produced by our method
+  are very close to the ideal thermal situations.}
 \label{vDist}
 \end{figure}
 
 \subsection{Bulk SPC/E water}
-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].
+We extend our applications of thermal conductivity and shear viscosity
+computations to a complex fluid model of SPC/E water. A simulation
+cell with 1000 molecules was set up in the similar manner as in
+\cite{kuang:164101}. For thermal conductivity simulations,
+measurements were taken to compare with previous RNEMD methods; for
+shear viscosity computations, simulations were run under a series of
+temperatures (with corresponding pressure relaxation using the
+isobaric-isothermal ensemble\cite{melchionna93}), and results were
+compared to available data from EMD
+methods\cite{10.1063/1.3330544,Medina2011}. Besides, a simulation with
+both thermal and momentum gradient was carried out to map out shear
+viscosity as a function of temperature to see the effectiveness and
+accuracy our method could reach.
 
 \subsubsection{Thermal conductivity}
 Table \ref{spceThermal} summarizes our thermal conductivity
@@ -549,7 +587,7 @@ energy even for systems involving electrostatic intera
 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.
+energy well for systems involving electrostatic interactions.
 
 Measurements using our method established similar temperature
 gradients to the previous NIVS method. Our simulation results are in
@@ -612,17 +650,19 @@ for shear viscosity computations.
         \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]\\
+        (K) & (10$^{10}$s$^{-1}$) & This work & Previous
+        simulations\cite{Medina2011} \\
         \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) &  \\
+        273 & 1.12 & 1.218(0.004) & 1.282(0.048) \\
+            & 1.79 & 1.140(0.012) &              \\
+        303 & 2.09 & 0.646(0.008) & 0.643(0.019) \\
+        318 & 2.50 & 0.536(0.007) &              \\
+            & 5.25 & 0.510(0.007) &              \\
+            & 2.82 & 0.474(0.003)\footnote{Simulation with $L_z$ twice
+                        as long.} &              \\
+        333 & 3.10 & 0.428(0.002) & 0.421(0.008) \\
+        363 & 2.34 & 0.279(0.014) & 0.291(0.005) \\
+            & 4.26 & 0.306(0.001) &              \\
         \hline\hline
       \end{tabular}
       \label{spceShear}
@@ -630,30 +670,70 @@ for shear viscosity computations.
   \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]
+A more effective way to map out $\eta$ vs $T$ is to combine a momentum
+flux with a thermal flux. Figure \ref{Tvxdvdz} shows the thermal and
+velocity gradient in one such simulation. At different positions with
+different temperatures, the velocity gradient is not a constant but
+can be computed locally. With the data provided in Figure
+\ref{Tvxdvdz}, a series of $\eta$ is calculated as in Figure
+\ref{etaT} and a linear fit was performed to $\partial v_x/\partial z$
+vs. $z$ so that the resulted $\eta$ can be present as a curve as
+well. For comparison, other results are also mapped in the figure.
+
+\begin{figure}
+\includegraphics[width=\linewidth]{tvxdvdz}
+\caption{With a combination of a thermal and a momentum flux, a
+  simulation can have both a temperature (top) and a velocity (middle)
+  gradient. Due to the thermal gradient, $\partial v_x/\partial z$ is
+  not constant but can be computed using finite difference
+  approximations (lower). These data can be used further to calculate
+  $\eta$ vs $T$ (Figure \ref{etaT}).}
+\label{Tvxdvdz}
+\end{figure}
+
+From Figure \ref{etaT}, one can see that the generated curve agrees
+well with the above RNEMD simulations at different temperatures, as
+well as results reported using EMD
+methods\cite{10.1063/1.3330544,Medina2011} in much of the temperature
+range simulated. However, this curve has relatively large error in
+lower temperature regions and has some difference in predicting $\eta$
+near 273K. Provided that this curve only takes one trajectory to
+generate, these results are of satisfactory efficiency and
+accuracy. Since previous work already pointed out that the SPC/E model
+tends to predict lower viscosity compared to experimental
+data,\cite{Medina2011} experimental comparison are not given here.
+
+\begin{figure}
+\includegraphics[width=\linewidth]{etaT}
+\caption{The curve generated by single simulation with thermal and
+  momentum gradient predicts satisfatory values in much of the
+  temperature range under test.}
+\label{etaT}
+\end{figure}
+
 \subsection{Interfacial frictions and slip lengths}
-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
+Another attractive aspect of our method is the ability to apply
+momentum and/or thermal flux in nonhomogeneous systems, where
+molecules of different identities (or phases) are segregated in
+different regions. We have previously studied the interfacial thermal
+transport of a series of metal gold-liquid
+surfaces\cite{kuang:164101,kuang:AuThl}, and would like to further
+investigate the relationship between this phenomenon and the
 interfacial frictions.
 
 Table \ref{etaKappaDelta} includes these computations and previous
 calculations of corresponding interfacial thermal conductance. For
 bare Au(111) surfaces, slip boundary conditions were observed for both
 organic and aqueous liquid phases, corresponding to previously
-computed low interfacial thermal conductance. 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]
+computed low interfacial thermal conductance. In comparison, the
+butanethiol covered Au(111) surface appeared to be sticky to the
+organic liquid layers in our simulations. We have reported conductance
+enhancement effect for this surface capping agent,\cite{kuang:AuThl}
+and these observations have a qualitative agreement with the thermal
+conductance results. This agreement also supports discussions on the
+relationship between surface wetting and slip effect and thermal
+conductance of the
+interface.\cite{PhysRevLett.82.4671,doi:10.1080/0026897031000068578,garde:PhysRevLett2009}
 
 \begin{table*}
   \begin{minipage}{\linewidth}
@@ -668,33 +748,33 @@ interface.[CITE BARRAT, GARDE]
         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 \\
+        surface & molecules & K & MPa & mPa$\cdot$s &
+        10$^4$Pa$\cdot$s/m & nm & MW/m$^2$/K \\
         \hline
-        Au(111) & hexane & 200 & 1.08 & 0.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 &
-            \\
+        Au(111) & hexane & 200 & 1.08 & 0.197(0.009) & 5.30(0.36) &
+        3.72 & 46.5 \\
+                &        &     & 2.15 & 0.141(0.002) & 5.31(0.26) &
+        2.76 &      \\
+        Au-SC$_4$H$_9$ & hexane & 200 & 2.16 & 0.286(0.019) & $\infty$
+        & 0 & 131 \\
+                       &        &     & 5.39 & 0.320(0.006) & $\infty$
+        & 0 &     \\
         \hline
-        Au(111) & toluene & 200 & 1.08 & 0.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
-        &     \\
+        Au(111) & toluene & 200 & 1.08 & 0.722(0.035) & 15.7(0.7) &
+        4.60 & 70.1 \\
+                &         &     & 2.16 & 0.544(0.030) & 11.2(0.5) &
+        4.86 &      \\
+        Au-SC$_4$H$_9$ & toluene & 200 & 5.39 & 0.980(0.057) &
+        $\infty$ & 0 & 187 \\
+                       &         &     & 10.8 & 0.995(0.005) &
+        $\infty$ & 0 &     \\
         \hline
-        Au(111) & water & 300 & 1.08 & 0.40() & 1.9$\times$10$^4$() &
+        Au(111) & water & 300 & 1.08 & 0.399(0.050) & 1.928(0.022) &
         20.7 & 1.65 \\
-                &       &     & 2.16 & 0.79() & 1.9$\times$10$^4$() &
+                &       &     & 2.16 & 0.794(0.255) & 1.895(0.003) &
         41.9 &      \\
         \hline
-        ice(basal) & water & 225 & 19.4 & 15.8() & $\infty$ & 0 & \\
+        ice(basal) & water & 225 & 19.4 & 15.8(0.2) & $\infty$ & 0 & \\
         \hline\hline
       \end{tabular}
       \label{etaKappaDelta}
@@ -704,24 +784,24 @@ solid surface. Note that $\eta$ measured near a ``slip
 
 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.
+solid surface. In our results, $\eta$ measured near a ``slip'' surface
+tends to be smaller than that near a ``stick'' surface. This may
+suggest the influence from an interface on the dynamic properties of
+liquid within its neighbor regions. It is known that diffusions of
+solid particles in liquid phase is affected by their surface
+conditions (stick or slip boundary).\cite{10.1063/1.1610442} Our
+observations could provide a support to this phenomenon.
 
 In addition to these previously studied interfaces, we attempt to
 construct ice-water interfaces and the basal plane of ice lattice was
-first studied. In contrast to the Au(111)/water interface, where the
-friction coefficient is relatively small and large slip effect
+studied here. In contrast to the Au(111)/water interface, where the
+friction coefficient is substantially small and large slip effect
 presents, the ice/liquid water interface demonstrates strong
-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.
+solid-liquid interactions and appears to be sticky. The supercooled
+liquid phase is an order of magnitude more viscous than measurements
+in previous section. It would be of interst to investigate the effect
+of different ice lattice planes (such as prism and other surfaces) on
+interfacial friction and the corresponding liquid viscosity.
 
 \section{Conclusions}
 Our simulations demonstrate the validity of our method in RNEMD
@@ -733,18 +813,20 @@ Furthermore, using this method, investigations can 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.
+Our method is capable of effectively imposing thermal and/or momentum
+flux accross an interface. This facilitates studies that relates
+dynamic property measurements to the chemical details of an
+interface. Therefore, investigations can be carried out to
+characterize interfacial interactions using the method.
 
 Another attractive feature of our method is the ability of
-simultaneously 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.
+simultaneously introducing thermal and momentum gradients in a
+system. This facilitates us to effectively map out the shear viscosity
+with respect to a range of temperature in single trajectory of
+simulation with satisafactory accuracy. Complex systems that involve
+thermal and momentum gradients might potentially benefit from
+this. For example, the Soret effects under a velocity gradient might
+be models of interest to purification and separation researches.
 
 \section{Acknowledgments}
 Support for this project was provided by the National Science
