--- trunk/nivsRnemd/nivsRnemd.tex	2010/03/10 02:00:06	3570
+++ trunk/nivsRnemd/nivsRnemd.tex	2010/07/14 15:20:31	3611
@@ -6,6 +6,7 @@
 \usepackage{caption}
 %\usepackage{tabularx}
 \usepackage{graphicx}
+\usepackage{multirow}
 %\usepackage{booktabs}
 %\usepackage{bibentry}
 %\usepackage{mathrsfs}
@@ -38,7 +39,22 @@ Notre Dame, Indiana 46556}
 \begin{doublespace}
 
 \begin{abstract}
-
+  We present a new method for introducing stable non-equilibrium
+  velocity and temperature distributions in molecular dynamics
+  simulations of heterogeneous systems.  This method extends earlier
+  Reverse Non-Equilibrium Molecular Dynamics (RNEMD) methods which use
+  momentum exchange swapping moves that can create non-thermal
+  velocity distributions and are difficult to use for interfacial
+  calculations.  By using non-isotropic velocity scaling (NIVS) on the
+  molecules in specific regions of a system, it is possible to impose
+  momentum or thermal flux between regions of a simulation and stable
+  thermal and momentum gradients can then be established.  The scaling
+  method we have developed conserves the total linear momentum and
+  total energy of the system.  To test the methods, we have computed
+  the thermal conductivity of model liquid and solid systems as well
+  as the interfacial thermal conductivity of a metal-water interface.
+  We find that the NIVS-RNEMD improves the problematic velocity
+  distributions that develop in other RNEMD methods.
 \end{abstract}
 
 \newpage
@@ -49,39 +65,51 @@ Notre Dame, Indiana 46556}
 %                          BODY OF TEXT
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
-
-
 \section{Introduction}
 The original formulation of Reverse Non-equilibrium Molecular Dynamics
 (RNEMD) obtains transport coefficients (thermal conductivity and shear
 viscosity) in a fluid by imposing an artificial momentum flux between
 two thin parallel slabs of material that are spatially separated in
 the simulation cell.\cite{MullerPlathe:1997xw,ISI:000080382700030} The
-artificial flux is typically created by periodically ``swapping'' either
-the entire momentum vector $\vec{p}$ or single components of this
-vector ($p_x$) between molecules in each of the two slabs.  If the two
-slabs are separated along the z coordinate, the imposed flux is either
-directional ($j_z(p_x)$) or isotropic ($J_z$), and the response of a
-simulated system to the imposed momentum flux will typically be a
-velocity or thermal gradient.  The transport coefficients (shear
-viscosity and thermal conductivity) are easily obtained by assuming
-linear response of the system,
+artificial flux is typically created by periodically ``swapping''
+either the entire momentum vector $\vec{p}$ or single components of
+this vector ($p_x$) between molecules in each of the two slabs.  If
+the two slabs are separated along the $z$ coordinate, the imposed flux
+is either directional ($j_z(p_x)$) or isotropic ($J_z$), and the
+response of a simulated system to the imposed momentum flux will
+typically be a velocity or thermal gradient (Fig. \ref{thermalDemo}).
+The shear viscosity ($\eta$) and thermal conductivity ($\lambda$) are
+easily obtained by assuming linear response of the system,
 \begin{eqnarray}
 j_z(p_x) & = & -\eta \frac{\partial v_x}{\partial z}\\
-J & = & \lambda \frac{\partial T}{\partial z}
+J_z & = & \lambda \frac{\partial T}{\partial z}
 \end{eqnarray}
-RNEMD has been widely used to provide computational estimates of thermal
-conductivities and shear viscosities in a wide range of materials,
-from liquid copper to monatomic liquids to molecular fluids
-(e.g. ionic liquids).\cite{ISI:000246190100032}
+RNEMD has been widely used to provide computational estimates of
+thermal conductivities and shear viscosities in a wide range of
+materials, from liquid copper to both monatomic and molecular fluids
+(e.g.  ionic
+liquids).\cite{Bedrov:2000-1,Bedrov:2000,Muller-Plathe:2002,ISI:000184808400018,ISI:000231042800044,Maginn:2007,Muller-Plathe:2008,ISI:000258460400020,ISI:000258840700015,ISI:000261835100054}
 
