--- interfacial/interfacial.tex	2011/06/29 18:14:23	3729
+++ interfacial/interfacial.tex	2011/07/12 22:11:11	3737
@@ -45,20 +45,22 @@ We have developed a Non-Isotropic Velocity Scaling alg
 
 \begin{abstract}
 
-We have developed a Non-Isotropic Velocity Scaling algorithm for
-setting up and maintaining stable thermal gradients in non-equilibrium
-molecular dynamics simulations. This approach effectively imposes
-unphysical thermal flux even between particles of different
-identities, conserves linear momentum and kinetic energy, and
-minimally perturbs the velocity profile of a system when compared with
-previous RNEMD methods. We have used this method to simulate thermal
-conductance at metal / organic solvent interfaces both with and
-without the presence of thiol-based capping agents.  We obtained
-values comparable with experimental values, and observed significant
-conductance enhancement with the presence of capping agents. Computed
-power spectra indicate the acoustic impedance mismatch between metal
-and liquid phase is greatly reduced by the capping agents and thus
-leads to higher interfacial thermal transfer efficiency.
+With the Non-Isotropic Velocity Scaling algorithm (NIVS) we have
+developed, an unphysical thermal flux can be effectively set up even
+for non-homogeneous systems like interfaces in non-equilibrium
+molecular dynamics simulations. In this work, this algorithm is
+applied for simulating thermal conductance at metal / organic solvent
+interfaces with various coverages of butanethiol capping
+agents. Different solvents and force field models were tested. Our
+results suggest that the United-Atom models are able to provide an
+estimate of the interfacial thermal conductivity comparable to
+experiments in our simulations with satisfactory computational
+efficiency. From our results, the acoustic impedance mismatch between
+metal and liquid phase is effectively reduced by the capping
+agents, and thus leads to interfacial thermal conductance
+enhancement. Furthermore, this effect is closely related to the
+capping agent coverage on the metal surfaces and the type of solvent
+molecules, and is affected by the models used in the simulations.
 
 \end{abstract}
 
@@ -71,66 +73,103 @@ leads to higher interfacial thermal transfer efficienc
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
 \section{Introduction}
-[BACKGROUND FOR INTERFACIAL THERMAL CONDUCTANCE PROBLEM]
 Interfacial thermal conductance is extensively studied both
-experimentally and computationally, and systems with interfaces
-present are generally heterogeneous. Although interfaces are commonly
-barriers to heat transfer, it has been
-reported\cite{doi:10.1021/la904855s} that under specific circustances,
-e.g. with certain capping agents present on the surface, interfacial
-conductance can be significantly enhanced. However, heat conductance
-of molecular and nano-scale interfaces will be affected by the
-chemical details of the surface and is challenging to
-experimentalist. The lower thermal flux through interfaces is even
-more difficult to measure with EMD and forward NEMD simulation
-methods. Therefore, developing good simulation methods will be
-desirable in order to investigate thermal transport across interfaces.
+experimentally and computationally\cite{cahill:793}, due to its
+importance in nanoscale science and technology. Reliability of
+nanoscale devices depends on their thermal transport
+properties. Unlike bulk homogeneous materials, nanoscale materials
+features significant presence of interfaces, and these interfaces
+could dominate the heat transfer behavior of these
+materials. Furthermore, these materials are generally heterogeneous,
+which challenges traditional research methods for homogeneous
+systems.
 
+Heat conductance of molecular and nano-scale interfaces will be
+affected by the chemical details of the surface. Experimentally,
+various interfaces have been investigated for their thermal
+conductance properties. Wang {\it et al.} studied heat transport
+through long-chain hydrocarbon monolayers on gold substrate at
+individual molecular level\cite{Wang10082007}; Schmidt {\it et al.}
+studied the role of CTAB on thermal transport between gold nanorods
+and solvent\cite{doi:10.1021/jp8051888}; Juv\'e {\it et al.} studied
+the cooling dynamics, which is controlled by thermal interface
+resistence of glass-embedded metal
+nanoparticles\cite{PhysRevB.80.195406}. Although interfaces are
+commonly barriers for heat transport, Alper {\it et al.} suggested
+that specific ligands (capping agents) could completely eliminate this
+barrier ($G\rightarrow\infty$)\cite{doi:10.1021/la904855s}.
+
+Theoretical and computational models have also been used to study the
+interfacial thermal transport in order to gain an understanding of
+this phenomena at the molecular level. Recently, Hase and coworkers
+employed Non-Equilibrium Molecular Dynamics (NEMD) simulations to
+study thermal transport from hot Au(111) substrate to a self-assembled
+monolayer of alkylthiolate with relatively long chain (8-20 carbon
+atoms)\cite{hase:2010,hase:2011}. However, ensemble averaged
+measurements for heat conductance of interfaces between the capping
+monolayer on Au and a solvent phase has yet to be studied. 
+The relatively low thermal flux through interfaces is
+difficult to measure with Equilibrium MD or forward NEMD simulation
+methods. Therefore, the Reverse NEMD (RNEMD) methods would have the
+advantage of having this difficult to measure flux known when studying
+the thermal transport across interfaces, given that the simulation
+methods being able to effectively apply an unphysical flux in
+non-homogeneous systems.
+
 Recently, we have developed the Non-Isotropic Velocity Scaling (NIVS)
 algorithm for RNEMD simulations\cite{kuang:164101}. This algorithm
 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. Furthermore, it allows more effective thermal exchange
-between particles of different identities, and thus enables extensive
-study of interfacial conductance.
+the two slabs. Furthermore, it allows effective thermal exchange
+between particles of different identities, and thus makes the study of
+interfacial conductance much simpler.
 
+The work presented here deals with the Au(111) surface covered to
+varying degrees by butanethiol, a capping agent with short carbon
+chain, and solvated with organic solvents of different molecular
+properties. Different models were used for both the capping agent and
+the solvent force field parameters. Using the NIVS algorithm, the
+thermal transport across these interfaces was studied and the
+underlying mechanism for this phenomena was investigated.
+
+[MAY ADD WHY STUDY AU-THIOL SURFACE; CITE SHAOYI JIANG]
+
 \section{Methodology}
