ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/trunk/nonperiodicVSS/nonperiodicVSS.tex
(Generate patch)

Comparing trunk/nonperiodicVSS/nonperiodicVSS.tex (file contents):
Revision 3978 by kstocke1, Tue Nov 26 20:33:01 2013 UTC vs.
Revision 4063 by gezelter, Thu Mar 13 15:44:27 2014 UTC

# Line 2 | Line 2
2   \setkeys{acs}{usetitle = true}
3  
4   \usepackage{caption}
5 \usepackage{endfloat}
5   \usepackage{geometry}
6   \usepackage{natbib}
7   \usepackage{setspace}
# Line 12 | Line 11
11   \usepackage{amssymb}
12   \usepackage{times}
13   \usepackage{mathptm}
15 \usepackage{setspace}
16 \usepackage{endfloat}
14   \usepackage{caption}
15   \usepackage{tabularx}
16   \usepackage{longtable}
17   \usepackage{graphicx}
21 \usepackage{multirow}
22 \usepackage{multicol}
18   \usepackage{achemso}
19   \usepackage[version=3]{mhchem}  % this is a great package for formatting chemical reactions
25 % \usepackage[square, comma, sort&compress]{natbib}
20   \usepackage{url}
27 \pagestyle{plain} \pagenumbering{arabic} \oddsidemargin 0.0cm
28 \evensidemargin 0.0cm \topmargin -21pt \headsep 10pt \textheight
29 9.0in \textwidth 6.5in \brokenpenalty=10000
21  
22 < % double space list of tables and figures
23 < % \AtBeginDelayedFloats{\renewcomand{\baselinestretch}{1.66}}
33 < \setlength{\abovecaptionskip}{20 pt}
34 < \setlength{\belowcaptionskip}{30 pt}
22 > \title{A method for creating thermal and angular momentum fluxes in
23 >  non-periodic simulations}
24  
36 % \bibpunct{}{}{,}{s}{}{;}
37
38 % \citestyle{nature}
39 % \bibliographystyle{achemso}
40
41 \title{A Method for Creating Thermal and Angular Momentum Fluxes in Non-Periodic Systems}
42
25   \author{Kelsey M. Stocker}
26   \author{J. Daniel Gezelter}
27   \email{gezelter@nd.edu}
# Line 47 | Line 29
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  
# Line 58 | Line 63
63   %       University of Notre Dame\\
64   %       Notre Dame, Indiana 46556}
65  
66 < \date{\today}
66 > %\date{\today}
67  
68 < \maketitle
68 > %\maketitle
69  
70 < \begin{doublespace}
70 > %\begin{doublespace}
71  
72   \begin{abstract}
73  
74 < We have adapted the Velocity Shearing and Scaling Reverse Non-Equilibium Molecular Dynamics (VSS-RNEMD) method for use with aperiodic system geometries. This new method is capable of creating stable temperature and angular velocity gradients in heterogeneous non-periodic systems while conserving total energy and angular momentum. To demonstrate the method, we have computed the thermal conductivities of a gold nanoparticle and water cluster, the shear viscosity of a water cluster, the interfacial thermal conductivity of a solvated gold nanoparticle and the interfacial friction of solvated gold nanostructures.
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  
# Line 79 | Line 99 | Non-equilibrium Molecular Dynamics (NEMD) methods impo
99   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
100   \section{Introduction}
101  
102 < Non-equilibrium Molecular Dynamics (NEMD) methods impose a temperature
102 > Non-equilibrium molecular dynamics (NEMD) methods impose a temperature
103   or velocity {\it gradient} on a
104 < system,\cite{ASHURST:1975tg,Evans:1982zk,ERPENBECK:1984sp,MAGINN:1993hc,Berthier:2002ij,Evans:2002ai,Schelling:2002dp,PhysRevA.34.1449,JiangHao_jp802942v}
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 flux to transport coefficients of bulk materials.  However,
107 < for heterogeneous systems, such as phase boundaries or interfaces, it
108 < is often unclear what shape of gradient should be imposed at the
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/VSS}
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
# Line 98 | Line 190 | Reverse Non-Equilibrium Molecular Dynamics (RNEMD) met
190   \label{fig:VSS}
191   \end{figure}
192  
101 Reverse Non-Equilibrium Molecular Dynamics (RNEMD) methods impose an
102 unphysical {\it flux} between different regions or ``slabs'' of the
103 simulation
104 box.\cite{MullerPlathe:1997xw,ISI:000080382700030,Kuang:2010uq} The
105 system responds by developing a temperature or velocity {\it gradient}
106 between the two regions.  The gradients which develop in response to
107 the applied flux are then related (via linear response theory) to the
108 transport coefficient of interest. Since the amount of the applied
109 flux is known exactly, and measurement of a gradient is generally less
110 complicated, imposed-flux methods typically take shorter simulation
111 times to obtain converged results. At interfaces, the observed
112 gradients often exhibit near-discontinuities at the boundaries between
113 dissimilar materials.  RNEMD methods do not need many trajectories to
114 provide information about transport properties, and they have become
115 widely used to compute thermal and mechanical transport in both
116 homogeneous liquids and
117 solids~\cite{MullerPlathe:1997xw,ISI:000080382700030,Maginn:2010} as
118 well as heterogeneous
119 interfaces.\cite{garde:nl2005,garde:PhysRevLett2009,kuang:AuThl}
120
121
122 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
123 % **METHODOLOGY**
124 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
125 \section{Velocity Shearing and Scaling (VSS) for non-periodic systems}
126 The VSS-RNEMD approach uses a series of simultaneous velocity shearing
127 and scaling exchanges between the two slabs.\cite{2012MolPh.110..691K}
128 This method imposes energy and momentum conservation constraints while
129 simultaneously creating a desired flux between the two slabs.  These
130 constraints ensure that all configurations are sampled from the same
131 microcanonical (NVE) ensemble.
132
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
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.
200 <
201 < We present here a new set of constraints that are more general than
202 < the VSS constraints.  For the non-periodic variant, the constraints
203 < fix both the total energy and total {\it angular} momentum of the
144 < system while simultaneously imposing a thermal and angular momentum
145 < flux between the two regions.
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 each $\Delta t$ time interval, the particle velocities
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$).
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}~~~~~~~~} & &
# Line 160 | Line 228 | Here $\langle\mathbf{\omega}_a\rangle$ and
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
232 < $\langle\mathbf{\omega}_b\rangle$ are the instantaneous angular
233 < velocities of each shell, and $\mathbf{r}_i$ is the position of
234 < particle $i$ relative to a fixed point in space (usually the center of
167 < mass of the cluster).  Particles in the shells also receive an
168 < additive ``angular shear'' to their velocities.  The amount of shear
169 < is governed by the imposed angular momentum flux,
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
# Line 174 | Line 239 | where $\overleftrightarrow{I}_{\{a,b\}}$ is the moment
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 tensor for
243 < each of the two shells.
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,
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
# Line 188 | Line 253 | Simultaneous solution of these quadratic formulae for
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
257 < coefficients, $a$ and $b$, will ensure that the simulation samples
258 < from the original microcanonical (NVE) ensemble.  Here $K_{\{a,b\}}$
259 < is the instantaneous translational kinetic energy of each shell.  At
260 < each time interval, we solve for $a$, $b$, $\mathbf{c}_a$, and
261 < $\mathbf{c}_b$, subject to the imposed angular momentum flux,
197 < $j_r(\mathbf{L})$, and thermal flux, $J_r$ values.  The new particle
198 < velocities are computed, and the simulation continues. System
199 < configurations after the transformations have exactly the same energy
200 < ({\it and} angular momentum) as before the moves.
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 quite simple to obtain both the thermal conductivity
268 < ($\lambda$) and shear viscosity ($\eta$),
269 < \begin{equation}
270 <  J_r = -\lambda \frac{\partial T}{\partial
210 <    r} \hspace{2in} j_r(\mathbf{L}_z) = -\eta \frac{\partial
211 <    \omega_z}{\partial r}
212 < \end{equation}
213 < of a liquid cluster.
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 < % NON-PERIODIC DYNAMICS
274 < %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
275 < \subsection{Dynamics for non-periodic systems}
272 > %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
273 > % **COMPUTATIONAL DETAILS**
274 > %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
275 > \section{Computational Details}
276  
277 < We have run all tests using the Langevin Hull nonperiodic simulation methodology.\cite{Vardeman2011} The Langevin Hull is especially useful for simulating heterogeneous mixtures of materials with different compressibilities, which are typically problematic for traditional affine transform methods. We have had success applying this method to
278 < several different systems including bare metal nanoparticles, liquid water clusters, and explicitly solvated nanoparticles. Calculated physical properties such as the isothermal compressibility of water and the bulk modulus of gold nanoparticles are in good agreement with previous theoretical and experimental results. The Langevin Hull uses a Delaunay tesselation to create a dynamic convex hull composed of triangular facets with vertices at atomic sites. Atomic sites included in the hull are coupled to an external bath defined by a temperature, pressure and viscosity. Atoms not included in the hull are subject to standard Newtonian dynamics. For these tests, thermal coupling to the bath was turned off to avoid interference with any imposed kinetic flux. Systems containing liquids were run under moderate pressure ($\sim$ 5 atm) to avoid the formation of a substantial vapor phase.
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 < % NON-PERIODIC RNEMD
288 > % FORCE FIELD PARAMETERS
289   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
290 < \subsection{VSS-RNEMD for non-periodic systems}
290 > \subsection{Force field}
291  
292 < The most useful RNEMD approach developed so far utilizes a series of
293 < simultaneous velocity shearing and scaling (VSS) exchanges between the two
294 < regions.\cite{Kuang2012} This method provides a set of conservation constraints
295 < while simultaneously creating a desired flux between the two regions. Satisfying
232 < the constraint equations ensures that the new configurations are sampled from the
233 < same NVE ensemble.
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 < We have implemented this new method in OpenMD, our group molecular dynamics code.\cite{openmd} The adaptation of VSS-RNEMD for non-periodic systems is relatively
298 < straightforward. The major modifications to the method are the addition of a rotational shearing term and the use of more versatile hot / cold regions instead
299 < of rectangular slabs. A temperature profile along the $r$ coordinate is created by recording the average temperature of concentric spherical shells.
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 < \begin{figure}
305 <        \center{\includegraphics[width=7in]{figures/npVSS2}}
306 <        \caption{Schematics of periodic (left) and nonperiodic (right) Velocity Shearing and Scaling RNEMD. A kinetic energy or momentum flux is applied from region B to region A. Thermal gradients are depicted by a color gradient. Linear or angular velocity gradients are shown as arrows.}
307 <        \label{fig:VSS}
308 < \end{figure}
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 < At each time interval, the particle velocities ($\mathbf{v}_i$ and
319 < $\mathbf{v}_j$) in the cold and hot shells ($C$ and $H$) are modified by a
320 < velocity scaling coefficient ($c$ and $h$), an additive linear velocity shearing
321 < term ($\mathbf{a}_c$ and $\mathbf{a}_h$), and an additive angular velocity
322 < shearing term ($\mathbf{b}_c$ and $\mathbf{b}_h$). Here, $\langle
323 < \mathbf{v}_{i,j} \rangle$ and $\langle \mathbf{\omega}_{i,j} \rangle$ are the
251 < average linear and angular velocities for each shell.
252 < \begin{displaymath}
253 < \begin{array}{rclcl}
254 < & \underline{~~~~~~~~~~\mathrm{scaling}~~~~~~~~~~} & &
255 < \underline{\mathrm{rotational \; shearing}} \\  \\
256 < \mathbf{v}_i $~~~$\leftarrow &
257 <  c \, \left(\mathbf{v}_i - \langle \omega_c
258 <  \rangle \times r_i\right) & + & \mathbf{b}_c \times r_i \\
259 < \mathbf{v}_j $~~~$\leftarrow &
260 <  h \, \left(\mathbf{v}_j - \langle \omega_h
261 <  \rangle \times r_j\right) & + & \mathbf{b}_h \times r_j
262 < \end{array}
263 < \end{displaymath}
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 < \begin{eqnarray}
326 < \mathbf{b}_c & = & - \mathbf{j}_r(\mathbf{L}) \; \Delta \, t \; I_c^{-1} + \langle \omega_c \rangle \label{eq:bc}\\
327 < \mathbf{b}_h & = & + \mathbf{j}_r(\mathbf{L}) \; \Delta \, t \; I_h^{-1} + \langle \omega_h \rangle \label{eq:bh}
328 < \end{eqnarray}
329 <
330 < The total energy is constrained via two quadratic formulae,
331 < \begin{eqnarray}
332 < K_c - J_r \; \Delta \, t & = & c^2 \, (K_c - \frac{1}{2}\langle \omega_c \rangle \cdot I_c\langle \omega_c \rangle) + \frac{1}{2} \mathbf{b}_c \cdot I_c \, \mathbf{b}_c \label{eq:Kc}\\
333 < K_h + J_r \; \Delta \, t & = & h^2 \, (K_h - \frac{1}{2}\langle \omega_h \rangle \cdot I_h\langle \omega_h \rangle) + \frac{1}{2} \mathbf{b}_h \cdot I_h \, \mathbf{b}_h \label{eq:Kh}
334 < \end{eqnarray}
335 <
336 < the simultaneous
337 < solution of which provide the velocity scaling coefficients $c$ and $h$. Given an
338 < imposed angular momentum flux, $\mathbf{j}_{r} \left( \mathbf{L} \right)$, and/or
339 < thermal flux, $J_r$, equations \ref{eq:bc} - \ref{eq:Kh} are sufficient to obtain
280 < the velocity scaling ($c$ and $h$) and shearing ($\mathbf{b}_c,\,$ and $\mathbf{b}_h$) at each time interval.
281 <
282 < %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
283 < % **COMPUTATIONAL DETAILS**
284 < %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
285 < \section{Computational Details}
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 < %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
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 < %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
356 < % FORCE FIELD PARAMETERS
357 < %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
358 < \subsection{Force field parameters}
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 < We have chosen the SPC/E water model for these simulations, which is particularly useful for method validation as there are many values for physical properties from previous simulations available for direct comparison.\cite{Bedrov:2000, Kuang2010}
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 < Gold -- gold interactions are described by the quantum Sutton-Chen (QSC) model.\cite{PhysRevB.59.3527} The QSC parameters are tuned to experimental properties such as density, cohesive energy, and elastic moduli and include zero-point quantum corrections.
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  
302 Hexane molecules are described by the TraPPE united atom model,\cite{TraPPE-UA.alkanes} which provides good computational efficiency and reasonable accuracy for bulk thermal conductivity values. In this model,
303 sites are located at the carbon centers for alkyl groups. Bonding
304 interactions, including bond stretches and bends and torsions, were
305 used for intra-molecular sites closer than 3 bonds. For non-bonded
306 interactions, Lennard-Jones potentials were used.  We have previously
307 utilized both united atom (UA) and all-atom (AA) force fields for
308 thermal conductivity,\cite{kuang:AuThl,Kuang2012} and since the united
309 atom force fields cannot populate the high-frequency modes that are
310 present in AA force fields, they appear to work better for modeling
311 thermal conductivity.
312
313 Gold -- hexane nonbonded interactions are governed by pairwise Lennard-Jones parameters derived from Vlugt \emph{et al}.\cite{vlugt:cpc2007154} They fitted parameters for the interaction between Au and CH$_{\emph x}$ pseudo-atoms based on the effective potential of Hautman and Klein for the Au(111) surface.\cite{hautman:4994}
314
386   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
387   % THERMAL CONDUCTIVITIES
388   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
389   \subsection{Thermal conductivities}
390  
391 < The thermal conductivity, $\lambda$, can be calculated using the imposed kinetic energy flux, $J_r$, and thermal gradient, $\frac{\partial T}{\partial r}$ measured from the resulting temperature profile
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 <        J_r = -\lambda \frac{\partial T}{\partial r}
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 < A thermal flux is created using VSS-RNEMD moves, and the resulting temperature
431 < profiles are analyzed to yield information about the interfacial thermal
432 < conductance. As the simulation progresses, the VSS moves are performed on a regular basis, and
433 < the system develops a thermal or velocity gradient in response to the applied
434 < flux. A definition of the interfacial conductance in terms of a temperature difference ($\Delta T$) measured at $r_0$,
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 <        G = \frac{J_r}{\Delta T_{r_0}}, \label{eq:G}
445 >        R_K = \frac{1}{q_r} \Delta T 4 \pi r^2
446 > \label{eq:RK}
447   \end{equation}
339 is useful once the RNEMD approach has generated a
340 stable temperature gap across the interface.
448  
449 < %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
343 < % INTERFACIAL FRICTION
344 < %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
345 < \subsection{Interfacial friction}
449 > making the total resistance of two neighboring shells
450  
347 The slip length, $\delta$, is defined as the ratio of the interfacial friction coefficient, $\kappa$, and dynamic solvent viscosity, $\eta$
348
451   \begin{equation}
452 <        \delta = \frac{\eta}{\kappa}
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 < and is measured from a radial angular momentum profile generated from a nonperiodic VSS-RNEMD simulation.
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 < Analytical solutions for the rotational friction coefficients for a solvated spherical body of radius $r$ under ``stick'' boundary conditions can be estimated using Stokes' law
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 = 8 \pi \eta r^3 \label{eq:Xi}.
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 and was determined for TraPPE-UA hexane under these condition by applying a traditional VSS-RNEMD linear momentum flux to a periodic box of solvent.
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 solutions for symmetric prolate $(a \geq b = c)$ and oblate $(a < b = c)$ ellipsoids. For simplicity, we define a Perrin Factor, $S$,
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]. \label{eq:S}
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}
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}.                                      \label{eq:Xibc}
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 < % However, the friction between hexane solvent and gold more likely falls within ``slip'' boundary conditions. Hu and Zwanzig\cite{Zwanzig} investigated the rotational friction coefficients for spheroids under slip boundary conditions and obtained numerial results for a scaling factor as a function of $\tau$, the ratio of the shorter semiaxes and the longer semiaxis of the spheroid. For the sphere and prolate ellipsoid shown here, the values of $\tau$ are $1$ and $0.3939$, respectively. According to the values tabulated by Hu and Zwanzig, the sphere friction coefficient approaches $0$, while the ellipsoidal friction coefficient must be scaled by a factor of $0.880$ to account for the reduced interfacial friction under ``slip'' boundary conditions.
506 > corresponding to rotation about the long axis ($a$), and each of the equivalent short axes ($b$ and $c$), respectively.
507  
508 < Another useful quantity is the friction factor $f$, or the friction coefficient of a non-spherical shape divided by the friction coefficient of a sphere of equivalent volume. The nanoparticle and ellipsoidal nanorod dimensions used here were specifically chosen to create structures with equal volumes.
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   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
# Line 389 | Line 542 | Calculated values for the thermal conductivity of a 40
542   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
543   \subsection{Thermal conductivities}
544  
545 < Calculated values for the thermal conductivity of a 40 \AA$~$ radius gold nanoparticle (15707 atoms) at different kinetic energy flux values are shown in Table \ref{table:goldconductivity}. For these calculations, the hot and cold slabs were excluded from the linear regression of the thermal gradient. The measured linear slope $\langle dT / dr \rangle$ is linearly dependent on the applied kinetic energy flux $J_r$. Calculated thermal conductivity values compare well with previous bulk QSC values of 1.08 -- 1.26 W m$^{-1}$ K$^{-1}$\cite{Kuang2010}, though still significantly lower than the experimental value of 320 W m$^{-1}$ K$^{-1}$, as the QSC force field neglects significant electronic contributions to heat conduction. The small increase relative to previous simulated bulk values is due to a slight increase in gold density -- as expected, an increase in density results in higher thermal conductivity values. The increased density is a result of nanoparticle curvature relative to an infinite bulk slab, which introduces surface tension that increases ambient density.
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 dT / dr \rangle$} & {$\boldsymbol \lambda$}\\
553 < {\small(kcal fs$^{-1}$ \AA$^{-2}$)} & {\small(K \AA$^{-1}$)} & {\small(W m$^{-1}$ K$^{-1}$)}\\ \hline
554 < 3.25$\times 10^{-6}$ & 0.11435 & 1.9753 \\
555 < 6.50$\times 10^{-6}$ & 0.2324 & 1.9438 \\
556 < 1.30$\times 10^{-5}$ & 0.44922 & 2.0113 \\
557 < 3.25$\times 10^{-5}$ & 1.1802 & 1.9139 \\
558 < 6.50$\times 10^{-5}$ & 2.339 & 1.9314
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:goldconductivity}
562 > \label{table:goldTC}
563   \end{longtable}
564  
565 < Calculated values for the thermal conductivity of a cluster of 6912 SPC/E water molecules are shown in Table \ref{table:waterconductivity}.
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 ambient density.}
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 dT / dr \rangle$} & {$\boldsymbol \lambda$}\\
578 < {\small(kcal fs$^{-1}$ \AA$^{-2}$)} & {\small(K \AA$^{-1}$)} & {\small(W m$^{-1}$ K$^{-1}$)}\\ \hline
579 < 1.00$\times 10^{-5}$ & 0.38683 & 1.79665486\\
580 < 3.00$\times 10^{-5}$ & 1.1643 & 1.79077557\\
581 < 6.00$\times 10^{-5}$ & 2.2262 & 1.87314707\\
582 < \hline \hline
583 < \label{table:waterconductivity}
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{Caption.}
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} & $J_r$ & {G}\\
608 < {\small(\AA)} & {\small(kcal fs$^{-1}$ \AA$^{-2}$)}  & {\small(W m$^{-2}$ K$^{-1}$)}\\ \hline
609 < 20 & & 59.66 \\
610 < 30 & & 57.88 \\
611 < 40 & & 37.48 \\
612 < $\infty$ & & 30.2 \\
613 < \hline \hline
614 < \label{table:interfacialconductance}
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:interfacialfrictionstick} shows the calculated interfacial friction coefficients $\kappa$ for spherical gold nanoparticles and a prolate ellipsoidal gold nanorod in TraPPE-UA hexane. An angular momentum flux was applied between the A and B regions defined as the gold structure and hexane molecules beyond a certain radius, respectively. The resulting angular velocity gradient resulted in the gold structure rotating about the prescribed axis within the solvent.
629 <        
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{Calculated ``stick'' interfacial friction coefficients ($\kappa$) and friction factors ($f$) of gold nanostructures solvated in TraPPE-UA hexane at 230 K. The ellipsoid is oriented with the long axis along the $z$ direction.}
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} & {$\kappa_{VSS}$} & {$\kappa_{calc}$} & {$f_{VSS}$} & {$f_{calc}$}\\
637 < {} & {} & {\small($10^{-29}$ Pa s m$^{3}$)} & {\small($10^{-29}$ Pa s m$^{3}$)} & {} & {}\\  \hline
638 < {Sphere (r = 20 \AA)} & {$x = y = z$} & {} & {6.13239} & {1} & {1}\\
639 < {Sphere (r = 30 \AA)} & {$x = y = z$} & {} & {20.6968} & {1} & {1}\\
640 < {Sphere (r = 40 \AA)} & {$x = y = z$} & {} & {49.0591} & {1} & {1}\\
641 < {Prolate Ellipsoid} & {$x = y$} & {} & {8.22846} & {} & {1.92477}\\
642 < {Prolate Ellipsoid} & {$z$} & {} & {3.28634} & {} & {0.768726}\\
643 <  \hline \hline
644 < \label{table:interfacialfrictionstick}
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 < % \begin{longtable}{lccc}
648 < % \caption{Calculated ``slip'' interfacial friction coefficients ($\kappa$) and friction factors ($f$) of gold nanostructures solvated in TraPPE-UA hexane. The ellipsoid is oriented with the long axis along the $z$ direction.}
649 < % \\ \hline \hline
650 < % {Structure} & {Axis of rotation} & {$\kappa_{VSS}$} & {$\kappa_{calc}$}\\
651 < % {} & {} & {\small($10^{-29}$ Pa s m$^{3}$)} & {\small($10^{-29}$ Pa s m$^{3}$)}\\  \hline
652 < % {Sphere} & {$x = y = z$} & {} & {0}\\
653 < % {Prolate Ellipsoid} & {$x = y$} & {} & {0}\\
654 < % {Prolate Ellipsoid} & {$z$} & {} & {7.9295392}\\  \hline \hline
655 < % \label{table:interfacialfrictionslip}
656 < % \end{longtable}
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 < \section*{Acknowledgments}
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  
483 We gratefully acknowledge conversations with Dr. Shenyu Kuang. Support for
484 this project was provided by the National Science Foundation under grant
485 CHE-0848243. Computational time was provided by the Center for Research
486 Computing (CRC) at the University of Notre Dame.
678  
679   \newpage
680  
681   \bibliography{nonperiodicVSS}
682  
683 < \end{doublespace}
683 > %\end{doublespace}
684   \end{document}

Diff Legend

Removed lines
+ Added lines
< Changed lines
> Changed lines