trunk/nonperiodicVSS/nonperiodicVSS.tex

\documentclass[journal = jctcce, manuscript = article]{achemso}
\setkeys{acs}{usetitle = true}

\usepackage{caption}
\usepackage{geometry}
\usepackage{natbib}
\usepackage{setspace}
\usepackage{xkeyval}
%%%%%%%%%%%%%%%%%%%%%%%
\usepackage{amsmath}
\usepackage{amssymb}
\usepackage{times}
\usepackage{mathptm}
\usepackage{caption}
\usepackage{tabularx}
\usepackage{longtable}
\usepackage{graphicx}
\usepackage{achemso}
\usepackage{wrapfig} 
\usepackage[version=3]{mhchem}  % this is a great package for formatting chemical reactions
\usepackage{url}

\title{A method for creating thermal and angular momentum fluxes in
  non-periodic simulations}

\author{Kelsey M. Stocker}
\author{J. Daniel Gezelter}
\email{gezelter@nd.edu}
\affiliation[University of Notre Dame]{251 Nieuwland Science Hall\\ Department of Chemistry and Biochemistry\\ University of Notre Dame\\ Notre Dame, Indiana 46556}

\begin{document}

\begin{tocentry}
% \includegraphics[width=3.2cm]{figures/npVSS} \includegraphics[width=3.7cm]{figures/NP20}
\includegraphics[width=9cm]{figures/TOC}
\end{tocentry}

\newcolumntype{A}{p{1.5in}}
\newcolumntype{B}{p{0.75in}}

% \author{Kelsey M. Stocker and J. Daniel
%   Gezelter\footnote{Corresponding author. \ Electronic mail:
%     gezelter@nd.edu} \\
%   251 Nieuwland Science Hall, \\
%       Department of Chemistry and Biochemistry,\\
%       University of Notre Dame\\
%       Notre Dame, Indiana 46556}

%\date{\today}

%\maketitle 

%\begin{doublespace}

\begin{abstract}

We present a new reverse non-equilibrium molecular dynamics (RNEMD)
method that can be used with non-periodic simulation cells. This
method applies thermal and/or angular momentum fluxes between two
arbitrary regions of the simulation, and is capable of creating
stable temperature and angular velocity gradients while conserving
total energy and angular momentum.  One particularly useful
application is the exchange of kinetic energy between two concentric
spherical regions, which can be used to generate thermal transport
between nanoparticles and the solvent that surrounds them.  The
rotational couple to the solvent (a measure of interfacial friction)
is also available via this method.  As tests of the new method, we
have computed the thermal conductivities of gold nanoparticles and
water clusters, the interfacial thermal conductivity ($G$) of a
solvated gold nanoparticle, and the interfacial friction of a
variety of solvated gold nanostructures.

\end{abstract}

\newpage

%\narrowtext

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% **INTRODUCTION**
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Introduction}

Non-equilibrium molecular dynamics (NEMD) methods impose a temperature
or velocity {\it gradient} on a
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}
and use linear response theory to connect the resulting thermal or
momentum {\it flux} to transport coefficients of bulk materials,
\begin{equation}
j_z(p_x) = -\eta \frac{\partial v_x}{\partial z},  \hspace{0.5in}
J_z = \lambda \frac{\partial T}{\partial z}.
\end{equation}
Here, $\frac{\partial T}{\partial z}$ and $\frac{\partial
  v_x}{\partial z}$ are the imposed thermal and momentum gradients,
and as long as the imposed gradients are relatively small, the
corresponding fluxes, $J_z$ and $j_z(p_x)$, have a linear relationship
to the gradients.  The coefficients that provide this relationship
correspond to physical properties of the bulk material, either the
shear viscosity $(\eta)$ or thermal conductivity $(\lambda)$.  For
systems which include phase boundaries or interfaces, however, it is
often unclear what gradient (or discontinuity) should be imposed at
the boundary between materials.

In contrast, reverse Non-Equilibrium Molecular Dynamics (RNEMD)
methods impose an unphysical {\it flux} between different regions or
``slabs'' of the simulation
box.\cite{Muller-Plathe:1997wq,Muller-Plathe:1999ao,Patel:2005zm,Shenogina:2009ix,Tenney:2010rp,Kuang:2010if,Kuang:2011ef,Kuang:2012fe,Stocker:2013cl}
The system responds by developing a temperature or velocity {\it
  gradient} between the two regions. The gradients which develop in
response to the applied flux have the same linear response
relationships to the transport coefficient of interest. Since the
amount of the applied flux is known exactly, and measurement of a
gradient is generally less complicated, imposed-flux methods typically
take shorter simulation times to obtain converged results. At
interfaces, the observed gradients often exhibit near-discontinuities
at the boundaries between dissimilar materials. RNEMD methods do not
need many trajectories to provide information about transport
properties, and they have become widely used to compute thermal and
mechanical transport in both homogeneous liquids and
solids~\cite{Muller-Plathe:1997wq,Muller-Plathe:1999ao,Tenney:2010rp}
as well as heterogeneous
interfaces.\cite{Patel:2005zm,Shenogina:2009ix,Kuang:2010if,Kuang:2011ef,Kuang:2012fe,Stocker:2013cl}