-\subsection{Algorithm}
-[BACKGROUND FOR MD METHODS]
-There have been many algorithms for computing thermal conductivity
-using molecular dynamics simulations. However, interfacial conductance
-is at least an order of magnitude smaller. This would make the
-calculation even more difficult for those slowly-converging
-equilibrium methods. Imposed-flux non-equilibrium
+\subsection{Imposd-Flux Methods in MD Simulations}
+For systems with low interfacial conductivity one must have a method
+capable of generating relatively small fluxes, compared to those
+required for bulk conductivity. This requirement makes the calculation
+even more difficult for those slowly-converging equilibrium
+methods\cite{Viscardy:2007lq}.
+Forward methods impose gradient, but in interfacail conditions it is
+not clear what behavior to impose at the boundary...
+ Imposed-flux reverse non-equilibrium
 methods\cite{MullerPlathe:1997xw} have the flux set {\it a priori} and
-the response of temperature or momentum gradients are easier to
-measure than the flux, if unknown, and thus, is a preferable way to
-the forward NEMD methods. Although the momentum swapping approach for
-flux-imposing can be used for exchanging energy between particles of
-different identity, the kinetic energy transfer efficiency is affected
-by the mass difference between the particles, which limits its
-application on heterogeneous interfacial systems.
+the thermal response becomes easier to
+measure than the flux. Although M\"{u}ller-Plathe's original momentum
+swapping approach can be used for exchanging energy between particles
+of different identity, the kinetic energy transfer efficiency is
+affected by the mass difference between the particles, which limits
+its application on heterogeneous interfacial systems.
 
-The non-isotropic velocity scaling (NIVS)\cite{kuang:164101} approach in
-non-equilibrium MD simulations is able to impose relatively large
-kinetic energy flux without obvious perturbation to the velocity
-distribution of the simulated systems. Furthermore, this approach has
+The non-isotropic velocity scaling (NIVS)\cite{kuang:164101} approach to
+non-equilibrium MD simulations is able to impose a wide range of
+kinetic energy fluxes without obvious perturbation to the velocity
+distributions of the simulated systems. Furthermore, this approach has
 the advantage in heterogeneous interfaces in that kinetic energy flux
 can be applied between regions of particles of arbitary identity, and
-the flux quantity is not restricted by particle mass difference.
+the flux will not be restricted by difference in particle mass.
 
 The NIVS algorithm scales the velocity vectors in two separate regions
 of a simulation system with respective diagonal scaling matricies. To
 determine these scaling factors in the matricies, a set of equations
 including linear momentum conservation and kinetic energy conservation
-constraints and target momentum/energy flux satisfaction is
-solved. With the scaling operation applied to the system in a set
-frequency, corresponding momentum/temperature gradients can be built,
-which can be used for computing transportation properties and other
-applications related to momentum/temperature gradients. The NIVS
-algorithm conserves momenta and energy and does not depend on an
-external thermostat. 
+constraints and target energy flux satisfaction is solved. With the
+scaling operation applied to the system in a set frequency, bulk
+temperature gradients can be easily established, and these can be used
+for computing thermal conductivities. The NIVS algorithm conserves
+momenta and energy and does not depend on an external thermostat.
 
 \subsection{Defining Interfacial Thermal Conductivity $G$}
 For interfaces with a relatively low interfacial conductance, the bulk
@@ -148,13 +187,14 @@ When the interfacial conductance is {\it not} small, t
   T_\mathrm{cold}\rangle}$ are the average observed temperature of the
 two separated phases.
 
-When the interfacial conductance is {\it not} small, two ways can be
-used to define $G$.
+When the interfacial conductance is {\it not} small, there are two
+ways to define $G$.
 
-One way is to assume the temperature is discretely different on two
-sides of the interface, $G$ can be calculated with the thermal flux
-applied $J$ and the maximum temperature difference measured along the
-thermal gradient max($\Delta T$), which occurs at the interface, as:
+One way is to assume the temperature is discrete on the two sides of
+the interface. $G$ can be calculated using the applied thermal flux
+$J$ and the maximum temperature difference measured along the thermal
+gradient max($\Delta T$), which occurs at the Gibbs deviding surface,
+as:
 \begin{equation}
 G=\frac{J}{\Delta T}
 \label{discreteG}
@@ -175,16 +215,16 @@ difference method and thus calculate $G^\prime$.
 
 With the temperature profile obtained from simulations, one is able to
 approximate the first and second derivatives of $T$ with finite
-difference method and thus calculate $G^\prime$.
+difference methods and thus calculate $G^\prime$.
 
-In what follows, both definitions are used for calculation and comparison.
+In what follows, both definitions have been used for calculation and
+are compared in the results.
 
-[IMPOSE G DEFINITION INTO OUR SYSTEMS]
-To facilitate the use of the above definitions in calculating $G$ and
-$G^\prime$, we have a metal slab with its (111) surfaces perpendicular
-to the $z$-axis of our simulation cells. With or withour capping
-agents on the surfaces, the metal slab is solvated with organic
-solvents, as illustrated in Figure \ref{demoPic}.
+To compare the above definitions ($G$ and $G^\prime$), we have modeled
+a metal slab with its (111) surfaces perpendicular to the $z$-axis of
+our simulation cells. Both with and withour capping agents on the
+surfaces, the metal slab is solvated with simple organic solvents, as
+illustrated in Figure \ref{demoPic}.
 
 \begin{figure}
 \includegraphics[width=\linewidth]{demoPic}
@@ -193,16 +233,14 @@ With a simulation cell setup following the above manne
 \label{demoPic}
 \end{figure}
 
-With a simulation cell setup following the above manner, one is able
-to equilibrate the system and impose an unphysical thermal flux
-between the liquid and the metal phase with the NIVS algorithm. Under
-a stablized thermal gradient induced by periodically applying the
-unphysical flux, one is able to obtain a temperature profile and the
-physical thermal flux corresponding to it, which equals to the
-unphysical flux applied by NIVS. These data enables the evaluation of
-the interfacial thermal conductance of a surface. Figure \ref{gradT}
-is an example how those stablized thermal gradient can be used to
-obtain the 1st and 2nd derivatives of the temperature profile.
+With the simulation cell described above, we are able to equilibrate
+the system and impose an unphysical thermal flux between the liquid
+and the metal phase using the NIVS algorithm. By periodically applying
+the unphysical flux, we are able to obtain a temperature profile and
+its spatial derivatives. These quantities enable the evaluation of the
+interfacial thermal conductance of a surface. Figure \ref{gradT} is an
+example how those applied thermal fluxes can be used to obtain the 1st
+and 2nd derivatives of the temperature profile.
 
 \begin{figure}
 \includegraphics[width=\linewidth]{gradT}
@@ -212,42 +250,45 @@ obtain the 1st and 2nd derivatives of the temperature 
 \end{figure}
 
 \section{Computational Details}
