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