--- trunk/nivsRnemd/nivsRnemd.tex	2010/07/30 21:13:43	3619
+++ trunk/nivsRnemd/nivsRnemd.tex	2010/08/12 14:48:03	3630
@@ -44,20 +44,19 @@ Notre Dame, Indiana 46556}
 
 \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
+  velocity and temperature gradients in molecular dynamics simulations
+  of heterogeneous systems.  This method extends earlier Reverse
+  Non-Equilibrium Molecular Dynamics (RNEMD) methods which use
+  momentum exchange swapping moves. The standard swapping moves 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 while conserving the 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}
 
@@ -123,14 +122,15 @@ Recently, Tenney and Maginn\cite{Maginn:2010} have dis
 typically samples from the same manifold of states in the
 microcanonical ensemble.
 
-Recently, Tenney and Maginn\cite{Maginn:2010} have discovered
-some problems with the original RNEMD swap technique.  Notably, large
+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.
+and proposed 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) which retains the desirable features of RNEMD
@@ -181,13 +181,13 @@ P_c^\alpha & = & \sum_{i = 1}^{N_c} m_i \left[\vec{v}_
 \end{equation}
 where
 \begin{eqnarray}
-P_c^\alpha & = & \sum_{i = 1}^{N_c} m_i \left[\vec{v}_i\right]_\alpha \\
-P_h^\alpha & = & \sum_{j = 1}^{N_h} m_j \left[\vec{v}_j\right]_\alpha
+P_c^\alpha & = & \sum_{i = 1}^{N_c} m_i v_{i\alpha} \\
+P_h^\alpha & = & \sum_{j = 1}^{N_h} m_j v_{j\alpha}
 \label{eq:momentumdef}
 \end{eqnarray}
 Therefore, for each of the three directions, the hot scaling
 parameters are a simple function of the cold scaling parameters and
-the instantaneous linear momentum in each of the two slabs.
+the instantaneous linear momenta in each of the two slabs.
 \begin{equation}
 \alpha^\prime = 1 + (1 - \alpha) p_\alpha
 \label{eq:hotcoldscaling}
@@ -200,8 +200,8 @@ Conservation of total energy also places constraints o
 
 Conservation of total energy also places constraints on the scaling:
 \begin{equation}
-\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
+\sum_{\alpha = x,y,z} \left\{ K_h^\alpha + K_c^\alpha \right\} = \sum_{\alpha = x,y,z}
+\left\{ \left(\alpha^\prime\right)^2 K_h^\alpha + \alpha^2 K_c^\alpha \right\}
 \end{equation}
 where the translational kinetic energies, $K_h^\alpha$ and
 $K_c^\alpha$, are computed for each of the three directions in a
@@ -269,10 +269,10 @@ problem.\cite{Hoffman:2001sf} One simplification is to
 degree 16.  There are a number of polynomial root-finding methods in
 the literature,\cite{Hoffman:2001sf,384119} but numerically finding
 the roots of high-degree polynomials is generally an ill-conditioned
-problem.\cite{Hoffman:2001sf} 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
+problem.\cite{Hoffman:2001sf} One simplification is to maintain
+velocity scalings that are {\it as isotropic as possible}.  To do
+this, we impose $x=y$, and 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.
 
@@ -330,15 +330,16 @@ after an MD step with a variable frequency.  We have t
 
 We have implemented this methodology in our molecular dynamics code,
 OpenMD,\cite{Meineke:2005gd,openmd} performing the NIVS scaling moves
-after an MD step with a variable frequency.  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 ({\sc q-sc} 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$).
+with a variable frequency after the molecular dynamics (MD) steps.  We
+have tested the method in a variety of different systems, including:
+homogeneous fluids (Lennard-Jones and SPC/E water), crystalline
+solids, using both the embedded atom method
+(EAM)~\cite{PhysRevB.33.7983} and quantum Sutton-Chen
+(QSC)~\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 the
+current method can easily provide estimates of the interfacial thermal
+conductivity ($G$).
 
 \subsection{Simulation Cells}
 
@@ -357,18 +358,37 @@ first performed simulations using the original techniq
 
 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.
+first performed simulations using the original technique. At fixed
+intervals, kinetic energy or momentum exchange moves were performed
+between the hot and the cold slabs.  The interval between exchange
+moves governs the effective momentum flux ($j_z(p_x)$) or energy flux
+($J_z$) between the two slabs so to vary this quantity, we performed
+simulations with a variety of delay intervals between the swapping moves.
 
+For thermal conductivity measurements, the particle with smallest
+speed in the hot slab and the one with largest speed in the cold slab
+had their entire momentum vectors swapped.  In the test cases run
+here, all particles had the same chemical identity and mass, so this
+move preserves both total linear momentum and total energy.  It is
+also possible to exchange energy by assuming an elastic collision
+between the two particles which are exchanging energy.
+
+For shear stress simulations, the particle with the most negative
+$p_x$ in the hot slab and the one with the most positive $p_x$ in the
+cold slab exchanged only this component of their momentum vectors.
+
 \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 produce the same flux experienced 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.
+To test the temperature homogeneity, directional momentum and
+temperature distributions were accumulated for molecules in each of
+the slabs.  Transport coefficients were computed using the temperature
+(and momentum) gradients across the $z$-axis as well as the total
+momentum flux the system experienced during data collection portion of
+the simulation.
 
 \subsection{Shear viscosities}
 
@@ -379,18 +399,18 @@ momentum transfer occurs in two directions due to our 
 \end{equation}
 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}}$).
