1 |
\documentclass[journal = jctcce, manuscript = article]{achemso} |
2 |
\setkeys{acs}{usetitle = true} |
3 |
|
4 |
\usepackage{caption} |
5 |
\usepackage{geometry} |
6 |
\usepackage{natbib} |
7 |
\usepackage{setspace} |
8 |
\usepackage{xkeyval} |
9 |
%%%%%%%%%%%%%%%%%%%%%%% |
10 |
\usepackage{amsmath} |
11 |
\usepackage{amssymb} |
12 |
\usepackage{times} |
13 |
\usepackage{mathptm} |
14 |
\usepackage{caption} |
15 |
\usepackage{tabularx} |
16 |
\usepackage{longtable} |
17 |
\usepackage{graphicx} |
18 |
\usepackage{achemso} |
19 |
\usepackage{wrapfig} |
20 |
\usepackage[version=3]{mhchem} % this is a great package for formatting chemical reactions |
21 |
\usepackage{url} |
22 |
|
23 |
\title{A method for creating thermal and angular momentum fluxes in |
24 |
non-periodic simulations} |
25 |
|
26 |
\author{Kelsey M. Stocker} |
27 |
\author{J. Daniel Gezelter} |
28 |
\email{gezelter@nd.edu} |
29 |
\affiliation[University of Notre Dame]{251 Nieuwland Science Hall\\ Department of Chemistry and Biochemistry\\ University of Notre Dame\\ Notre Dame, Indiana 46556} |
30 |
|
31 |
\begin{document} |
32 |
|
33 |
\begin{tocentry} |
34 |
% \includegraphics[width=3.2cm]{figures/npVSS} \includegraphics[width=3.7cm]{figures/NP20} |
35 |
\includegraphics[width=9cm]{figures/TOC} |
36 |
\end{tocentry} |
37 |
|
38 |
|
39 |
\newcolumntype{A}{p{1.5in}} |
40 |
\newcolumntype{B}{p{0.75in}} |
41 |
|
42 |
% \author{Kelsey M. Stocker and J. Daniel |
43 |
% Gezelter\footnote{Corresponding author. \ Electronic mail: |
44 |
% gezelter@nd.edu} \\ |
45 |
% 251 Nieuwland Science Hall, \\ |
46 |
% Department of Chemistry and Biochemistry,\\ |
47 |
% University of Notre Dame\\ |
48 |
% Notre Dame, Indiana 46556} |
49 |
|
50 |
%\date{\today} |
51 |
|
52 |
%\maketitle |
53 |
|
54 |
%\begin{doublespace} |
55 |
|
56 |
\begin{abstract} |
57 |
|
58 |
We present a new reverse non-equilibrium molecular dynamics (RNEMD) |
59 |
method that can be used with non-periodic simulation cells. This |
60 |
method applies thermal and/or angular momentum fluxes between two |
61 |
arbitrary regions of the simulation, and is capable of creating |
62 |
stable temperature and angular velocity gradients while conserving |
63 |
total energy and angular momentum. One particularly useful |
64 |
application is the exchange of kinetic energy between two concentric |
65 |
spherical regions, which can be used to generate thermal transport |
66 |
between nanoparticles and the solvent that surrounds them. The |
67 |
rotational couple to the solvent (a measure of interfacial friction) |
68 |
is also available via this method. As demonstrations and tests of |
69 |
the new method, we have computed the thermal conductivities of gold |
70 |
nanoparticles and water clusters, the shear viscosity of a water |
71 |
cluster, the interfacial thermal conductivity ($G$) of a solvated |
72 |
gold nanoparticle and the interfacial friction of a variety of |
73 |
solvated gold nanostructures. |
74 |
|
75 |
\end{abstract} |
76 |
|
77 |
\newpage |
78 |
|
79 |
%\narrowtext |
80 |
|
81 |
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% |
82 |
% **INTRODUCTION** |
83 |
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% |
84 |
\section{Introduction} |
85 |
|
86 |
Non-equilibrium molecular dynamics (NEMD) methods impose a temperature |
87 |
or velocity {\it gradient} on a |
88 |
system,\cite{Ashurst:1975eu,Evans:1982oq,Erpenbeck:1984qe,Evans:1986nx,Vogelsang:1988qv,Maginn:1993kl,Hess:2002nr,Schelling:2002dp,Berthier:2002ai,Evans:2002tg,Vasquez:2004ty,Backer:2005sf,Jiang:2008hc,Picalek:2009rz} |
89 |
and use linear response theory to connect the resulting thermal or |
90 |
momentum {\it flux} to transport coefficients of bulk materials, |
91 |
\begin{equation} |
92 |
j_z(p_x) = -\eta \frac{\partial v_x}{\partial z}, \hspace{0.5in} |
93 |
J_z = \lambda \frac{\partial T}{\partial z}. |
94 |
\end{equation} |
95 |
Here, $\frac{\partial T}{\partial z}$ and $\frac{\partial |
96 |
v_x}{\partial z}$ are the imposed thermal and momentum gradients, |
97 |
and as long as the imposed gradients are relatively small, the |
98 |
corresponding fluxes, $J_z$ and $j_z(p_x)$, have a linear relationship |
99 |
to the gradients. The coefficients that provide this relationship |
100 |
correspond to physical properties of the bulk material, either the |
101 |
shear viscosity $(\eta)$ or thermal conductivity $(\lambda)$. For |
102 |
systems which include phase boundaries or interfaces, it is often |
103 |
unclear what gradient (or discontinuity) should be imposed at the |
104 |
boundary between materials. |
105 |
|
106 |
In contrast, reverse Non-Equilibrium Molecular Dynamics (RNEMD) |
107 |
methods impose an unphysical {\it flux} between different regions or |
108 |
``slabs'' of the simulation |
109 |
box.\cite{Muller-Plathe:1997wq,Muller-Plathe:1999ao,Patel:2005zm,Shenogina:2009ix,Tenney:2010rp,Kuang:2010if,Kuang:2011ef,Kuang:2012fe,Stocker:2013cl} |
110 |
The system responds by developing a temperature or velocity {\it |
111 |
gradient} between the two regions. The gradients which develop in |
112 |
response to the applied flux have the same linear response |
113 |
relationships to the transport coefficient of interest. Since the |
114 |
amount of the applied flux is known exactly, and measurement of a |
115 |
gradient is generally less complicated, imposed-flux methods typically |
116 |
take shorter simulation times to obtain converged results. At |
117 |
interfaces, the observed gradients often exhibit near-discontinuities |
118 |
at the boundaries between dissimilar materials. RNEMD methods do not |
119 |
need many trajectories to provide information about transport |
120 |
properties, and they have become widely used to compute thermal and |
121 |
mechanical transport in both homogeneous liquids and |
122 |
solids~\cite{Muller-Plathe:1997wq,Muller-Plathe:1999ao,Tenney:2010rp} |
123 |
as well as heterogeneous |
124 |
interfaces.\cite{Patel:2005zm,Shenogina:2009ix,Kuang:2010if,Kuang:2011ef,Kuang:2012fe,Stocker:2013cl} |
125 |
|
126 |
The strengths of specific algorithms for imposing the flux between two |
127 |
different slabs of the simulation cell has been the subject of some |
128 |
renewed interest. The original RNEMD approach used kinetic energy or |
129 |
momentum exchange between particles in the two slabs, either through |
130 |
direct swapping of momentum vectors or via virtual elastic collisions |
131 |
between atoms in the two regions. There have been recent |
132 |
methodological advances which involve scaling all particle velocities |
133 |
in both slabs.\cite{Kuang:2010if,Kuang:2012fe} Constraint equations |
134 |
are simultaneously imposed to require the simulation to conserve both |
135 |
total energy and total linear momentum. The most recent and simplest |
136 |
of the velocity scaling approaches allows for simultaneous shearing |
137 |
(to provide viscosity estimates) as well as scaling (to provide |
138 |
information about thermal conductivity).\cite{Kuang:2012fe} |
139 |
|
140 |
To date, the primary method of studying heat transfer at {\it curved} |
141 |
nanoscale interfaes has involved transient non-equilibrium molecular |
142 |
dynamics and temperature relaxation.\cite{Lervik:2009fk} RNEMD methods |
143 |
have only been used in periodic simulation cells where the exchange |
144 |
regions are physically separated along one of the axes of the |
145 |
simulation cell. This limits the applicability to infinite planar |
146 |
interfaces which are perpendicular to the applied flux. In order to |
147 |
model steady-state non-equilibrium distributions for curved surfaces |
148 |
(e.g. hot nanoparticles in contact with colder solvent), or for |
149 |
regions that are not planar slabs, the method requires some |
150 |
generalization for non-parallel exchange regions. In the following |
151 |
sections, we present a new velocity shearing and scaling (VSS) RNEMD |
152 |
algorithm which has been explicitly designed for non-periodic |
153 |
simulations, and use the method to compute some thermal transport and |
154 |
solid-liquid friction at the surfaces of spherical and ellipsoidal |
155 |
nanoparticles. |
156 |
|
157 |
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% |
158 |
% **METHODOLOGY** |
159 |
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% |
160 |
\section{Velocity shearing and scaling (VSS) for non-periodic systems} |
161 |
|
162 |
The original periodic VSS-RNEMD approach uses a series of simultaneous |
163 |
velocity shearing and scaling exchanges between the two |
164 |
slabs.\cite{Kuang:2012fe} This method imposes energy and linear |
165 |
momentum conservation constraints while simultaneously creating a |
166 |
desired flux between the two slabs. These constraints ensure that all |
167 |
configurations are sampled from the same microcanonical (NVE) |
168 |
ensemble. |
169 |
|
170 |
\begin{figure} |
171 |
\includegraphics[width=\linewidth]{figures/npVSS2} |
172 |
\caption{Schematics of periodic (left) and non-periodic (right) |
173 |
Velocity Shearing and Scaling RNEMD. A kinetic energy or momentum |
174 |
flux is applied from region B to region A. Thermal gradients are |
175 |
depicted by a color gradient. Linear or angular velocity gradients |
176 |
are shown as arrows.} |
177 |
\label{fig:VSS} |
178 |
\end{figure} |
179 |
|
180 |
We have extended the VSS method for use in {\it non-periodic} |
181 |
simulations, in which the ``slabs'' have been generalized to two |
182 |
separated regions of space. These regions could be defined as |
183 |
concentric spheres (as in figure \ref{fig:VSS}), or one of the regions |
184 |
can be defined in terms of a dynamically changing ``hull'' comprising |
185 |
the surface atoms of the cluster. This latter definition is identical |
186 |
to the hull used in the Langevin Hull algorithm.\cite{Vardeman2011} |
187 |
For the non-periodic variant, the constraints fix both the total |
188 |
energy and total {\it angular} momentum of the system while |
189 |
simultaneously imposing a thermal and angular momentum flux between |
190 |
the two regions. |
191 |
|
192 |
After a time interval of $\Delta t$, the particle velocities |
193 |
($\mathbf{v}_i$ and $\mathbf{v}_j$) in the two shells ($A$ and $B$) |
194 |
are modified by a velocity scaling coefficient ($a$ and $b$) and by a |
195 |
rotational shearing term ($\mathbf{c}_a$ and $\mathbf{c}_b$). The |
196 |
scalars $a$ and $b$ collectively provide a thermal exchange between |
197 |
the two regions. One of the values is larger than 1, and the other |
198 |
smaller. To conserve total energy and angular momentum, the values of |
199 |
these two scalars are coupled. The vectors ($\mathbf{c}_a$ and |
200 |
$\mathbf{c}_b$) provide a relative rotational shear to the velocities |
201 |
of the particles within the two regions, and these vectors must also |
202 |
be coupled to constrain the total angular momentum. |
203 |
|
204 |
Once the values of the scaling and shearing factors are known, the |
205 |
velocity changes are applied, |
206 |
\begin{displaymath} |
207 |
\begin{array}{rclcl} |
208 |
& \underline{~~~~~~~~\mathrm{scaling}~~~~~~~~} & & |
209 |
\underline{\mathrm{rotational~shearing}} \\ \\ |
210 |
\mathbf{v}_i $~~~$\leftarrow & |
211 |
a \left(\mathbf{v}_i - \langle \omega_a |
212 |
\rangle \times \mathbf{r}_i\right) & + & \mathbf{c}_a \times \mathbf{r}_i \\ |
213 |
\mathbf{v}_j $~~~$\leftarrow & |
214 |
b \left(\mathbf{v}_j - \langle \omega_b |
215 |
\rangle \times \mathbf{r}_j\right) & + & \mathbf{c}_b \times \mathbf{r}_j |
216 |
\end{array} |
217 |
\end{displaymath} |
218 |
Here $\langle\mathbf{\omega}_a\rangle$ and $\langle\mathbf{\omega}_b\rangle$ are the instantaneous angular |
219 |
velocities of each shell, and $\mathbf{r}_i$ is the position of particle $i$ relative to a fixed point in space |
220 |
(usually the center of mass of the cluster). Particles in the shells also receive an additive ``angular shear'' |
221 |
to their velocities. The amount of shear is governed by the imposed angular momentum flux, |
222 |
$\mathbf{j}_r(\mathbf{L})$, |
223 |
\begin{eqnarray} |
224 |
\mathbf{c}_a & = & - \mathbf{j}_r(\mathbf{L}) \cdot |
225 |
\overleftrightarrow{I_a}^{-1} \Delta t + \langle \omega_a \rangle \label{eq:bc}\\ |
226 |
\mathbf{c}_b & = & + \mathbf{j}_r(\mathbf{L}) \cdot |
227 |
\overleftrightarrow{I_b}^{-1} \Delta t + \langle \omega_b \rangle \label{eq:bh} |
228 |
\end{eqnarray} |
229 |
where $\overleftrightarrow{I}_{\{a,b\}}$ is the moment of inertia |
230 |
tensor for each of the two shells. |
231 |
|
232 |
To simultaneously impose a thermal flux ($J_r$) between the shells we |
233 |
use energy conservation constraints, |
234 |
\begin{eqnarray} |
235 |
K_a - J_r \Delta t & = & a^2 \left(K_a - \frac{1}{2}\langle |
236 |
\omega_a \rangle \cdot \overleftrightarrow{I_a} \cdot \langle \omega_a |
237 |
\rangle \right) + \frac{1}{2} \mathbf{c}_a \cdot \overleftrightarrow{I_a} |
238 |
\cdot \mathbf{c}_a \label{eq:Kc}\\ |
239 |
K_b + J_r \Delta t & = & b^2 \left(K_b - \frac{1}{2}\langle |
240 |
\omega_b \rangle \cdot \overleftrightarrow{I_b} \cdot \langle \omega_b |
241 |
\rangle \right) + \frac{1}{2} \mathbf{c}_b \cdot \overleftrightarrow{I_b} \cdot \mathbf{c}_b \label{eq:Kh} |
242 |
\end{eqnarray} |
243 |
Simultaneous solution of these quadratic formulae for the scaling coefficients, $a$ and $b$, will ensure that |
244 |
the simulation samples from the original microcanonical (NVE) ensemble. Here $K_{\{a,b\}}$ is the instantaneous |
245 |
translational kinetic energy of each shell. At each time interval, we solve for $a$, $b$, $\mathbf{c}_a$, and |
246 |
$\mathbf{c}_b$, subject to the imposed angular momentum flux, $j_r(\mathbf{L})$, and thermal flux, $J_r$, |
247 |
values. The new particle velocities are computed, and the simulation continues. System configurations after the |
248 |
transformations have exactly the same energy ({\it and} angular momentum) as before the moves. |
249 |
|
250 |
As the simulation progresses, the velocity transformations can be |
251 |
performed on a regular basis, and the system will develop a |
252 |
temperature and/or angular velocity gradient in response to the |
253 |
applied flux. By fitting the radial temperature gradient, it is |
254 |
straightforward to obtain the bulk thermal conductivity, |
255 |
\begin{equation} |
256 |
J_r = -\lambda \left( \frac{\partial T}{\partial r}\right) |
257 |
\end{equation} |
258 |
from the radial kinetic energy flux $(J_r)$ that was applied during |
259 |
the simulation. In practice, it is significantly easier to use the |
260 |
integrated form of Fourier's law for spherical geometries. In |
261 |
sections \ref{sec:thermal} -- \ref{sec:rotation} we outline ways of |
262 |
obtaining interfacial transport coefficients from these RNEMD |
263 |
simulations. |
264 |
|
265 |
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% |
266 |
% **COMPUTATIONAL DETAILS** |
267 |
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% |
268 |
\section{Computational Details} |
269 |
|
270 |
The new VSS-RNEMD methodology for non-periodic system geometries has |
271 |
been implemented in our group molecular dynamics code, |
272 |
OpenMD.\cite{Meineke:2005gd,openmd} We have tested the new method to |
273 |
calculate the thermal conductance of a gold nanoparticle and SPC/E |
274 |
water cluster, and compared the results with previous bulk RNEMD |
275 |
values, as well as experiment. We have also investigated the |
276 |
interfacial thermal conductance and interfacial rotational friction |
277 |
for gold nanostructures solvated in hexane as a function of |
278 |
nanoparticle size and shape. |
279 |
|
280 |
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% |
281 |
% FORCE FIELD PARAMETERS |
282 |
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% |
283 |
\subsection{Force field} |
284 |
|
285 |
Gold -- gold interactions are described by the quantum Sutton-Chen |
286 |
(QSC) model.\cite{PhysRevB.59.3527} The QSC parameters are tuned to |
287 |
experimental properties such as density, cohesive energy, and elastic |
288 |
moduli and include zero-point quantum corrections. |
289 |
|
290 |
The SPC/E water model~\cite{Berendsen87} is particularly useful for |
291 |
validation of conductivities and shear viscosities. This model has |
292 |
been used to previously test other RNEMD and NEMD approaches, and |
293 |
there are reported values for thermal conductivies and shear |
294 |
viscosities at a wide range of thermodynamic conditions that are |
295 |
available for direct comparison.\cite{Bedrov:2000,Kuang:2010if} |
296 |
|
297 |
Hexane molecules are described by the TraPPE united atom |
298 |
model,\cite{TraPPE-UA.alkanes} which provides good computational |
299 |
efficiency and reasonable accuracy for bulk thermal conductivity |
300 |
values. In this model, sites are located at the carbon centers for |
301 |
alkyl groups. Bonding interactions, including bond stretches and bends |
302 |
and torsions, were used for intra-molecular sites closer than 3 |
303 |
bonds. For non-bonded interactions, Lennard-Jones potentials were |
304 |
used. We have previously utilized both united atom (UA) and all-atom |
305 |
(AA) force fields for thermal |
306 |
conductivity,\cite{Kuang:2011ef,Kuang:2012fe,Stocker:2013cl} and since |
307 |
the united atom force fields cannot populate the high-frequency modes |
308 |
that are present in AA force fields, they appear to work better for |
309 |
modeling thermal conductance at metal/ligand interfaces. |
310 |
|
311 |
Gold -- hexane nonbonded interactions are governed by pairwise |
312 |
Lennard-Jones parameters derived from Vlugt \emph{et |
313 |
al}.\cite{vlugt:cpc2007154} They fitted parameters for the |
314 |
interaction between Au and CH$_{\emph x}$ pseudo-atoms based on the |
315 |
effective potential of Hautman and Klein for the Au(111) |
316 |
surface.\cite{hautman:4994} |
317 |
|
318 |
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% |
319 |
% NON-PERIODIC DYNAMICS |
320 |
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% |
321 |
% \subsection{Dynamics for non-periodic systems} |
322 |
% |
323 |
% We have run all tests using the Langevin Hull non-periodic simulation methodology.\cite{Vardeman2011} The |
324 |
% Langevin Hull is especially useful for simulating heterogeneous mixtures of materials with different |
325 |
% compressibilities, which are typically problematic for traditional affine transform methods. We have had |
326 |
% success applying this method to several different systems including bare metal nanoparticles, liquid water |
327 |
% clusters, and explicitly solvated nanoparticles. Calculated physical properties such as the isothermal |
328 |
% compressibility of water and the bulk modulus of gold nanoparticles are in good agreement with previous |
329 |
% theoretical and experimental results. The Langevin Hull uses a Delaunay tesselation to create a dynamic convex |
330 |
% hull composed of triangular facets with vertices at atomic sites. Atomic sites included in the hull are coupled |
331 |
% to an external bath defined by a temperature, pressure and viscosity. Atoms not included in the hull are |
332 |
% subject to standard Newtonian dynamics. |
333 |
|
334 |
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% |
335 |
% SIMULATION PROTOCOL |
336 |
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% |
337 |
\subsection{Simulation protocol} |
338 |
|
339 |
In all cases, systems were equilibrated under non-periodic |
340 |
isobaric-isothermal (NPT) conditions -- using the Langevin Hull |
341 |
methodology\cite{Vardeman2011} -- before any non-equilibrium methods |
342 |
were introduced. For heterogeneous systems, the gold nanoparticles and |
343 |
ellipsoids were created from a bulk fcc lattice and were thermally |
344 |
equilibrated before being solvated in hexane. Packmol\cite{packmol} |
345 |
was used to solvate previously equilibrated gold nanostructures within |
346 |
a spherical droplet of hexane. |
347 |
|
348 |
Once equilibrated, thermal or angular momentum fluxes were applied for |
349 |
1 - 2 ns, until stable temperature or angular velocity gradients had |
350 |
developed. Systems containing liquids were run under moderate pressure |
351 |
(5 atm) and temperatures (230 K) to avoid the formation of a vapor |
352 |
layer at the boundary of the cluster. Pressure was applied to the |
353 |
system via the non-periodic Langevin Hull.\cite{Vardeman2011} However, |
354 |
thermal coupling to the external temperature and pressure bath was |
355 |
removed to avoid interference with the imposed RNEMD flux. |
356 |
|
357 |
Because the method conserves \emph{total} angular momentum, systems |
358 |
which contain a metal nanoparticle embedded in a significant volume of |
359 |
solvent will still experience nanoparticle diffusion inside the |
360 |
solvent droplet. To aid in computing the rotational friction in these |
361 |
systems, a single gold atom at the origin of the coordinate system was |
362 |
assigned a mass $10,000 \times$ its original mass. The bonded and |
363 |
nonbonded interactions for this atom remain unchanged and the heavy |
364 |
atom is excluded from the RNEMD exchanges. The only effect of this |
365 |
gold atom is to effectively pin the nanoparticle at the origin of the |
366 |
coordinate system, while still allowing for rotation. For rotation of |
367 |
the gold ellipsoids we added two of these heavy atoms along the axis |
368 |
of rotation, separated by an equal distance from the origin of the |
369 |
coordinate system. These heavy atoms prevent off-axis tumbling of the |
370 |
nanoparticle and allow for measurement of rotational friction relative |
371 |
to a particular axis of the ellipsoid. |
372 |
|
373 |
Angular velocity data was collected for the heterogeneous systems |
374 |
after a brief period of imposed flux to initialize rotation of the |
375 |
solvated nanostructure. Doing so ensures that we overcome the initial |
376 |
static friction and calculate only the \emph{dynamic} interfacial |
377 |
rotational friction. |
378 |
|
379 |
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% |
380 |
% THERMAL CONDUCTIVITIES |
381 |
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% |
382 |
\subsection{Thermal conductivities} |
383 |
\label{sec:thermal} |
384 |
|
385 |
To compute the thermal conductivities of bulk materials, the |
386 |
integrated form of Fourier's Law of heat conduction in radial |
387 |
coordinates can be used to obtain an expression for the heat transfer |
388 |
rate between concentric spherical shells: |
389 |
\begin{equation} |
390 |
q_r = - \frac{4 \pi \lambda (T_b - T_a)}{\frac{1}{r_a} - \frac{1}{r_b}} |
391 |
\label{eq:Q} |
392 |
\end{equation} |
393 |
The heat transfer rate, $q_r$, is constant in spherical geometries, |
394 |
while the heat flux, $J_r$ depends on the surface area of the two |
395 |
shells. $\lambda$ is the thermal conductivity, and $T_{a,b}$ and |
396 |
$r_{a,b}$ are the temperatures and radii of the two concentric RNEMD |
397 |
regions, respectively. |
398 |
|
399 |
A kinetic energy flux is created using VSS-RNEMD moves, and the |
400 |
temperature in each of the radial shells is recorded. The resulting |
401 |
temperature profiles are analyzed to yield information about the |
402 |
interfacial thermal conductance. As the simulation progresses, the |
403 |
VSS moves are performed on a regular basis, and the system develops a |
404 |
thermal or velocity gradient in response to the applied flux. Once a |
405 |
stable thermal gradient has been established between the two regions, |
406 |
the thermal conductivity, $\lambda$, can be calculated using a linear |
407 |
regression of the thermal gradient, $\langle \frac{dT}{dr} \rangle$: |
408 |
|
409 |
\begin{equation} |
410 |
\lambda = \frac{r_a}{r_b} \frac{q_r}{4 \pi \langle \frac{dT}{dr} \rangle} |
411 |
\label{eq:lambda} |
412 |
\end{equation} |
413 |
|
414 |
The rate of heat transfer, $q_r$, between the two RNEMD regions is |
415 |
easily obtained from either the applied kinetic energy flux and the |
416 |
area of the smaller of the two regions, or from the total amount of |
417 |
transferred kinetic energy and the run time of the simulation. |
418 |
\begin{equation} |
419 |
q_r = J_r A = \frac{KE}{t} |
420 |
\label{eq:heat} |
421 |
\end{equation} |
422 |
|
423 |
\subsubsection{Thermal conductivity of nanocrystalline gold} |
424 |
Calculated values for the thermal conductivity of a 40 \AA$~$ radius |
425 |
gold nanoparticle (15707 atoms) at a range of kinetic energy flux |
426 |
values are shown in Table \ref{table:goldTC}. For these calculations, |
427 |
the hot and cold slabs were excluded from the linear regression of the |
428 |
thermal gradient. |
429 |
|
430 |
\begin{table} |
431 |
\caption{Calculated thermal gradients of a crystalline gold |
432 |
nanoparticle of radius 40 \AA\ subject to a range of applied |
433 |
kinetic energy fluxes. Calculations were performed at an average |
434 |
temperature of 300 K. Gold-gold interactions are described by the |
435 |
Quantum Sutton-Chen potential.} |
436 |
\label{table:goldTC} |
437 |
\begin{tabular}{cc} |
438 |
\\ \hline \hline |
439 |
{$J_r$} & {$\langle dT/dr \rangle$} \\ |
440 |
{\footnotesize(kcal / fs $\cdot$ \AA$^{2}$)} & {\footnotesize(K / |
441 |
\AA)} \\ \hline \\ |
442 |
3.25$\times 10^{-6}$ & 0.114 \\ |
443 |
6.50$\times 10^{-6}$ & 0.232 \\ |
444 |
1.30$\times 10^{-5}$ & 0.449 \\ |
445 |
3.25$\times 10^{-5}$ & 1.180 \\ |
446 |
6.50$\times 10^{-5}$ & 2.339 \\ \hline \hline |
447 |
\end{tabular} |
448 |
\end{table} |
449 |
|
450 |
The measured thermal gradients $\langle dT/dr \rangle$ are linearly |
451 |
dependent on the applied kinetic energy flux $J_r$. The calculated |
452 |
thermal conductivity value, $\lambda = 1.004 \pm |
453 |
0.009$~W~/~m~$\cdot$~K compares well with previous {\it bulk} QSC |
454 |
values of 1.08~--~1.26~W~/~m~$\cdot$~K\cite{Kuang:2010if}, though |
455 |
still significantly lower than the experimental value of |
456 |
320~W~/~m~$\cdot$~K, as the QSC force field neglects the significant |
457 |
electronic contributions to thermal conductivity. We note that there |
458 |
is only a minimal effect on the computed thermal conductivity of gold |
459 |
when comparing periodic (bulk) simulations carried out at the same |
460 |
densities as those used here. |
461 |
|
462 |
\subsubsection{Thermal conductivity of a droplet of SPC/E water} |
463 |
|
464 |
Calculated values for the thermal conductivity of a cluster of 6912 |
465 |
SPC/E water molecules were computed with a range of applied kinetic |
466 |
energy fluxes. As with the gold nanoparticle, the measured slopes were |
467 |
linearly dependent on the applied kinetic energy flux $J_r$, and the |
468 |
RNEMD regions were excluded from the $\langle dT/dr \rangle$ fit (see |
469 |
Fig. \ref{fig:lambda_G}). |
470 |
|
471 |
The computed mean value of the thermal conductivity, $\lambda = 0.884 |
472 |
\pm 0.013$~W~/~m~$\cdot$~K, compares well with previous |
473 |
non-equilibrium molecular dynamics simulations of bulk SPC/E that were |
474 |
carried out in periodic geometries.\cite{Zhang2005,Romer2012} These |
475 |
simualtions gave conductivity values from 0.81--0.87~W~/~m~$\cdot$~K, |
476 |
and all of the simulated values are approximately 30-45\% higher than the |
477 |
experimental thermal conductivity of water, which has been measured at |
478 |
0.61~W~/~m~$\cdot$~K.\cite{WagnerKruse} |
479 |
|
480 |
\begin{figure} |
481 |
\includegraphics[width=\linewidth]{figures/lambda_G} |
482 |
\caption{Upper panel: temperature profile of a typical SPC/E water |
483 |
droplet. The RNEMD simulation transfers kinetic energy from region |
484 |
A to region B, and the system responds by creating a temperature |
485 |
gradient (circles). The fit to determine $\langle dT/dr \rangle$ is |
486 |
carried out between the two RNEMD exchange regions. Lower panel: |
487 |
temperature profile for a solvated Au nanoparticle ($r = 20$\AA). |
488 |
The temperature differences used to compute the total Kapitza |
489 |
resistance (and interfacial thermal conductance) are indicated with |
490 |
arrows.} |
491 |
\label{fig:lambda_G} |
492 |
\end{figure} |
493 |
|
494 |
|
495 |
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% |
496 |
% INTERFACIAL THERMAL CONDUCTANCE |
497 |
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% |
498 |
\subsection{Interfacial thermal conductance} |
499 |
\label{sec:interfacial} |
500 |
|
501 |
The interfacial thermal conductance, $G$, of a heterogeneous interface |
502 |
located at $r_0$ can be understood as the change in thermal |
503 |
conductivity in a direction normal $(\mathbf{\hat{n}})$ to the |
504 |
interface, |
505 |
\begin{equation} |
506 |
G = \left(\nabla\lambda \cdot \mathbf{\hat{n}}\right)_{r_0} |
507 |
\end{equation} |
508 |
For heterogeneous systems such as the solvated nanoparticle shown in |
509 |
Figure \ref{fig:NP20}, the interfacial thermal conductance at the |
510 |
surface of the nanoparticle can be determined using a kinetic energy |
511 |
flux applied using the RNEMD method developed above. |
512 |
|
513 |
In spherical geometries, it is most convenient to begin by finding the |
514 |
Kapitza or interfacial thermal resistance for a thin spherical shell, |
515 |
\begin{equation} |
516 |
R_K = \frac{1}{G} = \frac{\Delta |
517 |
T}{J_r} |
518 |
\end{equation} |
519 |
where $\Delta T$ is the temperature drop from the interior to the |
520 |
exterior of the shell. For two concentric shells, the kinetic energy |
521 |
flux ($J_r$) is not the same (as the surface areas are not the same), |
522 |
but the heat transfer rate, $q_r = J_r A$ is constant. The thermal |
523 |
resistance of a shell with interior radius $r$ is most conveniently |
524 |
written as |
525 |
\begin{equation} |
526 |
R_K = \frac{1}{q_r} \Delta T 4 \pi r^2. |
527 |
\label{eq:RK} |
528 |
\end{equation} |
529 |
|
530 |
\begin{figure} |
531 |
\includegraphics[width=\linewidth]{figures/NP20} |
532 |
\caption{A gold nanoparticle with a radius of 20 \AA$\,$ solvated in |
533 |
TraPPE-UA hexane. A kinetic energy flux is applied between the |
534 |
nanoparticle and an outer shell of solvent to obtain the interfacial |
535 |
thermal conductance, $G$, and the interfacial rotational resistance |
536 |
$\Xi^{rr}$ is determined using an angular momentum flux. } |
537 |
\label{fig:NP20} |
538 |
\end{figure} |
539 |
|
540 |
To model the thermal conductance across a wide interface (or multiple |
541 |
interfaces) it is useful to consider concentric shells as resistors |
542 |
wired in series. The resistance of the shells is then additive, and |
543 |
the interfacial thermal conductance is the inverse of the total |
544 |
Kapitza resistance: |
545 |
\begin{equation} |
546 |
\frac{1}{G} = R_\mathrm{total} = \frac{1}{q_r} \sum_i \left(T_{i+i} - |
547 |
T_i\right) 4 \pi r_i^2 |
548 |
\end{equation} |
549 |
This series can be extended for any number of adjacent shells, |
550 |
allowing for the calculation of the interfacial thermal conductance |
551 |
for interfaces of considerable thickness, such as self-assembled |
552 |
ligand monolayers on a metal surface. |
553 |
|
554 |
\subsubsection{Interfacial thermal conductance of solvated gold |
555 |
nanoparticles} |
556 |
Calculated interfacial thermal conductance ($G$) values for three |
557 |
sizes of gold nanoparticles and a flat Au(111) surface solvated in |
558 |
TraPPE-UA hexane~\cite{Stocker:2013cl} are shown in Table |
559 |
\ref{table:G}. |
560 |
|
561 |
\begin{table} |
562 |
\caption{Calculated interfacial thermal conductance ($G$) values for |
563 |
gold nanoparticles of varying radii solvated in TraPPE-UA |
564 |
hexane. The nanoparticle $G$ values are compared to previous |
565 |
simulation results for a Au(111) interface in TraPPE-UA hexane.} |
566 |
\label{table:G} |
567 |
\begin{tabular}{cc} |
568 |
\hline \hline |
569 |
{Nanoparticle Radius} & {$G$}\\ |
570 |
{\footnotesize(\AA)} & {\footnotesize(MW / m$^{2}$ $\cdot$ K)}\\ \hline |
571 |
20 & {47.1 $\pm$ 5.0} \\ |
572 |
30 & {45.4 $\pm$ 1.2} \\ |
573 |
40 & {46.5 $\pm$ 2.1} \\ |
574 |
\hline |
575 |
Au(111) & {30.2} \\ |
576 |
\hline \hline |
577 |
\end{tabular} |
578 |
\end{table} |
579 |
|
580 |
The introduction of surface curvature increases the interfacial |
581 |
thermal conductance by a factor of approximately $1.5$ relative to the |
582 |
flat interface, although there is no apparent size-dependence in the |
583 |
$G$ values obtained from our calculations. This is quite different |
584 |
behavior than the size-dependence observed by Lervik \textit{et |
585 |
al.}\cite{Lervik:2009fk} in their NEMD simulations of decane |
586 |
droplets in water. In the liquid/liquid case, the surface tension is |
587 |
strongly dependent on the droplet size, and the interfacial |
588 |
resistivity increases with surface tension. In our case, the surface |
589 |
tension of the solid / liquid interface does not have a strong |
590 |
dependence on the particle radius at these particle sizes. |
591 |
|
592 |
Simulations of larger nanoparticles may yet demonstrate limiting $G$ |
593 |
values close to the flat Au(111) slab, although any spherical particle |
594 |
will have a significant fraction of the surface atoms in non-111 |
595 |
facets. It is not clear at this point whether the increase in thermal |
596 |
conductance for the spherical particles is due to the increased |
597 |
population of undercoordinated surface atoms, or to increased solid |
598 |
angle space available for phonon scattering into the solvent. These |
599 |
two possibilities do open up some interesting avenues for |
600 |
investigation. |
601 |
|
602 |
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% |
603 |
% INTERFACIAL ROTATIONAL FRICTION |
604 |
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% |
605 |
\subsection{Interfacial rotational friction} |
606 |
\label{sec:rotation} |
607 |
The interfacial rotational resistance tensor, $\Xi^{rr}$, can be |
608 |
calculated for heterogeneous nanostructure/solvent systems by applying |
609 |
an angular momentum flux between the solvated nanostructure and a |
610 |
spherical shell of solvent at the outer edge of the cluster. An |
611 |
angular velocity gradient develops in response to the applied flux, |
612 |
causing the nanostructure and solvent shell to rotate in opposite |
613 |
directions about a given axis. |
614 |
|
615 |
\begin{figure} |
616 |
\includegraphics[width=\linewidth]{figures/E25-75} |
617 |
\caption{A gold prolate ellipsoid of length 65 \AA$\,$ and width 25 |
618 |
\AA$\,$ solvated by TraPPE-UA hexane. An angular momentum flux is |
619 |
applied between the ellipsoid and an outer shell of solvent.} |
620 |
\label{fig:E25-75} |
621 |
\end{figure} |
622 |
|
623 |
Analytical solutions for the diagonal elements of the rotational |
624 |
resistance tensor for solvated spherical body of radius $r$ under |
625 |
ideal stick boundary conditions can be estimated using Stokes' law |
626 |
\begin{equation} |
627 |
\Xi^{rr}_{stick} = 8 \pi \eta r^3, |
628 |
\label{eq:Xisphere} |
629 |
\end{equation} |
630 |
where $\eta$ is the dynamic viscosity of the surrounding solvent. |
631 |
|
632 |
For general ellipsoids with semiaxes $\alpha$, $\beta$, and $\gamma$, |
633 |
Perrin's extension of Stokes' law provides exact solutions for |
634 |
symmetric prolate $(\alpha \geq \beta = \gamma)$ and oblate $(\alpha < |
635 |
\beta = \gamma)$ ellipsoids under ideal stick conditions. The Perrin |
636 |
elliptic integral parameter, |
637 |
\begin{equation} |
638 |
S = \frac{2}{\sqrt{\alpha^2 - \beta^2}} \ln \left[ \frac{\alpha + \sqrt{\alpha^2 - \beta^2}}{\beta} \right]. |
639 |
\label{eq:S} |
640 |
\end{equation} |
641 |
|
642 |
For a prolate ellipsoidal nanoparticle (see Fig. \ref{fig:E25-75}), |
643 |
the rotational resistance tensor $\Xi^{rr}$ is a $3 \times 3$ diagonal |
644 |
matrix with elements |
645 |
\begin{align} |
646 |
\Xi^{rr}_\alpha =& \frac{32 \pi}{3} \eta \frac{ \left( |
647 |
\alpha^2 - \beta^2 \right) \beta^2}{2\alpha - \beta^2 S} |
648 |
\\ \nonumber \\ |
649 |
\Xi^{rr}_{\beta,\gamma} =& \frac{32 \pi}{3} \eta \frac{ \left( \alpha^4 - \beta^4 \right)}{ \left( 2\alpha^2 - \beta^2 \right)S - 2\alpha}, |
650 |
\label{eq:Xirr} |
651 |
\end{align} |
652 |
corresponding to rotation about the long axis ($\alpha$), and each of |
653 |
the equivalent short axes ($\beta$ and $\gamma$), respectively. |
654 |
|
655 |
Previous VSS-RNEMD simulations of the interfacial friction of the |
656 |
planar Au(111) / hexane interface have shown that the flat interface |
657 |
exists within slip boundary conditions.\cite{Kuang:2012fe} Hu and |
658 |
Zwanzig\cite{Zwanzig} investigated the rotational friction |
659 |
coefficients for spheroids under slip boundary conditions and obtained |
660 |
numerial results for a scaling factor to be applied to |
661 |
$\Xi^{rr}_{\mathit{stick}}$ as a function of ratio of the shorter |
662 |
semiaxes to the longer. Under slip conditions, rotation of any |
663 |
sphere, and rotation of a prolate ellipsoid about its long axis does |
664 |
not displace any solvent, so the resulting $\Xi^{rr}_{\mathit{slip}}$ |
665 |
approaches $0$. For rotation of a prolate ellipsoid (with the aspect |
666 |
ratio shown here) about its short axis, Hu and Zwanzig obtained |
667 |
$\Xi^{rr}_{\mathit{slip}}$ values of $35.9\%$ of the analytical |
668 |
$\Xi^{rr}_{\mathit{stick}}$ result. This reduced rotational |
669 |
resistance accounts for the reduced interfacial friction under slip |
670 |
boundary conditions. |
671 |
|
672 |
The measured rotational friction coefficient, |
673 |
$\Xi^{rr}_{\mathit{meas}}$ at the interface can be extracted from |
674 |
non-periodic VSS-RNEMD simulations quite easily using the total |
675 |
applied torque ($\tau$) and the observed angular velocity of the solid |
676 |
structure ($\omega_\mathrm{solid}$), |
677 |
\begin{equation} |
678 |
\Xi^{rr}_{\mathit{meas}} = \frac{\tau}{\omega_\mathrm{solid}} |
679 |
\label{eq:Xieff} |
680 |
\end{equation} |
681 |
|
682 |
The total applied torque required to overcome the interfacial friction |
683 |
and maintain constant rotation of the gold is |
684 |
\begin{equation} |
685 |
\tau = \frac{j_r(\mathbf{L}) \cdot A}{2} |
686 |
\label{eq:tau} |
687 |
\end{equation} |
688 |
where $j_r(\mathbf{L})$ is the applied angular momentum flux, $A$ is |
689 |
the surface area of the solid nanoparticle, and the factor of $2$ is |
690 |
present because the torque is exerted on both the particle and the |
691 |
counter-rotating solvent shell. |
692 |
|
693 |
\subsubsection{Rotational friction for gold nanostructures in hexane} |
694 |
|
695 |
Table \ref{table:couple} shows the calculated rotational friction |
696 |
coefficients $\Xi^{rr}$ for spherical gold nanoparticles and a prolate |
697 |
ellipsoidal gold nanorod solvated in TraPPE-UA hexane. An angular |
698 |
momentum flux was applied between an outer shell of solvent that was |
699 |
at least 20 \AA\ away from the surface of the particle (region $A$) |
700 |
and the gold nanostructure (region $B$). Dynamic viscosity estimates |
701 |
$(\eta)$ for TraPPE-UA hexane under these particular temperature and |
702 |
pressure conditions was determined by applying a traditional VSS-RNEMD |
703 |
linear momentum flux to a periodic box of this solvent. |
704 |
|
705 |
\begin{table} |
706 |
\caption{Comparison of rotational friction coefficients under ideal stick |
707 |
($\Xi^{rr}_{\mathit{stick}}$) and |
708 |
slip ($\Xi^{rr}_{\mathit{slip}}$) conditions, along with the measured |
709 |
($\Xi^{rr}_{\mathit{meas}}$) rotational friction coefficients of gold |
710 |
nanostructures solvated in TraPPE-UA hexane at 230 K. The ellipsoid |
711 |
is oriented with the long axis along the $z$ direction.} |
712 |
\label{table:couple} |
713 |
\begin{tabular}{lccccc} |
714 |
\\ \hline \hline |
715 |
{Structure} & {Axis of Rotation} & {$\Xi^{rr}_{\mathit{stick}}$} & {$\Xi^{rr}_{\mathit{slip}}$} & {$\Xi^{rr}_{\mathit{meas}}$} & {$\Xi^{rr}_{\mathit{meas}}$ / $\Xi^{rr}_{\mathit{stick}}$}\\ |
716 |
{} & {} & {\footnotesize(amu A$^2$ / fs)} & {\footnotesize(amu A$^2$ / fs)} & {\footnotesize(amu A$^2$ / fs)} & \\ \hline |
717 |
Sphere (r = 20 \AA) & {$x = y = z$} & {3314} & 0 & {2386 $\pm$ 14} & {0.720 $\pm$ 0.004}\\ |
718 |
Sphere (r = 30 \AA) & {$x = y = z$} & {11749}& 0 & {8415 $\pm$ 274} & {0.716 $\pm$ 0.023}\\ |
719 |
Sphere (r = 40 \AA) & {$x = y = z$} & {34464}& 0 & {47544 $\pm$ 3051} & {1.380 $\pm$ 0.089}\\ |
720 |
Prolate Ellipsoid & {$x = y$} & {4991} & {1792} & {3128 $\pm$ 166} & {0.627 $\pm$ 0.033}\\ |
721 |
Prolate Ellipsoid & {$z$} & {1993} & 0 & {1590 $\pm$ 30} & {0.798 $\pm$ 0.015} |
722 |
\\ \hline \hline |
723 |
\end{tabular} |
724 |
\end{table} |
725 |
|
726 |
The results for $\Xi^{rr}_{\mathit{meas}}$ show that, in contrast with |
727 |
the flat Au(111) / hexane interface, gold nanostructures solvated by |
728 |
hexane are closer to stick than slip boundary conditions. At these |
729 |
length scales, the nanostructures are not perfect spheroids due to |
730 |
atomic `roughening' of the surface and therefore experience increased |
731 |
interfacial friction, deviating significantly from the ideal slip |
732 |
case. The 20 and 30~\AA\ radius nanoparticles experience approximately |
733 |
70\% of the ideal stick boundary interfacial friction. Rotation of the |
734 |
ellipsoid about its long axis more closely approaches the stick limit |
735 |
than rotation about the short axis, which at first seems |
736 |
counterintuitive. However, the `propellor' motion caused by rotation |
737 |
around the short axis may exclude solvent from the rotation cavity or |
738 |
move a sufficient amount of solvent along with the gold that a smaller |
739 |
interfacial friction is actually experienced. |
740 |
|
741 |
The largest nanoparticle |
742 |
(40 \AA$\,$ radius) appears to experience interfacial friction in |
743 |
excess of ideal stick conditions. Although the solvent velocity |
744 |
field does not appear in the calculation of $\Xi^{rr}$, we do note |
745 |
that the smaller particles have solvent velocity fields that can be |
746 |
fit with the $v(r,\theta) = \left( A r + B/r^2 |
747 |
\right) \sin\theta $ form,\cite{Jeffery:1915fk} while the solvent velocity fields |
748 |
for the 40 \AA\ radius particle approaches the infinite fluid |
749 |
solution and can be fit well with only the second term. |
750 |
|
751 |
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% |
752 |
% **DISCUSSION** |
753 |
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% |
754 |
\section{Discussion} |
755 |
|
756 |
We have adapted the VSS-RNEMD methodology so that it can be used to |
757 |
apply kinetic energy and angular momentum fluxes in explicitly |
758 |
non-periodic simulation geometries. Non-periodic VSS-RNEMD preserves |
759 |
the strengths of the original periodic variant, specifically Boltzmann |
760 |
velocity distributions and minimal thermal anisotropy, while extending |
761 |
the constraints to conserve total energy and total \emph{angular} |
762 |
momentum. This method also maintains the ability to impose the kinetic |
763 |
energy and angular momentum fluxes either jointly or individually. |
764 |
|
765 |
This method enables steady-state calculation of interfacial thermal |
766 |
conductance and interfacial rotational friction in heterogeneous |
767 |
clusters. We have demonstrated the abilities of the method on some |
768 |
familiar systems (nanocrystalline gold, SPC/E water) as well as on |
769 |
interfaces that have been studied only in planar geometries. For the |
770 |
bulk systems that are well-understood in periodic geometries, the |
771 |
nonperiodic VSS-RNEMD provides very similar estimates of the thermal |
772 |
conductivity to other simulation methods. |
773 |
|
774 |
For nanoparticle-to-solvent heat transfer, we have observed that the |
775 |
interfacial conductance exhibits a $1.5\times$ increase over the |
776 |
planar Au(111) interface, although we do not see any dependence of the |
777 |
conductance on the particle size. Because the surface tension effects |
778 |
present in liquid/liquid heat transfer are not strong contributors |
779 |
here, we suggest two possible mechanisms for the $1.5\times$ increase |
780 |
over the planar surface: (1) the nanoparticles have an increased |
781 |
population of undercoordinated surface atoms that are more efficient |
782 |
at transferring vibrational energy to the solvent, or (2) there is |
783 |
increased solid angle space available for phonon scattering into the |
784 |
solvent. The second of these explanations would have a significant |
785 |
dependence on the radius of spherical particles, although crystalline |
786 |
faceting of the particles may be enough to reduce this dependence on |
787 |
particle radius. These two possibilities do open up some interesting |
788 |
avenues for further exploration. |
789 |
|
790 |
One area that this method opens up that was not available for periodic |
791 |
simulation cells is direct computation of the solid/liquid rotational |
792 |
friction coefficients of nanostructures. The systems we have |
793 |
investigated (gold nanospheres and prolate ellipsoids) have analytic |
794 |
stick and approximate slip solutions provided via hydrodynamics, and |
795 |
our molecular simulations indicate that the bare metallic |
796 |
nanoparticles are closer to the stick than they are to slip |
797 |
conditions. This is markedly different behavior than we see for the |
798 |
solid / liquid friction at planar Au(111) interfaces, where the |
799 |
solvents experience a nearly-laminar flow over the surface in previous |
800 |
simulations. Schmidt and Skinner previously computed the behavior of |
801 |
spherical tag particles in molecular dynamics simulations, and showed |
802 |
that {\it slip} boundary conditions for translational resistance may |
803 |
be more appropriate for molecule-sized spheres embedded in a sea of |
804 |
spherical solvent particles.\cite{Schmidt:2004fj,Schmidt:2003kx} The |
805 |
size of the gold nanoparticles studied here and the structuring of the |
806 |
surfaces of the particles appears to bring the behavior closer to |
807 |
stick boundaries. Exactly where the slip-stick crossover takes place |
808 |
is also an interesting avenue for exploration. |
809 |
|
810 |
The ability to interrogate explicitly non-periodic effects |
811 |
(e.g. surface curvature, edge effects between facets, and the splay of |
812 |
ligand surface groups) on interfacial transport means that this method |
813 |
can be applied to systems of broad experimental and theoretical |
814 |
interest. |
815 |
|
816 |
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% |
817 |
% **ACKNOWLEDGMENTS** |
818 |
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% |
819 |
\begin{acknowledgement} |
820 |
The authors thank Dr. Shenyu Kuang for helpful discussions. Support |
821 |
for this project was provided by the National Science Foundation |
822 |
under grant CHE-0848243. Computational time was provided by the |
823 |
Center for Research Computing (CRC) at the University of Notre Dame. |
824 |
\end{acknowledgement} |
825 |
|
826 |
|
827 |
\newpage |
828 |
|
829 |
\bibliography{nonperiodicVSS} |
830 |
|
831 |
%\end{doublespace} |
832 |
\end{document} |