The strengths of specific algorithms for imposing the flux between two
different slabs of the simulation cell has been the subject of some
renewed interest.  The original RNEMD approach used kinetic energy or
momentum exchange between particles in the two slabs, either through
direct swapping of momentum vectors or via virtual elastic collisions
between atoms in the two regions.  There have been recent
methodological advances which involve scaling all particle velocities
in both slabs.\cite{Kuang:2010if,Kuang:2012fe} Constraint equations
can be simultaneously imposed to require the simulation to conserve
both total energy and total linear momentum.  The most recent and
simplest of the velocity scaling approaches allows for simultaneous
shearing (to provide viscosity estimates) as well as scaling (to
provide information about thermal conductivity).\cite{Kuang:2012fe}

To date, the primary method of studying heat transfer at {\it curved}
nanoscale interfaces has involved transient non-equilibrium molecular
dynamics and temperature relaxation.\cite{Lervik:2009fk} RNEMD methods
have only been used in periodic simulation cells where the exchange
regions are physically separated along one of the axes of the
simulation cell. This limits the applicability to infinite planar
interfaces which are perpendicular to the applied flux. In order to
model steady-state non-equilibrium distributions for curved surfaces
(e.g. hot nanoparticles in contact with colder solvent), or for
regions that are not planar slabs, the method requires some
generalization for non-parallel exchange regions. In the following
sections, we present a new velocity shearing and scaling (VSS) RNEMD
algorithm which has been explicitly designed for non-periodic
simulations, and use the method to compute some thermal transport and
solid-liquid friction at the surfaces of spherical and ellipsoidal
nanoparticles.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% **METHODOLOGY**
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Velocity shearing and scaling (VSS) for non-periodic systems}

The original periodic VSS-RNEMD approach uses a series of simultaneous
velocity shearing and scaling exchanges between the two
slabs.\cite{Kuang:2012fe} This method imposes energy and linear
momentum conservation constraints while simultaneously creating a
desired flux between the two slabs. These constraints ensure that all
configurations are sampled from the same microcanonical (NVE)
ensemble.

\begin{figure}
\includegraphics[width=\linewidth]{figures/npVSS2}
\caption{Schematics of periodic (left) and non-periodic (right)
  Velocity Shearing and Scaling RNEMD. A kinetic energy or momentum
  flux is applied from region A to region B. Thermal gradients are
  depicted by a color gradient. Linear or angular velocity gradients
  are shown as arrows.}
\label{fig:VSS}
\end{figure}

We have extended the VSS method for use in {\it non-periodic}
simulations, in which the ``slabs'' have been generalized to two
separated regions of space. These regions could be defined as
concentric spheres (as in Figure \ref{fig:VSS}), or one of the regions
can be defined in terms of a dynamically changing ``hull'' comprising
the surface atoms of the cluster. This latter definition is identical
to the hull used in the Langevin Hull algorithm.\cite{Vardeman2011}
For the non-periodic variant, the constraints fix both the total
energy and total {\it angular} momentum of the system while
simultaneously imposing a thermal and angular momentum flux between
the two regions.

After a time interval of $\Delta t$, the particle velocities
($\mathbf{v}_i$ and $\mathbf{v}_j$) in the two shells ($A$ and $B$)
are modified by a velocity scaling coefficient ($a$ and $b$) and by a
rotational shearing term ($\mathbf{c}_a$ and $\mathbf{c}_b$).  The
scalars $a$ and $b$ collectively provide a thermal exchange between
the two regions.  One of the values is larger than 1, and the other
smaller. To conserve total energy and angular momentum, the values of
these two scalars are coupled.  The vectors ($\mathbf{c}_a$ and
$\mathbf{c}_b$) provide a relative rotational shear to the velocities
of the particles within the two regions, and these vectors must also
be coupled to constrain the total angular momentum.

Once the values of the scaling and shearing factors are known, the
velocity changes are applied,
\begin{displaymath}
\begin{array}{rclcl}
 & \underline{~~~~~~~~\mathrm{scaling}~~~~~~~~} & &
 \underline{\mathrm{rotational~shearing}} \\  \\
\mathbf{v}_i $~~~$\leftarrow & 
  a \left(\mathbf{v}_i - \langle \omega_a
  \rangle \times \mathbf{r}_i\right) & + & \mathbf{c}_a \times \mathbf{r}_i \\
\mathbf{v}_j $~~~$\leftarrow & 
  b  \left(\mathbf{v}_j - \langle \omega_b
  \rangle \times \mathbf{r}_j\right) & + & \mathbf{c}_b \times \mathbf{r}_j