-RNEMD is preferable in many ways to the forward NEMD methods because
-it imposes what is typically difficult to measure (a flux or stress)
-and it is typically much easier to compute momentum gradients or
-strains (the response).  For similar reasons, RNEMD is also preferable
-to slowly-converging equilibrium methods for measuring thermal
-conductivity and shear viscosity (using Green-Kubo relations or the
-Helfand moment approach of Viscardy {\it et
+\begin{figure}
+\includegraphics[width=\linewidth]{thermalDemo}
+\caption{RNEMD methods impose an unphysical transfer of momentum or
+  kinetic energy between a ``hot'' slab and a ``cold'' slab in the
+  simulation box.  The molecular system responds to this imposed flux
+  by generating a momentum or temperature gradient.  The slope of the
+  gradient can then be used to compute transport properties (e.g.
+  shear viscosity and thermal conductivity).}
+\label{thermalDemo}
+\end{figure}
+
+RNEMD is preferable in many ways to the forward NEMD
+methods\cite{ISI:A1988Q205300014,hess:209,Vasquez:2004fk,backer:154503,ISI:000266247600008}
+because it imposes what is typically difficult to measure (a flux or
+stress) and it is typically much easier to compute the response
+(momentum gradients or strains).  For similar reasons, RNEMD is also
+preferable to slowly-converging equilibrium methods for measuring
+thermal conductivity and shear viscosity (using Green-Kubo
+relations\cite{daivis:541,mondello:9327} or the Helfand moment
+approach of Viscardy {\it et
   al.}\cite{Viscardy:2007bh,Viscardy:2007lq}) because these rely on
 computing difficult to measure quantities.
 
@@ -91,29 +119,29 @@ Recently, Tenney and Maginn\cite{ISI:000273472300004} 
 typically samples from the same manifold of states in the
 microcanonical ensemble.
 
-Recently, Tenney and Maginn\cite{ISI:000273472300004} have discovered
+Recently, Tenney and Maginn\cite{Maginn:2010} have discovered
 some problems with the original RNEMD swap technique.  Notably, large
 momentum fluxes (equivalent to frequent momentum swaps between the
-slabs) can result in ``notched'', ``peaked'' and generally non-thermal momentum
-distributions in the two slabs, as well as non-linear thermal and
-velocity distributions along the direction of the imposed flux ($z$).
-Tenney and Maginn obtained reasonable limits on imposed flux and
-self-adjusting metrics for retaining the usability of the method.
+slabs) can result in ``notched'', ``peaked'' and generally non-thermal
+momentum distributions in the two slabs, as well as non-linear thermal
+and velocity distributions along the direction of the imposed flux
+($z$). Tenney and Maginn obtained reasonable limits on imposed flux
+and self-adjusting metrics for retaining the usability of the method.
 
 In this paper, we develop and test a method for non-isotropic velocity
-scaling (NIVS-RNEMD) which retains the desirable features of RNEMD
+scaling (NIVS) which retains the desirable features of RNEMD
 (conservation of linear momentum and total energy, compatibility with
 periodic boundary conditions) while establishing true thermal
-distributions in each of the two slabs.  In the next section, we
-develop the method for determining the scaling constraints.  We then
-test the method on both single component, multi-component, and
-non-isotropic mixtures and show that it is capable of providing
+distributions in each of the two slabs. In the next section, we
+present the method for determining the scaling constraints.  We then
+test the method on both liquids and solids as well as a non-isotropic
+liquid-solid interface and show that it is capable of providing
 reasonable estimates of the thermal conductivity and shear viscosity
-in these cases.
+in all of these cases.
 
 \section{Methodology}
-We retain the basic idea of Muller-Plathe's RNEMD method; the periodic
-system is partitioned into a series of thin slabs along a particular
+We retain the basic idea of M\"{u}ller-Plathe's RNEMD method; the
+periodic system is partitioned into a series of thin slabs along one
 axis ($z$).  One of the slabs at the end of the periodic box is
 designated the ``hot'' slab, while the slab in the center of the box
 is designated the ``cold'' slab.  The artificial momentum flux will be
@@ -121,7 +149,7 @@ moves.  For molecules $\{i\}$ located within the cold 
 hot slab.
 
 Rather than using momentum swaps, we use a series of velocity scaling
-moves.  For molecules $\{i\}$ located within the cold slab,
+moves.  For molecules $\{i\}$  located within the cold slab,
 \begin{equation}
 \vec{v}_i \leftarrow \left( \begin{array}{ccc}
 x & 0 & 0 \\
@@ -129,9 +157,10 @@ where ${x, y, z}$ are a set of 3 scaling variables for
 0 & 0 & z \\
 \end{array} \right) \cdot \vec{v}_i
 \end{equation}
-where ${x, y, z}$ are a set of 3 scaling variables for each of the
-three directions in the system.  Likewise, the molecules $\{j\}$
-located in the hot slab will see a concomitant scaling of velocities,
+where ${x, y, z}$ are a set of 3 velocity-scaling variables for each
+of the three directions in the system.  Likewise, the molecules
+$\{j\}$ located in the hot slab will see a concomitant scaling of
+velocities,
 \begin{equation}
 \vec{v}_j \leftarrow \left( \begin{array}{ccc}
 x^\prime & 0 & 0 \\
@@ -141,7 +170,7 @@ Conservation of linear momentum in each of the three d
 \end{equation}
 
 Conservation of linear momentum in each of the three directions
-($\alpha = x,y,z$) ties the values of the hot and cold bin scaling
+($\alpha = x,y,z$) ties the values of the hot and cold scaling
 parameters together:
 \begin{equation}
 P_h^\alpha + P_c^\alpha = \alpha^\prime P_h^\alpha + \alpha P_c^\alpha
@@ -170,13 +199,15 @@ where the kinetic energies, $K_h^\alpha$ and $K_c^\alp
 \sum_{\alpha = x,y,z} K_h^\alpha + K_c^\alpha = \sum_{\alpha = x,y,z}
 \left(\alpha^\prime\right)^2 K_h^\alpha + \alpha^2 K_c^\alpha
 \end{equation}
-where the kinetic energies, $K_h^\alpha$ and $K_c^\alpha$, are computed
-for each of the three directions in a similar manner to the linear momenta
-(Eq. \ref{eq:momentumdef}).  Substituting in the expressions for the
-hot scaling parameters ($\alpha^\prime$) from Eq. (\ref{eq:hotcoldscaling}),
-we obtain the {\it constraint ellipsoid equation}:
+where the translational kinetic energies, $K_h^\alpha$ and
+$K_c^\alpha$, are computed for each of the three directions in a
+similar manner to the linear momenta (Eq. \ref{eq:momentumdef}).
+Substituting in the expressions for the hot scaling parameters
+($\alpha^\prime$) from Eq. (\ref{eq:hotcoldscaling}), we obtain the
+{\it constraint ellipsoid}:
 \begin{equation}
-\sum_{\alpha = x,y,z} a_\alpha \alpha^2 + b_\alpha \alpha + c_\alpha = 0
+\sum_{\alpha = x,y,z} \left( a_\alpha \alpha^2 + b_\alpha \alpha +
+  c_\alpha \right) = 0
 \label{eq:constraintEllipsoid}
 \end{equation}
 where the constants are obtained from the instantaneous values of the
@@ -188,230 +219,359 @@ This ellipsoid equation defines the set of cold slab s
 c_\alpha & = & K_h^\alpha p_\alpha^2 + 2 K_h^\alpha p_\alpha - K_c^\alpha
 \label{eq:constraintEllipsoidConsts}
 \end{eqnarray}
-This ellipsoid equation defines the set of cold slab scaling
-parameters which can be applied while preserving both linear momentum
-in all three directions as well as kinetic energy.
+This ellipsoid (Eq. \ref{eq:constraintEllipsoid}) defines the set of
+cold slab scaling parameters which, when applied, preserve the linear
+momentum of the system in all three directions as well as total
+kinetic energy.
 
-The goal of using velocity scaling variables is to transfer linear
-momentum or kinetic energy from the cold slab to the hot slab.  If the
-hot and cold slabs are separated along the z-axis, the energy flux is
-given simply by the decrease in kinetic energy of the cold bin:
+The goal of using these velocity scaling variables is to transfer
+kinetic energy from the cold slab to the hot slab.  If the hot and
+cold slabs are separated along the z-axis, the energy flux is given
+simply by the decrease in kinetic energy of the cold bin:
 \begin{equation}
 (1-x^2) K_c^x  + (1-y^2) K_c^y + (1-z^2) K_c^z = J_z\Delta t
 \end{equation}
 The expression for the energy flux can be re-written as another
 ellipsoid centered on $(x,y,z) = 0$:
 \begin{equation}
-x^2 K_c^x + y^2 K_c^y + z^2 K_c^z = K_c^x + K_c^y + K_c^z - J_z\Delta t
+\sum_{\alpha = x,y,z} K_c^\alpha \alpha^2 = \sum_{\alpha = x,y,z}
+K_c^\alpha -J_z \Delta t
 \label{eq:fluxEllipsoid}
 \end{equation}
-The spatial extent of the {\it flux ellipsoid equation} is governed
-both by a targetted value, $J_z$ as well as the instantaneous values of the
-kinetic energy components in the cold bin.
+The spatial extent of the {\it thermal flux ellipsoid} is governed
+both by the target flux, $J_z$ as well as the instantaneous values of
+the kinetic energy components in the cold bin.
 
 To satisfy an energetic flux as well as the conservation constraints,
-it is sufficient to determine the points ${x,y,z}$ which lie on both
-the constraint ellipsoid (Eq. \ref{eq:constraintEllipsoid}) and the
-flux ellipsoid (Eq. \ref{eq:fluxEllipsoid}), i.e. the intersection of
-the two ellipsoids in 3-dimensional space.
+we must determine the points ${x,y,z}$ that lie on both the constraint
+ellipsoid (Eq. \ref{eq:constraintEllipsoid}) and the flux ellipsoid
+(Eq. \ref{eq:fluxEllipsoid}), i.e. the intersection of the two
+ellipsoids in 3-dimensional space.
 
 \begin{figure}
 \includegraphics[width=\linewidth]{ellipsoids}
-\caption{Scaling points which maintain both constant energy and
-  constant linear momentum of the system lie on the surface of the
-  {\it constraint ellipsoid} while points which generate the target
-  momentum flux lie on the surface of the {\it flux ellipsoid}.  The
-  velocity distributions in the hot bin are scaled by only those
+\caption{Velocity scaling coefficients which maintain both constant
+  energy and constant linear momentum of the system lie on the surface
+  of the {\it constraint ellipsoid} while points which generate the
+  target momentum flux lie on the surface of the {\it flux ellipsoid}.
+  The velocity distributions in the cold bin are scaled by only those
   points which lie on both ellipsoids.}
 \label{ellipsoids}
 \end{figure}
 
-One may also define momentum flux (say along the x-direction) as:
+Since ellipsoids can be expressed as polynomials up to second order in
+each of the three coordinates, finding the the intersection points of
+two ellipsoids is isomorphic to finding the roots a polynomial of
+degree 16.  There are a number of polynomial root-finding methods in
+the literature, [CITATIONS NEEDED] but numerically finding the roots
+of high-degree polynomials is generally an ill-conditioned
+problem.[CITATION NEEDED] One simplification is to maintain velocity
+scalings that are {\it as isotropic as possible}.  To do this, we
+impose $x=y$, and to treat both the constraint and flux ellipsoids as
+2-dimensional ellipses.  In reduced dimensionality, the
+intersecting-ellipse problem reduces to finding the roots of
+polynomials of degree 4.
+
+Depending on the target flux and current velocity distributions, the
+ellipsoids can have between 0 and 4 intersection points.  If there are
+no intersection points, it is not possible to satisfy the constraints
+while performing a non-equilibrium scaling move, and no change is made
+to the dynamics.  
+
+With multiple intersection points, any of the scaling points will
+conserve the linear momentum and kinetic energy of the system and will
+generate the correct target flux.  Although this method is inherently
+non-isotropic, the goal is still to maintain the system as close to an
+isotropic fluid as possible.  With this in mind, we would like the
+kinetic energies in the three different directions could become as
+close as each other as possible after each scaling.  Simultaneously,
+one would also like each scaling as gentle as possible, i.e. ${x,y,z
+  \rightarrow 1}$, in order to avoid large perturbation to the system.
+To do this, we pick the intersection point which maintains the three
+scaling variables ${x, y, z}$ as well as the ratio of kinetic energies
+${K_c^z/K_c^x, K_c^z/K_c^y}$ as close as possible to 1.
+
+After the valid scaling parameters are arrived at by solving geometric
+intersection problems in $x, y, z$ space in order to obtain cold slab
+scaling parameters, Eq. (\ref{eq:hotcoldscaling}) is used to
+determine the conjugate hot slab scaling variables.
+
+\subsection{Introducing shear stress via velocity scaling}
+It is also possible to use this method to magnify the random
+fluctuations of the average momentum in each of the bins to induce a
+momentum flux.  Doing this repeatedly will create a shear stress on
+the system which will respond with an easily-measured strain.  The
+momentum flux (say along the $x$-direction) may be defined as:
 \begin{equation}
 (1-x) P_c^x = j_z(p_x)\Delta t
 \label{eq:fluxPlane}
 \end{equation}
-The above {\it flux equation} is essentially a plane which is
-perpendicular to the x-axis, with its position governed both by a
-targetted value, $j_z(p_x)$ as well as the instantaneous value of the
-momentum along the x-direction.
+This {\it momentum flux plane} is perpendicular to the $x$-axis, with
+its position governed both by a target value, $j_z(p_x)$ as well as
+the instantaneous value of the momentum along the $x$-direction.
 
-Similarly, to satisfy a momentum flux as well as the conservation
-constraints, it is sufficient to determine the points ${x,y,z}$ which
-lie on both the constraint ellipsoid (Eq. \ref{eq:constraintEllipsoid})
-and the flux plane (Eq. \ref{eq:fluxPlane}), i.e. the intersection of
-an ellipsoid and a plane in 3-dimensional space.
+In order to satisfy a momentum flux as well as the conservation
+constraints, we must determine the points ${x,y,z}$ which lie on both
+the constraint ellipsoid (Eq. \ref{eq:constraintEllipsoid}) and the
+flux plane (Eq. \ref{eq:fluxPlane}), i.e. the intersection of an
+ellipsoid and a plane in 3-dimensional space.
 
-To summarize, by solving respective equation sets, one can determine
-possible sets of scaling variables for cold slab. And corresponding
-sets of scaling variables for hot slab can be determine as well.
+In the case of momentum flux transfer, we also impose another
+constraint to set the kinetic energy transfer as zero. In other
+words, we apply Eq. \ref{eq:fluxEllipsoid} and let ${J_z = 0}$.  With
+one variable fixed by Eq. \ref{eq:fluxPlane}, one now have a similar
+set of quartic equations to the above kinetic energy transfer problem.
 
-The following problem will be choosing an optimal set of scaling
-variables among the possible sets. Although this method is inherently
-non-isotropic, the goal is still to maintain the system as isotropic
-as possible. Under this consideration, one would like the kinetic
-energies in different directions could become as close as each other
-after each scaling. Simultaneously, one would also like each scaling
-as gentle as possible, i.e. ${x,y,z \rightarrow 1}$, in order to avoid
-large perturbation to the system. Therefore, one approach to obtain the
-scaling variables would be constructing an criteria function, with
-constraints as above equation sets, and solving the function's minimum
-by method like Lagrange multipliers.
-
-In order to save computation time, we have a different approach to a
-relatively good set of scaling variables with much less calculation
-than above. Here is the detail of our simplification of the problem.
+\section{Computational Details}
 
-In the case of kinetic energy transfer, we impose another constraint
-${x = y}$, into the equation sets. Consequently, there are two
-variables left. And now one only needs to solve a set of two {\it
-  ellipses equations}. This problem would be transformed into solving
-one quartic equation for one of the two variables. There are known
-generic methods that solve real roots of quartic equations. Then one
-can determine the other variable and obtain sets of scaling
-variables. Among these sets, one can apply the above criteria to
-choose the best set, while much faster with only a few sets to choose.
+We have implemented this methodology in our molecular dynamics code,
+OpenMD,\cite{Meineke:2005gd,openmd} performing the NIVS scaling moves
+after each MD step.  We have tested the method in a variety of
+different systems, including homogeneous fluids (Lennard-Jones and
+SPC/E water), crystalline solids ({\sc eam}~\cite{PhysRevB.33.7983} and
+quantum Sutton-Chen ({\sc q-sc})~\cite{PhysRevB.59.3527}
+models for Gold), and heterogeneous interfaces (QSC gold - SPC/E
+water). The last of these systems would have been difficult to study
+using previous RNEMD methods, but using velocity scaling moves, we can
+even obtain estimates of the interfacial thermal conductivities ($G$).
 
-In the case of momentum flux transfer, we impose another constraint to
-set the kinetic energy transfer as zero. In another word, we apply
-Eq. \ref{eq:fluxEllipsoid} and let ${J_z = 0}$. After that, with one
-variable fixed by Eq. \ref{eq:fluxPlane}, one now have a similar set
-of equations on the above kinetic energy transfer problem. Therefore,
-an approach similar to the above would be sufficient for this as well.
+\subsection{Simulation Cells}
 
-\section{Computational Details}
-Our simulation consists of a series of systems. All of these
-simulations were run with the OpenMD simulation software
-package\cite{Meineke:2005gd} integrated with RNEMD methods.
+In each of the systems studied, the dynamics was carried out in a
+rectangular simulation cell using periodic boundary conditions in all
+three dimensions.  The cells were longer along the $z$ axis and the
+space was divided into $N$ slabs along this axis (typically $N=20$).
+The top slab ($n=1$) was designated the ``cold'' slab, while the
+central slab ($n= N/2 + 1$) was designated the ``hot'' slab.  In all
+cases, simulations were first thermalized in canonical ensemble (NVT)
+using a Nos\'{e}-Hoover thermostat.\cite{Hoover85} then equilibrated in
+microcanonical ensemble (NVE) before introducing any non-equilibrium
+method.
 
-A Lennard-Jones fluid system was built and tested first. In order to
-compare our method with swapping RNEMD, a series of simulations were
-performed to calculate the shear viscosity and thermal conductivity of
-argon. 2592 atoms were in a orthorhombic cell, which was ${10.06\sigma
-  \times 10.06\sigma \times 30.18\sigma}$ by size. The reduced density
-${\rho^* = \rho\sigma^3}$ was thus 0.85, which enabled direct
-comparison between our results and others. These simulations used
-velocity Verlet algorithm with reduced timestep ${\tau^* =
-  4.6\times10^{-4}}$. 
+\subsection{RNEMD with M\"{u}ller-Plathe swaps}
 
-For shear viscosity calculation, the reduced temperature was ${T^* =
-  k_B T/\varepsilon = 0.72}$. Simulations were first equilibrated in canonical
-ensemble (NVT), then equilibrated in microcanonical ensemble
-(NVE). Establishing and stablizing momentum gradient were followed
-also in NVE ensemble. For the swapping part, Muller-Plathe's algorithm was
-adopted.\cite{ISI:000080382700030} The simulation box was under
-periodic boundary condition, and devided into ${N = 20}$ slabs. In each swap,
-the top slab ${(n = 1)}$ exchange the most negative $x$ momentum with the
-most positive $x$ momentum in the center slab ${(n = N/2 + 1)}$. Referred
-to Tenney {\it et al.}\cite{ISI:000273472300004}, a series of swapping
-frequency were chosen. According to each result from swapping
-RNEMD, scaling RNEMD simulations were run with the target momentum
-flux set to produce a similar momentum flux and shear
-rate. Furthermore, various scaling frequencies can be tested for one
-single swapping rate. To compare the performance between swapping and
-scaling method, temperatures of different dimensions in all the slabs
-were observed. Most of the simulations include $10^5$ steps of
-equilibration without imposing momentum flux, $10^5$ steps of
-stablization with imposing momentum transfer, and $10^6$ steps of data
-collection under RNEMD. For relatively high momentum flux simulations,
-${5\times10^5}$ step data collection is sufficient. For some low momentum
-flux simulations, ${2\times10^6}$ steps were necessary.
+In order to compare our new methodology with the original
+M\"{u}ller-Plathe swapping algorithm,\cite{ISI:000080382700030} we
+first performed simulations using the original technique.
 
-After each simulation, the shear viscosity was calculated in reduced
-unit. The momentum flux was calculated with total unphysical
-transferred momentum ${P_x}$ and data collection time $t$:
+\subsection{RNEMD with NIVS scaling}
+
+For each simulation utilizing the swapping method, a corresponding
+NIVS-RNEMD simulation was carried out using a target momentum flux set
+to produce a the same momentum or energy flux exhibited in the
+swapping simulation.
+
+To test the temperature homogeneity (and to compute transport
+coefficients), directional momentum and temperature distributions were
+accumulated for molecules in each of the slabs.
+
+\subsection{Shear viscosities}
+
+The momentum flux was calculated using the total non-physical momentum
+transferred (${P_x}$) and the data collection time ($t$):
 \begin{equation}
 j_z(p_x) = \frac{P_x}{2 t L_x L_y}
 \end{equation}
-And the velocity gradient ${\langle \partial v_x /\partial z \rangle}$
-can be obtained by a linear regression of the velocity profile. From
-the shear viscosity $\eta$ calculated with the above parameters, one
-can further convert it into reduced unit ${\eta^* = \eta \sigma^2
-  (\varepsilon  m)^{-1/2}}$.
+where $L_x$ and $L_y$ denote the $x$ and $y$ lengths of the simulation
+box.  The factor of two in the denominator is present because physical
+momentum transfer occurs in two directions due to our periodic
+boundary conditions.  The velocity gradient ${\langle \partial v_x
+  /\partial z \rangle}$ was obtained using linear regression of the
+velocity profiles in the bins.  For Lennard-Jones simulations, shear
+viscosities are reporte in reduced units (${\eta^* = \eta \sigma^2
+  (\varepsilon m)^{-1/2}}$).
 
-For thermal conductivity calculation, simulations were first run under
-reduced temperature ${T^* = 0.72}$ in NVE ensemble. Muller-Plathe's
-algorithm was adopted in the swapping method. Under identical
-simulation box parameters, in each swap, the top slab exchange the
-molecule with least kinetic energy with the molecule in the center
-slab with most kinetic energy, unless this ``coldest'' molecule in the
-``hot'' slab is still ``hotter'' than the ``hottest'' molecule in the ``cold''
-slab. According to swapping RNEMD results, target energy flux for
-scaling RNEMD simulations can be set. Also, various scaling
-frequencies can be tested for one target energy flux. To compare the
-performance between swapping and scaling method, distributions of
-velocity and speed in different slabs were observed.
+\subsection{Thermal Conductivities}
 
-For each swapping rate, thermal conductivity was calculated in reduced
-unit. The energy flux was calculated similarly to the momentum flux,
-with total unphysical transferred energy ${E_{total}}$ and data collection
-time $t$:
+The energy flux was calculated similarly to the momentum flux, using
+the total non-physical energy transferred (${E_{total}}$) and the data
+collection time $t$:
 \begin{equation}
 J_z = \frac{E_{total}}{2 t L_x L_y}
 \end{equation}
-And the temperature gradient ${\langle\partial T/\partial z\rangle}$
-can be obtained by a linear regression of the temperature
-profile. From the thermal conductivity $\lambda$ calculated, one can
-further convert it into reduced unit ${\lambda^*=\lambda \sigma^2
-  m^{1/2} k_B^{-1}\varepsilon^{-1/2}}$.
+The temperature gradient ${\langle\partial T/\partial z\rangle}$ was
+obtained by a linear regression of the temperature profile. For
+Lennard-Jones simulations, thermal conductivities are reported in
+reduced units (${\lambda^*=\lambda \sigma^2 m^{1/2}
+  k_B^{-1}\varepsilon^{-1/2}}$).
 
-Another series of our simulation is to calculate the interfacial
-thermal conductivity of a Au/H${_2}$O system. Respective calculations of
-water (SPC/E) and gold (QSC) thermal conductivity were performed and
-compared with current results to ensure the validity of
-NIVS-RNEMD. After that, the mixture system was simulated.
+\subsection{Interfacial Thermal Conductivities}
 
-\section{Results And Discussion}
-\subsection{Shear Viscosity}
-Our calculations (Table \ref{shearRate}) shows that scale RNEMD method
-produced comparable shear viscosity to swap RNEMD method. In Table
-\ref{shearRate}, the names of the calculated samples are devided into
-two parts. The first number refers to total slabs in one simulation
-box. The second number refers to the swapping interval in swap method, or
-in scale method the equilvalent swapping interval that the same
-momentum flux would theoretically result in swap method. All the scale
-method results were from simulations that had a scaling interval of 10
-time steps. The average molecular momentum gradients of these samples
-are shown in Figure \ref{shearGrad}.
+For materials with a relatively low interfacial conductance, and in
+cases where the flux between the materials is small, the bulk regions
+on either side of an interface rapidly come to a state in which the
+two phases have relatively homogeneous (but distinct) temperatures.
+In calculating the interfacial thermal conductivity $G$, this
+assumption was made, and the conductance can be approximated as:
 
+\begin{equation}
+G = \frac{E_{total}}{2 t L_x L_y \left( \langle T_{gold}\rangle -
+    \langle T_{water}\rangle \right)}
+\label{interfaceCalc}
+\end{equation}
+where ${E_{total}}$ is the imposed non-physical kinetic energy
+transfer and ${\langle T_{gold}\rangle}$ and ${\langle
+  T_{water}\rangle}$ are the average observed temperature of gold and
+water phases respectively.
+
+\section{Results}
+
+\subsection{Lennard-Jones Fluid}
+2592 Lennard-Jones atoms were placed in an orthorhombic cell
+${10.06\sigma \times 10.06\sigma \times 30.18\sigma}$ on a side.  The
+reduced density ${\rho^* = \rho\sigma^3}$ was thus 0.85, which enabled
+direct comparison between our results and previous methods.  These
+simulations were carried out with a reduced timestep ${\tau^* =
+  4.6\times10^{-4}}$.  For the shear viscosity calculations, the mean
+temperature was ${T^* = k_B T/\varepsilon = 0.72}$.  For thermal
+conductivity calculations, simulations were first run under reduced
+temperature ${\langle T^*\rangle = 0.72}$ in the microcanonical
+ensemble, but other temperatures ([XXX, YYY, and ZZZ]) were also
+sampled.  The simulations included $10^5$ steps of equilibration
+without any momentum flux, $10^5$ steps of stablization with an
+imposed momentum transfer to create a gradient, and $10^6$ steps of
+data collection under RNEMD.
+
+\subsubsection*{Thermal Conductivity}
+
+Our thermal conductivity calculations show that the NIVS method agrees
+well with the swapping method. Four different swap intervals were
+tested (Table \ref{thermalLJRes}). With a fixed 10 fs [WHY NOT REDUCED
+UNITS???]  scaling interval, the target exchange kinetic energy
+produced equivalent kinetic energy flux as in the swapping method.
+Similar thermal gradients were observed with similar thermal flux
+under the two different methods (Figure \ref{thermalGrad}).
+
 \begin{table*}
-\begin{minipage}{\linewidth}
-\begin{center}
+  \begin{minipage}{\linewidth}
+    \begin{center}
 
-\caption{Calculation results for shear viscosity of Lennard-Jones
-  fluid at ${T^* = 0.72}$ and ${\rho^* = 0.85}$, with swap and scale
-  methods at various momentum exchange rates. Results in reduced
-  unit. Errors of calculations in parentheses. }
+      \caption{Thermal conductivity (in reduced units) of a
+        Lennard-Jones fluid at ${\langle T^* \rangle = 0.72}$ and
+        ${\rho^* = 0.85}$ for the swapping and scaling methods at
+        various kinetic energy exchange rates. Uncertainties are
+        indicated in parentheses.}
+      
+      \begin{tabular}{|cc|cc|}
+        \hline
+        \multicolumn{2}{|c|}{Swapping RNEMD} &
+        \multicolumn{2}{|c|}{NIVS-RNEMD} \\
+        \hline
+        Swap Interval (fs) & $\lambda^*_{swap}$ & Equilvalent $J_z^*$ & $\lambda^*_{scale}$\\
+        \hline 
+        250  & 7.03(0.34) & 0.16  & 7.30(0.10)\\
+        500  & 7.03(0.14) & 0.09  & 6.95(0.09)\\
+        1000 & 6.91(0.42) & 0.047 & 7.19(0.07)\\
+        2000 & 7.52(0.15) & 0.024 & 7.19(0.28)\\
+        \hline
+      \end{tabular}
+      \label{thermalLJRes}
+    \end{center}
+  \end{minipage}
+\end{table*}
 
-\begin{tabular}{ccc}
-\hline
-Series & $\eta^*_{swap}$ & $\eta^*_{scale}$\\
-\hline
-20-500 & 3.64(0.05) & 3.76(0.09)\\
-20-1000 & 3.52(0.16) & 3.66(0.06)\\
-20-2000 & 3.72(0.05) & 3.32(0.18)\\
-20-2500 & 3.42(0.06) & 3.43(0.08)\\
-\hline
-\end{tabular}
-\label{shearRate}
-\end{center}
-\end{minipage}
-\end{table*}
+\begin{figure}
+\includegraphics[width=\linewidth]{thermalGrad}
+\caption{NIVS-RNEMD method creates similar temperature gradients
+  compared with the swapping method under a variety of imposed kinetic
+  energy flux values.}
+\label{thermalGrad}
+\end{figure}
 
+During these simulations, velocities were recorded every 1000 steps
+and was used to produce distributions of both velocity and speed in
+each of the slabs. From these distributions, we observed that under
+relatively high non-physical kinetic energy flux, the spee of
+molecules in slabs where swapping occured could deviate from the
+Maxwell-Boltzmann distribution. This behavior was also noted by Tenney
+and Maginn.\cite{Maginn:2010} Figure \ref{thermalHist} shows how these
+distributions deviate from an ideal distribution. In the ``hot'' slab,
+the probability density is notched at low speeds and has a substantial
+shoulder at higher speeds relative to the ideal MB distribution.  In
+the cold slab, the opposite notching and shouldering occurs.  This
+phenomenon is more obvious at higher swapping rates.  
+
+In the velocity distributions, the ideal Gaussian peak is
+substantially flattened in the hot slab, and is overly sharp (with
+truncated wings) in the cold slab. This problem is rooted in the
+mechanism of the swapping method. Continually depleting low (high)
+speed particles in the high (low) temperature slab is not complemented
+by diffusions of low (high) speed particles from neighboring slabs,
+unless the swapping rate is sufficiently small. Simutaneously, surplus
+low speed particles in the low temperature slab do not have sufficient
+time to diffuse to neighboring slabs.  Since the thermal exchange rate
+must reach a minimum level to produce an observable thermal gradient,
+the swapping-method RNEMD has a relatively narrow choice of exchange
+times that can be utilized.
+
+For comparison, NIVS-RNEMD produces a speed distribution closer to the
+Maxwell-Boltzmann curve (Figure \ref{thermalHist}).  The reason for
+this is simple; upon velocity scaling, a Gaussian distribution remains
+Gaussian.  Although a single scaling move is non-isotropic in three
+dimensions, our criteria for choosing a set of scaling coefficients
+helps maintain the distributions as close to isotropic as possible.
+Consequently, NIVS-RNEMD is able to impose an unphysical thermal flux
+as the previous RNEMD methods but without large perturbations to the
+velocity distributions in the two slabs.
+
 \begin{figure}
-\includegraphics[width=\linewidth]{shearGrad}
-\caption{Average momentum gradients of shear viscosity simulations}
-\label{shearGrad}
+\includegraphics[width=\linewidth]{thermalHist}
+\caption{Speed distribution for thermal conductivity using a)
+  ``swapping'' and b) NIVS- RNEMD methods. Shown is from the
+  simulations with an exchange or equilvalent exchange interval of 250
+  fs. In circled areas, distributions from ``swapping'' RNEMD
+  simulation have deviation from ideal Maxwell-Boltzmann distribution
+  (curves fit for each distribution).}
+\label{thermalHist}
 \end{figure}
 
+
+\subsubsection*{Shear Viscosity}
+Our calculations (Table \ref{shearRate}) show that velocity-scaling
+RNEMD predicted comparable shear viscosities to swap RNEMD method.  All
+the scale method results were from simulations that had a scaling
+interval of 10 time steps. The average molecular momentum gradients of
+these samples are shown in Figure \ref{shear} (a) and (b).
+
+\begin{table*}
+  \begin{minipage}{\linewidth}
+    \begin{center}
+      
+      \caption{Shear viscosities of Lennard-Jones fluid at ${T^* =
+          0.72}$ and ${\rho^* = 0.85}$ using swapping and NIVS methods
+        at various momentum exchange rates.  Uncertainties are
+        indicated in parentheses. }
+      
+      \begin{tabular}{ccccc}
+        Swapping method & & & NIVS-RNEMD & \\
+        \hline
+        Swap Interval (fs) & $\eta^*_{swap}$ & & Equilvalent $j_p^*(v_x)$ &
+        $\eta^*_{scale}$\\
+        \hline
+        500  & 3.64(0.05) & & 0.09  & 3.76(0.09)\\
+        1000 & 3.52(0.16) & & 0.046 & 3.66(0.06)\\
+        2000 & 3.72(0.05) & & 0.024 & 3.32(0.18)\\
+        2500 & 3.42(0.06) & & 0.019 & 3.43(0.08)\\
+        \hline
+      \end{tabular}
+      \label{shearRate}
+    \end{center}
+  \end{minipage}
+\end{table*}
+
 \begin{figure}
-\includegraphics[width=\linewidth]{shearTempScale}
-\caption{Temperature profile for scaling RNEMD simulation.}
-\label{shearTempScale}
+  \includegraphics[width=\linewidth]{shear}
+  \caption{Average momentum gradients in shear viscosity simulations,
+    using (a) ``swapping'' method and (b) NIVS-RNEMD method
+    respectively. (c) Temperature difference among x and y, z dimensions
+    observed when using NIVS-RNEMD with equivalent exchange interval of
+    500 fs.}
+  \label{shear}
 \end{figure}
+
 However, observations of temperatures along three dimensions show that
 inhomogeneity occurs in scaling RNEMD simulations, particularly in the
-two slabs which were scaled. Figure \ref{shearTempScale} indicate that with
+two slabs which were scaled. Figure \ref{shear} (c) indicate that with
 relatively large imposed momentum flux, the temperature difference among $x$
 and the other two dimensions was significant. This would result from the
 algorithm of scaling method. From Eq. \ref{eq:fluxPlane}, after
@@ -431,146 +591,285 @@ attractive than swapping RNEMD in shear viscosity calc
 exchange momentum flux simulations makes scaling RNEMD method less
 attractive than swapping RNEMD in shear viscosity calculation. 
 
-\subsection{Thermal Conductivity}
 
-Our thermal conductivity calculations also show that scaling method results
-agree with swapping method. Table \ref{thermal} lists our simulation
-results with similar manner we used in shear viscosity
-calculation. All the data reported from scaling method were obtained
-by simulations of 10-step exchange frequency, and the target exchange
-kinetic energy were set to produce equivalent kinetic energy flux as
-in swapping method. Figure \ref{thermalGrad} exhibits similar thermal
-gradients of respective similar kinetic energy flux.
+\subsection{Bulk SPC/E water}
 
-\begin{table*}
-\begin{minipage}{\linewidth}
-\begin{center}
+We compared the thermal conductivity of SPC/E water using NIVS-RNEMD
+to the work of Bedrov {\it et al.}\cite{Bedrov:2000}, which employed
+the original swapping RNEMD method.  Bedrov {\it et
+  al.}\cite{Bedrov:2000} argued that exchange of the molecule
+center-of-mass velocities instead of single atom velocities in a
+molecule conserves the total kinetic energy and linear momentum.  This
+principle is also adopted in our simulations. Scaling was applied to
+the center-of-mass velocities of the rigid bodies of SPC/E model water
+molecules.
 
-\caption{Calculation results for thermal conductivity of Lennard-Jones
-  fluid at ${\langle T^* \rangle = 0.72}$ and ${\rho^* = 0.85}$, with
-  swap and scale methods at various kinetic energy exchange rates. Results
-  in reduced unit. Errors of calculations in parentheses.} 
+To construct the simulations, a simulation box consisting of 1000
+molecules were first equilibrated under ambient pressure and
+temperature conditions using the isobaric-isothermal (NPT)
+ensemble.\cite{melchionna93} A fixed volume was chosen to match the
+average volume observed in the NPT simulations, and this was followed
+by equilibration, first in the canonical (NVT) ensemble, followed by a
+[XXX ps] period under constant-NVE conditions without any momentum
+flux.  [YYY ps] was allowed to stabilize the system with an imposed
+momentum transfer to create a gradient, and [ZZZ ps] was alotted for
+data collection under RNEMD.
 
-\begin{tabular}{ccc}
-\hline
-Series & $\lambda^*_{swap}$ & $\lambda^*_{scale}$\\
-\hline 
-20-250  & 7.03(0.34) & 7.30(0.10)\\
-20-500  & 7.03(0.14) & 6.95(0.09)\\
-20-1000 & 6.91(0.42) & 7.19(0.07)\\
-20-2000 & 7.52(0.15) & 7.19(0.28)\\
-\hline
-\end{tabular}
-\label{thermal}
-\end{center}
-\end{minipage}
-\end{table*}
+As shown in Figure \ref{spceGrad}, temperature gradients were
+established similar to the previous work.  However, the average
+temperature of our system is 300K, while that in Bedrov {\it et al.}
+is 318K, which would be attributed for part of the difference between
+the final calculation results (Table \ref{spceThermal}). [WHY DIDN'T
+WE DO 318 K?]  Both methods yield values in reasonable agreement with
+experiment [DONE AT WHAT TEMPERATURE?]
 
 \begin{figure}
-\includegraphics[width=\linewidth]{thermalGrad}
-\caption{Temperature gradients of thermal conductivity simulations}
-\label{thermalGrad}
+  \includegraphics[width=\linewidth]{spceGrad}
+  \caption{Temperature gradients in SPC/E water thermal conductivity
+    simulations.}
+  \label{spceGrad}
 \end{figure}
 
-During these simulations, molecule velocities were recorded in 1000 of
-all the snapshots. These velocity data were used to produce histograms
-of velocity and speed distribution in different slabs. From these
-histograms, it is observed that with increasing unphysical kinetic
-energy flux, speed and velocity distribution of molecules in slabs
-where swapping occured could deviate from Maxwell-Boltzmann
-distribution. Figure \ref{histSwap} indicates how these distributions
-deviate from ideal condition. In high temperature slabs, probability
-density in low speed is confidently smaller than ideal distribution;
-in low temperature slabs, probability density in high speed is smaller
-than ideal. This phenomenon is observable even in our relatively low
-swapping rate simulations. And this deviation could also leads to
-deviation of distribution of velocity in various dimensions. One
-feature of these deviated distribution is that in high temperature
-slab, the ideal Gaussian peak was changed into a relatively flat
-plateau; while in low temperature slab, that peak appears sharper.
+\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}{|c|ccc|}
+        \hline
+        $\langle dT/dz\rangle$(K/\AA) & \multicolumn{3}{|c|}{$\lambda
+          (\mathrm{W m}^{-1} \mathrm{K}^{-1})$} \\
+        & This work (300K) & Previous simulations (318K)\cite{Bedrov:2000} &
+        Experiment\cite{WagnerKruse}\\
+        \hline
+        0.38 & 0.816(0.044) & & 0.64\\
+        0.81 & 0.770(0.008) & 0.784 & \\
+        1.54 & 0.813(0.007) & 0.730 & \\
+        \hline
+      \end{tabular}
+      \label{spceThermal}
+    \end{center}
+  \end{minipage}
+\end{table*}
 
-\begin{figure}
-\includegraphics[width=\linewidth]{histSwap}
-\caption{Speed distribution for thermal conductivity using swapping RNEMD.}
-\label{histSwap}
-\end{figure}
+\subsection{Crystalline Gold} 
 
-\begin{figure}
-\includegraphics[width=\linewidth]{histScale}
-\caption{Speed distribution for thermal conductivity using scaling RNEMD.}
-\label{histScale}
-\end{figure}
+To see how the method performed in a solid, we calculated thermal
+conductivities using two atomistic models for gold.  Several different
+potential models have been developed that reasonably describe
+interactions in transition metals. In particular, the Embedded Atom
+Model ({\sc eam})~\cite{PhysRevB.33.7983} and Sutton-Chen ({\sc
+  sc})~\cite{Chen90} potential have been used to study a wide range of
+phenomena in both bulk materials and
+nanoparticles.\cite{ISI:000207079300006,Clancy:1992,Vardeman:2008fk,Vardeman-II:2001jn,ShibataT._ja026764r,Sankaranarayanan:2005lr,Chui:2003fk,Wang:2005qy,Medasani:2007uq}
+Both potentials are based on a model of a metal which treats the
+nuclei and core electrons as pseudo-atoms embedded in the electron
+density due to the valence electrons on all of the other atoms in the
+system. The {\sc sc} potential has a simple form that closely
+resembles the Lennard Jones potential,
+\begin{equation}
+\label{eq:SCP1} 
+U_{tot}=\sum _{i}\left[ \frac{1}{2}\sum _{j\neq i}D_{ij}V^{pair}_{ij}(r_{ij})-c_{i}D_{ii}\sqrt{\rho_{i}}\right] , 
+\end{equation}
+where $V^{pair}_{ij}$ and $\rho_{i}$ are given by 
+\begin{equation}
+\label{eq:SCP2} 
+V^{pair}_{ij}(r)=\left( \frac{\alpha_{ij}}{r_{ij}}\right)^{n_{ij}}, \rho_{i}=\sum_{j\neq i}\left( \frac{\alpha_{ij}}{r_{ij}}\right) ^{m_{ij}}. 
+\end{equation}
+$V^{pair}_{ij}$ is a repulsive pairwise potential that accounts for
+interactions between the pseudoatom cores. The $\sqrt{\rho_i}$ term in
+Eq. (\ref{eq:SCP1}) is an attractive many-body potential that models
+the interactions between the valence electrons and the cores of the
+pseudo-atoms. $D_{ij}$, $D_{ii}$ set the appropriate overall energy
+scale, $c_i$ scales the attractive portion of the potential relative
+to the repulsive interaction and $\alpha_{ij}$ is a length parameter
+that assures a dimensionless form for $\rho$. These parameters are
+tuned to various experimental properties such as the density, cohesive
+energy, and elastic moduli for FCC transition metals. The quantum
+Sutton-Chen ({\sc q-sc}) formulation matches these properties while
+including zero-point quantum corrections for different transition
+metals.\cite{PhysRevB.59.3527}  The {\sc eam} functional forms differ
+slightly from {\sc sc} but the overall method is very similar.
 
-\subsection{Interfaciel Thermal Conductivity}
+In this work, we have utilized both the {\sc eam} and the {\sc q-sc}
+potentials to test the behavior of scaling RNEMD. 
 
-\begin{figure}
-\includegraphics[width=\linewidth]{spceGrad}
-\caption{Temperature gradients for SPC/E water thermal conductivity.}
-\label{spceGrad}
-\end{figure}
+A face-centered-cubic (FCC) lattice was prepared containing 2880 Au
+atoms.  [LxMxN UNIT CELLS].  Simulations were run both with and
+without isobaric-isothermal (NPT)~\cite{melchionna93}
+pre-equilibration at a target pressure of 1 atm. When equilibrated
+under NPT conditions, our simulation box expanded by approximately 1\%
+[IN VOLUME OR LINEAR DIMENSIONS ?].  Following adjustment of the box
+volume, equilibrations in both the canonical and microcanonical
+ensembles were carried out. With the simulation cell divided evenly
+into 10 slabs, different thermal gradients were established by
+applying a set of target thermal transfer fluxes.
 
+The results for the thermal conductivity of gold are shown in Table
+\ref{AuThermal}.  In these calculations, the end and middle slabs were
+excluded in thermal gradient linear regession. {\sc eam} predicts
+slightly larger thermal conductivities than {\sc q-sc}.  However, both
+values are smaller than experimental value by a factor of more than
+200. This behavior has been observed previously by Richardson and
+Clancy, and has been attributed to the lack of electronic effects in
+these force fields.\cite{Clancy:1992} The non-equilibrium MD method
+employed in their simulations produced comparable results to ours.  It
+should be noted that the density of the metal being simulated also
+greatly affects the thermal conductivity.  With an expanded lattice,
+lower thermal conductance is expected (and observed). We also observed
+a decrease in thermal conductance at higher temperatures, a trend that
+agrees with experimental measurements [PAGE
+NUMBERS?].\cite{AshcroftMermin}
+
 \begin{table*}
-\begin{minipage}{\linewidth}
-\begin{center}
+  \begin{minipage}{\linewidth}
+    \begin{center}
+      
+      \caption{Calculated thermal conductivity of crystalline gold
+        using two related force fields. Calculations were done at both
+        experimental and equilibrated densities and at a range of
+        temperatures and thermal flux rates.  Uncertainties are
+        indicated in parentheses. [CLANCY COMPARISON? SWAPPING
+        COMPARISON?]}
+      
+      \begin{tabular}{|c|c|c|cc|}
+        \hline
+        Force Field & $\rho$ (g/cm$^3$) & ${\langle T\rangle}$ (K) &
+        $\langle dT/dz\rangle$ (K/\AA) & $\lambda (\mathrm{W m}^{-1} \mathrm{K}^{-1})$\\
+        \hline
+        \multirow{7}{*}{\sc q-sc} & \multirow{3}{*}{19.188} & \multirow{3}{*}{300} & 1.44 & 1.10(0.06)\\
+        &        &     & 2.86 & 1.08(0.05)\\
+        &        &     & 5.14 & 1.15(0.07)\\\cline{2-5}
+        & \multirow{4}{*}{19.263} & \multirow{2}{*}{300} & 2.31 & 1.25(0.06)\\
+        &        &     & 3.02 & 1.26(0.05)\\\cline{3-5}
+        &        & \multirow{2}{*}{575} & 3.02 & 1.02(0.07)\\
+        &        &     & 4.84 & 0.92(0.05)\\
+        \hline
+        \multirow{8}{*}{\sc eam} & \multirow{3}{*}{19.045} & \multirow{3}{*}{300} & 1.24 & 1.24(0.16)\\
+        &        &     & 2.06 & 1.37(0.04)\\
+        &        &     & 2.55 & 1.41(0.07)\\\cline{2-5}
+        & \multirow{5}{*}{19.263} & \multirow{3}{*}{300} & 1.06 & 1.45(0.13)\\
+        &        &     & 2.04 & 1.41(0.07)\\
+        &        &     & 2.41 & 1.53(0.10)\\\cline{3-5}
+        &        & \multirow{2}{*}{575} & 2.82 & 1.08(0.03)\\
+        &        &     & 4.14 & 1.08(0.05)\\
+        \hline
+      \end{tabular}
+      \label{AuThermal}
+    \end{center}
+  \end{minipage}
+\end{table*}
 
-\caption{Calculation results for thermal conductivity of SPC/E water
-  at ${\langle T\rangle}$ = 300K at various thermal exchange rates. Errors of
-  calculations in parentheses. }
+\subsection{Thermal Conductance at the Au/H$_2$O interface}
+The most attractive aspect of the scaling approach for RNEMD is the
+ability to use the method in non-homogeneous systems, where molecules
+of different identities are segregated in different slabs.  To test
+this application, we simulated a Gold (111) / water interface.  To
+construct the interface, a box containing a lattice of 1188 Au atoms
+(with the 111 surface in the +z and -z directions) was allowed to
+relax under ambient temperature and pressure.  A separate (but
+identically sized) box of SPC/E water was also equilibrated at ambient
+conditions.  The two boxes were combined by removing all water
+molecules withing 3 \AA radius of any gold atom.  The final
+configuration contained 1862 SPC/E water molecules.
 
-\begin{tabular}{cccc}
-\hline
-$\langle dT/dz\rangle$(K/\AA) & & $\lambda$(W/m/K) & \\
-& This work & Previous simulations$^a$ & Experiment$^b$\\
-\hline
-0.3 & 0.82()  & 0.784 & 0.64\\
-0.8 & 0.770() & 0.730\\
-1.5 & 0.813() & \\
-\hline
-\end{tabular}
-\label{spceThermal}
-\end{center}
-\end{minipage}
-\end{table*}
+After simulations of bulk water and crystal gold, a mixture system was
+constructed, consisting of 1188 Au atoms and 1862 H$_2$O
+molecules. Spohr potential was adopted in depicting the interaction
+between metal atom and water molecule.\cite{ISI:000167766600035} A
+similar protocol of equilibration was followed. Several thermal
+gradients was built under different target thermal flux. It was found
+out that compared to our previous simulation systems, the two phases
+could have large temperature difference even under a relatively low
+thermal flux.
 
 
+After simulations of homogeneous water and gold systems using
+NIVS-RNEMD method were proved valid, calculation of gold/water
+interfacial thermal conductivity was followed. It is found out that
+the low interfacial conductance is probably due to the hydrophobic
+surface in our system. Figure \ref{interface} (a) demonstrates mass
+density change along $z$-axis, which is perpendicular to the
+gold/water interface. It is observed that water density significantly
+decreases when approaching the surface. Under this low thermal
+conductance, both gold and water phase have sufficient time to
+eliminate temperature difference inside respectively (Figure
+\ref{interface} b). With indistinguishable temperature difference
+within respective phase, it is valid to assume that the temperature
+difference between gold and water on surface would be approximately
+the same as the difference between the gold and water phase. This
+assumption enables convenient calculation of $G$ using
+Eq. \ref{interfaceCalc} instead of measuring temperatures of thin
+layer of water and gold close enough to surface, which would have
+greater fluctuation and lower accuracy. Reported results (Table
+\ref{interfaceRes}) are of two orders of magnitude smaller than our
+calculations on homogeneous systems, and thus have larger relative
+errors than our calculation results on homogeneous systems.
+
 \begin{figure}
-\includegraphics[width=\linewidth]{AuGrad}
-\caption{Temperature gradients for crystal gold thermal conductivity.}
-\label{AuGrad}
+\includegraphics[width=\linewidth]{interface}
+\caption{Simulation results for Gold/Water interfacial thermal
+  conductivity: (a) Significant water density decrease is observed on
+  crystalline gold surface, which indicates low surface contact and
+  leads to low thermal conductance. (b) Temperature profiles for a
+  series of simulations. Temperatures of different slabs in the same
+  phase show no significant differences.}
+\label{interface}
 \end{figure}
 
 \begin{table*}
 \begin{minipage}{\linewidth}
 \begin{center}
 
-\caption{Calculation results for thermal conductivity of crystal gold
-  at ${\langle T\rangle}$ = 300K at various thermal exchange rates. Errors of
-  calculations in parentheses. }
+\caption{Calculation results for interfacial thermal conductivity
+  at ${\langle T\rangle \sim}$ 300K at various thermal exchange
+  rates. Errors of calculations in parentheses. }
 
-\begin{tabular}{ccc}
+\begin{tabular}{cccc}
 \hline
-$\langle dT/dz\rangle$(K/\AA) & $\lambda$(W/m/K)\\
-& This work & Previous simulations$^a$ \\
+$J_z$ (MW/m$^2$) & $T_{gold}$ (K) & $T_{water}$ (K) & $G$
+(MW/m$^2$/K)\\
 \hline
-1.4 & 1.10() & \\
-2.8 & 1.08() & \\
-5.1 & 1.15() & \\
+98.0 & 355.2 & 295.8 & 1.65(0.21) \\
+78.8 & 343.8 & 298.0 & 1.72(0.32) \\
+73.6 & 344.3 & 298.0 & 1.59(0.24) \\
+49.2 & 330.1 & 300.4 & 1.65(0.35) \\
 \hline
 \end{tabular}
-\label{AuThermal}
+\label{interfaceRes}
 \end{center}
 \end{minipage}
 \end{table*}
 
 
+\section{Conclusions}
+NIVS-RNEMD simulation method is developed and tested on various
+systems. Simulation results demonstrate its validity in thermal
+conductivity calculations, from Lennard-Jones fluid to multi-atom
+molecule like water and metal crystals. NIVS-RNEMD improves
+non-Boltzmann-Maxwell distributions, which exist in previous RNEMD
+methods. Furthermore, it develops a valid means for unphysical thermal
+transfer between different species of molecules, and thus extends its
+applicability to interfacial systems. Our calculation of gold/water
+interfacial thermal conductivity demonstrates this advantage over
+previous RNEMD methods. NIVS-RNEMD has also limited application on
+shear viscosity calculations, but could cause temperature difference
+among different dimensions under high momentum flux. Modification is
+necessary to extend the applicability of NIVS-RNEMD in shear viscosity
+calculations.
+
 \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
 
-\bibliographystyle{jcp2}
+\bibliographystyle{aip}
 \bibliography{nivsRnemd}
+
 \end{doublespace}
 \end{document}
 
