| 2 |
|
\setkeys{acs}{usetitle = true} |
| 3 |
|
|
| 4 |
|
\usepackage{caption} |
| 5 |
– |
% \usepackage{endfloat} |
| 5 |
|
\usepackage{geometry} |
| 6 |
|
\usepackage{natbib} |
| 7 |
|
\usepackage{setspace} |
| 15 |
|
\usepackage{tabularx} |
| 16 |
|
\usepackage{longtable} |
| 17 |
|
\usepackage{graphicx} |
| 19 |
– |
\usepackage{multirow} |
| 20 |
– |
\usepackage{multicol} |
| 18 |
|
\usepackage{achemso} |
| 19 |
|
\usepackage[version=3]{mhchem} % this is a great package for formatting chemical reactions |
| 23 |
– |
% \usepackage[square, comma, sort&compress]{natbib} |
| 20 |
|
\usepackage{url} |
| 25 |
– |
\pagestyle{plain} \pagenumbering{arabic} \oddsidemargin 0.0cm |
| 26 |
– |
\evensidemargin 0.0cm \topmargin -21pt \headsep 10pt \textheight |
| 27 |
– |
9.0in \textwidth 6.5in \brokenpenalty=10000 |
| 21 |
|
|
| 22 |
< |
% double space list of tables and figures |
| 23 |
< |
% \AtBeginDelayedFloats{\renewcomand{\baselinestretch}{1.66}} |
| 31 |
< |
\setlength{\abovecaptionskip}{20 pt} |
| 32 |
< |
\setlength{\belowcaptionskip}{30 pt} |
| 22 |
> |
\title{A method for creating thermal and angular momentum fluxes in |
| 23 |
> |
non-periodic simulations} |
| 24 |
|
|
| 34 |
– |
% \bibpunct{}{}{,}{s}{}{;} |
| 35 |
– |
|
| 36 |
– |
% \citestyle{nature} |
| 37 |
– |
% \bibliographystyle{achemso} |
| 38 |
– |
|
| 39 |
– |
\title{A Method for Creating Thermal and Angular Momentum Fluxes in Non-Periodic Systems} |
| 40 |
– |
|
| 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 |
| 75 |
< |
for use with non-periodic system geometries. This new method is capable of creating stable temperature and |
| 76 |
< |
angular velocity gradients in heterogeneous non-periodic systems while conserving total energy and angular |
| 77 |
< |
momentum. To demonstrate the method, we have computed the thermal conductivities of a gold nanoparticle and |
| 78 |
< |
water cluster, the shear viscosity of a water cluster, the interfacial thermal conductivity of a solvated gold |
| 79 |
< |
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 or velocity {\it gradient} on a |
| 103 |
< |
system,\cite{ASHURST:1975tg,Evans:1982zk,ERPENBECK:1984sp,MAGINN:1993hc,Berthier:2002ij,Evans:2002ai,Schelling:2 |
| 104 |
< |
002dp,PhysRevA.34.1449,JiangHao_jp802942v} and use linear response theory to connect the resulting thermal or |
| 105 |
< |
momentum flux to transport coefficients of bulk materials. However, for heterogeneous systems, such as phase |
| 106 |
< |
boundaries or interfaces, it is often unclear what shape of gradient should be imposed at the boundary between |
| 107 |
< |
materials. |
| 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 |
< |
Reverse Non-Equilibrium Molecular Dynamics (RNEMD) methods impose an unphysical {\it flux} between different |
| 123 |
< |
regions or ``slabs'' of the simulation box.\cite{MullerPlathe:1997xw,ISI:000080382700030,Kuang2010} The system |
| 124 |
< |
responds by developing a temperature or velocity {\it gradient} between the two regions. The gradients which |
| 125 |
< |
develop in response to the applied flux are then related (via linear response theory) to the transport |
| 126 |
< |
coefficient of interest. Since the amount of the applied flux is known exactly, and measurement of a gradient |
| 127 |
< |
is generally less complicated, imposed-flux methods typically take shorter simulation times to obtain converged |
| 128 |
< |
results. At interfaces, the observed gradients often exhibit near-discontinuities at the boundaries between |
| 129 |
< |
dissimilar materials. RNEMD methods do not need many trajectories to provide information about transport |
| 130 |
< |
properties, and they have become widely used to compute thermal and mechanical transport in both homogeneous |
| 131 |
< |
liquids and solids~\cite{MullerPlathe:1997xw,ISI:000080382700030,Maginn:2010} as well as heterogeneous |
| 132 |
< |
interfaces.\cite{garde:nl2005,garde:PhysRevLett2009,kuang:AuThl} |
| 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} |
| 173 |
> |
\section{Velocity shearing and scaling (VSS) for non-periodic systems} |
| 174 |
|
|
| 175 |
< |
The periodic VSS-RNEMD approach uses a series of simultaneous velocity shearing and scaling exchanges between the two |
| 176 |
< |
slabs.\cite{Kuang2012} This method imposes energy and momentum conservation constraints while simultaneously |
| 177 |
< |
creating a desired flux between the two slabs. These constraints ensure that all configurations are sampled |
| 178 |
< |
from the same microcanonical (NVE) ensemble. |
| 175 |
> |
The original periodic VSS-RNEMD approach uses a series of simultaneous |
| 176 |
> |
velocity shearing and scaling exchanges between the two |
| 177 |
> |
slabs.\cite{Kuang:2012fe} This method imposes energy and linear |
| 178 |
> |
momentum conservation constraints while simultaneously creating a |
| 179 |
> |
desired flux between the two slabs. These constraints ensure that all |
| 180 |
> |
configurations are sampled from the same microcanonical (NVE) |
| 181 |
> |
ensemble. |
| 182 |
|
|
| 183 |
|
\begin{figure} |
| 184 |
|
\includegraphics[width=\linewidth]{figures/npVSS} |
| 190 |
|
\label{fig:VSS} |
| 191 |
|
\end{figure} |
| 192 |
|
|
| 193 |
< |
We have extended the VSS method for use in {\it non-periodic} simulations, in which the ``slabs'' have been |
| 194 |
< |
generalized to two separated regions of space. These regions could be defined as concentric spheres (as in |
| 195 |
< |
figure \ref{fig:VSS}), or one of the regions can be defined in terms of a dynamically changing ``hull'' |
| 196 |
< |
comprising the surface atoms of the cluster. This latter definition is identical to the hull used in the |
| 197 |
< |
Langevin Hull algorithm. |
| 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 |
< |
We present here a new set of constraints that are more general than the VSS constraints. For the non-periodic |
| 206 |
< |
variant, the constraints fix both the total energy and total {\it angular} momentum of the system while |
| 207 |
< |
simultaneously imposing a thermal and angular momentum flux between the two regions. |
| 208 |
< |
|
| 209 |
< |
After each $\Delta t$ time interval, the particle velocities ($\mathbf{v}_i$ and $\mathbf{v}_j$) in the two |
| 210 |
< |
shells ($A$ and $B$) are modified by a velocity scaling coefficient ($a$ and $b$) and by a rotational shearing |
| 211 |
< |
term ($\mathbf{c}_a$ and $\mathbf{c}_b$). |
| 205 |
> |
After a time interval of $\Delta t$, the particle velocities |
| 206 |
> |
($\mathbf{v}_i$ and $\mathbf{v}_j$) in the two shells ($A$ and $B$) |
| 207 |
> |
are modified by a velocity scaling coefficient ($a$ and $b$) and by a |
| 208 |
> |
rotational shearing term ($\mathbf{c}_a$ and $\mathbf{c}_b$). The |
| 209 |
> |
scalars $a$ and $b$ collectively provide a thermal exchange between |
| 210 |
> |
the two regions. One of the values is larger than 1, and the other |
| 211 |
> |
smaller. To conserve total energy and angular momentum, the values of |
| 212 |
> |
these two scalars are coupled. The vectors ($\mathbf{c}_a$ and |
| 213 |
> |
$\mathbf{c}_b$) provide a relative rotational shear to the velocities |
| 214 |
> |
of the particles within the two regions, and these vectors must also |
| 215 |
> |
be coupled to constrain the total angular momentum. |
| 216 |
|
|
| 217 |
+ |
Once the values of the scaling and shearing factors are known, the |
| 218 |
+ |
velocity changes are applied, |
| 219 |
|
\begin{displaymath} |
| 220 |
|
\begin{array}{rclcl} |
| 221 |
|
& \underline{~~~~~~~~\mathrm{scaling}~~~~~~~~} & & |
| 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 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 use energy conservation constraints, |
| 245 |
> |
To simultaneously impose a thermal flux ($J_r$) between the shells we |
| 246 |
> |
use energy conservation constraints, |
| 247 |
|
\begin{eqnarray} |
| 248 |
|
K_a - J_r \Delta t & = & a^2 \left(K_a - \frac{1}{2}\langle |
| 249 |
|
\omega_a \rangle \cdot \overleftrightarrow{I_a} \cdot \langle \omega_a |
| 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 performed on a regular basis, and the system |
| 264 |
< |
will develop a temperature and/or angular velocity gradient in response to the applied flux. Using the slope of |
| 265 |
< |
the radial temperature or velocity gradients, it is quite simple to obtain both the thermal conductivity |
| 266 |
< |
($\lambda$), interfacial thermal conductance ($G$), or rotational friction coefficients ($\Xi^{rr}$) of any |
| 267 |
< |
non-periodic system. |
| 263 |
> |
As the simulation progresses, the velocity transformations can be |
| 264 |
> |
performed on a regular basis, and the system will develop a |
| 265 |
> |
temperature and/or angular velocity gradient in response to the |
| 266 |
> |
applied flux. Using the slope of the radial temperature or velocity |
| 267 |
> |
gradients, it is straightforward to obtain both the thermal |
| 268 |
> |
conductivity ($\lambda$), interfacial thermal conductance ($G$), or |
| 269 |
> |
rotational friction coefficients ($\Xi^{rr}$) of any non-periodic |
| 270 |
> |
system. |
| 271 |
|
|
| 272 |
|
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% |
| 273 |
|
% **COMPUTATIONAL DETAILS** |
| 274 |
|
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% |
| 275 |
|
\section{Computational Details} |
| 276 |
|
|
| 277 |
< |
The new VSS-RNEMD methodology for non-periodic system geometries has been implemented in our group molecular |
| 278 |
< |
dynamics code, OpenMD.\cite{openmd} We have used the new method to calculate the thermal conductance of a gold |
| 279 |
< |
nanoparticle and SPC/E water cluster, and compared the results with previous bulk RNEMD values, as well as |
| 280 |
< |
experiment. We have also investigated the interfacial thermal conductance and interfacial rotational friction |
| 281 |
< |
for gold nanostructures solvated in hexane as a function of nanoparticle size and shape. |
| 277 |
> |
The new VSS-RNEMD methodology for non-periodic system geometries has |
| 278 |
> |
been implemented in our group molecular dynamics code, |
| 279 |
> |
OpenMD.\cite{Meineke:2005gd,openmd} We have tested the new method to |
| 280 |
> |
calculate the thermal conductance of a gold nanoparticle and SPC/E |
| 281 |
> |
water cluster, and compared the results with previous bulk RNEMD |
| 282 |
> |
values, as well as experiment. We have also investigated the |
| 283 |
> |
interfacial thermal conductance and interfacial rotational friction |
| 284 |
> |
for gold nanostructures solvated in hexane as a function of |
| 285 |
> |
nanoparticle size and shape. |
| 286 |
|
|
| 287 |
|
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% |
| 288 |
|
% FORCE FIELD PARAMETERS |
| 289 |
|
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% |
| 290 |
< |
\subsection{Force field parameters} |
| 290 |
> |
\subsection{Force field} |
| 291 |
|
|
| 292 |
< |
Gold -- gold interactions are described by the quantum Sutton-Chen (QSC) model.\cite{PhysRevB.59.3527} The QSC |
| 293 |
< |
parameters are tuned to experimental properties such as density, cohesive energy, and elastic moduli and |
| 294 |
< |
include zero-point quantum corrections. |
| 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 chosen the SPC/E water model for these simulations, which is particularly useful for method validation |
| 298 |
< |
as there are many values for physical properties from previous simulations available for direct |
| 299 |
< |
comparison.\cite{Bedrov:2000, Kuang2010} |
| 297 |
> |
The SPC/E water model~\cite{Berendsen87} is particularly useful for |
| 298 |
> |
validation of conductivities and shear viscosities. This model has |
| 299 |
> |
been used to previously test other RNEMD and NEMD approaches, and |
| 300 |
> |
there are reported values for thermal conductivies and shear |
| 301 |
> |
viscosities at a wide range of thermodynamic conditions that are |
| 302 |
> |
available for direct comparison.\cite{Bedrov:2000,Kuang:2010if} |
| 303 |
|
|
| 304 |
< |
Hexane molecules are described by the TraPPE united atom model,\cite{TraPPE-UA.alkanes} which provides good |
| 305 |
< |
computational efficiency and reasonable accuracy for bulk thermal conductivity values. In this model, sites are |
| 306 |
< |
located at the carbon centers for alkyl groups. Bonding interactions, including bond stretches and bends and |
| 307 |
< |
torsions, were used for intra-molecular sites closer than 3 bonds. For non-bonded interactions, Lennard-Jones |
| 308 |
< |
potentials were used. We have previously utilized both united atom (UA) and all-atom (AA) force fields for |
| 309 |
< |
thermal conductivity,\cite{kuang:AuThl,Kuang2012} and since the united atom force fields cannot populate the |
| 310 |
< |
high-frequency modes that are present in AA force fields, they appear to work better for modeling thermal |
| 311 |
< |
conductivity. |
| 304 |
> |
Hexane molecules are described by the TraPPE united atom |
| 305 |
> |
model,\cite{TraPPE-UA.alkanes} which provides good computational |
| 306 |
> |
efficiency and reasonable accuracy for bulk thermal conductivity |
| 307 |
> |
values. In this model, sites are located at the carbon centers for |
| 308 |
> |
alkyl groups. Bonding interactions, including bond stretches and bends |
| 309 |
> |
and torsions, were used for intra-molecular sites closer than 3 |
| 310 |
> |
bonds. For non-bonded interactions, Lennard-Jones potentials were |
| 311 |
> |
used. We have previously utilized both united atom (UA) and all-atom |
| 312 |
> |
(AA) force fields for thermal |
| 313 |
> |
conductivity,\cite{Kuang:2011ef,Kuang:2012fe,Stocker:2013cl} and since |
| 314 |
> |
the united atom force fields cannot populate the high-frequency modes |
| 315 |
> |
that are present in AA force fields, they appear to work better for |
| 316 |
> |
modeling thermal conductance at metal/ligand interfaces. |
| 317 |
|
|
| 318 |
< |
Gold -- hexane nonbonded interactions are governed by pairwise Lennard-Jones parameters derived from Vlugt |
| 319 |
< |
\emph{et al}.\cite{vlugt:cpc2007154} They fitted parameters for the interaction between Au and CH$_{\emph x}$ |
| 320 |
< |
pseudo-atoms based on the effective potential of Hautman and Klein for the Au(111) surface.\cite{hautman:4994} |
| 318 |
> |
Gold -- hexane nonbonded interactions are governed by pairwise |
| 319 |
> |
Lennard-Jones parameters derived from Vlugt \emph{et |
| 320 |
> |
al}.\cite{vlugt:cpc2007154} They fitted parameters for the |
| 321 |
> |
interaction between Au and CH$_{\emph x}$ pseudo-atoms based on the |
| 322 |
> |
effective potential of Hautman and Klein for the Au(111) |
| 323 |
> |
surface.\cite{hautman:4994} |
| 324 |
|
|
| 325 |
|
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% |
| 326 |
|
% NON-PERIODIC DYNAMICS |
| 343 |
|
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% |
| 344 |
|
\subsection{Simulation protocol} |
| 345 |
|
|
| 346 |
< |
In all cases, systems were fully equilibrated under non-periodic isobaric-isothermal (NPT) conditions -- using |
| 347 |
< |
the Langevin Hull methodology\cite{Vardeman2011} -- before any non-equilibrium methods were introduced. For |
| 348 |
< |
heterogeneous systems, the gold nanoparticles and ellipsoid were first created from a bulk lattice and |
| 349 |
< |
thermally equilibrated before being solvated in hexane. Packmol\cite{packmol} was used to solvate previously |
| 350 |
< |
equilibrated gold nanostructures within a spherical droplet of hexane. |
| 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 fully equilibrated, a thermal or angular momentum flux was applied for 1 - 2 |
| 356 |
< |
ns, until a stable temperature or angular velocity gradient had developed. Systems containing liquids were run |
| 357 |
< |
under moderate pressure (5 atm) and temperatures (230 K) to avoid the formation of a substantial vapor phase at the boundary of the cluster. Pressure was applied to the system via the non-periodic convex Langevin Hull. Thermal coupling to the external temperature and pressure bath was removed to avoid interference with any |
| 358 |
< |
imposed flux. |
| 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 |
< |
To stabilize the gold nanoparticle under the imposed angular momentum flux we altered the gold atom at the |
| 365 |
< |
designated coordinate origin to have $10,000$ times its original mass. The nonbonded interactions remain |
| 366 |
< |
unchanged and the heavy atom is excluded from the angular momentum exchange. For rotation of the ellipsoid about |
| 367 |
< |
its long axis we have added two heavy atoms along the axis of rotation, one at each end of the rod. We collected angular velocity data for the heterogeneous systems after a brief VSS-RNEMD simulation to initialize rotation of the solvated nanostructure. Doing so ensures that we overcome the initial static friction and calculate only the \emph{dynamic} interfacial rotational friction. |
| 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 |
< |
Fourier's Law of heat conduction in radial coordinates yields an expression for the heat flow between the |
| 392 |
< |
concentric spherical RNEMD shells: |
| 393 |
< |
|
| 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 |
< |
where $\lambda$ is the thermal conductivity, and $T_{a,b}$ and $r_{a,b}$ are the temperatures and radii of the |
| 403 |
< |
two RNEMD regions, respectively. |
| 404 |
< |
|
| 405 |
< |
A thermal flux is created using VSS-RNEMD moves, and the temperature in each of the radial shells is recorded. |
| 406 |
< |
The resulting temperature profiles are analyzed to yield information about the interfacial thermal conductance. |
| 407 |
< |
As the simulation progresses, the VSS moves are performed on a regular basis, and the system develops a thermal |
| 408 |
< |
or velocity gradient in response to the applied flux. Once a stable thermal gradient has been established |
| 409 |
< |
between the two regions, the thermal conductivity, $\lambda$, can be calculated using a linear regression of |
| 410 |
< |
the thermal gradient, $\langle \frac{dT}{dr} \rangle$: |
| 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} |
| 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{Kuang2012} Hu and Zwanzig\cite{Zwanzig} |
| 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 |
| 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{Kuang2010}, though still significantly lower than the experimental value |
| 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 |
|
|
| 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 |
|
|
| 545 |
– |
We gratefully acknowledge conversations with Dr. Shenyu Kuang. Support for this project was provided by the |
| 546 |
– |
National Science Foundation under grant CHE-0848243. Computational time was provided by the Center for Research |
| 547 |
– |
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} |