\end{array}
\end{displaymath}
Here $\langle\mathbf{\omega}_a\rangle$ and $\langle\mathbf{\omega}_b\rangle$ are the instantaneous angular
velocities of each shell, and $\mathbf{r}_i$ is the position of particle $i$ relative to a fixed point in space
(usually the center of mass of the cluster). Particles in the shells also receive an additive ``angular shear''
to their velocities. The amount of shear is governed by the imposed angular momentum flux,
$\mathbf{j}_r(\mathbf{L})$,
\begin{eqnarray}
\mathbf{c}_a & = & - \mathbf{j}_r(\mathbf{L}) \cdot
\overleftrightarrow{I_a}^{-1} \Delta t   + \langle \omega_a \rangle \label{eq:bc}\\
\mathbf{c}_b & = & + \mathbf{j}_r(\mathbf{L}) \cdot
\overleftrightarrow{I_b}^{-1}  \Delta t  + \langle \omega_b \rangle \label{eq:bh}
\end{eqnarray}
where $\overleftrightarrow{I}_{\{a,b\}}$ is the moment of inertia
tensor for each of the two shells. 

To simultaneously impose a thermal flux ($J_r$) between the shells we
use energy conservation constraints, 
\begin{eqnarray}
K_a - J_r \Delta t & = & a^2 \left(K_a - \frac{1}{2}\langle
\omega_a \rangle \cdot \overleftrightarrow{I_a} \cdot \langle \omega_a
\rangle \right) + \frac{1}{2} \mathbf{c}_a \cdot \overleftrightarrow{I_a}
\cdot \mathbf{c}_a \label{eq:Kc}\\
K_b + J_r \Delta t & = & b^2 \left(K_b - \frac{1}{2}\langle
\omega_b \rangle \cdot \overleftrightarrow{I_b} \cdot \langle \omega_b
\rangle \right) + \frac{1}{2} \mathbf{c}_b \cdot \overleftrightarrow{I_b} \cdot \mathbf{c}_b \label{eq:Kh}
\end{eqnarray}
Simultaneous solution of these quadratic formulae for the scaling
coefficients, $a$ and $b$, will ensure that the simulation samples
from the original microcanonical (NVE) ensemble. Here $K_{\{a,b\}}$ is
the instantaneous translational kinetic energy of each shell. At each
time interval, we solve for $a$, $b$, $\mathbf{c}_a$, and
$\mathbf{c}_b$, subject to the imposed angular momentum flux,
$j_r(\mathbf{L})$, and kinetic energy flux, $J_r$, values. The new
particle velocities are computed, and the simulation continues. System
configurations after the transformations have exactly the same energy
({\it and} angular momentum) as before the moves.

As the simulation progresses, the velocity transformations can be
performed on a regular basis, and the system will develop a
temperature and/or angular velocity gradient in response to the
applied flux.  By fitting the radial temperature gradient, it is
straightforward to obtain the bulk thermal conductivity,
\begin{equation}
J_r = -\lambda \left( \frac{\partial T}{\partial r}\right)
\end{equation}
from the radial kinetic energy flux $(J_r)$ that was applied during
the simulation.  In practice, it is significantly easier to use the
integrated form of Fourier's law for spherical geometries.  In
sections \ref{sec:thermal} -- \ref{sec:rotation} we outline ways of
obtaining interfacial transport coefficients from these RNEMD
simulations.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% **COMPUTATIONAL DETAILS**
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Computational Details}

The new VSS-RNEMD methodology for non-periodic system geometries has
been implemented in our group molecular dynamics code,
OpenMD.\cite{Meineke:2005gd,openmd} We have tested the new method to
calculate the thermal conductance of a gold nanoparticle and SPC/E
water cluster, and compared the results with previous bulk RNEMD and
experimental values. We have also investigated the interfacial thermal
conductance and interfacial rotational friction for gold
nanostructures solvated in hexane as a function of nanoparticle size
and shape.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% FORCE FIELD PARAMETERS
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{Force fields}

Gold -- gold interactions are described by the quantum Sutton-Chen
(QSC) model.\cite{PhysRevB.59.3527} The QSC parameters are tuned to
experimental properties such as density, cohesive energy, and elastic
moduli and include zero-point quantum corrections.

The SPC/E water model~\cite{Berendsen87} is particularly useful for
validation of conductivities and shear viscosities.  This model has
been used to previously test other RNEMD and NEMD approaches, and
there are reported values for thermal conductivies and shear
viscosities at a wide range of thermodynamic conditions that are
available for direct comparison.\cite{Bedrov:2000,Kuang:2010if}

Hexane molecules are described by the TraPPE united atom
model,\cite{TraPPE-UA.alkanes} which provides good computational
efficiency and reasonable accuracy for bulk thermal conductivity
values. In this model, sites are located at the carbon centers for
alkyl groups. Bonding interactions, including bond stretches and bends
and torsions, were used for intra-molecular sites closer than 3
bonds. For non-bonded interactions, Lennard-Jones potentials were
used. We have previously utilized both united atom (UA) and all-atom
(AA) force fields for thermal
conductivity,\cite{Kuang:2011ef,Kuang:2012fe,Stocker:2013cl} and since
the united atom force fields cannot populate the high-frequency modes
that are present in AA force fields, they appear to work better for
modeling thermal conductance at metal/ligand interfaces.

Gold -- hexane nonbonded interactions are governed by pairwise
Lennard-Jones parameters derived from Vlugt \emph{et
  al}.\cite{vlugt:cpc2007154} They fitted parameters for the