+momentum transfer between the slabs occurs in two directions ($+z$ and
+$-z$).  The velocity gradient ${\langle \partial v_x /\partial z
+  \rangle}$ was obtained using linear regression of the mean $x$
+component of the velocity, $\langle v_x \rangle$, in each of the bins.
+For Lennard-Jones simulations, shear viscosities are reported in
+reduced units (${\eta^* = \eta \sigma^2 (\varepsilon m)^{-1/2}}$).
 
 \subsection{Thermal Conductivities}
 
-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$:
+The energy flux was calculated in a similar manner 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}
@@ -402,12 +422,11 @@ For materials with a relatively low interfacial conduc
 
 \subsection{Interfacial Thermal Conductivities}
 
-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:
+For interfaces with a relatively low interfacial conductance, 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.  The interfacial thermal conductivity $G$ can therefore
+be approximated as:
 
 \begin{equation}
 G = \frac{E_{total}}{2 t L_x L_y \left( \langle T_{gold}\rangle -
@@ -417,7 +436,11 @@ water phases respectively.
 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.
+water phases respectively.  If the interfacial conductance is {\it
+  not} small, it is also be possible to compute the interfacial
+thermal conductivity using this method utilizing the change in the
+slope of the thermal gradient ($\partial^2 \langle T \rangle / \partial
+z^2$) at the interface.
 
 \section{Results}
 
@@ -440,14 +463,13 @@ tested (Table \ref{LJ}). With a fixed scaling interval
 
 Our thermal conductivity calculations show that the NIVS method agrees
 well with the swapping method. Five different swap intervals were
-tested (Table \ref{LJ}). With a fixed scaling interval of 10 time steps,
-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}). Furthermore, with appropriate choice of
-scaling variables, temperature along $x$, $y$ and $z$ axis has no
-observable difference(Figure TO BE ADDED). The system is able
-to maintain temperature homogeneity even under high flux.
+tested (Table \ref{LJ}). Similar thermal gradients were observed with
+similar thermal flux under the two different methods (Figure
+\ref{thermalGrad}). Furthermore, the 1-d temperature profiles showed
+no observable differences between the $x$, $y$ and $z$ axes (Figure
+\ref{thermalGrad} c), so even though we are using a non-isotropic
+scaling method, none of the three directions are experience
+disproportionate heating due to the velocity scaling.
 
 \begin{table*}
   \begin{minipage}{\linewidth}
