--- interfacial/interfacial.tex	2011/02/05 00:02:11	3721
+++ interfacial/interfacial.tex	2011/07/22 03:35:19	3746
@@ -22,9 +22,9 @@
 \setlength{\abovecaptionskip}{20 pt}
 \setlength{\belowcaptionskip}{30 pt}
 
-%\renewcommand\citemid{\ } % no comma in optional referenc note
-\bibpunct{[}{]}{,}{s}{}{;}
-\bibliographystyle{aip}
+%\renewcommand\citemid{\ } % no comma in optional reference note
+\bibpunct{[}{]}{,}{n}{}{;}
+\bibliographystyle{achemso}
 
 \begin{document}
 
@@ -44,7 +44,24 @@ The abstract version 2
 \begin{doublespace}
 
 \begin{abstract}
-The abstract version 2
+
+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}
 
 \newpage
@@ -56,96 +73,847 @@ The abstract version 2
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
 \section{Introduction}
+Interfacial thermal conductance is extensively studied both
+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.
 
-The intro.
+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 alkylthiol 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 comparatively 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 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}
-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 interfacial 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.
 
-(wondering how much detail of algorithm should be put here...)
+\subsection{Defining Interfacial Thermal Conductivity $G$}
+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_\mathrm{hot}\rangle -
+    \langle T_\mathrm{cold}\rangle \right)}
+\label{lowG}
+\end{equation}
+where ${E_{total}}$ is the imposed non-physical kinetic energy
+transfer and ${\langle T_\mathrm{hot}\rangle}$ and ${\langle
+  T_\mathrm{cold}\rangle}$ are the average observed temperature of the
+two separated phases.
 
-\subsection{Simulation Parameters}
-Our simulation systems consists of metal gold lattice slab solvated by
-organic solvents. In order to study the role of capping agents in
-interfacial thermal conductance, butanethiol is chosen to cover gold
-surfaces in comparison to no capping agent present.
+When the interfacial conductance is {\it not} small, there are two
+ways to define $G$.
 
-The Au-Au interactions in metal lattice slab is described by the
-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}. 
+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
+(Figure \ref{demoPic}):
+\begin{equation}
+G=\frac{J}{\Delta T}
+\label{discreteG}
+\end{equation}
 
-Straight chain {\it n}-hexane and aromatic toluene are respectively
-used as solvents. For hexane, both United-Atom\cite{TraPPE-UA.alkanes}
-and All-Atom\cite{OPLSAA} force fields are used for comparison; for
-toluene, United-Atom\cite{TraPPE-UA.alkylbenzenes} force fields are
-used with rigid body constraints applied. (maybe needs more details
-about rigid body)
+\begin{figure}
+\includegraphics[width=\linewidth]{method}
+\caption{Interfacial conductance can be calculated by applying an
+  (unphysical) kinetic energy flux between two slabs, one located
+  within the metal and another on the edge of the periodic box.  The
+  system responds by forming a thermal response or a gradient.  In
+  bulk liquids, this gradient typically has a single slope, but in
+  interfacial systems, there are distinct thermal conductivity
+  domains.  The interfacial conductance, $G$ is found by measuring the
+  temperature gap at the Gibbs dividing surface, or by using second
+  derivatives of the thermal profile.}
+\label{demoPic}
+\end{figure}
 
-Buatnethiol molecules are used as capping agent for some of our
-simulations. United-Atom\cite{TraPPE-UA.thiols} and All-Atom models
-are respectively used corresponding to the force field type of
-solvent.
+The other approach is to assume a continuous temperature profile along
+the thermal gradient axis (e.g. $z$) and define $G$ at the point where
+the magnitude of thermal conductivity $\lambda$ change reach its
+maximum, given that $\lambda$ is well-defined throughout the space:
+\begin{equation}
+G^\prime = \Big|\frac{\partial\lambda}{\partial z}\Big|
+         = \Big|\frac{\partial}{\partial z}\left(-J_z\Big/
+           \left(\frac{\partial T}{\partial z}\right)\right)\Big|
+         = |J_z|\Big|\frac{\partial^2 T}{\partial z^2}\Big|
+         \Big/\left(\frac{\partial T}{\partial z}\right)^2
+\label{derivativeG}
+\end{equation}
 
-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 parametrized in their work. (can add
-hautman and klein's paper here and more discussion; need to put
-aromatic-metal interaction approximation here)
+With the temperature profile obtained from simulations, one is able to
+approximate the first and second derivatives of $T$ with finite
+difference methods and thus calculate $G^\prime$.
 