interaction between Au and CH$_{\emph x}$ pseudo-atoms based on the
effective potential of Hautman and Klein for the Au(111)
surface.\cite{hautman:4994}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% NON-PERIODIC DYNAMICS
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% \subsection{Dynamics for non-periodic systems}
% 
% We have run all tests using the Langevin Hull non-periodic simulation methodology.\cite{Vardeman2011} The
% Langevin Hull is especially useful for simulating heterogeneous mixtures of materials with different
% compressibilities, which are typically problematic for traditional affine transform methods. We have had
% success applying this method to several different systems including bare metal nanoparticles, liquid water
% clusters, and explicitly solvated nanoparticles. Calculated physical properties such as the isothermal
% compressibility of water and the bulk modulus of gold nanoparticles are in good agreement with previous
% theoretical and experimental results. The Langevin Hull uses a Delaunay tesselation to create a dynamic convex
% hull composed of triangular facets with vertices at atomic sites. Atomic sites included in the hull are coupled
% to an external bath defined by a temperature, pressure and viscosity. Atoms not included in the hull are
% subject to standard Newtonian dynamics.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% SIMULATION PROTOCOL
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
\subsection{Simulation protocol}

In all cases, systems were equilibrated under non-periodic
isobaric-isothermal (NPT) conditions -- using the Langevin Hull
methodology\cite{Vardeman2011} -- before any non-equilibrium methods
were introduced. For heterogeneous systems, the gold nanoparticles and
ellipsoids were created from a bulk fcc lattice and were thermally
equilibrated before being solvated in hexane.  Packmol\cite{packmol}
was used to solvate previously equilibrated gold nanostructures within
a spherical droplet of liquid hexane.

Once equilibrated, thermal or angular momentum fluxes were applied for
1 - 2 ns, until stable temperature or angular velocity gradients had
developed. Systems containing liquids were run under moderate pressure
(5 atm) and temperatures (230 K) to avoid the formation of a vapor
layer at the boundary of the cluster.  Pressure was applied to the
system via the non-periodic Langevin Hull.\cite{Vardeman2011} However,
 coupling to the external temperature and pressure bath was
removed to avoid interference with the imposed RNEMD flux.

Because the method conserves \emph{total} angular momentum, systems
which contain a metal nanoparticle embedded in a significant volume of
solvent will still experience nanoparticle diffusion inside the
solvent droplet.  To aid in computing the interfacial transport properties in these
systems, a single gold atom at the origin of the coordinate system was
assigned a mass $10,000 \times$ its original mass. The bonded and
nonbonded interactions for this atom remain unchanged and the heavy
atom is excluded from the RNEMD exchanges.  The only effect of this
gold atom is to effectively pin the nanoparticle at the origin of the
coordinate system, while still allowing for any desired rotation about a given axis. For rotation of
the gold ellipsoids we added two of these heavy atoms along the axis
of rotation, separated by an equal distance from the origin of the
coordinate system.  These heavy atoms prevent off-axis tumbling of the
nanoparticle and allow for measurement of rotational friction relative
to a particular axis of the ellipsoid.

Angular velocity data was collected for the heterogeneous systems
after a brief period of imposed flux to initialize rotation of the
solvated nanostructure. Doing so ensures that we overcome the initial
static friction and calculate only the dynamic interfacial rotational
friction.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% THERMAL CONDUCTIVITIES
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{Thermal conductivities}
\label{sec:thermal}

To compute the thermal conductivities of bulk materials, the
integrated form of Fourier's Law of heat conduction in radial
coordinates can be used to obtain an expression for the heat transfer
rate between concentric spherical shells:
\begin{equation}
  q_r = -  \frac{4 \pi \lambda (T_b - T_a)}{\frac{1}{r_a} - \frac{1}{r_b}}
\label{eq:Q}
\end{equation}
The heat transfer rate, $q_r$, is constant in spherical geometries,
while the heat flux, $J_r$ depends on the surface area of the two
shells.  $\lambda$ is the thermal conductivity, and $T_{a,b}$ and
$r_{a,b}$ are the temperatures and radii of the two concentric RNEMD
regions, respectively.

A kinetic energy flux is created using VSS-RNEMD moves, and the
temperature in a set of 2~\AA\ thick radial shells is recorded.  The
resulting temperature profiles are analyzed to yield information about
the interfacial thermal conductance.  As the simulation progresses,
the VSS moves are performed on a regular basis, and the system
develops a thermal or velocity gradient in response to the applied
flux. Once a stable thermal gradient has been established between the
two regions, the thermal conductivity, $\lambda$, can be calculated
using a linear regression of the thermal gradient, $\langle
\frac{dT}{dr} \rangle$:

\begin{equation}
\lambda = \frac{r_a}{r_b} \frac{q_r}{4 \pi \langle \frac{dT}{dr} \rangle}
\label{eq:lambda}
\end{equation}

The rate of heat transfer, $q_r$, between the two RNEMD regions is
easily obtained from either the applied kinetic energy flux and the
area of the smaller of the two regions, or from the total amount of
transferred kinetic energy and the run time of the simulation.
\begin{equation}
        q_r = J_r \cdot A = \frac{\Delta KE}{\Delta t}
\label{eq:heat}
\end{equation}