@@ -487,42 +509,41 @@ to maintain temperature homogeneity even under high fl
 
 \begin{figure}
   \includegraphics[width=\linewidth]{thermalGrad}
-  \caption{NIVS-RNEMD method (b) creates similar temperature gradients
-    compared with the swapping method (a) under a variety of imposed
+  \caption{The NIVS-RNEMD method creates similar temperature gradients
+    compared with the swapping method under a variety of imposed
     kinetic energy flux values. Furthermore, the implementation of
     Non-Isotropic Velocity Scaling does not cause temperature
-    difference among the three dimensions (c).}
+    anisotropy to develop in thermal conductivity calculations.}
   \label{thermalGrad}
 \end{figure}
 
 \subsubsection*{Velocity Distributions}
 
 During these simulations, velocities were recorded every 1000 steps
-and was used to produce distributions of both velocity and speed in
+and were 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 speed 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} in their simulations for shear viscosity
-calculations. 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
+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.
+The peak of the velocity distribution 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
@@ -536,55 +557,55 @@ velocity distributions in the two slabs.
 
 \begin{figure}
 \includegraphics[width=\linewidth]{thermalHist}
-\caption{Speed distribution for thermal conductivity using
-  ``swapping'' and NIVS-RNEMD methods. Shown is from simulations under
-  ${\langle T^* \rangle = 0.8}$ with an swapping interval of 200
-  time steps (equivalent ${J_z^*\sim 0.2}$). In circled areas,
-  distributions from ``swapping'' RNEMD simulation have deviations
-  from ideal Maxwell-Boltzmann distributions.}
+\caption{Velocity and speed distributions that develop under the
+  swapping and NIVS-RNEMD methods at high flux.  The distributions for
+  the hot bins (upper panels) and cold bins (lower panels) were
+  obtained from Lennard-Jones simulations with $\langle T^* \rangle =
+  4.5$ with a flux of $J_z^* \sim 5$ (equivalent to a swapping interval
+  of 10 time steps).  This is a relatively large flux which shows the
+  non-thermal distributions that develop under the swapping method.
+  NIVS does a better job of producing near-thermal distributions in
+  the bins.}
 \label{thermalHist}
 \end{figure}
 
 
 \subsubsection*{Shear Viscosity}
-Our calculations (Table \ref{LJ}) 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).
+Our calculations (Table \ref{LJ}) show that velocity-scaling RNEMD
+predicted comparable shear viscosities to swap RNEMD method. The
+average molecular momentum gradients of these samples are shown in
+Figure \ref{shear} (a) and (b).
 
 \begin{figure}
   \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 ${j_z^*(p_x)\sim 0.09}$.}