-\subsection{System Geometry}
-In our simulations, Au is used to construct a metal slab with bare
-(111) surface perpendicular to the $z$-axis. Different slab thickness
-(layer numbers of Au) are simulated. This metal slab is first
-equilibrated under normal pressure (1 atm) and a desired
-temperature. After equilibration, butanethiol is used as the capping
-agent molecule to cover the bare Au (111) surfaces evenly. The sulfur
-atoms in the butanethiol molecules would occupy the three-fold sites
-of the surfaces, and the maximal butanethiol capacity on Au surface is
-$1/3$ of the total number of surface Au atoms[CITATION]. A series of
-different coverage surfaces is investigated in order to study the
-relation between coverage and conductance.
+\subsection{Simulation Protocol}
+The NIVS algorithm has been implemented in our MD simulation code,
+OpenMD\cite{Meineke:2005gd,openmd}, and was used for our
+simulations. Different slab thickness (layer numbers of Au) were
+simulated. Metal slabs were first equilibrated under atmospheric
+pressure (1 atm) and a desired temperature (e.g. 200K). After
+equilibration, butanethiol capping agents were placed at three-fold
+sites on the Au(111) surfaces. The maximum butanethiol capacity on Au
+surface is $1/3$ of the total number of surface Au
+atoms\cite{vlugt:cpc2007154}. A series of different coverages was
+investigated in order to study the relation between coverage and
+interfacial conductance.
 
-[COVERAGE DISCRIPTION] However, since the interactions between surface
-Au and butanethiol is non-bonded, the capping agent molecules are
-allowed to migrate to an empty neighbor three-fold site during a
-simulation. Therefore, the initial configuration would not severely
-affect the sampling of a variety of configurations of the same
-coverage, and the final conductance measurement would be an average
-effect of these configurations explored in the simulations. [MAY NEED FIGURES]
+The capping agent molecules were allowed to migrate during the
+simulations. They distributed themselves uniformly and sampled a
+number of three-fold sites throughout out study. Therefore, the
+initial configuration would not noticeably affect the sampling of a
+variety of configurations of the same coverage, and the final
+conductance measurement would be an average effect of these
+configurations explored in the simulations. [MAY NEED FIGURES]
 
-After the modified Au-butanethiol surface systems are equilibrated
-under canonical ensemble, Packmol\cite{packmol} is used to pack
-organic solvent molecules in the previously vacuum part of the
-simulation cells, which guarantees that short range repulsive
-interactions do not disrupt the simulations. Two solvents are
-investigated, one which has little vibrational overlap with the
-alkanethiol and plane-like shape (toluene), and one which has similar
-vibrational frequencies and chain-like shape ({\it n}-hexane). The
-spacing filled by solvent molecules, i.e. the gap between periodically
-repeated Au-butanethiol surfaces should be carefully chosen so that it
-would not be too short to affect the liquid phase structure, nor too
-long, leading to over cooling (freezing) or heating (boiling) when a
-thermal flux is applied. In our simulations, this spacing is usually
-$35 \sim 60$\AA.
+After the modified Au-butanethiol surface systems were equilibrated
+under canonical ensemble, organic solvent molecules were packed in the
+previously empty part of the simulation cells\cite{packmol}. Two
+solvents were investigated, one which has little vibrational overlap
+with the alkanethiol and a planar shape (toluene), and one which has
+similar vibrational frequencies and chain-like shape ({\it n}-hexane).
 
+The space filled by solvent molecules, i.e. the gap between
+periodically repeated Au-butanethiol surfaces should be carefully
+chosen. A very long length scale for the thermal gradient axis ($z$)
+may cause excessively hot or cold temperatures in the middle of the
+solvent region and lead to undesired phenomena such as solvent boiling
+or freezing when a thermal flux is applied. Conversely, too few
+solvent molecules would change the normal behavior of the liquid
+phase. Therefore, our $N_{solvent}$ values were chosen to ensure that
+these extreme cases did not happen to our simulations. And the
+corresponding spacing is usually $35 \sim 60$\AA.
+
 The initial configurations generated by Packmol are further
 equilibrated with the $x$ and $y$ dimensions fixed, only allowing
 length scale change in $z$ dimension. This is to ensure that the
@@ -282,11 +323,21 @@ quantum Sutton-Chen (QSC) formulation.\cite{PhysRevB.5
 same type of particles and between particles of different species.
 
 The Au-Au interactions in metal lattice slab is described by the
-quantum Sutton-Chen (QSC) formulation.\cite{PhysRevB.59.3527} The QSC
+quantum Sutton-Chen (QSC) formulation\cite{PhysRevB.59.3527}. The QSC
 potentials include zero-point quantum corrections and are
 reparametrized for accurate surface energies compared to the
 Sutton-Chen potentials\cite{Chen90}. 
 
+Figure \ref{demoMol} demonstrates how we name our pseudo-atoms of the
+organic solvent molecules in our simulations.
+
+\begin{figure}
+\includegraphics[width=\linewidth]{demoMol}
+\caption{Denomination of atoms or pseudo-atoms in our simulations: a)
+  UA-hexane; b) AA-hexane; c) UA-toluene; d) AA-toluene.}
+\label{demoMol}
+\end{figure}
+
 For both solvent molecules, straight chain {\it n}-hexane and aromatic
 toluene, United-Atom (UA) and All-Atom (AA) models are used
 respectively. The TraPPE-UA
@@ -300,73 +351,220 @@ simulations should be much lower than the liquid boili
 to a system, the temperature of ``hot'' area in the liquid phase would be
 significantly higher than the average, to prevent over heating and
 boiling of the liquid phase, the average temperature in our
-simulations should be much lower than the liquid boiling point. [NEED MORE DISCUSSION]
+simulations should be much lower than the liquid boiling point. [MORE DISCUSSION]
 For UA-toluene model, rigid body constraints are applied, so that the
-benzene ring and the methyl-C(aromatic) bond are kept rigid. This
-would save computational time.[MORE DETAILS NEEDED]
+benzene ring and the methyl-CRar bond are kept rigid. This would save
+computational time.[MORE DETAILS]
 
 Besides the TraPPE-UA models, AA models for both organic solvents are
-included in our studies as well. For hexane, the OPLS
-all-atom\cite{OPLSAA} force field is used. [MORE DETAILS] 
+included in our studies as well. For hexane, the OPLS-AA\cite{OPLSAA}
+force field is used. [MORE DETAILS]
 For toluene, the United Force Field developed by Rapp\'{e} {\it et
   al.}\cite{doi:10.1021/ja00051a040} is adopted.[MORE DETAILS]
 
 The capping agent in our simulations, the butanethiol molecules can
 either use UA or AA model. The TraPPE-UA force fields includes