\subsubsection{Thermal conductivity of nanocrystalline gold}
Calculated values for the thermal conductivity of a 40 \AA\ radius
gold nanoparticle (15707 atoms) at a range of kinetic energy flux
values are shown in Table \ref{table:goldTC}. For these calculations,
the hot and cold slabs were excluded from the linear regression of the
thermal gradient.

\begin{table}
  \caption{Calculated thermal gradients of a crystalline gold
    nanoparticle of radius 40 \AA\  subject to a range of applied
    kinetic energy fluxes. Calculations were performed at an average
    temperature of 300 K. Gold-gold interactions are described by the
    Quantum Sutton-Chen potential.}
  \label{table:goldTC}
  \begin{tabular}{cc}
    \\ \hline \hline
    {$J_r$} & {$\langle dT/dr \rangle$} \\ 
    {\footnotesize(kcal / fs $\cdot$ \AA$^{2}$)} & {\footnotesize(K /
      \AA)} \\ \hline \\
    3.25$\times 10^{-6}$ & 0.114  \\
    6.50$\times 10^{-6}$ & 0.232   \\
    1.30$\times 10^{-5}$ & 0.449  \\
    3.25$\times 10^{-5}$ & 1.180   \\
    6.50$\times 10^{-5}$ & 2.339    \\ \hline \hline
  \end{tabular}
\end{table}

The measured thermal gradients $\langle dT/dr \rangle$ are linearly
dependent on the applied kinetic energy flux $J_r$. The calculated
thermal conductivity value, $\lambda = 1.004 \pm
0.009$~W~/~m~$\cdot$~K compares well with previous {\it bulk} QSC
values of 1.08~--~1.26~W~/~m~$\cdot$~K\cite{Kuang:2010if}, though
still significantly lower than the experimental value of
320~W~/~m~$\cdot$~K, as the QSC force field neglects the significant
electronic contributions to thermal conductivity.  We note that there
is only a minimal effect on the computed thermal conductivity of gold
when comparing periodic (bulk) simulations carried out at the same
densities as those used here.

\subsubsection{Thermal conductivity of an SPC/E water droplet}

Calculated values for the thermal conductivity of a cluster of 6912
SPC/E water molecules were computed with a range of applied kinetic
energy fluxes. As with the gold nanoparticle, the measured slopes were
linearly dependent on the applied kinetic energy flux $J_r$, and the
RNEMD regions were excluded from the $\langle dT/dr \rangle$ fit (see
Figure \ref{fig:lambda_G}).

The computed mean value of the thermal conductivity, $\lambda = 0.884
\pm 0.013$~W~/~m~$\cdot$~K, compares well with previous
non-equilibrium molecular dynamics simulations of bulk SPC/E that were
carried out in periodic geometries.\cite{Zhang2005,Romer2012} These
simualtions gave conductivity values from 0.81--0.87~W~/~m~$\cdot$~K,
and all of the simulated values are approximately 30-45\% higher than the
experimental thermal conductivity of water, which has been measured at
0.61~W~/~m~$\cdot$~K.\cite{WagnerKruse}

\begin{figure}
  \includegraphics[width=\linewidth]{figures/lambda_G}
  \caption{Upper panel: temperature profile of a typical SPC/E water
    droplet.  The RNEMD simulation transfers kinetic energy from region
    A to region B, and the system responds by creating a temperature
    gradient (circles).  The fit to determine $\langle dT/dr \rangle$ is
    carried out between the two RNEMD exchange regions.  Lower panel:
    temperature profile for a solvated Au nanoparticle ($r = 20$ \AA).
    The temperature differences used to compute the total Kapitza
    resistance (and interfacial thermal conductance) are indicated with
    arrows.}
  \label{fig:lambda_G}
\end{figure}


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
% INTERFACIAL THERMAL CONDUCTANCE
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
\subsection{Interfacial thermal conductance}
\label{sec:interfacial}

The interfacial thermal conductance, $G$, of a heterogeneous interface
located at $r_0$ can be understood as the change in thermal
conductivity in a direction normal $(\mathbf{\hat{n}})$ to the
interface,
\begin{equation}
  G = \left(\nabla\lambda \cdot \mathbf{\hat{n}}\right)_{r_0}
\end{equation}
For heterogeneous systems such as the solvated nanoparticle shown in
Figure \ref{fig:NP20}, the interfacial thermal conductance at the
surface of the nanoparticle can be determined using a kinetic energy
flux applied using the RNEMD method developed above.  

In spherical geometries, it is most convenient to begin by finding the
Kapitza or interfacial thermal resistance for a thin spherical shell,
\begin{equation}
  R_K = \frac{1}{G} = \frac{\Delta
    T}{J_r}
\end{equation}
where $\Delta T$ is the temperature drop from the interior to the
exterior of the shell.  For two concentric shells, the kinetic energy
flux ($J_r$) is not the same (as the surface areas are not the same),
but the heat transfer rate, $q_r = J_r \cdot A$ is constant.  The thermal
resistance of a shell with interior radius $r$ is most conveniently
written as
\begin{equation}
  R_K = \frac{1}{q_r} \Delta T 4 \pi r^2.
  \label{eq:RK}