+    using ``swapping'' method (top panel) and NIVS-RNEMD method
+    (middle panel). NIVS-RNEMD produces a thermal anisotropy artifact
+    in the hot and cold bins when used for shear viscosity.  This
+    artifact does not appear in thermal conductivity calculations.}
   \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{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
-momentum gradient is set up, $P_c^x$ would be roughly stable
-($<0$). Consequently, scaling factor $x$ would most probably larger
-than 1. Therefore, the kinetic energy in $x$-dimension $K_c^x$ would
-keep increase after most scaling steps. And if there is not enough time
-for the kinetic energy to exchange among different dimensions and
-different slabs, the system would finally build up temperature
-(kinetic energy) difference among the three dimensions. Also, between
-$y$ and $z$ dimensions in the scaled slabs, temperatures of $z$-axis
-are closer to neighbor slabs. This is due to momentum transfer along
-$z$ dimension between slabs.
+Observations of the three one-dimensional temperatures in each of the
+slabs shows that NIVS-RNEMD does produce some thermal anisotropy,
+particularly in the hot and cold slabs.  Figure \ref{shear} (c)
+indicates that with a relatively large imposed momentum flux,
+$j_z(p_x)$, the $x$ direction approaches a different temperature from
+the $y$ and $z$ directions in both the hot and cold bins.  This is an
+artifact of the scaling constraints.  After the momentum gradient has
+been established, $P_c^x < 0$.  Consequently, the scaling factor $x$
+is nearly always greater than one in order to satisfy the constraints.
+This will continually increase the kinetic energy in $x$-dimension,
+$K_c^x$.  If there is not enough time for the kinetic energy to
+exchange among different directions and different slabs, the system
+will exhibit the observed thermal anisotropy in the hot and cold bins.
 
 Although results between scaling and swapping methods are comparable,
-the inherent temperature inhomogeneity even in relatively low imposed
-exchange momentum flux simulations makes scaling RNEMD method less
-attractive than swapping RNEMD in shear viscosity calculation. 
+the inherent temperature anisotropy does make NIVS-RNEMD method less
+attractive than swapping RNEMD for shear viscosity calculations.  We
+note that this problem appears only when momentum flux is applied, and
+does not appear in thermal flux calculations.
 
-
 \subsection{Bulk SPC/E water}
 
 We compared the thermal conductivity of SPC/E water using NIVS-RNEMD
@@ -593,7 +614,7 @@ principle is also adopted in our simulations. Scaling 
   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
+principle is also adopted Fin our simulations. Scaling was applied to
 the center-of-mass velocities of the rigid bodies of SPC/E model water
 molecules.
 
@@ -603,26 +624,17 @@ by equilibration, first in the canonical (NVT) ensembl
 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
-100ps period under constant-NVE conditions without any momentum
-flux. 100ps was allowed to stabilize the system with an imposed
-momentum transfer to create a gradient, and 1ns was alotted for
-data collection under RNEMD.
+100~ps period under constant-NVE conditions without any momentum flux.
+Another 100~ps was allowed to stabilize the system with an imposed
+momentum transfer to create a gradient, and 1~ns was allotted for data
+collection under RNEMD.
 
-In our simulations, temperature gradients were established similar to
-the previous work. Our simulation results under 318K are in well
+In our simulations, the established temperature gradients were similar
+to the previous work.  Our simulation results at 318K are in good
 agreement with those from Bedrov {\it et al.} (Table
 \ref{spceThermal}). And both methods yield values in reasonable
-agreement with experimental value. A larger difference between
-simulation result and experiment is found under 300K. This could
-result from the force field that is used in simulation.
+agreement with experimental values.  
 
-\begin{figure}
-  \includegraphics[width=\linewidth]{spceGrad}
-  \caption{Temperature gradients in SPC/E water thermal conductivity
-    simulations.}
-  \label{spceGrad}
-\end{figure}
-
 \begin{table*}
   \begin{minipage}{\linewidth}
     \begin{center}
@@ -661,15 +673,15 @@ Model ({\sc eam})~\cite{PhysRevB.33.7983} and Sutton-C
 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
+Model (EAM)~\cite{PhysRevB.33.7983} and Sutton-Chen (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,
+system. The 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] , 
@@ -689,13 +701,13 @@ Sutton-Chen ({\sc q-sc}) formulation matches these pro
 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.
+Sutton-Chen (QSC) formulation matches these properties while including
+zero-point quantum corrections for different transition
+metals.\cite{PhysRevB.59.3527} The EAM functional forms differ
+slightly from SC but the overall method is very similar.
 
-In this work, we have utilized both the {\sc eam} and the {\sc q-sc}
-potentials to test the behavior of scaling RNEMD. 
+In this work, we have utilized both the EAM and the QSC potentials to
+test the behavior of scaling RNEMD.
 
 A face-centered-cubic (FCC) lattice was prepared containing 2880 Au
 atoms (i.e. ${6\times 6\times 20}$ unit cells). Simulations were run
@@ -710,22 +722,17 @@ excluded in thermal gradient linear regession. {\sc ea
 
 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 contribution
-in these force fields.\cite{Clancy:1992} The non-equilibrium MD method
-employed in their simulations was only able to give a rough estimation
-of thermal conductance for {\sc eam} gold, and the result was an
-average over a wide temperature range (300-800K). Comparatively, our
-results were based on measurements with linear temperature gradients,
-and thus of higher reliability and accuracy. It should be noted that
-the density of the metal being simulated also has an observable effect
-on thermal conductance. 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.\cite{AshcroftMermin}
+excluded in thermal gradient linear regession. EAM predicts slightly
+larger thermal conductivities than QSC.  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 contribution in these
+force fields.\cite{Clancy:1992} It should be noted that the density of
+the metal being simulated has an effect on thermal conductance. 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.\cite{AshcroftMermin}
 
 \begin{table*}
   \begin{minipage}{\linewidth}
@@ -736,15 +743,16 @@ experimental measurements.\cite{AshcroftMermin}
         experimental and equilibrated densities and at a range of
         temperatures and thermal flux rates. Uncertainties are
         indicated in parentheses. Richardson {\it et
-          al.}\cite{Clancy:1992} gave an estimatioin for {\sc eam} gold
-        of 1.74$\mathrm{W m}^{-1}\mathrm{K}^{-1}$.}
+          al.}\cite{Clancy:1992} give an estimate of 1.74 $\mathrm{W
+          m}^{-1}\mathrm{K}^{-1}$ for EAM gold
+         at a density of 19.263 g / cm$^3$.}
       
       \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)\\
