| 2 |
|
\setkeys{acs}{usetitle = true} |
| 3 |
|
|
| 4 |
|
\usepackage{caption} |
| 5 |
– |
\usepackage{float} |
| 5 |
|
\usepackage{geometry} |
| 6 |
|
\usepackage{natbib} |
| 7 |
|
\usepackage{setspace} |
| 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} |
| 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 |
|
|
| 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 |
|
|
| 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{Methodology} |
| 173 |
> |
\section{Velocity shearing and scaling (VSS) for non-periodic systems} |
| 174 |
|
|
| 175 |
< |
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% |
| 176 |
< |
% FORCE FIELD PARAMETERS |
| 177 |
< |
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% |
| 178 |
< |
\subsection{Force field parameters} |
| 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 |
< |
We have chosen the SPC/E water model for these simulations. There are many values for physical properties from previous simulations available for direct comparison. |
| 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 |
< |
Gold-gold interactions are described by the quantum Sutton-Chen (QSC) model. The QSC parameters are tuned to experimental properties such as density, cohesive energy, and elastic moduli and include zero-point quantum corrections. |
| 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 |
< |
Hexane molecules are described by the TraPPE united atom model. This model provides good computational efficiency and reasonable accuracy for bulk thermal conductivity values. |
| 206 |
< |
|
| 207 |
< |
Metal-nonmetal interactions are governed by parameters derived from Luedtke and Landman. |
| 208 |
< |
|
| 209 |
< |
|
| 210 |
< |
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% |
| 211 |
< |
% NON-PERIODIC DYNAMICS |
| 212 |
< |
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% |
| 213 |
< |
\subsection{Dynamics for non-periodic systems} |
| 214 |
< |
|
| 215 |
< |
We have run all tests using the Langevin Hull nonperiodic simulation methodology. 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. 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, which would have a hull created out of a relatively small number of molecules. |
| 108 |
< |
|
| 109 |
< |
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% |
| 110 |
< |
% NON-PERIODIC RNEMD |
| 111 |
< |
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% |
| 112 |
< |
\subsection{VSS-RNEMD for non-periodic systems} |
| 113 |
< |
|
| 114 |
< |
The adaptation of VSS-RNEMD for non-periodic systems is relatively |
| 115 |
< |
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 |
| 116 |
< |
of rectangular slabs. A temperature profile along the $r$ coordinate is created by recording the average temperature of concentric spherical shells. |
| 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 |
< |
At each time interval, the particle velocities ($\mathbf{v}_i$ and |
| 218 |
< |
$\mathbf{v}_j$) in the cold and hot shells ($C$ and $H$) are modified by a |
| 120 |
< |
velocity scaling coefficient ($c$ and $h$), an additive linear velocity shearing |
| 121 |
< |
term ($\mathbf{a}_c$ and $\mathbf{a}_h$), and an additive angular velocity |
| 122 |
< |
shearing term ($\mathbf{b}_c$ and $\mathbf{b}_h$). Here, $\langle |
| 123 |
< |
\mathbf{v}_{i,j} \rangle$ and $\langle \mathbf{\omega}_{i,j} \rangle$ are the |
| 124 |
< |
average linear and angular velocities for each shell. |
| 125 |
< |
|
| 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}} \\ \\ |
| 221 |
> |
& \underline{~~~~~~~~\mathrm{scaling}~~~~~~~~} & & |
| 222 |
> |
\underline{\mathrm{rotational~shearing}} \\ \\ |
| 223 |
|
\mathbf{v}_i $~~~$\leftarrow & |
| 224 |
< |
c \, \left(\mathbf{v}_i - \langle \omega_c |
| 225 |
< |
\rangle \times r_i\right) & + & \mathbf{b}_c \times r_i \\ |
| 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 |
< |
h \, \left(\mathbf{v}_j - \langle \omega_h |
| 228 |
< |
\rangle \times r_j\right) & + & \mathbf{b}_h \times r_j |
| 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 |
< |
|
| 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{b}_c & = & - \mathbf{j}_r(\mathbf{L}) \; \Delta \, t \; I_c^{-1} + \langle \omega_c \rangle \label{eq:bc}\\ |
| 238 |
< |
\mathbf{b}_h & = & + \mathbf{j}_r(\mathbf{L}) \; \Delta \, t \; I_h^{-1} + \langle \omega_h \rangle \label{eq:bh} |
| 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 |
< |
The total energy is constrained via two quadratic formulae, |
| 246 |
< |
|
| 245 |
> |
To simultaneously impose a thermal flux ($J_r$) between the shells we |
| 246 |
> |
use energy conservation constraints, |
| 247 |
|
\begin{eqnarray} |
| 248 |
< |
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}\\ |
| 249 |
< |
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} |
| 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 |
< |
the simultaneous |
| 264 |
< |
solution of which provide the velocity scaling coefficients $c$ and $h$. Given an |
| 265 |
< |
imposed angular momentum flux, $\mathbf{j}_{r} \left( \mathbf{L} \right)$, and/or |
| 266 |
< |
thermal flux, $J_r$, equations \ref{eq:bc} - \ref{eq:Kh} are sufficient to obtain |
| 267 |
< |
the velocity scaling ($c$ and $h$) and shearing ($\mathbf{b}_c,\,$ and $\mathbf{b}_h$) at each time interval. |
| 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 |
< |
% **TESTS AND APPLICATIONS** |
| 273 |
> |
% **COMPUTATIONAL DETAILS** |
| 274 |
|
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% |
| 275 |
< |
\section{Tests and Applications} |
| 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 |
< |
% THERMAL CONDUCTIVITIES |
| 288 |
> |
% FORCE FIELD PARAMETERS |
| 289 |
|
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% |
| 290 |
< |
\subsection{Thermal conductivities} |
| 290 |
> |
\subsection{Force field} |
| 291 |
|
|
| 292 |
< |
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 NIVS paper], 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. |
| 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 |
< |
\begin{longtable}{p{2.7cm} p{2.5cm} p{2.5cm}} |
| 298 |
< |
\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.} |
| 299 |
< |
\\ \hline \hline |
| 300 |
< |
\centering {$J_r$} & \centering\arraybackslash {$\langle dT / dr \rangle$} & \centering\arraybackslash {$\boldsymbol \lambda$}\\ |
| 301 |
< |
\centering {\small(kcal fs$^{-1}$ \AA$^{-2}$)} & \centering\arraybackslash {\small(K \AA$^{-1}$)} & \centering\arraybackslash {\small(W m$^{-1}$ K$^{-1}$)}\\ \hline |
| 302 |
< |
\endhead |
| 175 |
< |
\hline |
| 176 |
< |
% \endfoot |
| 177 |
< |
\centering3.25$\times 10^{-6}$ & \centering\arraybackslash 0.11435 & \centering\arraybackslash 1.9753 \\ |
| 178 |
< |
\centering6.50$\times 10^{-6}$ & \centering\arraybackslash 0.2324 & \centering\arraybackslash 1.9438 \\ |
| 179 |
< |
\centering1.30$\times 10^{-5}$ & \centering\arraybackslash 0.44922 & \centering\arraybackslash 2.0113 \\ |
| 180 |
< |
\centering3.25$\times 10^{-5}$ & \centering\arraybackslash 1.1802 & \centering\arraybackslash 1.9139 \\ |
| 181 |
< |
\centering6.50$\times 10^{-5}$ & \centering\arraybackslash 2.339 & \centering\arraybackslash 1.9314 |
| 182 |
< |
\\ \hline \hline |
| 183 |
< |
\label{table:goldconductivity} |
| 184 |
< |
\end{longtable} |
| 185 |
< |
|
| 186 |
< |
SPC/E Water Cluster: |
| 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{longtable}{p{2.7cm} p{2.5cm} p{2.5cm}} |
| 305 |
< |
\caption{Calculated thermal conductivity of a cluster of 6912 SPC/E water molecules. Calculations were performed at 300 K and ambient density.} |
| 306 |
< |
\\ \hline \hline |
| 307 |
< |
\centering {$J_r$} & \centering\arraybackslash {$\langle dT / dr \rangle$} & \centering\arraybackslash {$\boldsymbol \lambda$}\\ |
| 308 |
< |
\centering {\small(kcal fs$^{-1}$ \AA$^{-2}$)} & \centering\arraybackslash {\small(K \AA$^{-1}$)} & \centering\arraybackslash {\small(W m$^{-1}$ K$^{-1}$)}\\ \hline |
| 309 |
< |
\endhead |
| 310 |
< |
\hline |
| 311 |
< |
% \endfoot |
| 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 |
< |
\\ \hline \hline |
| 319 |
< |
\label{table:waterconductivity} |
| 320 |
< |
\end{longtable} |
| 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 |
< |
% SHEAR VISCOSITY |
| 326 |
> |
% NON-PERIODIC DYNAMICS |
| 327 |
|
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% |
| 328 |
< |
\subsection{Shear viscosity} |
| 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 |
< |
SPC/E Water Cluster: |
| 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 |
< |
Gold Nanoparticle in Hexane: |
| 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 FRICTION |
| 461 |
> |
% INTERFACIAL ROTATIONAL FRICTION |
| 462 |
|
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% |
| 463 |
< |
\subsection{Interfacial friction} |
| 463 |
> |
\subsection{Interfacial rotational friction} |
| 464 |
|
|
| 465 |
< |
Table \ref{table:interfacialfriction} gives the calculated interfacial friction coefficients $\kappa$ for a spherical gold nanoparticle and 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 included in the convex hull, respectively. The resulting angular velocity gradient resulted in the gold structure rotating about the prescribed axis within the solvent. |
| 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 |
< |
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 |
| 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 |
< |
f_r = 8 \pi \eta r^3 \label{eq:fr}. |
| 479 |
> |
\Xi^{rr}_{stick} = 8 \pi \eta r^3 |
| 480 |
> |
\label{eq:Xisphere}. |
| 481 |
|
\end{equation} |
| 482 |
|
|
| 483 |
< |
For general ellipsoids with semiaxes $a$, $b$, and $c$, Perrin's extension of Stokes' law provides exact solutions for prolate $(a \geq b = c)$ and oblate $(a < b = c)$ ellipsoids. For a prolate ellipsoidal rod, demonstrated here, |
| 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 |
< |
\begin{eqnarray} |
| 488 |
< |
f_a = \label{eq:fa}\\ |
| 489 |
< |
f_b = f_c = \label{eq:fb} |
| 233 |
< |
\end{eqnarray} |
| 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 |
< |
The dynamic viscosity of the solvent, $\eta$, was calculated by applying a linear momentum flux to a periodic box of TraPPE-UA hexane. |
| 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 |
< |
\begin{longtable}{p{3.8cm} p{3cm} p{2.8cm} p{2.5cm} p{2.5cm}} |
| 497 |
< |
\caption{Calculated interfacial friction coefficients ($\kappa$) and slip length ($\delta$) of gold nanostructures solvated in TraPPE-UA hexane. The ellipsoid is oriented with the long axis along the $z$ direction.} |
| 498 |
< |
\\ \hline \hline |
| 499 |
< |
{Structure} & \centering{Axis of rotation} & \centering\arraybackslash {$\kappa$} & \centering\arraybackslash {$\delta$} & \centering\arraybackslash Stokes' Law $F$\\ |
| 500 |
< |
\centering {} & {} & \centering\arraybackslash {\small($10^4$ Pa s m$^{-1}$)} & \centering\arraybackslash {\small(nm)} & \centering\arraybackslash{\small()}\\ \hline |
| 501 |
< |
\endhead |
| 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 |
< |
% \endfoot |
| 245 |
< |
Nanoparticle & \centering$x = y = z$ & & & \\ |
| 246 |
< |
Prolate Ellipsoidal rod & \centering$x = y$ & & & \\ |
| 247 |
< |
Prolate Ellipsoidal rod & \centering$z$ & & & |
| 560 |
> |
This work & & 1.0040 |
| 561 |
|
\\ \hline \hline |
| 562 |
< |
\label{table:interfacialfriction} |
| 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 |
< |
\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 |
|
|
| 264 |
– |
We gratefully acknowledge conversations with Dr. Shenyu Kuang. Support for |
| 265 |
– |
this project was provided by the National Science Foundation under grant |
| 266 |
– |
CHE-0848243. Computational time was provided by the Center for Research |
| 267 |
– |
Computing (CRC) at the University of Notre Dame. |
| 678 |
|
|
| 679 |
|
\newpage |
| 680 |
|
|
| 681 |
|
\bibliography{nonperiodicVSS} |
| 682 |
|
|
| 683 |
< |
\end{doublespace} |
| 684 |
< |
\end{document} |
| 683 |
> |
%\end{doublespace} |
| 684 |
> |
\end{document} |