\end{equation}

\begin{figure}
  \includegraphics[width=\linewidth]{figures/NP20}
  \caption{A gold nanoparticle with a radius of 20 \AA\ solvated in
    TraPPE-UA hexane. A kinetic energy flux is applied between the
    nanoparticle and an outer shell of solvent to obtain the interfacial
    thermal conductance, $G$, and the interfacial rotational resistance
    $\Xi^{rr}$ is determined using an angular momentum flux. }
  \label{fig:NP20}
\end{figure}

To model the thermal conductance across a wide interface (or multiple
interfaces) it is useful to consider concentric shells as resistors
wired in series.  The resistance of the shells is then additive, and
the interfacial thermal conductance is the inverse of the total
Kapitza resistance:
\begin{equation}
  \frac{1}{G} = R_\mathrm{total} = \frac{1}{q_r} \sum_i \left(T_{i+i} -
    T_i\right) 4 \pi r_i^2 
\end{equation}
This series can be extended for any number of adjacent shells,
allowing for the calculation of the interfacial thermal conductance
for interfaces of considerable thickness, such as self-assembled
ligand monolayers on a metal surface.

\subsubsection{Interfacial thermal conductance of solvated gold
  nanoparticles}
Calculated interfacial thermal conductance ($G$) values for three
sizes of gold nanoparticles and a flat Au(111) surface solvated in
TraPPE-UA hexane~\cite{Stocker:2013cl} are shown in Table
\ref{table:G}.

\begin{table}
  \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.}
  \label{table:G}
  \begin{tabular}{cc}
    \hline \hline 
    {Nanoparticle Radius} & {$G$}\\ 
    {\footnotesize(\AA)} & {\footnotesize(MW / m$^{2}$ $\cdot$ K)}\\ \hline
    20 & {47.1 $\pm$ 5.0} \\
    30 & {45.4 $\pm$ 1.2} \\
    40 & {46.5 $\pm$ 2.1} \\
    \hline
    Au(111) & {30.2} \\ 
    \hline  \hline
  \end{tabular}
\end{table}

The introduction of surface curvature increases the interfacial
thermal conductance by a factor of approximately $1.5$ relative to the
flat interface, although there is no apparent size-dependence in the
$G$ values obtained from our calculations.  This is quite different
behavior than the size-dependence observed by Lervik \textit{et
  al.}\cite{Lervik:2009fk} in their NEMD simulations of decane
droplets in water.  In the liquid / liquid case, the surface tension is
strongly dependent on the droplet size, and the interfacial
resistivity increases with surface tension.  In our case, the surface
tension of the solid / liquid interface does not have a strong
dependence on the particle radius at these particle sizes.

Simulations of larger nanoparticles may yet demonstrate limiting $G$
values close to the flat Au(111) slab, although any spherical particle
will have a significant fraction of the surface atoms in non-111
facets.  It is not clear at this point whether the increase in thermal
conductance for the spherical particles is due to the increased
population of undercoordinated surface atoms, or to increased solid
angle space available for phonon scattering into the solvent. These
two possibilities do open up some interesting avenues for
investigation.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% INTERFACIAL ROTATIONAL FRICTION
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{Interfacial rotational friction}
\label{sec:rotation}
The interfacial rotational resistance tensor, $\Xi^{rr}$, can be
calculated for heterogeneous nanostructure / solvent systems by applying
an angular momentum flux between the solvated nanostructure and a
spherical shell of solvent at the outer edge of the cluster. An
angular velocity gradient develops in response to the applied flux,
causing the nanostructure and solvent shell to rotate in opposite
directions about a given axis.

\begin{figure}
\includegraphics[width=\linewidth]{figures/E25-75}
\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.}
\label{fig:E25-75}
\end{figure}

Analytical solutions for the diagonal elements of the rotational
resistance tensor for a solvated spherical body of radius $r$ under
ideal stick boundary conditions can be estimated using Stokes' law
\begin{equation}
        \Xi^{rr}_{stick} = 8 \pi \eta r^3
\label{eq:Xisphere}
\end{equation}
where $\eta$ is the dynamic viscosity of the surrounding solvent. 

For general ellipsoids with semiaxes $\alpha$, $\beta$, and $\gamma$,
Perrin's extension of Stokes' law provides exact solutions for
symmetric prolate $(\alpha \geq \beta = \gamma)$ and oblate $(\alpha <
\beta = \gamma)$ ellipsoids under ideal stick conditions.  The Perrin
elliptic integral parameter,
\begin{equation}
  S = \frac{2}{\sqrt{\alpha^2 - \beta^2}} \ln \left[ \frac{\alpha + \sqrt{\alpha^2 - \beta^2}}{\beta} \right].
\label{eq:S}
\end{equation}

For a prolate ellipsoidal nanoparticle (see Figure \ref{fig:E25-75}),
the rotational resistance tensor $\Xi^{rr}$ is a $3 \times 3$ diagonal
matrix with elements
\begin{align}
        \Xi^{rr}_\alpha =& \frac{32 \pi}{3} \eta \frac{ \left(
            \alpha^2 - \beta^2 \right) \beta^2}{2\alpha - \beta^2 S}
        \\ \nonumber \\
        \Xi^{rr}_{\beta,\gamma} =& \frac{32 \pi}{3} \eta \frac{ \left( \alpha^4 - \beta^4 \right)}{ \left( 2\alpha^2 - \beta^2 \right)S - 2\alpha},
