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