-parameters for thiol molecules\cite{TraPPE-UA.thiols} and are used in
-our simulations corresponding to our TraPPE-UA models for solvent. 
- and All-Atom models [NEED CITATIONS]
-However, the model choice (UA or AA) of capping agent can be different
-from the solvent. Regardless of model choice, the force field
-parameters for interactions between capping agent and solvent can be
-derived using Lorentz-Berthelot Mixing Rule.
+parameters for thiol molecules\cite{TraPPE-UA.thiols} and are used for
+UA butanethiol model in our simulations. The OPLS-AA also provides
+parameters for alkyl thiols. However, alkyl thiols adsorbed on Au(111)
+surfaces do not have the hydrogen atom bonded to sulfur. To adapt this
+change and derive suitable parameters for butanethiol adsorbed on
+Au(111) surfaces, we adopt the S parameters from Luedtke and
+Landman\cite{landman:1998} and modify parameters for its neighbor C
+atom for charge balance in the molecule. Note that the model choice
+(UA or AA) of capping agent can be different from the
+solvent. Regardless of model choice, the force field parameters for
+interactions between capping agent and solvent can be derived using
+Lorentz-Berthelot Mixing Rule:[EQN'S]
 
+
 To describe the interactions between metal Au and non-metal capping
-agent and solvent, we refer to Vlugt\cite{vlugt:cpc2007154} and derive
-other interactions which are not yet finely parametrized. [can add
-hautman and klein's paper here and more discussion; need to put
-aromatic-metal interaction approximation here]\cite{doi:10.1021/jp034405s}
+agent and solvent particles, we refer to an adsorption study of alkyl
+thiols on gold surfaces by Vlugt {\it et
+  al.}\cite{vlugt:cpc2007154} They fitted an effective Lennard-Jones
+form of potential parameters for the interaction between Au and
+pseudo-atoms CH$_x$ and S based on a well-established and widely-used
+effective potential of Hautman and Klein\cite{hautman:4994} for the
+Au(111) surface. As our simulations require the gold lattice slab to
+be non-rigid so that it could accommodate kinetic energy for thermal
+transport study purpose, the pair-wise form of potentials is
+preferred.
 
-[TABULATED FORCE FIELD PARAMETERS NEEDED]
+Besides, the potentials developed from {\it ab initio} calculations by
+Leng {\it et al.}\cite{doi:10.1021/jp034405s} are adopted for the
+interactions between Au and aromatic C/H atoms in toluene.[MORE DETAILS]
 
+However, the Lennard-Jones parameters between Au and other types of
+particles in our simulations are not yet well-established. For these
+interactions, we attempt to derive their parameters using the Mixing
+Rule. To do this, the ``Metal-non-Metal'' (MnM) interaction parameters
+for Au is first extracted from the Au-CH$_x$ parameters by applying
+the Mixing Rule reversely. Table \ref{MnM} summarizes these ``MnM''
+parameters in our simulations.
 
-[SURFACE RECONSTRUCTION PREVENTS SIMULATION TEMP TO GO HIGHER]
+\begin{table*}
+  \begin{minipage}{\linewidth}
+    \begin{center}
+      \caption{Lennard-Jones parameters for Au-non-Metal
+        interactions in our simulations.}
+      
+      \begin{tabular}{ccc}
+        \hline\hline
+        Non-metal atom   & $\sigma$ & $\epsilon$ \\
+        (or pseudo-atom) & \AA      & kcal/mol  \\
+        \hline
+        S    & 2.40   & 8.465   \\
+        CH3  & 3.54   & 0.2146  \\
+        CH2  & 3.54   & 0.1749  \\
+        CT3  & 3.365  & 0.1373  \\
+        CT2  & 3.365  & 0.1373  \\
+        CTT  & 3.365  & 0.1373  \\
+        HC   & 2.865  & 0.09256 \\
+        CHar & 3.4625 & 0.1680  \\
+        CRar & 3.555  & 0.1604  \\
+        CA   & 3.173  & 0.0640  \\
+        HA   & 2.746  & 0.0414  \\
+        \hline\hline
+      \end{tabular}
+      \label{MnM}
+    \end{center}
+  \end{minipage}
+\end{table*}
 
 
-\section{Results}
-[REARRANGEMENT NEEDED]
-\subsection{Toluene Solvent}
+\section{Results and Discussions}
+[MAY HAVE A BRIEF SUMMARY]
+\subsection{How Simulation Parameters Affects $G$}
+[MAY NOT PUT AT FIRST]
+We have varied our protocol or other parameters of the simulations in
+order to investigate how these factors would affect the measurement of
+$G$'s. It turned out that while some of these parameters would not
+affect the results substantially, some other changes to the
+simulations would have a significant impact on the measurement
+results.
 
-The results (Table \ref{AuThiolToluene}) show a
-significant conductance enhancement compared to the gold/water
-interface without capping agent and agree with available experimental
-data. This indicates that the metal-metal potential, though not
-predicting an accurate bulk metal thermal conductivity, does not
-greatly interfere with the simulation of the thermal conductance
-behavior across a non-metal interface. The solvent model is not
-particularly volatile, so the simulation cell does not expand
-significantly under higher temperature. We did not observe a
-significant conductance decrease when the temperature was increased to
-300K. The results show that the two definitions used for $G$ yield
-comparable values, though $G^\prime$ tends to be smaller.
+In some of our simulations, we allowed $L_x$ and $L_y$ to change
+during equilibrating the liquid phase. Due to the stiffness of the Au
+slab, $L_x$ and $L_y$ would not change noticeably after
+equilibration. Although $L_z$ could fluctuate $\sim$1\% after a system
+is fully equilibrated in the NPT ensemble, this fluctuation, as well
+as those comparably smaller to $L_x$ and $L_y$, would not be magnified
+on the calculated $G$'s, as shown in Table \ref{AuThiolHexaneUA}. This
+insensivity to $L_i$ fluctuations allows reliable measurement of $G$'s
+without the necessity of extremely cautious equilibration process. 
 
+As stated in our computational details, the spacing filled with
+solvent molecules can be chosen within a range. This allows some
+change of solvent molecule numbers for the same Au-butanethiol
+surfaces. We did this study on our Au-butanethiol/hexane
+simulations. Nevertheless, the results obtained from systems of
+different $N_{hexane}$ did not indicate that the measurement of $G$ is
+susceptible to this parameter. For computational efficiency concern,
+smaller system size would be preferable, given that the liquid phase
+structure is not affected.
+
+Our NIVS algorithm allows change of unphysical thermal flux both in
+direction and in quantity. This feature extends our investigation of
+interfacial thermal conductance. However, the magnitude of this
+thermal flux is not arbitary if one aims to obtain a stable and
+reliable thermal gradient. A temperature profile would be
+substantially affected by noise when $|J_z|$ has a much too low
+magnitude; while an excessively large $|J_z|$ that overwhelms the
+conductance capacity of the interface would prevent a thermal gradient
+to reach a stablized steady state. NIVS has the advantage of allowing
+$J$ to vary in a wide range such that the optimal flux range for $G$
+measurement can generally be simulated by the algorithm. Within the
+optimal range, we were able to study how $G$ would change according to
+the thermal flux across the interface. For our simulations, we denote
+$J_z$ to be positive when the physical thermal flux is from the liquid
+to metal, and negative vice versa. The $G$'s measured under different
+$J_z$ is listed in Table \ref{AuThiolHexaneUA} and [REF]. These
+results do not suggest that $G$ is dependent on $J_z$ within this flux
+range. The linear response of flux to thermal gradient simplifies our
+investigations in that we can rely on $G$ measurement with only a
+couple $J_z$'s and do not need to test a large series of fluxes.
+
+%ADD MORE TO TABLE
 \begin{table*}
   \begin{minipage}{\linewidth}
     \begin{center}
       \caption{Computed interfacial thermal conductivity ($G$ and
-        $G^\prime$) values for the Au/butanethiol/toluene interface at
-        different temperatures using a range of energy fluxes.}
+        $G^\prime$) values for the 100\% covered Au-butanethiol/hexane
+        interfaces with UA model and different hexane molecule numbers
+        at different temperatures using a range of energy fluxes.}
       