-\section{Computational Details}
-Our simulation systems consists of a lattice Au slab with the (111)
-surface perpendicular to the $z$-axis, and a solvent layer between the
-periodic Au slabs along the $z$-axis. To set up the interfacial
-system, the Au slab is first equilibrated without solvent under room
-pressure and a desired temperature. After the metal slab is
-equilibrated, United-Atom or All-Atom butanethiols are replicated on
-the Au surface, each occupying the (??) among three Au atoms, and is
-equilibrated under NVT ensemble. According to (CITATION), the maximal
-thiol capacity on Au surface is $1/3$ of the total number of surface
-Au atoms. 
+In what follows, both definitions have been used for calculation and
+are compared in the results.
+
+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 without capping agents on the
+surfaces, the metal slab is solvated with simple organic solvents, as
+illustrated in Figure \ref{gradT}.
 
-\cite{packmol}
+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}
+\caption{A sample of Au-butanethiol/hexane interfacial system and the
+  temperature profile after a kinetic energy flux is imposed to
+  it. The 1st and 2nd derivatives of the temperature profile can be
+  obtained with finite difference approximation (lower panel).}
+\label{gradT}
+\end{figure}
+
+\section{Computational Details}
+\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 metal 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}[CITE CHEM REV].
+ A series of different coverages was
+investigated in order to study the relation between coverage and
+interfacial conductance.
+
+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 SNAPSHOTS]
+
+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 75$\AA.
+
+The initial configurations generated are further equilibrated with the
+$x$ and $y$ dimensions fixed, only allowing length scale change in $z$
+dimension. This is to ensure that the equilibration of liquid phase
+does not affect the metal crystal structure in $x$ and $y$ dimensions.
+To investigate this effect, comparisons were made with simulations
+allowed to change $L_x$ and $L_y$ during NPT equilibration, and the
+results are shown in later sections. After ensuring the liquid phase
+reaches equilibrium at atmospheric pressure (1 atm), further
+equilibration are followed under NVT and then NVE ensembles.
+
+After the systems reach equilibrium, NIVS is implemented to impose a
+periodic unphysical thermal flux between the metal and the liquid
+phase. Most of our simulations are under an average temperature of
+$\sim$200K. Therefore, this flux usually comes from the metal to the
+liquid so that the liquid has a higher temperature and would not
+freeze due to excessively low temperature. This induced temperature
+gradient is stablized and the simulation cell is devided evenly into
+N slabs along the $z$-axis and the temperatures of each slab are
+recorded. When the slab width $d$ of each slab is the same, the
+derivatives of $T$ with respect to slab number $n$ can be directly
+used for $G^\prime$ calculations:
+\begin{equation}
+G^\prime = |J_z|\Big|\frac{\partial^2 T}{\partial z^2}\Big|
+         \Big/\left(\frac{\partial T}{\partial z}\right)^2
+         = |J_z|\Big|\frac{1}{d^2}\frac{\partial^2 T}{\partial n^2}\Big|
+         \Big/\left(\frac{1}{d}\frac{\partial T}{\partial n}\right)^2
+         = |J_z|\Big|\frac{\partial^2 T}{\partial n^2}\Big|
+         \Big/\left(\frac{\partial T}{\partial n}\right)^2
+\label{derivativeG2}
+\end{equation}
+
+\subsection{Force Field Parameters}
+Our simulations include various components. Figure \ref{demoMol}
+demonstrates the sites defined for both United-Atom and All-Atom
+models of the organic solvent and capping agent molecules in our
+simulations. Force field parameter descriptions are needed for
+interactions both between the same type of particles and between
+particles of different species.
+
+\begin{figure}
+\includegraphics[width=\linewidth]{structures}
+\caption{Structures of the capping agent and solvents utilized in
+  these simulations. The chemically-distinct sites (a-e) are expanded
+  in terms of constituent atoms for both United Atom (UA) and All Atom
+  (AA) force fields.  Most parameters are from
+  Refs. \protect\cite{TraPPE-UA.alkanes,TraPPE-UA.alkylbenzenes} (UA) and
+  \protect\cite{OPLSAA} (AA).  Cross-interactions with the Au atoms are given
+  in Table \ref{MnM}.}
+\label{demoMol}
+\end{figure}
+
+The Au-Au interactions in metal lattice slab is described by the
+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}. 
+
+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
+parameters\cite{TraPPE-UA.alkanes,TraPPE-UA.alkylbenzenes} are used
+for our UA solvent molecules. In these models, sites are located at
+the carbon centers for alkyl groups. Bonding interactions, including
+bond stretches and bends and torsions, were used for intra-molecular
+sites not separated by more than 3 bonds. Otherwise, for non-bonded
+interactions, Lennard-Jones potentials are used. [MORE CITATION?]
+
+By eliminating explicit hydrogen atoms, these models are simple and
+computationally efficient, while maintains good accuracy. However, the
+TraPPE-UA for alkanes is known to predict a lower boiling point than
+experimental values. Considering that after an unphysical thermal flux
+is applied 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.
+
+For UA-toluene model, the non-bonded potentials between
+inter-molecular sites have a similar Lennard-Jones formulation. For
+intra-molecular interactions, considering the stiffness of the benzene
+ring, rigid body constraints are applied for further computational
+efficiency. All bonds in the benzene ring and between the ring and the
+methyl group remain rigid during the progress of simulations.
+
+Besides the TraPPE-UA models, AA models for both organic solvents are
+included in our studies as well. For hexane, the OPLS-AA\cite{OPLSAA}
+force field is used. Additional explicit hydrogen sites were
+included. Besides bonding and non-bonded site-site interactions,
+partial charges and the electrostatic interactions were added to each
+CT and HC site. For toluene, the United Force Field developed by
+Rapp\'{e} {\it et al.}\cite{doi:10.1021/ja00051a040} is
+adopted. Without the rigid body constraints, bonding interactions were
+included. For the aromatic ring, improper torsions (inversions) were
+added as an extra potential for maintaining the planar shape.
+[MORE CITATION?]
+
+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 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:
+\begin{eqnarray}
+\sigma_{ij} & = & \frac{1}{2} \left(\sigma_{ii} + \sigma_{jj}\right) \\
+\epsilon_{ij} & = & \sqrt{\epsilon_{ii}\epsilon_{jj}}
+\end{eqnarray}
+
+To describe the interactions between metal Au and non-metal capping
+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.
+
+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. A set of
+pseudo Lennard-Jones parameters were provided for Au in their force
+fields. By using the Mixing Rule, this can be used to derive pair-wise
+potentials for non-bonded interactions between Au and non-metal sites.
+
+However, the Lennard-Jones parameters between Au and other types of
+particles, such as All-Atom normal alkanes in our simulations are not
+yet well-established. For these interactions, we attempt to derive
+their parameters using the Mixing Rule. To do this, Au pseudo
+Lennard-Jones parameters for ``Metal-non-Metal'' (MnM) interactions
+were first extracted from the Au-CH$_x$ parameters by applying the
+Mixing Rule reversely. Table \ref{MnM} summarizes these ``MnM''
+parameters in our simulations.
+
+\begin{table*}
+  \begin{minipage}{\linewidth}
+    \begin{center}
+      \caption{Non-bonded interaction parameters (including cross
+        interactions with Au atoms) for both force fields used in this
+        work.}      
+      \begin{tabular}{lllllll}
+        \hline\hline
+        & Site  & $\sigma_{ii}$ & $\epsilon_{ii}$ & $q_i$ &
+        $\sigma_{Au-i}$ & $\epsilon_{Au-i}$ \\
+        & & (\AA) & (kcal/mol) & ($e$) & (\AA) & (kcal/mol) \\
+        \hline
+        United Atom (UA)
+        &CH3  & 3.75  & 0.1947  & -      & 3.54   & 0.2146  \\
+        &CH2  & 3.95  & 0.0914  & -      & 3.54   & 0.1749  \\
+        &CHar & 3.695 & 0.1003  & -      & 3.4625 & 0.1680  \\
+        &CRar & 3.88  & 0.04173 & -      & 3.555  & 0.1604  \\
+        \hline
+        All Atom (AA) 
+        &CT3  & 3.50  & 0.066   & -0.18  & 3.365  & 0.1373  \\
+        &CT2  & 3.50  & 0.066   & -0.12  & 3.365  & 0.1373  \\
+        &CTT  & 3.50  & 0.066   & -0.065 & 3.365  & 0.1373  \\
+        &HC   & 2.50  & 0.030   &  0.06  & 2.865  & 0.09256 \\
+        &CA   & 3.55  & 0.070   & -0.115 & 3.173  & 0.0640  \\
+        &HA   & 2.42  & 0.030   &  0.115 & 2.746  & 0.0414  \\
+        \hline
+        Both UA and AA 
+        & S   & 4.45  & 0.25    & -      & 2.40   & 8.465   \\
+        \hline\hline
+      \end{tabular}
+      \label{MnM}
+    \end{center}
+  \end{minipage}
+\end{table*}
 
