ViewVC Help

View File | Revision Log | Show Annotations | View Changeset | Root Listing

root/group/stokes/stokes.tex

Comparing stokes/stokes.tex (file contents):
Revision 3770 by skuang, Fri Dec 2 20:14:03 2011 UTC vs.
Revision 3779 by skuang, Sat Dec 10 23:46:08 2011 UTC

#	Line 43 | Line 43 | Notre Dame, Indiana 46556}
43		\begin{doublespace}
44
45		\begin{abstract}
46	<	REPLACE ABSTRACT HERE
47	<	With the Non-Isotropic Velocity Scaling (NIVS) approach to Reverse
48	<	Non-Equilibrium Molecular Dynamics (RNEMD) it is possible to impose
49	<	an unphysical thermal flux between different regions of
50	<	inhomogeneous systems such as solid / liquid interfaces.  We have
51	<	applied NIVS to compute the interfacial thermal conductance at a
52	<	metal / organic solvent interface that has been chemically capped by
53	<	butanethiol molecules.  Our calculations suggest that coupling
54	<	between the metal and liquid phases is enhanced by the capping
55	<	agents, leading to a greatly enhanced conductivity at the interface.
56	<	Specifically, the chemical bond between the metal and the capping
57	<	agent introduces a vibrational overlap that is not present without
58	<	the capping agent, and the overlap between the vibrational spectra
59	<	(metal to cap, cap to solvent) provides a mechanism for rapid
60	<	thermal transport across the interface. Our calculations also
61	<	suggest that this is a non-monotonic function of the fractional
62	<	coverage of the surface, as moderate coverages allow diffusive heat
63	<	transport of solvent molecules that have been in close contact with
64	<	the capping agent.
46	>	We present a new method for introducing stable nonequilibrium
47	>	velocity and temperature gradients in molecular dynamics simulations
48	>	of heterogeneous systems. This method conserves the linear momentum
49	>	and total energy of the system and improves previous Reverse
50	>	Non-Equilibrium Molecular Dynamics (RNEMD) methods and maintains
51	>	thermal velocity distributions. It also avoid thermal anisotropy
52	>	occured in NIVS simulations by using isotropic velocity scaling on
53	>	the molecules in specific regions of a system. To test the method,
54	>	we have computed the thermal conductivity and shear viscosity of
55	>	model liquid systems as well as the interfacial frictions of a
56	>	series of  metal/liquid interfaces.
57
58		\end{abstract}
59
#	Line 74 | Line 66 | Notre Dame, Indiana 46556}
66		%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
67
68		\section{Introduction}
69	<	[DO THIS LATER]
69	>	[REFINE LATER, ADD MORE REF.S]
70	>	Imposed-flux methods in Molecular Dynamics (MD)
71	>	simulations\cite{MullerPlathe:1997xw} can establish steady state
72	>	systems with a set applied flux vs a corresponding gradient that can
73	>	be measured. These methods does not need many trajectories to provide
74	>	information of transport properties of a given system. Thus, they are
75	>	utilized in computing thermal and mechanical transfer of homogeneous
76	>	or bulk systems as well as heterogeneous systems such as liquid-solid
77	>	interfaces.\cite{kuang:AuThl}
78
79	<	[IMPORTANCE OF NANOSCALE TRANSPORT PROPERTIES STUDIES]
79	>	The Reverse Non-Equilibrium MD (RNEMD) methods adopt constraints that
80	>	satisfy linear momentum and total energy conservation of a system when
81	>	imposing fluxes in a simulation. Thus they are compatible with various
82	>	ensembles, including the micro-canonical (NVE) ensemble, without the
83	>	need of an external thermostat. The original approaches by
84	>	M\"{u}ller-Plathe {\it et
85	>	al.}\cite{MullerPlathe:1997xw,ISI:000080382700030} utilize simple
86	>	momentum swapping for generating energy/momentum fluxes, which is also
87	>	compatible with particles of different identities. Although simple to
88	>	implement in a simulation, this approach can create nonthermal
89	>	velocity distributions, as discovered by Tenney and
90	>	Maginn\cite{Maginn:2010}. Furthermore, this approach to kinetic energy
91	>	transfer between particles of different identities is less efficient
92	>	when the mass difference between the particles becomes significant,
93	>	which also limits its application on heterogeneous interfacial
94	>	systems.
95
96	<	Due to the importance of heat flow (and heat removal) in
97	<	nanotechnology, interfacial thermal conductance has been studied
98	<	extensively both experimentally and computationally.\cite{cahill:793}
99	<	Nanoscale materials have a significant fraction of their atoms at
100	<	interfaces, and the chemical details of these interfaces govern the
101	<	thermal transport properties.  Furthermore, the interfaces are often
102	<	heterogeneous (e.g. solid - liquid), which provides a challenge to
103	<	computational methods which have been developed for homogeneous or
104	<	bulk systems.
96	>	Recently, we developed a different approach, using Non-Isotropic
97	>	Velocity Scaling (NIVS) \cite{kuang:164101} algorithm to impose
98	>	fluxes. Compared to the momentum swapping move, it scales the velocity
99	>	vectors in two separate regions of a simulated system with respective
100	>	diagonal scaling matrices. These matrices are determined by solving a
101	>	set of equations including linear momentum and kinetic energy
102	>	conservation constraints and target flux satisfaction. This method is
103	>	able to effectively impose a wide range of kinetic energy fluxes
104	>	without obvious perturbation to the velocity distributions of the
105	>	simulated systems, regardless of the presence of heterogeneous
106	>	interfaces. We have successfully applied this approach in studying the
107	>	interfacial thermal conductance at metal-solvent
108	>	interfaces.\cite{kuang:AuThl}
109
110	<	Experimentally, the thermal properties of a number of interfaces have
111	<	been investigated.  Cahill and coworkers studied nanoscale thermal
112	<	transport from metal nanoparticle/fluid interfaces, to epitaxial
113	<	TiN/single crystal oxides interfaces, and hydrophilic and hydrophobic
114	<	interfaces between water and solids with different self-assembled
115	<	monolayers.\cite{Wilson:2002uq,PhysRevB.67.054302,doi:10.1021/jp048375k,PhysRevLett.96.186101}
116	<	Wang {\it et al.} studied heat transport through long-chain
117	<	hydrocarbon monolayers on gold substrate at individual molecular
99	<	level,\cite{Wang10082007} Schmidt {\it et al.} studied the role of
100	<	cetyltrimethylammonium bromide (CTAB) on the thermal transport between
101	<	gold nanorods and solvent,\cite{doi:10.1021/jp8051888} and Juv\'e {\it
102	<	et al.} studied the cooling dynamics, which is controlled by thermal
103	<	interface resistance of glass-embedded metal
104	<	nanoparticles.\cite{PhysRevB.80.195406} Although interfaces are
105	<	normally considered barriers for heat transport, Alper {\it et al.}
106	<	suggested that specific ligands (capping agents) could completely
107	<	eliminate this barrier
108	<	($G\rightarrow\infty$).\cite{doi:10.1021/la904855s}
110	>	However, the NIVS approach limits its application in imposing momentum
111	>	fluxes. Temperature anisotropy can happen under high momentum fluxes,
112	>	due to the nature of the algorithm. Thus, combining thermal and
113	>	momentum flux is also difficult to implement with this
114	>	approach. However, such combination may provide a means to simulate
115	>	thermal/momentum gradient coupled processes such as freeze
116	>	desalination. Therefore, developing novel approaches to extend the
117	>	application of imposed-flux method is desired.
118
119	<	The acoustic mismatch model for interfacial conductance utilizes the
120	<	acoustic impedance ($Z_a = \rho_a v^s_a$) on both sides of the
121	<	interface.\cite{swartz1989} Here, $\rho_a$ and $v^s_a$ are the density
122	<	and speed of sound in material $a$.  The phonon transmission
123	<	probability at the $a-b$ interface is
124	<	\begin{equation}
125	<	t_{ab} = \frac{4 Z_a Z_b}{\left(Z_a + Z_b \right)^2},
126	<	\end{equation}
127	<	and the interfacial conductance can then be approximated as
128	<	\begin{equation}
120	<	G_{ab} \approx \frac{1}{4} C_D v_D t_{ab}
121	<	\end{equation}
122	<	where $C_D$ is the Debye heat capacity of the hot side, and $v_D$ is
123	<	the Debye phonon velocity ($1/v_D^3 = 1/3v_L^3 + 2/3 v_T^3$) where
124	<	$v_L$ and $v_T$ are the longitudinal and transverse speeds of sound,
125	<	respectively.  For the Au/hexane and Au/toluene interfaces, the
126	<	acoustic mismatch model predicts room-temperature $G \approx 87 \mbox{
127	<	and } 129$ MW m$^{-2}$ K$^{-1}$, respectively.  However, it is not
128	<	clear how to apply the acoustic mismatch model to a
129	<	chemically-modified surface, particularly when the acoustic properties
130	<	of a monolayer film may not be well characterized.
119	>	In this paper, we improve the NIVS method and propose a novel approach
120	>	to impose fluxes. This approach separate the means of applying
121	>	momentum and thermal flux with operations in one time step and thus is
122	>	able to simutaneously impose thermal and momentum flux. Furthermore,
123	>	the approach retains desirable features of previous RNEMD approaches
124	>	and is simpler to implement compared to the NIVS method. In what
125	>	follows, we first present the method to implement the method in a
126	>	simulation. Then we compare the method on bulk fluids to previous
127	>	methods. Also, interfacial frictions are computed for a series of
128	>	interfaces.
129
132	–	[PREVIOUS METHODS INCLUDING NIVS AND THEIR LIMITATIONS]
133	–	[DIFFICULTY TO GENERATE JZKE AND JZP SIMUTANEOUSLY]
134	–
135	–	More precise computational models have also been used to study the
136	–	interfacial thermal transport in order to gain an understanding of
137	–	this phenomena at the molecular level. Recently, Hase and coworkers
138	–	employed Non-Equilibrium Molecular Dynamics (NEMD) simulations to
139	–	study thermal transport from hot Au(111) substrate to a self-assembled
140	–	monolayer of alkylthiol with relatively long chain (8-20 carbon
141	–	atoms).\cite{hase:2010,hase:2011} However, ensemble averaged
142	–	measurements for heat conductance of interfaces between the capping
143	–	monolayer on Au and a solvent phase have yet to be studied with their
144	–	approach. The comparatively low thermal flux through interfaces is
145	–	difficult to measure with Equilibrium
146	–	MD\cite{doi:10.1080/0026897031000068578} or forward NEMD simulation
147	–	methods. Therefore, the Reverse NEMD (RNEMD)
148	–	methods\cite{MullerPlathe:1997xw,kuang:164101} would be advantageous
149	–	in that they {\it apply} the difficult to measure quantity (flux),
150	–	while {\it measuring} the easily-computed quantity (the thermal
151	–	gradient).  This is particularly true for inhomogeneous interfaces
152	–	where it would not be clear how to apply a gradient {\it a priori}.
153	–	Garde and coworkers\cite{garde:nl2005,garde:PhysRevLett2009} applied
154	–	this approach to various liquid interfaces and studied how thermal
155	–	conductance (or resistance) is dependent on chemical details of a
156	–	number of hydrophobic and hydrophilic aqueous interfaces. And
157	–	recently, Luo {\it et al.} studied the thermal conductance of
158	–	Au-SAM-Au junctions using the same approach, comparing to a constant
159	–	temperature difference method.\cite{Luo20101} While this latter
160	–	approach establishes more ideal Maxwell-Boltzmann distributions than
161	–	previous RNEMD methods, it does not guarantee momentum or kinetic
162	–	energy conservation.
163	–
164	–	Recently, we have developed a Non-Isotropic Velocity Scaling (NIVS)
165	–	algorithm for RNEMD simulations\cite{kuang:164101}. This algorithm
166	–	retains the desirable features of RNEMD (conservation of linear
167	–	momentum and total energy, compatibility with periodic boundary
168	–	conditions) while establishing true thermal distributions in each of
169	–	the two slabs. Furthermore, it allows effective thermal exchange
170	–	between particles of different identities, and thus makes the study of
171	–	interfacial conductance much simpler.
172	–
173	–	[WHAT IS COVERED IN THIS MANUSCRIPT]
174	–	[MAY PUT FIGURE 1 HERE]
175	–	The work presented here deals with the Au(111) surface covered to
176	–	varying degrees by butanethiol, a capping agent with short carbon
177	–	chain, and solvated with organic solvents of different molecular
178	–	properties. Different models were used for both the capping agent and
179	–	the solvent force field parameters. Using the NIVS algorithm, the
180	–	thermal transport across these interfaces was studied and the
181	–	underlying mechanism for the phenomena was investigated.
182	–
130		\section{Methodology}
131		Similar to the NIVS methodology,\cite{kuang:164101} we consider a
132		periodic system divided into a series of slabs along a certain axis
#	Line 218 | Line 165 | thermal flux $J_z$, one would have
165		The above operations conserve the linear momentum of a periodic
166		system. To satisfy total energy conservation as well as to impose a
167		thermal flux $J_z$, one would have
168	<	%SUPPORT INFO MIGHT BE NECESSARY TO PUT EXTRA MATH IN
168	>	[SUPPORT INFO MIGHT BE NECESSARY TO PUT EXTRA MATH IN]
169		\begin{eqnarray}
170		K_c - J_z\Delta t & = & c^2 (K_c - \frac{1}{2}M_c \langle\vec{v}_c
171		\rangle^2) + \frac{1}{2}M_c (\langle \vec{v}_c \rangle + \vec{a}_c)^2 \\
#	Line 236 | Line 183 | the positive roots (which are closer to 1) are chosen.
183		scaling coefficients ($c$ and $h$) as well as $\vec{a}_c$ and
184		$\vec{a}_h$. Note that two roots of $c$ and $h$ exist
185		respectively. However, to avoid dramatic perturbations to a system,
186	<	the positive roots (which are closer to 1) are chosen.
187	<
186	>	the positive roots (which are closer to 1) are chosen. Figure
187	>	\ref{method} illustrates the implementation of this algorithm in an
188	>	individual step.
189	>
190	>	\begin{figure}
191	>	\includegraphics[width=\linewidth]{method}
192	>	\caption{Illustration of the implementation of the algorithm in a
193	>	single step. Starting from an ideal velocity distribution, the
194	>	transformation is used to apply both thermal and momentum flux from
195	>	the ``c'' slab to the ``h'' slab. As the figure shows, the thermal
196	>	distributions preserve after this operation.}
197	>	\label{method}
198	>	\end{figure}
199	>
200		By implementing these operations at a certain frequency, a steady
201		thermal and/or momentum flux can be applied and the corresponding
202		temperature and/or momentum gradients can be established.
244	–	[REFER TO NIVS PAPER]
245	–	[ADVANTAGES]
203
204	<	Steady state MD simulations have an advantage in that not many
205	<	trajectories are needed to study the relationship between thermal flux
206	<	and thermal gradients. For systems with low interfacial conductance,
207	<	one must have a method capable of generating or measuring relatively
208	<	small fluxes, compared to those required for bulk conductivity. This
209	<	requirement makes the calculation even more difficult for
210	<	slowly-converging equilibrium methods.\cite{Viscardy:2007lq} Forward
211	<	NEMD methods impose a gradient (and measure a flux), but at interfaces
212	<	it is not clear what behavior should be imposed at the boundaries
256	<	between materials.  Imposed-flux reverse non-equilibrium
257	<	methods\cite{MullerPlathe:1997xw} set the flux {\it a priori} and
258	<	the thermal response becomes an easy-to-measure quantity.  Although
259	<	M\"{u}ller-Plathe's original momentum swapping approach can be used
260	<	for exchanging energy between particles of different identity, the
261	<	kinetic energy transfer efficiency is affected by the mass difference
262	<	between the particles, which limits its application on heterogeneous
263	<	interfacial systems.
204	>	This approach is more computationaly efficient compared to the
205	>	previous NIVS method, in that only quadratic equations are involved,
206	>	while the NIVS method needs to solve a quartic equations. Furthermore,
207	>	the method implements isotropic scaling of velocities in respective
208	>	slabs, unlike the NIVS, where an extra criteria function is necessary
209	>	to choose a set of coefficients that performs the most isotropic
210	>	scaling. More importantly, separating the momentum flux imposing from
211	>	velocity scaling avoids the underlying cause that NIVS produced
212	>	thermal anisotropy when applying a momentum flux.
213
214	<	The non-isotropic velocity scaling (NIVS) \cite{kuang:164101} approach
215	<	to non-equilibrium MD simulations is able to impose a wide range of
216	<	kinetic energy fluxes without obvious perturbation to the velocity
217	<	distributions of the simulated systems. Furthermore, this approach has
218	<	the advantage in heterogeneous interfaces in that kinetic energy flux
219	<	can be applied between regions of particles of arbitrary identity, and
220	<	the flux will not be restricted by difference in particle mass.
214	>	The advantages of the approach over the original momentum swapping
215	>	approach lies in its nature to preserve a Gaussian
216	>	distribution. Because the momentum swapping tends to render a
217	>	nonthermal distribution, when the imposed flux is relatively large,
218	>	diffusion of the neighboring slabs could no longer remedy this effect,
219	>	and nonthermal distributions would be observed. Results in later
220	>	section will illustrate this effect.
221
222	<	The NIVS algorithm scales the velocity vectors in two separate regions
223	<	of a simulation system with respective diagonal scaling matrices. To
224	<	determine these scaling factors in the matrices, a set of equations
225	<	including linear momentum conservation and kinetic energy conservation
226	<	constraints and target energy flux satisfaction is solved. With the
227	<	scaling operation applied to the system in a set frequency, bulk
228	<	temperature gradients can be easily established, and these can be used
229	<	for computing thermal conductivities. The NIVS algorithm conserves
281	<	momenta and energy and does not depend on an external thermostat.
222	>	\section{Computational Details}
223	>	The algorithm has been implemented in our MD simulation code,
224	>	OpenMD\cite{Meineke:2005gd,openmd}. We compare the method with
225	>	previous RNEMD methods or equilibrium MD methods in homogeneous fluids
226	>	(Lennard-Jones and SPC/E water). And taking advantage of the method,
227	>	we simulate the interfacial friction of different heterogeneous
228	>	interfaces (gold-organic solvent and gold-SPC/E water and ice-liquid
229	>	water).
230
231	<	\subsection{Defining Interfacial Thermal Conductivity ($G$)}
231	>	\subsection{Simulation Protocols}
232	>	The systems to be investigated are set up in a orthorhombic simulation
233	>	cell with periodic boundary conditions in all three dimensions. The
234	>	$z$ axis of these cells were longer and was set as the gradient axis
235	>	of temperature and/or momentum. Thus the cells were divided into $N$
236	>	slabs along this axis, with various $N$ depending on individual
237	>	system. The $x$ and $y$ axis were usually of the same length in
238	>	homogeneous systems or close to each other where interfaces
239	>	presents. In all cases, before introducing a nonequilibrium method to
240	>	establish steady thermal and/or momentum gradients for further
241	>	measurements and calculations, canonical ensemble with a Nos\'e-Hoover
242	>	thermostat\cite{hoover85} and microcanonical ensemble equilibrations
243	>	were used to prepare systems ready for data
244	>	collections. Isobaric-isothermal equilibrations are performed before
245	>	this for SPC/E water systems to reach normal pressure (1 bar), while
246	>	similar equilibrations are used for interfacial systems to relax the
247	>	surface tensions.
248
249	<	For an interface with relatively low interfacial conductance, and a
250	<	thermal flux between two distinct bulk regions, the regions on either
251	<	side of the interface rapidly come to a state in which the two phases
252	<	have relatively homogeneous (but distinct) temperatures. The
253	<	interfacial thermal conductivity $G$ can therefore be approximated as:
254	<	\begin{equation}
291	<	G = \frac{E_{total}}{2 t L_x L_y \left( \langle T_\mathrm{hot}\rangle -
292	<	\langle T_\mathrm{cold}\rangle \right)}
293	<	\label{lowG}
294	<	\end{equation}
295	<	where ${E_{total}}$ is the total imposed non-physical kinetic energy
296	<	transfer during the simulation and ${\langle T_\mathrm{hot}\rangle}$
297	<	and ${\langle T_\mathrm{cold}\rangle}$ are the average observed
298	<	temperature of the two separated phases.  For an applied flux $J_z$
299	<	operating over a simulation time $t$ on a periodically-replicated slab
300	<	of dimensions $L_x \times L_y$, $E_{total} = 2 J_z t L_x L_y$.
249	>	While homogeneous fluid systems can be set up with random
250	>	configurations, our interfacial systems needs extra steps to ensure
251	>	the interfaces be established properly for computations. The
252	>	preparation and equilibration of butanethiol covered gold (111)
253	>	surface and further solvation and equilibration process is described
254	>	as in reference \cite{kuang:AuThl}.
255
256	<	When the interfacial conductance is {\it not} small, there are two
257	<	ways to define $G$. One common way is to assume the temperature is
258	<	discrete on the two sides of the interface. $G$ can be calculated
259	<	using the applied thermal flux $J$ and the maximum temperature
260	<	difference measured along the thermal gradient max($\Delta T$), which
261	<	occurs at the Gibbs dividing surface (Figure \ref{demoPic}). This is
262	<	known as the Kapitza conductance, which is the inverse of the Kapitza
263	<	resistance.
264	<	\begin{equation}
265	<	G=\frac{J}{\Delta T}
266	<	\label{discreteG}
267	<	\end{equation}
256	>	As for the ice/liquid water interfaces, the basal surface of ice
257	>	lattice was first constructed. Hirsch {\it et
258	>	al.}\cite{doi:10.1021/jp048434u} explored the energetics of ice
259	>	lattices with different proton orders. We refer to their results and
260	>	choose the configuration of the lowest energy after geometry
261	>	optimization as the unit cells of our ice lattices. Although
262	>	experimental solid/liquid coexistant temperature near normal pressure
263	>	is 273K, Bryk and Haymet's simulations of ice/liquid water interfaces
264	>	with different models suggest that for SPC/E, the most stable
265	>	interface is observed at 225$\pm$5K. Therefore, all our ice/liquid
266	>	water simulations were carried out under 225K. To have extra
267	>	protection of the ice lattice during initial equilibration (when the
268	>	randomly generated liquid phase configuration could release large
269	>	amount of energy in relaxation), a constraint method (REF?) was
270	>	adopted until the high energy configuration was relaxed.
271	>	[MAY ADD A FIGURE HERE FOR BASAL PLANE, MAY INCLUDE PRISM IF POSSIBLE]
272
273	<	\begin{figure}
274	<	\includegraphics[width=\linewidth]{method}
275	<	\caption{Interfacial conductance can be calculated by applying an
276	<	(unphysical) kinetic energy flux between two slabs, one located
277	<	within the metal and another on the edge of the periodic box.  The
320	<	system responds by forming a thermal gradient.  In bulk liquids,
321	<	this gradient typically has a single slope, but in interfacial
322	<	systems, there are distinct thermal conductivity domains.  The
323	<	interfacial conductance, $G$ is found by measuring the temperature
324	<	gap at the Gibbs dividing surface, or by using second derivatives of
325	<	the thermal profile.}
326	<	\label{demoPic}
327	<	\end{figure}
273	>	\subsection{Force Field Parameters}
274	>	For comparison of our new method with previous work, we retain our
275	>	force field parameters consistent with the results we will compare
276	>	with. The Lennard-Jones fluid used here for argon , and reduced unit
277	>	results are reported for direct comparison purpose.
278
279	<	Another approach is to assume that the temperature is continuous and
280	<	differentiable throughout the space. Given that $\lambda$ is also
281	<	differentiable, $G$ can be defined as its gradient ($\nabla\lambda$)
282	<	projected along a vector normal to the interface ($\mathbf{\hat{n}}$)
333	<	and evaluated at the interface location ($z_0$). This quantity,
334	<	\begin{align}
335	<	G^\prime &= \left(\nabla\lambda \cdot \mathbf{\hat{n}}\right)_{z_0} \\
336	<	&= \frac{\partial}{\partial z}\left(-\frac{J_z}{
337	<	\left(\frac{\partial T}{\partial z}\right)}\right)_{z_0} \\
338	<	&= J_z\left(\frac{\partial^2 T}{\partial z^2}\right)_{z_0} \Big/
339	<	\left(\frac{\partial T}{\partial z}\right)_{z_0}^2 \label{derivativeG}
340	<	\end{align}
341	<	has the same units as the common definition for $G$, and the maximum
342	<	of its magnitude denotes where thermal conductivity has the largest
343	<	change, i.e. the interface.  In the geometry used in this study, the
344	<	vector normal to the interface points along the $z$ axis, as do
345	<	$\vec{J}$ and the thermal gradient.  This yields the simplified
346	<	expressions in Eq. \ref{derivativeG}.
279	>	As for our water simulations, SPC/E model is used throughout this work
280	>	for consistency. Previous work for transport properties of SPC/E water
281	>	model is available\cite{Bedrov:2000,10.1063/1.3330544,Medina2011} so
282	>	that unnecessary repetition of previous methods can be avoided.
283
284	<	With temperature profiles obtained from simulation, one is able to
285	<	approximate the first and second derivatives of $T$ with finite
286	<	difference methods and calculate $G^\prime$. In what follows, both
287	<	definitions have been used, and are compared in the results.
284	>	The Au-Au interaction parameters in all simulations are described by
285	>	the quantum Sutton-Chen (QSC) formulation.\cite{PhysRevB.59.3527} The
286	>	QSC potentials include zero-point quantum corrections and are
287	>	reparametrized for accurate surface energies compared to the
288	>	Sutton-Chen potentials.\cite{Chen90} For gold/water interfaces, the
289	>	Spohr potential was adopted\cite{ISI:000167766600035} to depict
290	>	Au-H$_2$O interactions.
291
292	<	To investigate the interfacial conductivity at metal / solvent
293	<	interfaces, we have modeled a metal slab with its (111) surfaces
294	<	perpendicular to the $z$-axis of our simulation cells. The metal slab
295	<	has been prepared both with and without capping agents on the exposed
296	<	surface, and has been solvated with simple organic solvents, as
297	<	illustrated in Figure \ref{gradT}.
292	>	The small organic molecules included in our simulations are the Au
293	>	surface capping agent butanethiol and liquid hexane and toluene. The
294	>	United-Atom
295	>	models\cite{TraPPE-UA.thiols,TraPPE-UA.alkanes,TraPPE-UA.alkylbenzenes}
296	>	for these components were used in this work for better computational
297	>	efficiency, while maintaining good accuracy. We refer readers to our
298	>	previous work\cite{kuang:AuThl} for further details of these models,
299	>	as well as the interactions between Au and the above organic molecule
300	>	components.
301
302	<	With the simulation cell described above, we are able to equilibrate
303	<	the system and impose an unphysical thermal flux between the liquid
304	<	and the metal phase using the NIVS algorithm. By periodically applying
305	<	the unphysical flux, we obtained a temperature profile and its spatial
306	<	derivatives. Figure \ref{gradT} shows how an applied thermal flux can
307	<	be used to obtain the 1st and 2nd derivatives of the temperature
308	<	profile.
302	>	\subsection{Thermal conductivities}
303	>	When $\vec{j}_z(\vec{p})$ is set to zero and a target $J_z$ is set to
304	>	impose kinetic energy transfer, the method can be used for thermal
305	>	conductivity computations. Similar to previous RNEMD methods, we
306	>	assume linear response of the temperature gradient with respect to the
307	>	thermal flux in general case. And the thermal conductivity ($\lambda$)
308	>	can be obtained with the imposed kinetic energy flux and the measured
309	>	thermal gradient:
310	>	\begin{equation}
311	>	J_z = -\lambda \frac{\partial T}{\partial z}
312	>	\end{equation}
313	>	Like other imposed-flux methods, the energy flux was calculated using
314	>	the total non-physical energy transferred (${E_{total}}$) from slab
315	>	``c'' to slab ``h'', which is recorded throughout a simulation, and
316	>	the time for data collection $t$:
317	>	\begin{equation}
318	>	J_z = \frac{E_{total}}{2 t L_x L_y}
319	>	\end{equation}
320	>	where $L_x$ and $L_y$ denotes the dimensions of the plane in a
321	>	simulation cell perpendicular to the thermal gradient, and a factor of
322	>	two in the denominator is present for the heat transport occurs in
323	>	both $+z$ and $-z$ directions. The temperature gradient
324	>	${\langle\partial T/\partial z\rangle}$ can be obtained by a linear
325	>	regression of the temperature profile, which is recorded during a
326	>	simulation for each slab in a cell. For Lennard-Jones simulations,
327	>	thermal conductivities are reported in reduced units
328	>	(${\lambda^*=\lambda \sigma^2 m^{1/2} k_B^{-1}\varepsilon^{-1/2}}$).
329
330	<	\begin{figure}
331	<	\includegraphics[width=\linewidth]{gradT}
332	<	\caption{A sample of Au (111) / butanethiol / hexane interfacial
333	<	system with the temperature profile after a kinetic energy flux has
334	<	been imposed.  Note that the largest temperature jump in the thermal
335	<	profile (corresponding to the lowest interfacial conductance) is at
336	<	the interface between the butanethiol molecules (blue) and the
337	<	solvent (grey).  First and second derivatives of the temperature
338	<	profile are obtained using a finite difference approximation (lower
377	<	panel).}
378	<	\label{gradT}
379	<	\end{figure}
330	>	\subsection{Shear viscosities}
331	>	Alternatively, the method can carry out shear viscosity calculations
332	>	by switching off $J_z$. One can specify the vector
333	>	$\vec{j}_z(\vec{p})$ by choosing the three components
334	>	respectively. For shear viscosity simulations, $j_z(p_z)$ is usually
335	>	set to zero. Although for isotropic systems, the direction of
336	>	$\vec{j}_z(\vec{p})$ on the $xy$ plane should not matter, the ability
337	>	of arbitarily specifying the vector direction in our method provides
338	>	convenience in anisotropic simulations.
339
340	<	\section{Computational Details}
341	<	\subsection{Simulation Protocol}
342	<	The NIVS algorithm has been implemented in our MD simulation code,
343	<	OpenMD\cite{Meineke:2005gd,openmd}, and was used for our simulations.
344	<	Metal slabs of 6 or 11 layers of Au atoms were first equilibrated
345	<	under atmospheric pressure (1 atm) and 200K. After equilibration,
346	<	butanethiol capping agents were placed at three-fold hollow sites on
347	<	the Au(111) surfaces. These sites are either {\it fcc} or {\it
348	<	hcp} sites, although Hase {\it et al.} found that they are
349	<	equivalent in a heat transfer process,\cite{hase:2010} so we did not
350	<	distinguish between these sites in our study. The maximum butanethiol
351	<	capacity on Au surface is $1/3$ of the total number of surface Au
352	<	atoms, and the packing forms a $(\sqrt{3}\times\sqrt{3})R30^\circ$
353	<	structure\cite{doi:10.1021/ja00008a001,doi:10.1021/cr9801317}. A
354	<	series of lower coverages was also prepared by eliminating
355	<	butanethiols from the higher coverage surface in a regular manner. The
356	<	lower coverages were prepared in order to study the relation between
357	<	coverage and interfacial conductance.
340	>	Similar to thermal conductivity computations, linear response of the
341	>	momentum gradient with respect to the shear stress is assumed, and the
342	>	shear viscosity ($\eta$) can be obtained with the imposed momentum
343	>	flux (e.g. in $x$ direction) and the measured gradient:
344	>	\begin{equation}
345	>	j_z(p_x) = -\eta \frac{\partial v_x}{\partial z}
346	>	\end{equation}
347	>	where the flux is similarly defined:
348	>	\begin{equation}
349	>	j_z(p_x) = \frac{P_x}{2 t L_x L_y}
350	>	\end{equation}
351	>	with $P_x$ being the total non-physical momentum transferred within
352	>	the data collection time. Also, the velocity gradient
353	>	${\langle\partial v_x/\partial z\rangle}$ can be obtained using linear
354	>	regression of the $x$ component of the mean velocity, $\langle
355	>	v_x\rangle$, in each of the bins. For Lennard-Jones simulations, shear
356	>	viscosities are reported in reduced units
357	>	(${\eta^*=\eta\sigma^2(\varepsilon m)^{-1/2}}$).
358
359	<	The capping agent molecules were allowed to migrate during the
360	<	simulations. They distributed themselves uniformly and sampled a
361	<	number of three-fold sites throughout out study. Therefore, the
362	<	initial configuration does not noticeably affect the sampling of a
363	<	variety of configurations of the same coverage, and the final
364	<	conductance measurement would be an average effect of these
365	<	configurations explored in the simulations.
366	<
367	<	After the modified Au-butanethiol surface systems were equilibrated in
409	<	the canonical (NVT) ensemble, organic solvent molecules were packed in
410	<	the previously empty part of the simulation cells.\cite{packmol} Two
411	<	solvents were investigated, one which has little vibrational overlap
412	<	with the alkanethiol and which has a planar shape (toluene), and one
413	<	which has similar vibrational frequencies to the capping agent and
414	<	chain-like shape ({\it n}-hexane).
415	<
416	<	The simulation cells were not particularly extensive along the
417	<	$z$-axis, as a very long length scale for the thermal gradient may
418	<	cause excessively hot or cold temperatures in the middle of the
419	<	solvent region and lead to undesired phenomena such as solvent boiling
420	<	or freezing when a thermal flux is applied. Conversely, too few
421	<	solvent molecules would change the normal behavior of the liquid
422	<	phase. Therefore, our $N_{solvent}$ values were chosen to ensure that
423	<	these extreme cases did not happen to our simulations. The spacing
424	<	between periodic images of the gold interfaces is $45 \sim 75$\AA in
425	<	our simulations.
426	<
427	<	The initial configurations generated are further equilibrated with the
428	<	$x$ and $y$ dimensions fixed, only allowing the $z$-length scale to
429	<	change. This is to ensure that the equilibration of liquid phase does
430	<	not affect the metal's crystalline structure. Comparisons were made
431	<	with simulations that allowed changes of $L_x$ and $L_y$ during NPT
432	<	equilibration. No substantial changes in the box geometry were noticed
433	<	in these simulations. After ensuring the liquid phase reaches
434	<	equilibrium at atmospheric pressure (1 atm), further equilibration was
435	<	carried out under canonical (NVT) and microcanonical (NVE) ensembles.
436	<
437	<	After the systems reach equilibrium, NIVS was used to impose an
438	<	unphysical thermal flux between the metal and the liquid phases. Most
439	<	of our simulations were done under an average temperature of
440	<	$\sim$200K. Therefore, thermal flux usually came from the metal to the
441	<	liquid so that the liquid has a higher temperature and would not
442	<	freeze due to lowered temperatures. After this induced temperature
443	<	gradient had stabilized, the temperature profile of the simulation cell
444	<	was recorded. To do this, the simulation cell is divided evenly into
445	<	$N$ slabs along the $z$-axis. The average temperatures of each slab
446	<	are recorded for 1$\sim$2 ns. When the slab width $d$ of each slab is
447	<	the same, the derivatives of $T$ with respect to slab number $n$ can
448	<	be directly used for $G^\prime$ calculations: \begin{equation}
449	<	G^\prime = |J_z|\Big|\frac{\partial^2 T}{\partial z^2}\Big|
450	<	\Big/\left(\frac{\partial T}{\partial z}\right)^2
451	<	= |J_z|\Big|\frac{1}{d^2}\frac{\partial^2 T}{\partial n^2}\Big|
452	<	\Big/\left(\frac{1}{d}\frac{\partial T}{\partial n}\right)^2
453	<	= |J_z|\Big|\frac{\partial^2 T}{\partial n^2}\Big|
454	<	\Big/\left(\frac{\partial T}{\partial n}\right)^2
455	<	\label{derivativeG2}
359	>	\subsection{Interfacial friction and Slip length}
360	>	While the shear stress results in a velocity gradient within bulk
361	>	fluid phase, its effect at a solid-liquid interface could vary due to
362	>	the interaction strength between the two phases. The interfacial
363	>	friction coefficient $\kappa$ is defined to relate the shear stress
364	>	(e.g. along $x$-axis) and the relative fluid velocity tangent to the
365	>	interface:
366	>	\begin{equation}
367	>	j_z(p_x)|_{interface} = \kappa\Delta v_x|_{interface}
368		\end{equation}
369	<	The absolute values in Eq. \ref{derivativeG2} appear because the
370	<	direction of the flux $\vec{J}$ is in an opposing direction on either
371	<	side of the metal slab.
369	>	Under ``stick'' boundary condition, $\Delta v_x|_{interface}
370	>	\rightarrow 0$, which leads to $\kappa\rightarrow\infty$. However, for
371	>	``slip'' boundary condition at the solid-liquid interface, $\kappa$
372	>	becomes finite. To characterize the interfacial boundary conditions,
373	>	slip length ($\delta$) is defined using $\kappa$ and the shear
374	>	viscocity of liquid phase ($\eta$):
375	>	\begin{equation}
376	>	\delta = \frac{\eta}{\kappa}
377	>	\end{equation}
378	>	so that $\delta\rightarrow 0$ in the ``no-slip'' boundary condition,
379	>	and depicts how ``slippery'' an interface is. Figure \ref{slipLength}
380	>	illustrates how this quantity is defined and computed for a
381	>	solid-liquid interface. [MAY INCLUDE A SNAPSHOT IN FIGURE]
382
461	–	All of the above simulation procedures use a time step of 1 fs. Each
462	–	equilibration stage took a minimum of 100 ps, although in some cases,
463	–	longer equilibration stages were utilized.
464	–
465	–	\subsection{Force Field Parameters}
466	–	Our simulations include a number of chemically distinct components.
467	–	Figure \ref{demoMol} demonstrates the sites defined for both
468	–	United-Atom and All-Atom models of the organic solvent and capping
469	–	agents in our simulations. Force field parameters are needed for
470	–	interactions both between the same type of particles and between
471	–	particles of different species.
472	–
383		\begin{figure}
384	<	\includegraphics[width=\linewidth]{structures}
385	<	\caption{Structures of the capping agent and solvents utilized in
386	<	these simulations. The chemically-distinct sites (a-e) are expanded
387	<	in terms of constituent atoms for both United Atom (UA) and All Atom
388	<	(AA) force fields.  Most parameters are from References
479	<	\protect\cite{TraPPE-UA.alkanes,TraPPE-UA.alkylbenzenes,TraPPE-UA.thiols}
480	<	(UA) and \protect\cite{OPLSAA} (AA). Cross-interactions with the Au
481	<	atoms are given in Table 1 in the supporting information.}
482	<	\label{demoMol}
384	>	\includegraphics[width=\linewidth]{defDelta}
385	>	\caption{The slip length $\delta$ can be obtained from a velocity
386	>	profile of a solid-liquid interface. An example of Au/hexane
387	>	interfaces is shown.}
388	>	\label{slipLength}
389		\end{figure}
390
391	<	The Au-Au interactions in metal lattice slab is described by the
392	<	quantum Sutton-Chen (QSC) formulation.\cite{PhysRevB.59.3527} The QSC
393	<	potentials include zero-point quantum corrections and are
394	<	reparametrized for accurate surface energies compared to the
395	<	Sutton-Chen potentials.\cite{Chen90}
396	<
397	<	For the two solvent molecules, {\it n}-hexane and toluene, two
398	<	different atomistic models were utilized. Both solvents were modeled
493	<	using United-Atom (UA) and All-Atom (AA) models. The TraPPE-UA
494	<	parameters\cite{TraPPE-UA.alkanes,TraPPE-UA.alkylbenzenes} are used
495	<	for our UA solvent molecules. In these models, sites are located at
496	<	the carbon centers for alkyl groups. Bonding interactions, including
497	<	bond stretches and bends and torsions, were used for intra-molecular
498	<	sites closer than 3 bonds. For non-bonded interactions, Lennard-Jones
499	<	potentials are used.
500	<
501	<	By eliminating explicit hydrogen atoms, the TraPPE-UA models are
502	<	simple and computationally efficient, while maintaining good accuracy.
503	<	However, the TraPPE-UA model for alkanes is known to predict a slightly
504	<	lower boiling point than experimental values. This is one of the
505	<	reasons we used a lower average temperature (200K) for our
506	<	simulations. If heat is transferred to the liquid phase during the
507	<	NIVS simulation, the liquid in the hot slab can actually be
508	<	substantially warmer than the mean temperature in the simulation. The
509	<	lower mean temperatures therefore prevent solvent boiling.
510	<
511	<	For UA-toluene, the non-bonded potentials between intermolecular sites
512	<	have a similar Lennard-Jones formulation. The toluene molecules were
513	<	treated as a single rigid body, so there was no need for
514	<	intramolecular interactions (including bonds, bends, or torsions) in
515	<	this solvent model.
516	<
517	<	Besides the TraPPE-UA models, AA models for both organic solvents are
518	<	included in our studies as well. The OPLS-AA\cite{OPLSAA} force fields
519	<	were used. For hexane, additional explicit hydrogen sites were
520	<	included. Besides bonding and non-bonded site-site interactions,
521	<	partial charges and the electrostatic interactions were added to each
522	<	CT and HC site. For toluene, a flexible model for the toluene molecule
523	<	was utilized which included bond, bend, torsion, and inversion
524	<	potentials to enforce ring planarity.
525	<
526	<	The butanethiol capping agent in our simulations, were also modeled
527	<	with both UA and AA model. The TraPPE-UA force field includes
528	<	parameters for thiol molecules\cite{TraPPE-UA.thiols} and are used for
529	<	UA butanethiol model in our simulations. The OPLS-AA also provides
530	<	parameters for alkyl thiols. However, alkyl thiols adsorbed on Au(111)
531	<	surfaces do not have the hydrogen atom bonded to sulfur. To derive
532	<	suitable parameters for butanethiol adsorbed on Au(111) surfaces, we
533	<	adopt the S parameters from Luedtke and Landman\cite{landman:1998} and
534	<	modify the parameters for the CTS atom to maintain charge neutrality
535	<	in the molecule.  Note that the model choice (UA or AA) for the capping
536	<	agent can be different from the solvent. Regardless of model choice,
537	<	the force field parameters for interactions between capping agent and
538	<	solvent can be derived using Lorentz-Berthelot Mixing Rule:
539	<	\begin{eqnarray}
540	<	\sigma_{ij} & = & \frac{1}{2} \left(\sigma_{ii} + \sigma_{jj}\right) \\
541	<	\epsilon_{ij} & = & \sqrt{\epsilon_{ii}\epsilon_{jj}}
542	<	\end{eqnarray}
543	<
544	<	To describe the interactions between metal (Au) and non-metal atoms,
545	<	we refer to an adsorption study of alkyl thiols on gold surfaces by
546	<	Vlugt {\it et al.}\cite{vlugt:cpc2007154} They fitted an effective
547	<	Lennard-Jones form of potential parameters for the interaction between
548	<	Au and pseudo-atoms CH$_x$ and S based on a well-established and
549	<	widely-used effective potential of Hautman and Klein for the Au(111)
550	<	surface.\cite{hautman:4994} As our simulations require the gold slab
551	<	to be flexible to accommodate thermal excitation, the pair-wise form
552	<	of potentials they developed was used for our study.
553	<
554	<	The potentials developed from {\it ab initio} calculations by Leng
555	<	{\it et al.}\cite{doi:10.1021/jp034405s} are adopted for the
556	<	interactions between Au and aromatic C/H atoms in toluene. However,
557	<	the Lennard-Jones parameters between Au and other types of particles,
558	<	(e.g. AA alkanes) have not yet been established. For these
559	<	interactions, the Lorentz-Berthelot mixing rule can be used to derive
560	<	effective single-atom LJ parameters for the metal using the fit values
561	<	for toluene. These are then used to construct reasonable mixing
562	<	parameters for the interactions between the gold and other atoms.
563	<	Table 1 in the supporting information summarizes the
564	<	``metal/non-metal'' parameters utilized in our simulations.
565	<
566	<	\section{Results}
567	<	[L-J COMPARED TO RENMD NIVS; WATER COMPARED TO RNEMD NIVS;
568	<	SLIP BOUNDARY VS STICK BOUNDARY; ICE-WATER INTERFACES]
569	<
570	<	There are many factors contributing to the measured interfacial
571	<	conductance; some of these factors are physically motivated
572	<	(e.g. coverage of the surface by the capping agent coverage and
573	<	solvent identity), while some are governed by parameters of the
574	<	methodology (e.g. applied flux and the formulas used to obtain the
575	<	conductance). In this section we discuss the major physical and
576	<	calculational effects on the computed conductivity.
577	<
578	<	\subsection{Effects due to capping agent coverage}
579	<
580	<	A series of different initial conditions with a range of surface
581	<	coverages was prepared and solvated with various with both of the
582	<	solvent molecules. These systems were then equilibrated and their
583	<	interfacial thermal conductivity was measured with the NIVS
584	<	algorithm. Figure \ref{coverage} demonstrates the trend of conductance
585	<	with respect to surface coverage.
586	<
587	<	\begin{figure}
588	<	\includegraphics[width=\linewidth]{coverage}
589	<	\caption{The interfacial thermal conductivity ($G$) has a
590	<	non-monotonic dependence on the degree of surface capping.  This
591	<	data is for the Au(111) / butanethiol / solvent interface with
592	<	various UA force fields at $\langle T\rangle \sim $200K.}
593	<	\label{coverage}
594	<	\end{figure}
595	<
596	<	In partially covered surfaces, the derivative definition for
597	<	$G^\prime$ (Eq. \ref{derivativeG}) becomes difficult to apply, as the
598	<	location of maximum change of $\lambda$ becomes washed out.  The
599	<	discrete definition (Eq. \ref{discreteG}) is easier to apply, as the
600	<	Gibbs dividing surface is still well-defined. Therefore, $G$ (not
601	<	$G^\prime$) was used in this section.
602	<
603	<	From Figure \ref{coverage}, one can see the significance of the
604	<	presence of capping agents. When even a small fraction of the Au(111)
605	<	surface sites are covered with butanethiols, the conductivity exhibits
606	<	an enhancement by at least a factor of 3.  Capping agents are clearly
607	<	playing a major role in thermal transport at metal / organic solvent
608	<	surfaces.
391	>	In our method, a shear stress can be applied similar to shear
392	>	viscosity computations by applying an unphysical momentum flux
393	>	(e.g. $j_z(p_x)$). A corresponding velocity profile can be obtained as
394	>	shown in Figure \ref{slipLength}, in which the velocity gradients
395	>	within liquid phase and velocity difference at the liquid-solid
396	>	interface can be measured respectively. Further calculations and
397	>	characterizations of the interface can be carried out using these
398	>	data.
399
400	<	We note a non-monotonic behavior in the interfacial conductance as a
401	<	function of surface coverage. The maximum conductance (largest $G$)
402	<	happens when the surfaces are about 75\% covered with butanethiol
403	<	caps.  The reason for this behavior is not entirely clear.  One
404	<	explanation is that incomplete butanethiol coverage allows small gaps
405	<	between butanethiols to form. These gaps can be filled by transient
406	<	solvent molecules.  These solvent molecules couple very strongly with
407	<	the hot capping agent molecules near the surface, and can then carry
408	<	away (diffusively) the excess thermal energy from the surface.
400	>	\section{Results and Discussions}
401	>	\subsection{Lennard-Jones fluid}
402	>	Our orthorhombic simulation cell of Lennard-Jones fluid has identical
403	>	parameters to our previous work\cite{kuang:164101} to facilitate
404	>	comparison. Thermal conductivitis and shear viscosities were computed
405	>	with the algorithm applied to the simulations. The results of thermal
406	>	conductivity are compared with our previous NIVS algorithm. However,
407	>	since the NIVS algorithm could produce temperature anisotropy for
408	>	shear viscocity computations, these results are instead compared to
409	>	the momentum swapping approaches. Table \ref{LJ} lists these
410	>	calculations with various fluxes in reduced units.
411
620	–	There appears to be a competition between the conduction of the
621	–	thermal energy away from the surface by the capping agents (enhanced
622	–	by greater coverage) and the coupling of the capping agents with the
623	–	solvent (enhanced by interdigitation at lower coverages).  This
624	–	competition would lead to the non-monotonic coverage behavior observed
625	–	here.
626	–
627	–	Results for rigid body toluene solvent, as well as the UA hexane, are
628	–	within the ranges expected from prior experimental
629	–	work.\cite{Wilson:2002uq,cahill:793,PhysRevB.80.195406} This suggests
630	–	that explicit hydrogen atoms might not be required for modeling
631	–	thermal transport in these systems.  C-H vibrational modes do not see
632	–	significant excited state population at low temperatures, and are not
633	–	likely to carry lower frequency excitations from the solid layer into
634	–	the bulk liquid.
635	–
636	–	The toluene solvent does not exhibit the same behavior as hexane in
637	–	that $G$ remains at approximately the same magnitude when the capping
638	–	coverage increases from 25\% to 75\%.  Toluene, as a rigid planar
639	–	molecule, cannot occupy the relatively small gaps between the capping
640	–	agents as easily as the chain-like {\it n}-hexane.  The effect of
641	–	solvent coupling to the capping agent is therefore weaker in toluene
642	–	except at the very lowest coverage levels.  This effect counters the
643	–	coverage-dependent conduction of heat away from the metal surface,
644	–	leading to a much flatter $G$ vs. coverage trend than is observed in
645	–	{\it n}-hexane.
646	–
647	–	\subsection{Effects due to Solvent \& Solvent Models}
648	–	In addition to UA solvent and capping agent models, AA models have
649	–	also been included in our simulations.  In most of this work, the same
650	–	(UA or AA) model for solvent and capping agent was used, but it is
651	–	also possible to utilize different models for different components.
652	–	We have also included isotopic substitutions (Hydrogen to Deuterium)
653	–	to decrease the explicit vibrational overlap between solvent and
654	–	capping agent. Table \ref{modelTest} summarizes the results of these
655	–	studies.
656	–
412		\begin{table*}
413		\begin{minipage}{\linewidth}
414		\begin{center}
415	+
416	+	\caption{Thermal conductivity ($\lambda^*$) and shear viscosity
417	+	($\eta^*$) (in reduced units) of a Lennard-Jones fluid at
418	+	${\langle T^* \rangle = 0.72}$ and ${\rho^* = 0.85}$ computed
419	+	at various momentum fluxes.  The new method yields similar
420	+	results to previous RNEMD methods. All results are reported in
421	+	reduced unit. Uncertainties are indicated in parentheses.}
422
423	<	\caption{Computed interfacial thermal conductance ($G$ and
662	<	$G^\prime$) values for interfaces using various models for
663	<	solvent and capping agent (or without capping agent) at
664	<	$\langle T\rangle\sim$200K.  Here ``D'' stands for deuterated
665	<	solvent or capping agent molecules. Error estimates are
666	<	indicated in parentheses.}
667	<
668	<	\begin{tabular}{llccc}
423	>	\begin{tabular}{cccccc}
424		\hline\hline
425	<	Butanethiol model & Solvent & $G$ & $G^\prime$ \\
426	<	(or bare surface) & model & \multicolumn{2}{c}{(MW/m$^2$/K)} \\
425	>	\multicolumn{2}{c}{Momentum Exchange} &
426	>	\multicolumn{2}{c}{$\lambda^*$} &
427	>	\multicolumn{2}{c}{$\eta^*$} \\
428		\hline
429	<	UA    & UA hexane    & 131(9)    & 87(10)    \\
430	<	& UA hexane(D) & 153(5)    & 136(13)   \\
675	<	& AA hexane    & 131(6)    & 122(10)   \\
676	<	& UA toluene   & 187(16)   & 151(11)   \\
677	<	& AA toluene   & 200(36)   & 149(53)   \\
429	>	Swap Interval $t^*$ & Equivalent $J_z^*$ or $j_z^*(p_x)$ &
430	>	NIVS & This work & Swapping & This work \\
431		\hline
432	<	AA    & UA hexane    & 116(9)    & 129(8)    \\
433	<	& AA hexane    & 442(14)   & 356(31)   \\
434	<	& AA hexane(D) & 222(12)   & 234(54)   \\
435	<	& UA toluene   & 125(25)   & 97(60)    \\
436	<	& AA toluene   & 487(56)   & 290(42)   \\
684	<	\hline
685	<	AA(D) & UA hexane    & 158(25)   & 172(4)    \\
686	<	& AA hexane    & 243(29)   & 191(11)   \\
687	<	& AA toluene   & 364(36)   & 322(67)   \\
688	<	\hline
689	<	bare  & UA hexane    & 46.5(3.2) & 49.4(4.5) \\
690	<	& UA hexane(D) & 43.9(4.6) & 43.0(2.0) \\
691	<	& AA hexane    & 31.0(1.4) & 29.4(1.3) \\
692	<	& UA toluene   & 70.1(1.3) & 65.8(0.5) \\
432	>	0.116 & 0.16  & 7.30(0.10) & 6.25(0.23) & 3.57(0.06) & 3.53(0.16)\\
433	>	0.232 & 0.09  & 6.95(0.09) & 8.02(0.56) & 3.64(0.05) & 3.43(0.17)\\
434	>	0.463 & 0.047 & 7.19(0.07) & 6.48(0.15) & 3.52(0.16) & 3.51(0.08)\\
435	>	0.926 & 0.024 & 7.19(0.28) & 7.02(0.13) & 3.72(0.05) & 3.79(0.11)\\
436	>	1.158 & 0.019 & 7.98(0.33) & 7.37(0.23) & 3.42(0.06) & 3.53(0.06)\\
437		\hline\hline
438		\end{tabular}
439	<	\label{modelTest}
439	>	\label{LJ}
440		\end{center}
441		\end{minipage}
442		\end{table*}
443
444	<	To facilitate direct comparison between force fields, systems with the
445	<	same capping agent and solvent were prepared with the same length
446	<	scales for the simulation cells.
444	>	\subsubsection{Thermal conductivity}
445	>	Our thermal conductivity calculations with this method yields
446	>	comparable results to the previous NIVS algorithm. This indicates that
447	>	the thermal gradients rendered using this method are also close to
448	>	previous RNEMD methods. Simulations with moderately higher thermal
449	>	fluxes tend to yield more reliable thermal gradients and thus avoid
450	>	large errors, while overly high thermal fluxes could introduce side
451	>	effects such as non-linear temperature gradient response or
452	>	inadvertent phase transitions.
453
454	<	On bare metal / solvent surfaces, different force field models for
455	<	hexane yield similar results for both $G$ and $G^\prime$, and these
456	<	two definitions agree with each other very well. This is primarily an
457	<	indicator of weak interactions between the metal and the solvent.
454	>	Since the scaling operation is isotropic in this method, one does not
455	>	need extra care to ensure temperature isotropy between the $x$, $y$
456	>	and $z$ axes, while thermal anisotropy might happen if the criteria
457	>	function for choosing scaling coefficients does not perform as
458	>	expected. Furthermore, this method avoids inadvertent concomitant
459	>	momentum flux when only thermal flux is imposed, which could not be
460	>	achieved with swapping or NIVS approaches. The thermal energy exchange
461	>	in swapping ($\vec{p}_i$ in slab ``c'' with $\vec{p}_j$ in slab ``h'')
462	>	or NIVS (total slab momentum conponemts $P^\alpha$ scaled to $\alpha
463	>	P^\alpha$) would not obtain this result unless thermal flux vanishes
464	>	(i.e. $\vec{p}_i=\vec{p}_j$ or $\alpha=1$ which are trivial to apply a
465	>	thermal flux). In this sense, this method contributes to having
466	>	minimal perturbation to a simulation while imposing thermal flux.
467
468	<	For the fully-covered surfaces, the choice of force field for the
469	<	capping agent and solvent has a large impact on the calculated values
470	<	of $G$ and $G^\prime$.  The AA thiol to AA solvent conductivities are
471	<	much larger than their UA to UA counterparts, and these values exceed
472	<	the experimental estimates by a large measure.  The AA force field
473	<	allows significant energy to go into C-H (or C-D) stretching modes,
474	<	and since these modes are high frequency, this non-quantum behavior is
475	<	likely responsible for the overestimate of the conductivity.  Compared
476	<	to the AA model, the UA model yields more reasonable conductivity
718	<	values with much higher computational efficiency.
719	<
720	<	\subsubsection{Are electronic excitations in the metal important?}
721	<	Because they lack electronic excitations, the QSC and related embedded
722	<	atom method (EAM) models for gold are known to predict unreasonably
723	<	low values for bulk conductivity
724	<	($\lambda$).\cite{kuang:164101,ISI:000207079300006,Clancy:1992} If the
725	<	conductance between the phases ($G$) is governed primarily by phonon
726	<	excitation (and not electronic degrees of freedom), one would expect a
727	<	classical model to capture most of the interfacial thermal
728	<	conductance.  Our results for $G$ and $G^\prime$ indicate that this is
729	<	indeed the case, and suggest that the modeling of interfacial thermal
730	<	transport depends primarily on the description of the interactions
731	<	between the various components at the interface.  When the metal is
732	<	chemically capped, the primary barrier to thermal conductivity appears
733	<	to be the interface between the capping agent and the surrounding
734	<	solvent, so the excitations in the metal have little impact on the
735	<	value of $G$.
736	<
737	<	\subsection{Effects due to methodology and simulation parameters}
738	<
739	<	We have varied the parameters of the simulations in order to
740	<	investigate how these factors would affect the computation of $G$.  Of
741	<	particular interest are: 1) the length scale for the applied thermal
742	<	gradient (modified by increasing the amount of solvent in the system),
743	<	2) the sign and magnitude of the applied thermal flux, 3) the average
744	<	temperature of the simulation (which alters the solvent density during
745	<	equilibration), and 4) the definition of the interfacial conductance
746	<	(Eqs. (\ref{discreteG}) or (\ref{derivativeG})) used in the
747	<	calculation.
748	<
749	<	Systems of different lengths were prepared by altering the number of
750	<	solvent molecules and extending the length of the box along the $z$
751	<	axis to accomodate the extra solvent.  Equilibration at the same
752	<	temperature and pressure conditions led to nearly identical surface
753	<	areas ($L_x$ and $L_y$) available to the metal and capping agent,
754	<	while the extra solvent served mainly to lengthen the axis that was
755	<	used to apply the thermal flux.  For a given value of the applied
756	<	flux, the different $z$ length scale has only a weak effect on the
757	<	computed conductivities.
758	<
759	<	\subsubsection{Effects of applied flux}
760	<	The NIVS algorithm allows changes in both the sign and magnitude of
761	<	the applied flux.  It is possible to reverse the direction of heat
762	<	flow simply by changing the sign of the flux, and thermal gradients
763	<	which would be difficult to obtain experimentally ($5$ K/\AA) can be
764	<	easily simulated.  However, the magnitude of the applied flux is not
765	<	arbitrary if one aims to obtain a stable and reliable thermal gradient.
766	<	A temperature gradient can be lost in the noise if $|J_z|$ is too
767	<	small, and excessive $|J_z|$ values can cause phase transitions if the
768	<	extremes of the simulation cell become widely separated in
769	<	temperature.  Also, if $|J_z|$ is too large for the bulk conductivity
770	<	of the materials, the thermal gradient will never reach a stable
771	<	state.
772	<
773	<	Within a reasonable range of $J_z$ values, we were able to study how
774	<	$G$ changes as a function of this flux.  In what follows, we use
775	<	positive $J_z$ values to denote the case where energy is being
776	<	transferred by the method from the metal phase and into the liquid.
777	<	The resulting gradient therefore has a higher temperature in the
778	<	liquid phase.  Negative flux values reverse this transfer, and result
779	<	in higher temperature metal phases.  The conductance measured under
780	<	different applied $J_z$ values is listed in Tables 2 and 3 in the
781	<	supporting information. These results do not indicate that $G$ depends
782	<	strongly on $J_z$ within this flux range. The linear response of flux
783	<	to thermal gradient simplifies our investigations in that we can rely
784	<	on $G$ measurement with only a small number $J_z$ values.
785	<
786	<	The sign of $J_z$ is a different matter, however, as this can alter
787	<	the temperature on the two sides of the interface. The average
788	<	temperature values reported are for the entire system, and not for the
789	<	liquid phase, so at a given $\langle T \rangle$, the system with
790	<	positive $J_z$ has a warmer liquid phase.  This means that if the
791	<	liquid carries thermal energy via diffusive transport, {\it positive}
792	<	$J_z$ values will result in increased molecular motion on the liquid
793	<	side of the interface, and this will increase the measured
794	<	conductivity.
795	<
796	<	\subsubsection{Effects due to average temperature}
797	<
798	<	We also studied the effect of average system temperature on the
799	<	interfacial conductance.  The simulations are first equilibrated in
800	<	the NPT ensemble at 1 atm.  The TraPPE-UA model for hexane tends to
801	<	predict a lower boiling point (and liquid state density) than
802	<	experiments.  This lower-density liquid phase leads to reduced contact
803	<	between the hexane and butanethiol, and this accounts for our
804	<	observation of lower conductance at higher temperatures.  In raising
805	<	the average temperature from 200K to 250K, the density drop of
806	<	$\sim$20\% in the solvent phase leads to a $\sim$40\% drop in the
807	<	conductance.
808	<
809	<	Similar behavior is observed in the TraPPE-UA model for toluene,
810	<	although this model has better agreement with the experimental
811	<	densities of toluene.  The expansion of the toluene liquid phase is
812	<	not as significant as that of the hexane (8.3\% over 100K), and this
813	<	limits the effect to $\sim$20\% drop in thermal conductivity.
814	<
815	<	Although we have not mapped out the behavior at a large number of
816	<	temperatures, is clear that there will be a strong temperature
817	<	dependence in the interfacial conductance when the physical properties
818	<	of one side of the interface (notably the density) change rapidly as a
819	<	function of temperature.
820	<
821	<	Besides the lower interfacial thermal conductance, surfaces at
822	<	relatively high temperatures are susceptible to reconstructions,
823	<	particularly when butanethiols fully cover the Au(111) surface. These
824	<	reconstructions include surface Au atoms which migrate outward to the
825	<	S atom layer, and butanethiol molecules which embed into the surface
826	<	Au layer. The driving force for this behavior is the strong Au-S
827	<	interactions which are modeled here with a deep Lennard-Jones
828	<	potential. This phenomenon agrees with reconstructions that have been
829	<	experimentally
830	<	observed.\cite{doi:10.1021/j100035a033,doi:10.1021/la026493y}.  Vlugt
831	<	{\it et al.} kept their Au(111) slab rigid so that their simulations
832	<	could reach 300K without surface
833	<	reconstructions.\cite{vlugt:cpc2007154} Since surface reconstructions
834	<	blur the interface, the measurement of $G$ becomes more difficult to
835	<	conduct at higher temperatures.  For this reason, most of our
836	<	measurements are undertaken at $\langle T\rangle\sim$200K where
837	<	reconstruction is minimized.
838	<
839	<	However, when the surface is not completely covered by butanethiols,
840	<	the simulated system appears to be more resistent to the
841	<	reconstruction. Our Au / butanethiol / toluene system had the Au(111)
842	<	surfaces 90\% covered by butanethiols, but did not see this above
843	<	phenomena even at $\langle T\rangle\sim$300K.  That said, we did
844	<	observe butanethiols migrating to neighboring three-fold sites during
845	<	a simulation.  Since the interface persisted in these simulations, we
846	<	were able to obtain $G$'s for these interfaces even at a relatively
847	<	high temperature without being affected by surface reconstructions.
848	<
849	<	\section{Discussion}
850	<	[COMBINE W. RESULTS]
851	<	The primary result of this work is that the capping agent acts as an
852	<	efficient thermal coupler between solid and solvent phases.  One of
853	<	the ways the capping agent can carry out this role is to down-shift
854	<	between the phonon vibrations in the solid (which carry the heat from
855	<	the gold) and the molecular vibrations in the liquid (which carry some
856	<	of the heat in the solvent).
857	<
858	<	To investigate the mechanism of interfacial thermal conductance, the
859	<	vibrational power spectrum was computed. Power spectra were taken for
860	<	individual components in different simulations. To obtain these
861	<	spectra, simulations were run after equilibration in the
862	<	microcanonical (NVE) ensemble and without a thermal
863	<	gradient. Snapshots of configurations were collected at a frequency
864	<	that is higher than that of the fastest vibrations occurring in the
865	<	simulations. With these configurations, the velocity auto-correlation
866	<	functions can be computed:
867	<	\begin{equation}
868	<	C_A (t) = \langle\vec{v}_A (t)\cdot\vec{v}_A (0)\rangle
869	<	\label{vCorr}
870	<	\end{equation}
871	<	The power spectrum is constructed via a Fourier transform of the
872	<	symmetrized velocity autocorrelation function,
873	<	\begin{equation}
874	<	\hat{f}(\omega) =
875	<	\int_{-\infty}^{\infty} C_A (t) e^{-2\pi it\omega}\,dt
876	<	\label{fourier}
877	<	\end{equation}
878	<
879	<	\subsection{The role of specific vibrations}
880	<	The vibrational spectra for gold slabs in different environments are
881	<	shown as in Figure \ref{specAu}. Regardless of the presence of
882	<	solvent, the gold surfaces which are covered by butanethiol molecules
883	<	exhibit an additional peak observed at a frequency of
884	<	$\sim$165cm$^{-1}$.  We attribute this peak to the S-Au bonding
885	<	vibration. This vibration enables efficient thermal coupling of the
886	<	surface Au layer to the capping agents. Therefore, in our simulations,
887	<	the Au / S interfaces do not appear to be the primary barrier to
888	<	thermal transport when compared with the butanethiol / solvent
889	<	interfaces.  This supports the results of Luo {\it et
890	<	al.}\cite{Luo20101}, who reported $G$ for Au-SAM junctions roughly
891	<	twice as large as what we have computed for the thiol-liquid
892	<	interfaces.
893	<
894	<	\begin{figure}
895	<	\includegraphics[width=\linewidth]{vibration}
896	<	\caption{The vibrational power spectrum for thiol-capped gold has an
897	<	additional vibrational peak at $\sim $165cm$^{-1}$.  Bare gold
898	<	surfaces (both with and without a solvent over-layer) are missing
899	<	this peak.   A similar peak at  $\sim $165cm$^{-1}$ also appears in
900	<	the vibrational power spectrum for the butanethiol capping agents.}
901	<	\label{specAu}
902	<	\end{figure}
903	<
904	<	Also in this figure, we show the vibrational power spectrum for the
905	<	bound butanethiol molecules, which also exhibits the same
906	<	$\sim$165cm$^{-1}$ peak.
907	<
908	<	\subsection{Overlap of power spectra}
909	<	A comparison of the results obtained from the two different organic
910	<	solvents can also provide useful information of the interfacial
911	<	thermal transport process.  In particular, the vibrational overlap
912	<	between the butanethiol and the organic solvents suggests a highly
913	<	efficient thermal exchange between these components.  Very high
914	<	thermal conductivity was observed when AA models were used and C-H
915	<	vibrations were treated classically. The presence of extra degrees of
916	<	freedom in the AA force field yields higher heat exchange rates
917	<	between the two phases and results in a much higher conductivity than
918	<	in the UA force field. The all-atom classical models include high
919	<	frequency modes which should be unpopulated at our relatively low
920	<	temperatures.  This artifact is likely the cause of the high thermal
921	<	conductance in all-atom MD simulations.
468	>	\subsubsection{Shear viscosity}
469	>	Table \ref{LJ} also compares our shear viscosity results with momentum
470	>	swapping approach. Our calculations show that our method predicted
471	>	similar values for shear viscosities to the momentum swapping
472	>	approach, as well as the velocity gradient profiles. Moderately larger
473	>	momentum fluxes are helpful to reduce the errors of measured velocity
474	>	gradients and thus the final result. However, it is pointed out that
475	>	the momentum swapping approach tends to produce nonthermal velocity
476	>	distributions.\cite{Maginn:2010}
477
478	<	The similarity in the vibrational modes available to solvent and
479	<	capping agent can be reduced by deuterating one of the two components
480	<	(Fig. \ref{aahxntln}).  Once either the hexanes or the butanethiols
481	<	are deuterated, one can observe a significantly lower $G$ and
482	<	$G^\prime$ values (Table \ref{modelTest}).
478	>	To examine that temperature isotropy holds in simulations using our
479	>	method, we measured the three one-dimensional temperatures in each of
480	>	the slabs (Figure \ref{tempXyz}). Note that the $x$-dimensional
481	>	temperatures were calculated after subtracting the effects from bulk
482	>	velocities of the slabs. The one-dimensional temperature profiles
483	>	showed no observable difference between the three dimensions. This
484	>	ensures that isotropic scaling automatically preserves temperature
485	>	isotropy and that our method is useful in shear viscosity
486	>	computations.
487
488		\begin{figure}
489	<	\includegraphics[width=\linewidth]{aahxntln}
490	<	\caption{Spectra obtained for all-atom (AA) Au / butanethiol / solvent
491	<	systems. When butanethiol is deuterated (lower left), its
492	<	vibrational overlap with hexane decreases significantly.  Since
493	<	aromatic molecules and the butanethiol are vibrationally dissimilar,
494	<	the change is not as dramatic when toluene is the solvent (right).}
495	<	\label{aahxntln}
489	>	\includegraphics[width=\linewidth]{tempXyz}
490	>	\caption{Unlike the previous NIVS algorithm, the new method does not
491	>	produce a thermal anisotropy. No temperature difference between
492	>	different dimensions were observed beyond the magnitude of the error
493	>	bars. Note that the two ``hotter'' regions are caused by the shear
494	>	stress, as reported by Tenney and Maginn\cite{Maginn:2010}, but not
495	>	an effect that only observed in our methods.}
496	>	\label{tempXyz}
497		\end{figure}
498
499	<	For the Au / butanethiol / toluene interfaces, having the AA
500	<	butanethiol deuterated did not yield a significant change in the
501	<	measured conductance. Compared to the C-H vibrational overlap between
502	<	hexane and butanethiol, both of which have alkyl chains, the overlap
503	<	between toluene and butanethiol is not as significant and thus does
504	<	not contribute as much to the heat exchange process.
499	>	Furthermore, the velocity distribution profiles are tested by imposing
500	>	a large shear stress into the simulations. Figure \ref{vDist}
501	>	demonstrates how our method is able to maintain thermal velocity
502	>	distributions against the momentum swapping approach even under large
503	>	imposed fluxes. Previous swapping methods tend to deplete particles of
504	>	positive velocities in the negative velocity slab (``c'') and vice
505	>	versa in slab ``h'', where the distributions leave a notch. This
506	>	problematic profiles become significant when the imposed-flux becomes
507	>	larger and diffusions from neighboring slabs could not offset the
508	>	depletion. Simutaneously, abnormal peaks appear corresponding to
509	>	excessive velocity swapped from the other slab. This nonthermal
510	>	distributions limit applications of the swapping approach in shear
511	>	stress simulations. Our method avoids the above problematic
512	>	distributions by altering the means of applying momentum
513	>	fluxes. Comparatively, velocity distributions recorded from
514	>	simulations with our method is so close to the ideal thermal
515	>	prediction that no observable difference is shown in Figure
516	>	\ref{vDist}. Conclusively, our method avoids problems happened in
517	>	previous RNEMD methods and provides a useful means for shear viscosity
518	>	computations.
519
946	–	Previous observations of Zhang {\it et al.}\cite{hase:2010} indicate
947	–	that the {\it intra}molecular heat transport due to alkylthiols is
948	–	highly efficient.  Combining our observations with those of Zhang {\it
949	–	et al.}, it appears that butanethiol acts as a channel to expedite
950	–	heat flow from the gold surface and into the alkyl chain.  The
951	–	vibrational coupling between the metal and the liquid phase can
952	–	therefore be enhanced with the presence of suitable capping agents.
953	–
954	–	Deuterated models in the UA force field did not decouple the thermal
955	–	transport as well as in the AA force field.  The UA models, even
956	–	though they have eliminated the high frequency C-H vibrational
957	–	overlap, still have significant overlap in the lower-frequency
958	–	portions of the infrared spectra (Figure \ref{uahxnua}).  Deuterating
959	–	the UA models did not decouple the low frequency region enough to
960	–	produce an observable difference for the results of $G$ (Table
961	–	\ref{modelTest}).
962	–
520		\begin{figure}
521	<	\includegraphics[width=\linewidth]{uahxnua}
522	<	\caption{Vibrational power spectra for UA models for the butanethiol
523	<	and hexane solvent (upper panel) show the high degree of overlap
524	<	between these two molecules, particularly at lower frequencies.
525	<	Deuterating a UA model for the solvent (lower panel) does not
526	<	decouple the two spectra to the same degree as in the AA force
527	<	field (see Fig \ref{aahxntln}).}
528	<	\label{uahxnua}
521	>	\includegraphics[width=\linewidth]{velDist}
522	>	\caption{Velocity distributions that develop under the swapping and
523	>	our methods at high flux. These distributions were obtained from
524	>	Lennard-Jones simulations with $j_z^*(p_x)\sim 0.4$ (equivalent to a
525	>	swapping interval of 20 time steps). This is a relatively large flux
526	>	to demonstrate the nonthermal distributions that develop under the
527	>	swapping method. Distributions produced by our method are very close
528	>	to the ideal thermal situations.}
529	>	\label{vDist}
530		\end{figure}
531
532	<	\section{Conclusions}
533	<	The NIVS algorithm has been applied to simulations of
534	<	butanethiol-capped Au(111) surfaces in the presence of organic
535	<	solvents. This algorithm allows the application of unphysical thermal
536	<	flux to transfer heat between the metal and the liquid phase. With the
537	<	flux applied, we were able to measure the corresponding thermal
538	<	gradients and to obtain interfacial thermal conductivities. Under
539	<	steady states, 2-3 ns trajectory simulations are sufficient for
540	<	computation of this quantity.
532	>	\subsection{Bulk SPC/E water}
533	>	Since our method was in good performance of thermal conductivity and
534	>	shear viscosity computations for simple Lennard-Jones fluid, we extend
535	>	our applications of these simulations to complex fluid like SPC/E
536	>	water model. A simulation cell with 1000 molecules was set up in the
537	>	same manner as in \cite{kuang:164101}. For thermal conductivity
538	>	simulations, measurements were taken to compare with previous RNEMD
539	>	methods; for shear viscosity computations, simulations were run under
540	>	a series of temperatures (with corresponding pressure relaxation using
541	>	the isobaric-isothermal ensemble[CITE NIVS REF 32]), and results were
542	>	compared to available data from Equilibrium MD methods[CITATIONS].
543
544	<	Our simulations have seen significant conductance enhancement in the
545	<	presence of capping agent, compared with the bare gold / liquid
546	<	interfaces. The vibrational coupling between the metal and the liquid
547	<	phase is enhanced by a chemically-bonded capping agent. Furthermore,
548	<	the coverage percentage of the capping agent plays an important role
549	<	in the interfacial thermal transport process. Moderately low coverages
550	<	allow higher contact between capping agent and solvent, and thus could
551	<	further enhance the heat transfer process, giving a non-monotonic
992	<	behavior of conductance with increasing coverage.
544	>	\subsubsection{Thermal conductivity}
545	>	Table \ref{spceThermal} summarizes our thermal conductivity
546	>	computations under different temperatures and thermal gradients, in
547	>	comparison to the previous NIVS results\cite{kuang:164101} and
548	>	experimental measurements\cite{WagnerKruse}. Note that no appreciable
549	>	drift of total system energy or temperature was observed when our
550	>	method is applied, which indicates that our algorithm conserves total
551	>	energy even for systems involving electrostatic interactions.
552
553	<	Our results, particularly using the UA models, agree well with
554	<	available experimental data.  The AA models tend to overestimate the
555	<	interfacial thermal conductance in that the classically treated C-H
556	<	vibrations become too easily populated. Compared to the AA models, the
557	<	UA models have higher computational efficiency with satisfactory
558	<	accuracy, and thus are preferable in modeling interfacial thermal
559	<	transport.
553	>	Measurements using our method established similar temperature
554	>	gradients to the previous NIVS method. Our simulation results are in
555	>	good agreement with those from previous simulations. And both methods
556	>	yield values in reasonable agreement with experimental
557	>	values. Simulations using moderately higher thermal gradient or those
558	>	with longer gradient axis ($z$) for measurement seem to have better
559	>	accuracy, from our results.
560
561	<	Of the two definitions for $G$, the discrete form
562	<	(Eq. \ref{discreteG}) was easier to use and gives out relatively
563	<	consistent results, while the derivative form (Eq. \ref{derivativeG})
564	<	is not as versatile. Although $G^\prime$ gives out comparable results
565	<	and follows similar trend with $G$ when measuring close to fully
566	<	covered or bare surfaces, the spatial resolution of $T$ profile
567	<	required for the use of a derivative form is limited by the number of
568	<	bins and the sampling required to obtain thermal gradient information.
561	>	\begin{table*}
562	>	\begin{minipage}{\linewidth}
563	>	\begin{center}
564	>
565	>	\caption{Thermal conductivity of SPC/E water under various
566	>	imposed thermal gradients. Uncertainties are indicated in
567	>	parentheses.}
568	>
569	>	\begin{tabular}{ccccc}
570	>	\hline\hline
571	>	$\langle T\rangle$ & $\langle dT/dz\rangle$ & \multicolumn{3}{c}
572	>	{$\lambda (\mathrm{W m}^{-1} \mathrm{K}^{-1})$} \\
573	>	(K) & (K/\AA) & This work & Previous NIVS\cite{kuang:164101} &
574	>	Experiment\cite{WagnerKruse} \\
575	>	\hline
576	>	300 & 0.8 & 0.815(0.027)     & 0.770(0.008) & 0.61 \\
577	>	318 & 0.8 & 0.801(0.024)     & 0.750(0.032) & 0.64 \\
578	>	& 1.6 & 0.766(0.007)     & 0.778(0.019) &      \\
579	>	& 0.8 & 0.786(0.009)\footnote{Simulation with $L_z$
580	>	twice as long.} &              &      \\
581	>	\hline\hline
582	>	\end{tabular}
583	>	\label{spceThermal}
584	>	\end{center}
585	>	\end{minipage}
586	>	\end{table*}
587
588	<	Vlugt {\it et al.} have investigated the surface thiol structures for
589	<	nanocrystalline gold and pointed out that they differ from those of
590	<	the Au(111) surface.\cite{landman:1998,vlugt:cpc2007154} This
591	<	difference could also cause differences in the interfacial thermal
592	<	transport behavior. To investigate this problem, one would need an
593	<	effective method for applying thermal gradients in non-planar
594	<	(i.e. spherical) geometries.
588	>	\subsubsection{Shear viscosity}
589	>	The improvement our method achieves for shear viscosity computations
590	>	enables us to apply it on SPC/E water models. The series of
591	>	temperatures under which our shear viscosity calculations were carried
592	>	out covers the liquid range under normal pressure. Our simulations
593	>	predict a similar trend of $\eta$ vs. $T$ to EMD results we refer to
594	>	(Table \ref{spceShear}). Considering subtlties such as temperature or
595	>	pressure/density errors in these two series of measurements, our
596	>	results show no significant difference from those with EMD
597	>	methods. Since each value reported using our method takes only one
598	>	single trajectory of simulation, instead of average from many
599	>	trajectories when using EMD, our method provides an effective means
600	>	for shear viscosity computations.
601
602	+	\begin{table*}
603	+	\begin{minipage}{\linewidth}
604	+	\begin{center}
605	+
606	+	\caption{Computed shear viscosity of SPC/E water under different
607	+	temperatures. Results are compared to those obtained with EMD
608	+	method[CITATION]. Uncertainties are indicated in parentheses.}
609	+
610	+	\begin{tabular}{cccc}
611	+	\hline\hline
612	+	$T$ & $\partial v_x/\partial z$ & \multicolumn{2}{c}
613	+	{$\eta (\mathrm{mPa}\cdot\mathrm{s})$} \\
614	+	(K) & (10$^{10}$s$^{-1}$) & This work & Previous simulations[CITATION]\\
615	+	\hline
616	+	273 &  & 1.218(0.004) &  \\
617	+	&  & 1.140(0.012) &  \\
618	+	303 &  & 0.646(0.008) &  \\
619	+	318 &  & 0.536(0.007) &  \\
620	+	&  & 0.510(0.007) &  \\
621	+	&  &  &  \\
622	+	333 &  & 0.428(0.002) &  \\
623	+	363 &  & 0.279(0.014) &  \\
624	+	&  & 0.306(0.001) &  \\
625	+	\hline\hline
626	+	\end{tabular}
627	+	\label{spceShear}
628	+	\end{center}
629	+	\end{minipage}
630	+	\end{table*}
631	+
632	+	[MAY COMBINE JZPX AND JZKE TO RUN ONE JOB BUT GET ETA=F(T)]
633	+	[PUT RESULTS AND FIGURE HERE IF IT WORKS]
634	+	\subsection{Interfacial frictions and slip lengths}
635	+	An attractive aspect of our method is the ability to apply momentum
636	+	and/or thermal flux in nonhomogeneous systems, where molecules of
637	+	different identities (or phases) are segregated in different
638	+	regions. We have previously studied the interfacial thermal transport
639	+	of a series of metal gold-liquid
640	+	surfaces\cite{kuang:164101,kuang:AuThl}, and attemptions have been
641	+	made to investigate the relationship between this phenomenon and the
642	+	interfacial frictions.
643	+
644	+	Table \ref{etaKappaDelta} includes these computations and previous
645	+	calculations of corresponding interfacial thermal conductance. For
646	+	bare Au(111) surfaces, slip boundary conditions were observed for both
647	+	organic and aqueous liquid phases, corresponding to previously
648	+	computed low interfacial thermal conductance. Instead, the butanethiol
649	+	covered Au(111) surface appeared to be sticky to the organic liquid
650	+	molecules in our simulations. We have reported conductance enhancement
651	+	effect for this surface capping agent,\cite{kuang:AuThl} and these
652	+	observations have a qualitative agreement with the thermal conductance
653	+	results. This agreement also supports discussions on the relationship
654	+	between surface wetting and slip effect and thermal conductance of the
655	+	interface.[CITE BARRAT, GARDE]
656	+
657	+	\begin{table*}
658	+	\begin{minipage}{\linewidth}
659	+	\begin{center}
660	+
661	+	\caption{Computed interfacial friction coefficient values for
662	+	interfaces with various components for liquid and solid
663	+	phase. Error estimates are indicated in parentheses.}
664	+
665	+	\begin{tabular}{llcccccc}
666	+	\hline\hline
667	+	Solid & Liquid & $T$ & $j_z(p_x)$ & $\eta_{liquid}$ & $\kappa$
668	+	& $\delta$ & $G$\footnote{References \cite{kuang:AuThl} and
669	+	\cite{kuang:164101}.} \\
670	+	surface & model & K & MPa & mPa$\cdot$s & Pa$\cdot$s/m & nm &
671	+	MW/m$^2$/K \\
672	+	\hline
673	+	Au(111) & hexane & 200 & 1.08 & 0.20() & 5.3$\times$10$^4$() &
674	+	3.7 & 46.5 \\
675	+	&        &     & 2.15 & 0.14() & 5.3$\times$10$^4$() &
676	+	2.7 &      \\
677	+	Au-SC$_4$H$_9$ & hexane & 200 & 2.16 & 0.29() & $\infty$ & 0 &
678	+	131 \\
679	+	&        &     & 5.39 & 0.32() & $\infty$ & 0 &
680	+	\\
681	+	\hline
682	+	Au(111) & toluene & 200 & 1.08 & 0.72() & 1.?$\times$10$^5$() &
683	+	4.6 & 70.1 \\
684	+	&         &     & 2.16 & 0.54() & 1.?$\times$10$^5$() &
685	+	4.9 &      \\
686	+	Au-SC$_4$H$_9$ & toluene & 200 & 5.39 & 0.98() & $\infty$ & 0
687	+	& 187 \\
688	+	&         &     & 10.8 & 0.99() & $\infty$ & 0
689	+	&     \\
690	+	\hline
691	+	Au(111) & water & 300 & 1.08 & 0.40() & 1.9$\times$10$^4$() &
692	+	20.7 & 1.65 \\
693	+	&       &     & 2.16 & 0.79() & 1.9$\times$10$^4$() &
694	+	41.9 &      \\
695	+	\hline
696	+	ice(basal) & water & 225 & 19.4 & 15.8() & $\infty$ & 0 & \\
697	+	\hline\hline
698	+	\end{tabular}
699	+	\label{etaKappaDelta}
700	+	\end{center}
701	+	\end{minipage}
702	+	\end{table*}
703	+
704	+	An interesting effect alongside the surface friction change is
705	+	observed on the shear viscosity of liquids in the regions close to the
706	+	solid surface. Note that $\eta$ measured near a ``slip'' surface tends
707	+	to be smaller than that near a ``stick'' surface. This suggests that
708	+	an interface could affect the dynamic properties on its neighbor
709	+	regions. It is known that diffusions of solid particles in liquid
710	+	phase is affected by their surface conditions (stick or slip
711	+	boundary).[CITE SCHMIDT AND SKINNER] Our observations could provide
712	+	support to this phenomenon.
713	+
714	+	In addition to these previously studied interfaces, we attempt to
715	+	construct ice-water interfaces and the basal plane of ice lattice was
716	+	first studied. In contrast to the Au(111)/water interface, where the
717	+	friction coefficient is relatively small and large slip effect
718	+	presents, the ice/liquid water interface demonstrates strong
719	+	interactions and appears to be sticky. The supercooled liquid phase is
720	+	an order of magnitude viscous than measurements in previous
721	+	section. It would be of interst to investigate the effect of different
722	+	ice lattice planes (such as prism surface) on interfacial friction and
723	+	corresponding liquid viscosity.
724	+
725	+	\section{Conclusions}
726	+	Our simulations demonstrate the validity of our method in RNEMD
727	+	computations of thermal conductivity and shear viscosity in atomic and
728	+	molecular liquids. Our method maintains thermal velocity distributions
729	+	and avoids thermal anisotropy in previous NIVS shear stress
730	+	simulations, as well as retains attractive features of previous RNEMD
731	+	methods. There is no {\it a priori} restrictions to the method to be
732	+	applied in various ensembles, so prospective applications to
733	+	extended-system methods are possible.
734	+
735	+	Furthermore, using this method, investigations can be carried out to
736	+	characterize interfacial interactions. Our method is capable of
737	+	effectively imposing both thermal and momentum flux accross an
738	+	interface and thus facilitates studies that relates dynamic property
739	+	measurements to the chemical details of an interface.
740	+
741	+	Another attractive feature of our method is the ability of
742	+	simultaneously imposing thermal and momentum flux in a
743	+	system. potential researches that might be benefit include complex
744	+	systems that involve thermal and momentum gradients. For example, the
745	+	Soret effects under a velocity gradient would be of interest to
746	+	purification and separation researches.
747	+
748		\section{Acknowledgments}
749		Support for this project was provided by the National Science
750		Foundation under grant CHE-0848243. Computational time was provided by
#	Line 1028 | Line 757 | Dame.
757
758		\end{doublespace}
759		\end{document}
1031	–

Diff Legend

–	Removed lines
+	Added lines
<	Changed lines
>	Changed lines