+      \begin{tabular}{cccccccc}
+        \hline\hline
+        $\langle T\rangle$ & & $L_x$ & $L_y$ & $L_z$ & $J_z$ &
+        $G$ & $G^\prime$ \\
+        (K) & $N_{hexane}$ & \multicolumn{3}{c}{(\AA)} & (GW/m$^2$) &
+        \multicolumn{2}{c}{(MW/m$^2$/K)} \\
+        \hline
+        200 & 266 & 29.86 & 25.80 & 113.1 & -0.96 &
+        102()  & 80.0() \\
+            & 200 & 29.84 & 25.81 &  93.9 &  1.92 &
+        129()  & 87.3() \\
+            &     & 29.84 & 25.81 &  95.3 &  1.93 &
+        131()  & 77.5() \\
+            & 166 & 29.84 & 25.81 &  85.7 &  0.97 &
+        115()  & 69.3() \\
+            &     &       &       &       &  1.94 &
+        125()  & 87.1() \\
+        250 & 200 & 29.84 & 25.87 & 106.8 &  0.96 &
+        81.8() & 67.0() \\
+            & 166 & 29.87 & 25.84 &  94.8 &  0.98 &
+        79.0() & 62.9() \\
+            &     & 29.84 & 25.85 &  95.0 &  1.44 &
+        76.2() & 64.8() \\
+        \hline\hline
+      \end{tabular}
+      \label{AuThiolHexaneUA}
+    \end{center}
+  \end{minipage}
+\end{table*}
+
+Furthermore, we also attempted to increase system average temperatures
+to above 200K. These simulations are first equilibrated in the NPT
+ensemble under normal pressure. As stated above, the TraPPE-UA model
+for hexane tends to predict a lower boiling point. In our simulations,
+hexane had diffculty to remain in liquid phase when NPT equilibration
+temperature is higher than 250K. Additionally, the equilibrated liquid
+hexane density under 250K becomes lower than experimental value. This
+expanded liquid phase leads to lower contact between hexane and
+butanethiol as well.[MAY NEED FIGURE] And this reduced contact would
+probably be accountable for a lower interfacial thermal conductance,
+as shown in Table \ref{AuThiolHexaneUA}. 
+
+A similar study for TraPPE-UA toluene agrees with the above result as
+well. Having a higher boiling point, toluene tends to remain liquid in
+our simulations even equilibrated under 300K in NPT
+ensembles. Furthermore, the expansion of the toluene liquid phase is
+not as significant as that of the hexane. This prevents severe
+decrease of liquid-capping agent contact and the results (Table
+\ref{AuThiolToluene}) show only a slightly decreased interface
+conductance. Therefore, solvent-capping agent contact should play an
+important role in the thermal transport process across the interface
+in that higher degree of contact could yield increased conductance.
+
+[ADD Lxyz AND ERROR ESTIMATE TO TABLE]
+\begin{table*}
+  \begin{minipage}{\linewidth}
+    \begin{center}
+      \caption{Computed interfacial thermal conductivity ($G$ and
+        $G^\prime$) values for a 90\% coverage Au-butanethiol/toluene
+        interface at different temperatures using a range of energy
+        fluxes.}
+      
       \begin{tabular}{cccc}
         \hline\hline
         $\langle T\rangle$ & $J_z$ & $G$ & $G^\prime$ \\
         (K) & (GW/m$^2$) & \multicolumn{2}{c}{(MW/m$^2$/K)} \\
         \hline
-        200 & 1.86 & 180 & 135 \\
-            & 2.15 & 204 & 113 \\
-            & 3.93 & 175 & 114 \\
-        300 & 1.91 & 143 & 125 \\
-            & 4.19 & 134 & 113 \\
+        200 & -1.86 & 180() & 135() \\
+            &  2.15 & 204() & 113() \\
+            & -3.93 & 175() & 114() \\
+        300 & -1.91 & 143() & 125() \\
+            & -4.19 & 134() & 113() \\
         \hline\hline
       \end{tabular}
       \label{AuThiolToluene}
@@ -374,145 +572,252 @@ comparable values, though $G^\prime$ tends to be small
   \end{minipage}
 \end{table*}
 
-\subsection{Hexane Solvent}
+Besides lower interfacial thermal conductance, surfaces in relatively
+high temperatures are susceptible to reconstructions, when
+butanethiols have a full coverage on the Au(111) surface. These
+reconstructions include surface Au atoms migrated outward to the S
+atom layer, and butanethiol molecules embedded into the original
+surface Au layer. The driving force for this behavior is the strong
+Au-S interactions in our simulations. And these reconstructions lead
+to higher ratio of Au-S attraction and thus is energetically
+favorable. Furthermore, this phenomenon agrees with experimental
+results\cite{doi:10.1021/j100035a033,doi:10.1021/la026493y}. Vlugt
+{\it et al.} had kept their Au(111) slab rigid so that their
+simulations can reach 300K without surface reconstructions. Without
+this practice, simulating 100\% thiol covered interfaces under higher
+temperatures could hardly avoid surface reconstructions. However, our
+measurement is based on assuming homogeneity on $x$ and $y$ dimensions
+so that measurement of $T$ at particular $z$ would be an effective
+average of the particles of the same type. Since surface
+reconstructions could eliminate the original $x$ and $y$ dimensional
+homogeneity, measurement of $G$ is more difficult to conduct under
+higher temperatures. Therefore, most of our measurements are
+undertaken at $\langle T\rangle\sim$200K.
 