+\subsection{Vibrational Spectrum}
+
+[MAY ADD EQN'S]
+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.
+
+\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.
+
+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
+crystalline Au structure, $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 of $L_x$ and $L_y$ (which is significantly smaller),
+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{AuThiolToluene}. 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.
+
+\begin{table*}
+  \begin{minipage}{\linewidth}
+    \begin{center}
+      \caption{Computed interfacial thermal conductivity ($G$ and
+        $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. Error estimates indicated in parenthesis.}
+      
+      \begin{tabular}{ccccccc}
+        \hline\hline
+        $\langle T\rangle$ & $N_{hexane}$ & Fixed & $\rho_{hexane}$ &
+        $J_z$ & $G$ & $G^\prime$ \\
+        (K) & & $L_x$ \& $L_y$? & (g/cm$^3$) & (GW/m$^2$) &
+        \multicolumn{2}{c}{(MW/m$^2$/K)} \\
+        \hline
+        200 & 266 & No  & 0.672 & -0.96 & 102(3)    & 80.0(0.8) \\
+            & 200 & Yes & 0.694 &  1.92 & 129(11)   & 87.3(0.3) \\
+            &     & Yes & 0.672 &  1.93 & 131(16)   & 78(13)    \\
+            &     & No  & 0.688 &  0.96 & 125(16)   & 90.2(15)  \\
+            &     &     &       &  1.91 & 139(10)   & 101(10)   \\
+            &     &     &       &  2.83 & 141(6)    & 89.9(9.8) \\
+            & 166 & Yes & 0.679 &  0.97 & 115(19)   & 69(18)    \\
+            &     &     &       &  1.94 & 125(9)    & 87.1(0.2) \\
+            &     & No  & 0.681 &  0.97 & 141(30)   & 78(22)    \\
+            &     &     &       &  1.92 & 138(4)    & 98.9(9.5) \\
+        \hline
+        250 & 200 & No  & 0.560 &  0.96 & 75(10)    & 61.8(7.3) \\
+            &     &     &       & -0.95 & 49.4(0.3) & 45.7(2.1) \\
+            & 166 & Yes & 0.570 &  0.98 & 79.0(3.5) & 62.9(3.0) \\
+            &     & No  & 0.569 &  0.97 & 80.3(0.6) & 67(11)    \\
+            &     &     &       &  1.44 & 76.2(5.0) & 64.8(3.8) \\
+            &     &     &       & -0.95 & 56.4(2.5) & 54.4(1.1) \\
+            &     &     &       & -1.85 & 47.8(1.1) & 53.5(1.5) \\
+        \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 SLAB DENSITY 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.
+
+\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. Error estimates indicated in parenthesis.}
+      
+      \begin{tabular}{ccccc}
+        \hline\hline
+        $\langle T\rangle$ & $\rho_{toluene}$ & $J_z$ & $G$ & $G^\prime$ \\
+        (K) & (g/cm$^3$) & (GW/m$^2$) & \multicolumn{2}{c}{(MW/m$^2$/K)} \\
+        \hline
+        200 & 0.933 &  2.15 & 204(12) & 113(12) \\
+            &       & -1.86 & 180(3)  & 135(21) \\
+            &       & -3.93 & 176(5)  & 113(12) \\
+        \hline
+        300 & 0.855 & -1.91 & 143(5)  & 125(2)  \\
+            &       & -4.19 & 135(9)  & 113(12) \\
+        \hline\hline
+      \end{tabular}
+      \label{AuThiolToluene}
+    \end{center}
+  \end{minipage}
+\end{table*}
+
+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.
+
+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 had the Au(111) surfaces 90\%
+covered by butanethiols, but did not see this above phenomena even at
+$\langle T\rangle\sim$300K. The empty three-fold sites not occupied by
+capping agents 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.
+
+\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. Figure
+\ref{coverage} demonstrates the trend of conductance change with
+respect to different coverages of butanethiol. To study the isotope
+effect in interfacial thermal conductance, deuterated UA-hexane is
+included as well.
+
+It turned out that with partial covered butanethiol on the Au(111)
+surface, the derivative definition for $G^\prime$
+(Eq. \ref{derivativeG}) was difficult 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 the Gibbs
+deviding surface can still be well-defined. Therefore, $G$ (not
+$G^\prime$) was used for this section.
+
+From Figure \ref{coverage}, 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.
+
+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. [NEED
+SNAPSHOT SHOWING THE PHENOMENA / SLAB DENSITY ANALYSIS]
+
+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 SPECTRA 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$ remains at
+approximately the same magnitue when butanethiol coverage differs from
+25\% to 75\%. This might be rooted in the molecule shape difference
+for planar 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.
+
+\begin{figure}
+\includegraphics[width=\linewidth]{coverage}
+\caption{Comparison of 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.}
+\label{coverage}
+\end{figure}
+
+\subsection{Influence of Chosen Molecule Model on $G$}
+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.
+
+\begin{table*}
+  \begin{minipage}{\linewidth}
+    \begin{center}
+      
+      \caption{Computed interfacial thermal conductivity ($G$ and
+        $G^\prime$) values for interfaces using various models for
+        solvent and capping agent (or without capping agent) at
+        $\langle T\rangle\sim$200K. (D stands for deuterated solvent
+        or capping agent molecules; ``Avg.'' denotes results that are
+        averages of simulations under different $J_z$'s. Error
+        estimates indicated in parenthesis.)}
+      
+      \begin{tabular}{llccc}
+        \hline\hline
+        Butanethiol model & Solvent & $J_z$ & $G$ & $G^\prime$ \\
+        (or bare surface) & model & (GW/m$^2$) &
+        \multicolumn{2}{c}{(MW/m$^2$/K)} \\
+        \hline
+        UA    & UA hexane    & Avg. & 131(9)    & 87(10)    \\
+              & UA hexane(D) & 1.95 & 153(5)    & 136(13)   \\
+              & AA hexane    & Avg. & 131(6)    & 122(10)   \\
+              & UA toluene   & 1.96 & 187(16)   & 151(11)   \\
+              & AA toluene   & 1.89 & 200(36)   & 149(53)   \\
+        \hline
+        AA    & UA hexane    & 1.94 & 116(9)    & 129(8)    \\
+              & AA hexane    & Avg. & 442(14)   & 356(31)   \\
+              & AA hexane(D) & 1.93 & 222(12)   & 234(54)   \\
+              & UA toluene   & 1.98 & 125(25)   & 97(60)    \\
+              & AA toluene   & 3.79 & 487(56)   & 290(42)   \\
+        \hline
+        AA(D) & UA hexane    & 1.94 & 158(25)   & 172(4)    \\
+              & AA hexane    & 1.92 & 243(29)   & 191(11)   \\
+              & AA toluene   & 1.93 & 364(36)   & 322(67)   \\
+        \hline
+        bare  & UA hexane    & Avg. & 46.5(3.2) & 49.4(4.5) \\
+              & UA hexane(D) & 0.98 & 43.9(4.6) & 43.0(2.0) \\
+              & AA hexane    & 0.96 & 31.0(1.4) & 29.4(1.3) \\
+              & UA toluene   & 1.99 & 70.1(1.3) & 65.8(0.5) \\
+        \hline\hline
+      \end{tabular}
+      \label{modelTest}
+    \end{center}
+  \end{minipage}
+\end{table*}
+
+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 general range of experimental measurement results. 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 SPECTRA 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. Compared to the C-H vibrational overlap between
+hexane and butanethiol, both of which have alkyl chains, that overlap
+between toluene and butanethiol is not so significant and thus does
+not have as much contribution to the ``Intramolecular Vibration
+Redistribution''[CITE HASE]. Conversely, extra degrees of freedom such
+as the C-H vibrations could yield higher heat exchange rate 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{kuang:164101}, our computational
+results for $G$ and $G^\prime$ do not seem to be affected by this
+drawback of the model for metal. Instead, our results suggest that the
+modeling of interfacial thermal transport behavior relies mainly on
+the accuracy of the interaction descriptions between components
+occupying the interfaces.
+
+\subsection{Role of Capping Agent in Interfacial Thermal Conductance}
+To investigate the mechanism of this interfacial thermal conductance,
+the vibrational spectra of various gold systems were obtained and are
+shown as in Fig. \ref{vibration}. 
+[MAY RELATE TO HASE'S]
+The gold surfaces covered by butanethiol molecules, compared to bare
+gold surfaces, exhibit an additional peak observed at the frequency of
+$\sim$170cm$^{-1}$, which is attributed to the S-Au bonding
+vibration. This vibration enables efficient thermal transport from
+surface Au layer to the capping agents.
+[MAY PUT IN OTHER SECTION] 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.
+
+[NEED SEPARATE FIGURE. MAY NEED TO CONVERT TO JPEG]
+\begin{figure}
+\includegraphics[width=\linewidth]{vibration}
+\caption{Vibrational spectra obtained for gold in different
+  environments.}
+\label{vibration}
+\end{figure}
+
+[MAY ADD COMPARISON OF G AND G', 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}
 