+        \multirow{7}{*}{QSC} & \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)\\
@@ -752,7 +760,7 @@ experimental measurements.\cite{AshcroftMermin}
         &        & \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)\\
+        \multirow{8}{*}{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)\\
@@ -780,34 +788,35 @@ After simulations of bulk water and crystal gold, a mi
 molecules within 3 \AA radius of any gold atom.  The final
 configuration contained 1862 SPC/E water molecules.
 
-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.
+The Spohr potential was adopted in depicting the interaction between
+metal atoms and water molecules.\cite{ISI:000167766600035} A similar
+protocol of equilibration to our water simulations was followed.  We
+observed that the two phases developed large temperature differences
+even under a relatively low thermal flux.
 
+The low interfacial conductance is probably due to an acoustic
+impedance mismatch between the solid and the liquid
+phase.\cite{Cahill:793,RevModPhys.61.605} Experiments on the thermal
+conductivity of gold nanoparticles and nanorods in solvent and in
+glass cages have predicted values for $G$ between 100 and 350
+(MW/m$^2$/K).  The experiments typically have multiple gold surfaces
+that have been protected by a capping agent (citrate or CTAB) or which
+are in direct contact with various glassy solids.  In these cases, the
+acoustic impedance mismatch would be substantially reduced, leading to
+much higher interfacial conductances.  It is also possible, however,
+that the lack of electronic effects that gives rise to the low thermal
+conductivity of EAM gold is also causing a low reading for this
+particular interface.
 
-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
+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
@@ -815,12 +824,11 @@ errors than our calculation results on homogeneous sys
 
 \begin{figure}
 \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.}
+\caption{Temperature profiles of the Gold / Water interface at four
+  different values for the thermal flux.  Temperatures for slabs
+  either in the gold or in the water show no significant differences,
+  although there is a large discontinuity between the materials
+  because of the relatively low interfacial thermal conductivity.}
 \label{interface}
 \end{figure}
 
@@ -868,10 +876,12 @@ Support for this project was provided by the National 
 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
+The authors would like to thank Craig Tenney and Ed Maginn for many
+helpful discussions.  Support for this project was provided by the
+National Science Foundation under grant CHE-0848243. Computational
+time was provided by the Center for Research Computing (CRC) at the
+University of Notre Dame.  
+\newpage
 
 \bibliography{nivsRnemd}