-Using the united-atom model, different coverages of capping agent,
-temperatures of simulations and numbers of solvent molecules were all
-investigated and Table \ref{AuThiolHexaneUA} shows the results of
-these computations. The number of hexane molecules in our simulations
-does not affect the calculations significantly. However, a very long
-length scale for the thermal gradient axis ($z$) may cause excessively
-hot or cold temperatures in the middle of the solvent region and lead
-to undesired phenomena such as solvent boiling or freezing, while too
-few solvent molecules would change the normal behavior of the liquid
-phase. Our $N_{hexane}$ values were chosen to ensure that these
-extreme cases did not happen to our simulations.
+However, when the surface is not completely covered by butanethiols,
+the simulated system is more resistent to the reconstruction
+above. Our Au-butanethiol/toluene system did not see this phenomena
+even at $\langle T\rangle\sim$300K. The Au(111) surfaces have a 90\% coverage of
+butanethiols and have empty three-fold sites. These empty sites could
+help prevent surface reconstruction in that they provide other means
+of capping agent relaxation. It is observed that butanethiols can
+migrate to their neighbor empty sites during a simulation. Therefore,
+we were able to obtain $G$'s for these interfaces even at a relatively
+high temperature without being affected by surface reconstructions.
 
-Table \ref{AuThiolHexaneUA} enables direct comparison between
-different coverages of capping agent, when other system parameters are
-held constant. With high coverage of butanethiol on the gold surface,
-the interfacial thermal conductance is enhanced
-significantly. Interestingly, a slightly lower butanethiol coverage
-leads to a moderately higher conductivity. This is probably due to
-more solvent/capping agent contact when butanethiol molecules are
-not densely packed, which enhances the interactions between the two
-phases and lowers the thermal transfer barrier of this interface.
-% [COMPARE TO AU/WATER IN PAPER]
+\subsection{Influence of Capping Agent Coverage on $G$}
+To investigate the influence of butanethiol coverage on interfacial
+thermal conductance, a series of different coverage Au-butanethiol
+surfaces is prepared and solvated with various organic
+molecules. These systems are then equilibrated and their interfacial
+thermal conductivity are measured with our NIVS algorithm. Table
+\ref{tlnUhxnUhxnD} lists these results for direct comparison between
+different coverages of butanethiol. To study the isotope effect in
+interfacial thermal conductance, deuterated UA-hexane is included as
+well.
 
-It is also noted that the overall simulation temperature is another
-factor that affects the interfacial thermal conductance. One
-possibility of this effect may be rooted in the decrease in density of
-the liquid phase. We observed that when the average temperature
-increases from 200K to 250K, the bulk hexane density becomes lower
-than experimental value, as the system is equilibrated under NPT
-ensemble. This leads to lower contact between solvent and capping
-agent, and thus lower conductivity.
+It turned out that with partial covered butanethiol on the Au(111)
+surface, the derivative definition for $G$ (Eq. \ref{derivativeG}) has
+difficulty to apply, due to the difficulty in locating the maximum of
+change of $\lambda$. Instead, the discrete definition
+(Eq. \ref{discreteG}) is easier to apply, as max($\Delta T$) can still
+be well-defined. Therefore, $G$'s (not $G^\prime$) are used for this
+section. 
 
-Conductivity values are more difficult to obtain under higher
-temperatures. This is because the Au surface tends to undergo
-reconstructions in relatively high temperatures. Surface Au atoms can
-migrate outward to reach higher Au-S contact; and capping agent
-molecules can be embedded into the surface Au layer due to the same
-driving force. This phenomenon agrees with experimental
-results\cite{doi:10.1021/j100035a033,doi:10.1021/la026493y}. A surface
-fully covered in capping agent is more susceptible to reconstruction,
-possibly because fully coverage prevents other means of capping agent
-relaxation, such as migration to an empty neighbor three-fold site.
+From Table \ref{tlnUhxnUhxnD}, one can see the significance of the
+presence of capping agents. Even when a fraction of the Au(111)
+surface sites are covered with butanethiols, the conductivity would
+see an enhancement by at least a factor of 3. This indicates the
+important role cappping agent is playing for thermal transport
+phenomena on metal/organic solvent surfaces.
 
-%MAY ADD MORE DATA TO TABLE
+Interestingly, as one could observe from our results, the maximum
+conductance enhancement (largest $G$) happens while the surfaces are
+about 75\% covered with butanethiols. This again indicates that
+solvent-capping agent contact has an important role of the thermal
+transport process. Slightly lower butanethiol coverage allows small
+gaps between butanethiols to form. And these gaps could be filled with
+solvent molecules, which acts like ``heat conductors'' on the
+surface. The higher degree of interaction between these solvent
+molecules and capping agents increases the enhancement effect and thus
+produces a higher $G$ than densely packed butanethiol arrays. However,
+once this maximum conductance enhancement is reached, $G$ decreases
+when butanethiol coverage continues to decrease. Each capping agent
+molecule reaches its maximum capacity for thermal
+conductance. Therefore, even higher solvent-capping agent contact
+would not offset this effect. Eventually, when butanethiol coverage
+continues to decrease, solvent-capping agent contact actually
+decreases with the disappearing of butanethiol molecules. In this
+case, $G$ decrease could not be offset but instead accelerated.
+
+A comparison of the results obtained from differenet organic solvents
+can also provide useful information of the interfacial thermal
+transport process. The deuterated hexane (UA) results do not appear to
+be much different from those of normal hexane (UA), given that
+butanethiol (UA) is non-deuterated for both solvents. These UA model
+studies, even though eliminating C-H vibration samplings, still have
+C-C vibrational frequencies different from each other. However, these
+differences in the infrared range do not seem to produce an observable
+difference for the results of $G$. [MAY NEED FIGURE]
+
+Furthermore, results for rigid body toluene solvent, as well as other
+UA-hexane solvents, are reasonable within the general experimental
+ranges[CITATIONS]. This suggests that explicit hydrogen might not be a
+required factor for modeling thermal transport phenomena of systems
+such as Au-thiol/organic solvent.
+
+However, results for Au-butanethiol/toluene do not show an identical
+trend with those for Au-butanethiol/hexane in that $G$'s remain at
+approximately the same magnitue when butanethiol coverage differs from
+25\% to 75\%. This might be rooted in the molecule shape difference
+for plane-like toluene and chain-like {\it n}-hexane. Due to this
+difference, toluene molecules have more difficulty in occupying
+relatively small gaps among capping agents when their coverage is not
+too low. Therefore, the solvent-capping agent contact may keep
+increasing until the capping agent coverage reaches a relatively low
+level. This becomes an offset for decreasing butanethiol molecules on
+its effect to the process of interfacial thermal transport. Thus, one
+can see a plateau of $G$ vs. butanethiol coverage in our results.
+
+[NEED ERROR ESTIMATE, MAY ALSO PUT J HERE]
 \begin{table*}
   \begin{minipage}{\linewidth}
     \begin{center}