\label{eq:Xirr}
\end{align}
corresponding to rotation about the long axis ($\alpha$), and each of
the equivalent short axes ($\beta$ and $\gamma$), respectively.

Previous VSS-RNEMD simulations of the interfacial friction of the
planar Au(111) / hexane interface have shown that the flat interface
exists within slip boundary conditions.\cite{Kuang:2012fe} Hu and
Zwanzig\cite{Zwanzig} investigated the rotational friction
coefficients for spheroids under slip boundary conditions and obtained
numerial results for a scaling factor to be applied to
$\Xi^{rr}_{\mathit{stick}}$ as a function of ratio of the shorter
semiaxes to the longer.  Under slip conditions, rotation of any
sphere and rotation of a prolate ellipsoid about its long axis does
not displace any solvent, so the resulting $\Xi^{rr}_{\mathit{slip}}$
approaches $0$. For rotation of a prolate ellipsoid (with the aspect
ratio shown here) about its short axis, Hu and Zwanzig obtained
$\Xi^{rr}_{\mathit{slip}}$ values of $35.9\%$ of the analytical
$\Xi^{rr}_{\mathit{stick}}$ result.  This reduced rotational
resistance accounts for the reduced interfacial friction under slip
boundary conditions.

The simulated rotational friction coefficient,
$\Xi^{rr}_{\mathit{sim}}$, at the interface can be extracted from
non-periodic VSS-RNEMD simulations quite easily using the total
applied torque ($\tau$) and the observed angular velocity of the solid
structure ($\omega_\mathrm{solid}$),
\begin{equation}
        \Xi^{rr}_{\mathit{sim}} = \frac{\tau}{\omega_\mathrm{solid}}
\label{eq:Xieff}
\end{equation}

The total applied torque required to overcome the interfacial friction
and maintain constant angular rotation of the solid particle is
\begin{equation}
        \tau = \frac{j_r(\mathbf{L}) \cdot A}{2} 
\label{eq:tau}  
\end{equation}
where $j_r(\mathbf{L})$ is the applied angular momentum flux, $A$ is
the surface area of the solid nanoparticle, and the factor of $2$ is
present because the torque is exerted on both the particle and the
counter-rotating solvent shell.

\subsubsection{Rotational friction for gold nanostructures in hexane}

Table \ref{table:couple} shows the calculated rotational friction
coefficients $\Xi^{rr}$ for spherical gold nanoparticles and a prolate
ellipsoidal gold nanorod solvated in TraPPE-UA hexane. An angular
momentum flux was applied between an outer shell of solvent that was
at least 20 \AA\ away from the surface of the particle (region $A$)
and the gold nanostructure (region $B$).  Dynamic viscosity estimates
$(\eta)$ for TraPPE-UA hexane under these particular temperature and
pressure conditions were determined by applying a traditional VSS-RNEMD
linear momentum flux to a periodic box of this solvent.

\begin{table}
  \caption{Comparison of rotational friction coefficients under ideal stick
    ($\Xi^{rr}_{\mathit{stick}}$) and 
    slip ($\Xi^{rr}_{\mathit{slip}}$) conditions, along with the simulated
    ($\Xi^{rr}_{\mathit{sim}}$) 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.} 
\label{table:couple}
\begin{tabular}{lccccc}
\\ \hline \hline
{Structure} & {Axis of Rotation} & {$\Xi^{rr}_{\mathit{stick}}$} & {$\Xi^{rr}_{\mathit{slip}}$} & {$\Xi^{rr}_{\mathit{sim}}$}  & {$\Xi^{rr}_{\mathit{sim}}$ / $\Xi^{rr}_{\mathit{stick}}$}\\ 
{} & {} & {\footnotesize(amu A$^2$ / fs)} & {\footnotesize(amu A$^2$ / fs)} & {\footnotesize(amu A$^2$ / fs)} & \\  \hline
Sphere (r = 20 \AA) & {$x = y = z$} & {3314} & 0      & {2386 $\pm$ 14}     & {0.720 $\pm$ 0.004}\\
Sphere (r = 30 \AA) & {$x = y = z$} & {11749}& 0      & {8415 $\pm$ 274}    & {0.716 $\pm$ 0.023}\\
Sphere (r = 40 \AA) & {$x = y = z$} & {34464}& 0      & {47544 $\pm$ 3051}  & {1.380 $\pm$ 0.089}\\ 
Prolate Ellipsoid   & {$x = y$}     & {4991} & {1792} & {3128 $\pm$ 166}    & {0.627 $\pm$ 0.033}\\
Prolate Ellipsoid   & {$z$}         & {1993} & 0      & {1590 $\pm$ 30}     & {0.798 $\pm$ 0.015}
\\ \hline \hline
\end{tabular}
\end{table}