-      \caption{Computed interfacial thermal conductivity ($G$ and
-        $G^\prime$) values for the Au/butanethiol/hexane interface
-        with united-atom model and different capping agent coverage
-        and solvent molecule numbers at different temperatures using a
-        range of energy fluxes.}
+      \caption{Computed interfacial thermal conductivity ($G$) values
+        for the Au-butanethiol/solvent interface with various UA
+        models and different capping agent coverages at $\langle
+        T\rangle\sim$200K using certain energy flux respectively.}
       
-      \begin{tabular}{cccccc}
+      \begin{tabular}{cccc}
         \hline\hline
-        Thiol & $\langle T\rangle$ & & $J_z$ & $G$ & $G^\prime$ \\
-        coverage (\%) & (K) & $N_{hexane}$ & (GW/m$^2$) &
-        \multicolumn{2}{c}{(MW/m$^2$/K)} \\
+        Thiol & \multicolumn{3}{c}{$G$(MW/m$^2$/K)} \\
+        coverage (\%) & hexane & hexane(D) & toluene \\
         \hline
-        0.0   & 200 & 200 & 0.96 & 43.3 & 42.7 \\
-              &     &     & 1.91 & 45.7 & 42.9 \\
-              &     & 166 & 0.96 & 43.1 & 53.4 \\
-        88.9  & 200 & 166 & 1.94 & 172  & 108  \\
-        100.0 & 250 & 200 & 0.96 & 81.8 & 67.0 \\
-              &     & 166 & 0.98 & 79.0 & 62.9 \\
-              &     &     & 1.44 & 76.2 & 64.8 \\
-              & 200 & 200 & 1.92 & 129  & 87.3 \\
-              &     &     & 1.93 & 131  & 77.5 \\
-              &     & 166 & 0.97 & 115  & 69.3 \\
-              &     &     & 1.94 & 125  & 87.1 \\
+        0.0   & 46.5() & 43.9() & 70.1() \\
+        25.0  & 151()  & 153()  & 249()  \\
+        50.0  & 172()  & 182()  & 214()  \\
+        75.0  & 242()  & 229()  & 244()  \\
+        88.9  & 178()  & -      & -      \\
+        100.0 & 137()  & 153()  & 187()  \\
         \hline\hline
       \end{tabular}
-      \label{AuThiolHexaneUA}
+      \label{tlnUhxnUhxnD}
     \end{center}
   \end{minipage}
 \end{table*}
 
-For the all-atom model, the liquid hexane phase was not stable under NPT
-conditions. Therefore, the simulation length scale parameters are
-adopted from previous equilibration results of the united-atom model
-at 200K. Table \ref{AuThiolHexaneAA} shows the results of these
-simulations. The conductivity values calculated with full capping
-agent coverage are substantially larger than observed in the
-united-atom model, and is even higher than predicted by
-experiments. It is possible that our parameters for metal-non-metal
-particle interactions lead to an overestimate of the interfacial
-thermal conductivity, although the active C-H vibrations in the
-all-atom model (which should not be appreciably populated at normal
-temperatures) could also account for this high conductivity. The major
-thermal transfer barrier of Au/butanethiol/hexane interface is between
-the liquid phase and the capping agent, so extra degrees of freedom
-such as the C-H vibrations could enhance heat exchange between these
-two phases and result in a much higher conductivity.
+\subsection{Influence of Chosen Molecule Model on $G$}
+[MAY COMBINE W MECHANISM STUDY]
 
+In addition to UA solvent/capping agent models, AA models are included
+in our simulations as well. Besides simulations of the same (UA or AA)
+model for solvent and capping agent, different models can be applied
+to different components. Furthermore, regardless of models chosen,
+either the solvent or the capping agent can be deuterated, similar to
+the previous section. Table \ref{modelTest} summarizes the results of
+these studies.
+
+[MORE DATA; ERROR ESTIMATE]
 \begin{table*}
   \begin{minipage}{\linewidth}
     \begin{center}
       
       \caption{Computed interfacial thermal conductivity ($G$ and
-        $G^\prime$) values for the Au/butanethiol/hexane interface
-        with all-atom model and different capping agent coverage at
-        200K using a range of energy fluxes.}
+        $G^\prime$) values for interfaces using various models for
+        solvent and capping agent (or without capping agent) at
+        $\langle T\rangle\sim$200K.}
       
-      \begin{tabular}{cccc}
+      \begin{tabular}{ccccc}
         \hline\hline
-        Thiol & $J_z$ & $G$ & $G^\prime$ \\
-        coverage (\%) & (GW/m$^2$) & \multicolumn{2}{c}{(MW/m$^2$/K)} \\
+        Butanethiol model & Solvent & $J_z$ & $G$ & $G^\prime$ \\
+        (or bare surface) & model & (GW/m$^2$) &
+        \multicolumn{2}{c}{(MW/m$^2$/K)} \\
         \hline
-        0.0   & 0.95 & 28.5 & 27.2 \\
-              & 1.88 & 30.3 & 28.9 \\
-        100.0 & 2.87 & 551  & 294  \\
-              & 3.81 & 494  & 193  \\
+        UA    & AA hexane  & 1.94 & 135()  & 129()  \\
+              &            & 2.86 & 126()  & 115()  \\
+              & AA toluene & 1.89 & 200()  & 149()  \\
+        AA    & UA hexane  & 1.94 & 116()  & 129()  \\
+              & AA hexane  & 3.76 & 451()  & 378()  \\
+              &            & 4.71 & 432()  & 334()  \\
+              & AA toluene & 3.79 & 487()  & 290()  \\
+        AA(D) & UA hexane  & 1.94 & 158()  & 172()  \\
+        bare  & AA hexane  & 0.96 & 31.0() & 29.4() \\
         \hline\hline
       \end{tabular}
-      \label{AuThiolHexaneAA}
+      \label{modelTest}
     \end{center}
   \end{minipage}
 \end{table*}
 
-%subsubsection{Vibrational spectrum study on conductance mechanism}
+To facilitate direct comparison, the same system with differnt models
+for different components uses the same length scale for their
+simulation cells. Without the presence of capping agent, using
+different models for hexane yields similar results for both $G$ and
+$G^\prime$, and these two definitions agree with eath other very
+well. This indicates very weak interaction between the metal and the
+solvent, and is a typical case for acoustic impedance mismatch between
+these two phases.
+
+As for Au(111) surfaces completely covered by butanethiols, the choice
+of models for capping agent and solvent could impact the measurement
+of $G$ and $G^\prime$ quite significantly. For Au-butanethiol/hexane
+interfaces, using AA model for both butanethiol and hexane yields
+substantially higher conductivity values than using UA model for at
+least one component of the solvent and capping agent, which exceeds
+the upper bond of experimental value range. This is probably due to
+the classically treated C-H vibrations in the AA model, which should
+not be appreciably populated at normal temperatures. In comparison,
+once either the hexanes or the butanethiols are deuterated, one can
+see a significantly lower $G$ and $G^\prime$. In either of these
+cases, the C-H(D) vibrational overlap between the solvent and the
+capping agent is removed. [MAY NEED FIGURE] Conclusively, the
+improperly treated C-H vibration in the AA model produced
+over-predicted results accordingly. Compared to the AA model, the UA
+model yields more reasonable results with higher computational
+efficiency.
+
+However, for Au-butanethiol/toluene interfaces, having the AA
+butanethiol deuterated did not yield a significant change in the
+measurement results.
+. , so extra degrees of freedom
+such as the C-H vibrations could enhance heat exchange between these
+two phases and result in a much higher conductivity.
+
+
+Although the QSC model for Au is known to predict an overly low value
+for bulk metal gold conductivity[CITE NIVSRNEMD], our computational
+results for $G$ and $G^\prime$ do not seem to be affected by this
+drawback of the model for metal. Instead, the modeling of interfacial
+thermal transport behavior relies mainly on an accurate description of
+the interactions between components occupying the interfaces.
+
+\subsection{Mechanism of Interfacial Thermal Conductance Enhancement
+  by Capping Agent}
+%OR\subsection{Vibrational spectrum study on conductance mechanism}
+
+[MAY INTRODUCE PROTOCOL IN METHODOLOGY/COMPUTATIONAL DETAIL, EQN'S]
+
 To investigate the mechanism of this interfacial thermal conductance,
 the vibrational spectra of various gold systems were obtained and are
 shown as in the upper panel of Fig. \ref{vibration}. To obtain these
 spectra, one first runs a simulation in the NVE ensemble and collects
 snapshots of configurations; these configurations are used to compute
 the velocity auto-correlation functions, which is used to construct a
-power spectrum via a Fourier transform. The gold surfaces covered by
-butanethiol molecules exhibit an additional peak observed at a
-frequency of $\sim$170cm$^{-1}$, which is attributed to the vibration
-of the S-Au bond. This vibration enables efficient thermal transport
-from surface Au atoms to the capping agents. Simultaneously, as shown
-in the lower panel of Fig. \ref{vibration}, the large overlap of the
-vibration spectra of butanethiol and hexane in the all-atom model,
-including the C-H vibration, also suggests high thermal exchange
-efficiency. The combination of these two effects produces the drastic
-interfacial thermal conductance enhancement in the all-atom model.
+power spectrum via a Fourier transform.
 
+ The gold surfaces covered by
+butanethiol molecules, compared to bare gold surfaces, exhibit an
+additional peak observed at a frequency of $\sim$170cm$^{-1}$, which
+is attributed to the vibration of the S-Au bond. This vibration
+enables efficient thermal transport from surface Au atoms to the
+capping agents. Simultaneously, as shown in the lower panel of
+Fig. \ref{vibration}, the large overlap of the vibration spectra of
+butanethiol and hexane in the all-atom model, including the C-H
+vibration, also suggests high thermal exchange efficiency. The
+combination of these two effects produces the drastic interfacial
+thermal conductance enhancement in the all-atom model.
+
+[MAY NEED TO CONVERT TO JPEG]
 \begin{figure}
 \includegraphics[width=\linewidth]{vibration}
 \caption{Vibrational spectra obtained for gold in different
@@ -520,14 +825,48 @@ interfacial thermal conductance enhancement in the all
   all-atom model (lower panel).}
 \label{vibration}
 \end{figure}
-% 600dpi, letter size. too large?
 
+[COMPARISON OF TWO G'S; AU SLAB WIDTHS; ETC]
+% The results show that the two definitions used for $G$ yield
+% comparable values, though $G^\prime$ tends to be smaller.
 
+\section{Conclusions}
+The NIVS algorithm we developed has been applied to simulations of
+Au-butanethiol surfaces with organic solvents. This algorithm allows
+effective unphysical thermal flux transferred between the metal and
+the liquid phase. With the flux applied, we were able to measure the
+corresponding thermal gradient and to obtain interfacial thermal
+conductivities. Our simulations have seen significant conductance
+enhancement with the presence of capping agent, compared to the bare
+gold/liquid interfaces. The acoustic impedance mismatch between the
+metal and the liquid phase is effectively eliminated by proper capping
+agent. Furthermore, the coverage precentage of the capping agent plays
+an important role in the interfacial thermal transport process.
+
+Our measurement results, particularly of the UA models, agree with
+available experimental data. This indicates that our force field
+parameters have a nice description of the interactions between the
+particles at the interfaces. AA models tend to overestimate the
+interfacial thermal conductance in that the classically treated C-H
+vibration would be overly sampled. Compared to the AA models, the UA
+models have higher computational efficiency with satisfactory
+accuracy, and thus are preferable in interfacial thermal transport
+modelings.
+
+Vlugt {\it et al.} has investigated the surface thiol structures for
+nanocrystal gold and pointed out that they differs from those of the
+Au(111) surface\cite{vlugt:cpc2007154}. This difference might lead to
+change of interfacial thermal transport behavior as well. To
+investigate this problem, an effective means to introduce thermal flux
+and measure the corresponding thermal gradient is desirable for
+simulating structures with spherical symmetry.
+
+
 \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
+Dame. \newpage
 
 \bibliography{interfacial}
 