The results for $\Xi^{rr}_{\mathit{sim}}$ show that, in contrast with
the flat Au(111) / hexane interface, gold nanostructures solvated by
hexane are closer to stick than slip boundary conditions.  At these
length scales, the nanostructures are not perfect spheroids due to
atomic `roughening' of the surface and therefore experience increased
interfacial friction, deviating significantly from the ideal slip
case. The 20 and 30~\AA\ radius nanoparticles experience approximately
70\% of the ideal stick boundary interfacial friction. Rotation of the
ellipsoid about its long axis more closely approaches the stick limit
than rotation about the short axis, which at first seems
counterintuitive. However, the `propellor' motion caused by rotation
around the short axis may exclude solvent from the rotation cavity or
move a sufficient amount of solvent along with the gold that a smaller
interfacial friction is actually experienced. 

The largest nanoparticle (40 \AA\ radius) appears to experience
interfacial friction in excess of ideal stick conditions.  Although
the solvent velocity field does not appear in the calculation of
$\Xi^{rr}$, we do note that the smaller particles exhibit solvent
velocity fields that can be fit with the $v(r,\theta) = \sin\theta
\left( A r + B/r^2 \right)$ form first discovered by
Jeffery,\cite{Jeffery:1915fk} while the solvent velocity fields for
the $r=40$ \AA\ particle approaches the infinite fluid solution and
can be fit well with only the second term.  The larger particle (and
the resulting larger solvent volume) appears to be transitioning to a
regime that approaches the infinite solvent solutions.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% **DISCUSSION**
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Discussion} 

We have adapted the VSS-RNEMD methodology so that it can be used to
apply kinetic energy and angular momentum fluxes in explicitly
non-periodic simulation geometries. Non-periodic VSS-RNEMD preserves
the strengths of the original periodic variant, specifically Boltzmann
velocity distributions and minimal thermal anisotropy, while extending
the constraints to conserve total energy and total \emph{angular}
momentum. This method also maintains the ability to impose the kinetic
energy and angular momentum fluxes either jointly or individually.

This method enables steady-state calculation of interfacial thermal
conductance and interfacial rotational friction in heterogeneous
clusters.  We have demonstrated the abilities of the method on some
familiar systems (nanocrystalline gold, SPC/E water) as well as on
interfaces that have been studied only in planar geometries.  For the
bulk systems that are well-understood in periodic geometries, the
non-periodic VSS-RNEMD provides very similar estimates of the thermal
conductivity to other simulation methods.

For nanoparticle-to-solvent heat transfer, we have observed that the
interfacial conductance exhibits a $1.5\times$ increase over the
planar Au(111) interface, although we do not see any dependence of the
conductance on the particle size.  Because the surface tension effects
present in liquid/liquid heat transfer are not strong contributors
here, we suggest two possible mechanisms for the $1.5\times$ increase
over the planar surface: (1) the nanoparticles have an increased
population of undercoordinated surface atoms that are more efficient
at transferring vibrational energy to the solvent, or (2) there is
increased solid angle space available for phonon scattering into the
solvent. The second of these explanations would have a significant
dependence on the radius of spherical particles, although crystalline
faceting of the particles may be enough to reduce this dependence on
particle radius.  These two possibilities open up some interesting
avenues for further exploration.

One area that this method opens up that was not available for periodic
simulation cells is direct computation of the solid / liquid
rotational friction coefficients of nanostructures.  The systems we
have investigated (gold nanospheres and prolate ellipsoids) have
analytic stick and numerical slip solutions provided via
hydrodynamics, and our molecular simulations indicate that the bare
metallic nanoparticles are closer to stick than they are to slip
conditions. This is markedly different behavior than we have seen for
the solid / liquid friction at planar Au(111) interfaces in previous
simulations, where the solvents experience a nearly-laminar flow over
the surface.  Schmidt and Skinner previously computed the behavior of
spherical tag particles in molecular dynamics simulations, and showed
that {\it slip} boundary conditions for translational resistance may
be more appropriate for molecule-sized spheres embedded in a sea of
spherical solvent particles.\cite{Schmidt:2004fj,Schmidt:2003kx} The
size of the gold nanoparticles studied here and the structuring of the
particle surfaces appear to bring the behavior closer to stick
boundaries.  Exactly where the slip-stick crossover takes place is
also an interesting avenue for future exploration.

The ability to interrogate non-periodic effects (e.g. surface
curvature, edge effects between facets, and the splay of ligand
surface groups) on interfacial transport means that this method can be
applied to systems of broad experimental and theoretical interest.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% **ACKNOWLEDGMENTS**
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{acknowledgement}
  The authors thank Dr. Shenyu Kuang for helpful discussions. Support
  for this project was provided by the National Science Foundation
  under grant CHE-0848243. Computational time was provided by the
  Center for Research Computing (CRC) at the University of Notre Dame.
\end{acknowledgement}


\newpage

\bibliography{nonperiodicVSS}

%\end{doublespace}
\end{document}
Revision:	4083
Committed:	Sat Mar 15 23:00:32 2014 UTC (11 years, 4 months ago) by gezelter
Content type:	application/x-tex
File size:	41737 byte(s)
Log Message:	Some minor edits