| 29 |
|
9.0in \textwidth 6.5in \brokenpenalty=10000 |
| 30 |
|
|
| 31 |
|
% double space list of tables and figures |
| 32 |
< |
% \AtBeginDelayedFloats{\renewcomand{\baselinestretch}{1.66}} |
| 32 |
> |
% \AtBeginDelayedFloats{\renewcomand{\baselinestretch}{1.66}} |
| 33 |
|
\setlength{\abovecaptionskip}{20 pt} |
| 34 |
|
\setlength{\belowcaptionskip}{30 pt} |
| 35 |
|
|
| 79 |
|
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% |
| 80 |
|
\section{Introduction} |
| 81 |
|
|
| 82 |
+ |
Non-equilibrium Molecular Dynamics (NEMD) methods impose a temperature |
| 83 |
+ |
or velocity {\it gradient} on a |
| 84 |
+ |
system,\cite{ASHURST:1975tg,Evans:1982zk,ERPENBECK:1984sp,MAGINN:1993hc,Berthier:2002ij,Evans:2002ai,Schelling:2002dp,PhysRevA.34.1449,JiangHao_jp802942v} |
| 85 |
+ |
and use linear response theory to connect the resulting thermal or |
| 86 |
+ |
momentum flux to transport coefficients of bulk materials. However, |
| 87 |
+ |
for heterogeneous systems, such as phase boundaries or interfaces, it |
| 88 |
+ |
is often unclear what shape of gradient should be imposed at the |
| 89 |
+ |
boundary between materials. |
| 90 |
|
|
| 91 |
+ |
\begin{figure} |
| 92 |
+ |
\includegraphics[width=\linewidth]{figures/npVSS} |
| 93 |
+ |
\caption{Schematics of periodic (left) and non-periodic (right) |
| 94 |
+ |
Velocity Shearing and Scaling RNEMD. A kinetic energy or momentum |
| 95 |
+ |
flux is applied from region B to region A. Thermal gradients are |
| 96 |
+ |
depicted by a color gradient. Linear or angular velocity gradients |
| 97 |
+ |
are shown as arrows.} |
| 98 |
+ |
\label{fig:VSS} |
| 99 |
+ |
\end{figure} |
| 100 |
+ |
|
| 101 |
+ |
Reverse Non-Equilibrium Molecular Dynamics (RNEMD) methods impose an |
| 102 |
+ |
unphysical {\it flux} between different regions or ``slabs'' of the |
| 103 |
+ |
simulation |
| 104 |
+ |
box.\cite{MullerPlathe:1997xw,ISI:000080382700030,Kuang2010} The |
| 105 |
+ |
system responds by developing a temperature or velocity {\it gradient} |
| 106 |
+ |
between the two regions. The gradients which develop in response to |
| 107 |
+ |
the applied flux are then related (via linear response theory) to the |
| 108 |
+ |
transport coefficient of interest. Since the amount of the applied |
| 109 |
+ |
flux is known exactly, and measurement of a gradient is generally less |
| 110 |
+ |
complicated, imposed-flux methods typically take shorter simulation |
| 111 |
+ |
times to obtain converged results. At interfaces, the observed |
| 112 |
+ |
gradients often exhibit near-discontinuities at the boundaries between |
| 113 |
+ |
dissimilar materials. RNEMD methods do not need many trajectories to |
| 114 |
+ |
provide information about transport properties, and they have become |
| 115 |
+ |
widely used to compute thermal and mechanical transport in both |
| 116 |
+ |
homogeneous liquids and |
| 117 |
+ |
solids~\cite{MullerPlathe:1997xw,ISI:000080382700030,Maginn:2010} as |
| 118 |
+ |
well as heterogeneous |
| 119 |
+ |
interfaces.\cite{garde:nl2005,garde:PhysRevLett2009,kuang:AuThl} |
| 120 |
+ |
|
| 121 |
|
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% |
| 122 |
|
% **METHODOLOGY** |
| 123 |
|
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% |
| 124 |
< |
\section{Methodology} |
| 124 |
> |
\section{Velocity Shearing and Scaling (VSS) for non-periodic systems} |
| 125 |
|
|
| 126 |
< |
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% |
| 127 |
< |
% NON-PERIODIC DYNAMICS |
| 128 |
< |
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% |
| 129 |
< |
\subsection{Dynamics for non-periodic systems} |
| 126 |
> |
The VSS-RNEMD approach uses a series of simultaneous velocity shearing |
| 127 |
> |
and scaling exchanges between the two slabs.\cite{Kuang2012} |
| 128 |
> |
This method imposes energy and momentum conservation constraints while |
| 129 |
> |
simultaneously creating a desired flux between the two slabs. These |
| 130 |
> |
constraints ensure that all configurations are sampled from the same |
| 131 |
> |
microcanonical (NVE) ensemble. |
| 132 |
|
|
| 133 |
< |
We have run all tests using the Langevin Hull nonperiodic simulation methodology.\cite{Vardeman2011} The Langevin Hull is especially useful for simulating heterogeneous mixtures of materials with different compressibilities, which are typically problematic for traditional affine transform methods. We have had success applying this method to |
| 134 |
< |
several different systems including bare metal nanoparticles, liquid water clusters, and explicitly solvated nanoparticles. Calculated physical properties such as the isothermal compressibility of water and the bulk modulus of gold nanoparticles are in good agreement with previous theoretical and experimental results. The Langevin Hull uses a Delaunay tesselation to create a dynamic convex hull composed of triangular facets with vertices at atomic sites. Atomic sites included in the hull are coupled to an external bath defined by a temperature, pressure and viscosity. Atoms not included in the hull are subject to standard Newtonian dynamics. For these tests, thermal coupling to the bath was turned off to avoid interference with any imposed kinetic flux. Systems containing liquids were run under moderate pressure ($\sim$ 5 atm) to avoid the formation of a substantial vapor phase. |
| 133 |
> |
We have extended the VSS method for use in {\it non-periodic} |
| 134 |
> |
simulations, in which the ``slabs'' have been generalized to two |
| 135 |
> |
separated regions of space. These regions could be defined as |
| 136 |
> |
concentric spheres (as in figure \ref{fig:VSS}), or one of the regions |
| 137 |
> |
can be defined in terms of a dynamically changing ``hull'' comprising |
| 138 |
> |
the surface atoms of the cluster. This latter definition is identical |
| 139 |
> |
to the hull used in the Langevin Hull algorithm. |
| 140 |
|
|
| 141 |
< |
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% |
| 142 |
< |
% NON-PERIODIC RNEMD |
| 143 |
< |
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% |
| 144 |
< |
\subsection{VSS-RNEMD for non-periodic systems} |
| 141 |
> |
We present here a new set of constraints that are more general than |
| 142 |
> |
the VSS constraints. For the non-periodic variant, the constraints |
| 143 |
> |
fix both the total energy and total {\it angular} momentum of the |
| 144 |
> |
system while simultaneously imposing a thermal and angular momentum |
| 145 |
> |
flux between the two regions. |
| 146 |
|
|
| 147 |
< |
The most useful RNEMD approach developed so far utilizes a series of |
| 148 |
< |
simultaneous velocity shearing and scaling (VSS) exchanges between the two |
| 149 |
< |
regions.\cite{Kuang2012} This method provides a set of conservation constraints |
| 150 |
< |
while simultaneously creating a desired flux between the two regions. Satisfying |
| 105 |
< |
the constraint equations ensures that the new configurations are sampled from the |
| 106 |
< |
same NVE ensemble. |
| 107 |
< |
|
| 108 |
< |
We have implemented this new method in OpenMD, our group molecular dynamics code.\cite{openmd} The adaptation of VSS-RNEMD for non-periodic systems is relatively |
| 109 |
< |
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 |
| 110 |
< |
of rectangular slabs. A temperature profile along the $r$ coordinate is created by recording the average temperature of concentric spherical shells. |
| 111 |
< |
|
| 112 |
< |
\begin{figure} |
| 113 |
< |
\center{\includegraphics[width=7in]{figures/VSS}} |
| 114 |
< |
\caption{Schematics of periodic (left) and nonperiodic (right) Velocity Shearing and Scaling RNEMD. A kinetic energy or momentum flux is applied from region B to region A. Thermal gradients are depicted by a color gradient. Linear or angular velocity gradients are shown as arrows.} |
| 115 |
< |
\label{fig:VSS} |
| 116 |
< |
\end{figure} |
| 117 |
< |
|
| 118 |
< |
At each time interval, the particle velocities ($\mathbf{v}_i$ and |
| 119 |
< |
$\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. |
| 147 |
> |
After each $\Delta t$ time interval, the particle velocities |
| 148 |
> |
($\mathbf{v}_i$ and $\mathbf{v}_j$) in the two shells ($A$ and $B$) |
| 149 |
> |
are modified by a velocity scaling coefficient ($a$ and $b$) and by a |
| 150 |
> |
rotational shearing term ($\mathbf{c}_a$ and $\mathbf{c}_b$). |
| 151 |
|
\begin{displaymath} |
| 152 |
|
\begin{array}{rclcl} |
| 153 |
< |
& \underline{~~~~~~~~~~\mathrm{scaling}~~~~~~~~~~} & & |
| 154 |
< |
\underline{\mathrm{rotational \; shearing}} \\ \\ |
| 153 |
> |
& \underline{~~~~~~~~\mathrm{scaling}~~~~~~~~} & & |
| 154 |
> |
\underline{\mathrm{rotational~shearing}} \\ \\ |
| 155 |
|
\mathbf{v}_i $~~~$\leftarrow & |
| 156 |
< |
c \, \left(\mathbf{v}_i - \langle \omega_c |
| 157 |
< |
\rangle \times r_i\right) & + & \mathbf{b}_c \times r_i \\ |
| 156 |
> |
a \left(\mathbf{v}_i - \langle \omega_a |
| 157 |
> |
\rangle \times \mathbf{r}_i\right) & + & \mathbf{c}_a \times \mathbf{r}_i \\ |
| 158 |
|
\mathbf{v}_j $~~~$\leftarrow & |
| 159 |
< |
h \, \left(\mathbf{v}_j - \langle \omega_h |
| 160 |
< |
\rangle \times r_j\right) & + & \mathbf{b}_h \times r_j |
| 159 |
> |
b \left(\mathbf{v}_j - \langle \omega_b |
| 160 |
> |
\rangle \times \mathbf{r}_j\right) & + & \mathbf{c}_b \times \mathbf{r}_j |
| 161 |
|
\end{array} |
| 162 |
|
\end{displaymath} |
| 163 |
< |
|
| 163 |
> |
Here $\langle\mathbf{\omega}_a\rangle$ and |
| 164 |
> |
$\langle\mathbf{\omega}_b\rangle$ are the instantaneous angular |
| 165 |
> |
velocities of each shell, and $\mathbf{r}_i$ is the position of |
| 166 |
> |
particle $i$ relative to a fixed point in space (usually the center of |
| 167 |
> |
mass of the cluster). Particles in the shells also receive an |
| 168 |
> |
additive ``angular shear'' to their velocities. The amount of shear |
| 169 |
> |
is governed by the imposed angular momentum flux, |
| 170 |
> |
$\mathbf{j}_r(\mathbf{L})$, |
| 171 |
|
\begin{eqnarray} |
| 172 |
< |
\mathbf{b}_c & = & - \mathbf{j}_r(\mathbf{L}) \; \Delta \, t \; I_c^{-1} + \langle \omega_c \rangle \label{eq:bc}\\ |
| 173 |
< |
\mathbf{b}_h & = & + \mathbf{j}_r(\mathbf{L}) \; \Delta \, t \; I_h^{-1} + \langle \omega_h \rangle \label{eq:bh} |
| 172 |
> |
\mathbf{c}_a & = & - \mathbf{j}_r(\mathbf{L}) \cdot |
| 173 |
> |
\overleftrightarrow{I_a}^{-1} \Delta t + \langle \omega_a \rangle \label{eq:bc}\\ |
| 174 |
> |
\mathbf{c}_b & = & + \mathbf{j}_r(\mathbf{L}) \cdot |
| 175 |
> |
\overleftrightarrow{I_b}^{-1} \Delta t + \langle \omega_b \rangle \label{eq:bh} |
| 176 |
|
\end{eqnarray} |
| 177 |
+ |
where $\overleftrightarrow{I}_{\{a,b\}}$ is the moment of inertia tensor for |
| 178 |
+ |
each of the two shells. |
| 179 |
|
|
| 180 |
< |
The total energy is constrained via two quadratic formulae, |
| 180 |
> |
To simultaneously impose a thermal flux ($J_r$) between the shells we |
| 181 |
> |
use energy conservation constraints, |
| 182 |
|
\begin{eqnarray} |
| 183 |
< |
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}\\ |
| 184 |
< |
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} |
| 183 |
> |
K_a - J_r \Delta t & = & a^2 \left(K_a - \frac{1}{2}\langle |
| 184 |
> |
\omega_a \rangle \cdot \overleftrightarrow{I_a} \cdot \langle \omega_a |
| 185 |
> |
\rangle \right) + \frac{1}{2} \mathbf{c}_a \cdot \overleftrightarrow{I_a} |
| 186 |
> |
\cdot \mathbf{c}_a \label{eq:Kc}\\ |
| 187 |
> |
K_b + J_r \Delta t & = & b^2 \left(K_b - \frac{1}{2}\langle |
| 188 |
> |
\omega_b \rangle \cdot \overleftrightarrow{I_b} \cdot \langle \omega_b |
| 189 |
> |
\rangle \right) + \frac{1}{2} \mathbf{c}_b \cdot \overleftrightarrow{I_b} \cdot \mathbf{c}_b \label{eq:Kh} |
| 190 |
|
\end{eqnarray} |
| 191 |
+ |
Simultaneous solution of these quadratic formulae for the scaling |
| 192 |
+ |
coefficients, $a$ and $b$, will ensure that the simulation samples |
| 193 |
+ |
from the original microcanonical (NVE) ensemble. Here $K_{\{a,b\}}$ |
| 194 |
+ |
is the instantaneous translational kinetic energy of each shell. At |
| 195 |
+ |
each time interval, we solve for $a$, $b$, $\mathbf{c}_a$, and |
| 196 |
+ |
$\mathbf{c}_b$, subject to the imposed angular momentum flux, |
| 197 |
+ |
$j_r(\mathbf{L})$, and thermal flux, $J_r$, values. The new particle |
| 198 |
+ |
velocities are computed, and the simulation continues. System |
| 199 |
+ |
configurations after the transformations have exactly the same energy |
| 200 |
+ |
({\it and} angular momentum) as before the moves. |
| 201 |
|
|
| 202 |
< |
the simultaneous |
| 203 |
< |
solution of which provide the velocity scaling coefficients $c$ and $h$. Given an |
| 204 |
< |
imposed angular momentum flux, $\mathbf{j}_{r} \left( \mathbf{L} \right)$, and/or |
| 205 |
< |
thermal flux, $J_r$, equations \ref{eq:bc} - \ref{eq:Kh} are sufficient to obtain |
| 206 |
< |
the velocity scaling ($c$ and $h$) and shearing ($\mathbf{b}_c,\,$ and $\mathbf{b}_h$) at each time interval. |
| 202 |
> |
As the simulation progresses, the velocity transformations can be |
| 203 |
> |
performed on a regular basis, and the system will develop a |
| 204 |
> |
temperature and/or angular velocity gradient in response to the |
| 205 |
> |
applied flux. Using the slope of the radial temperature or velocity |
| 206 |
> |
gradients, it is quite simple to obtain both the thermal conductivity |
| 207 |
> |
($\lambda$), interfacial thermal conductance ($G$), or rotational friction coefficients ($\Xi^{rr}$) of any nonperiodic system. |
| 208 |
|
|
| 209 |
|
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% |
| 210 |
|
% **COMPUTATIONAL DETAILS** |
| 212 |
|
\section{Computational Details} |
| 213 |
|
|
| 214 |
|
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% |
| 215 |
< |
% SIMULATION PROTOCOL |
| 215 |
> |
% NON-PERIODIC DYNAMICS |
| 216 |
|
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% |
| 217 |
+ |
\subsection{Dynamics for non-periodic systems} |
| 218 |
+ |
|
| 219 |
+ |
We have run all tests using the Langevin Hull nonperiodic simulation methodology.\cite{Vardeman2011} The Langevin Hull is especially useful for simulating heterogeneous mixtures of materials with different compressibilities, which are typically problematic for traditional affine transform methods. We have had success applying this method to |
| 220 |
+ |
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. |
| 221 |
+ |
|
| 222 |
+ |
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% |
| 223 |
+ |
% SIMULATION PROTOCOL |
| 224 |
+ |
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% |
| 225 |
|
\subsection{Simulation protocol} |
| 226 |
|
|
| 227 |
+ |
Systems containing liquids were run under moderate pressure ($\sim$ 5 atm) to avoid the formation of a substantial vapor phase. Thermal coupling to the Langevin Hull external bath was turned off to avoid interference with any imposed flux. |
| 228 |
|
|
| 229 |
|
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% |
| 230 |
|
% FORCE FIELD PARAMETERS |
| 253 |
|
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% |
| 254 |
|
\subsection{Thermal conductivities} |
| 255 |
|
|
| 256 |
< |
The thermal conductivity, $\lambda$, can be calculated using the imposed kinetic energy flux, $J_r$, and thermal gradient, $\frac{\partial T}{\partial r}$ measured from the resulting temperature profile |
| 256 |
> |
Fourier's Law of heat conduction in radial coordinates yields an expression for the heat flow between the concentric spherical RNEMD shells: |
| 257 |
|
|
| 258 |
|
\begin{equation} |
| 259 |
< |
J_r = -\lambda \frac{\partial T}{\partial r} |
| 259 |
> |
q_r = - \frac{4 \pi \lambda (T_b - T_a)}{\frac{1}{r_a} - \frac{1}{r_b}} |
| 260 |
> |
\label{eq:Q} |
| 261 |
> |
\end{equation} |
| 262 |
> |
|
| 263 |
> |
where $\lambda$ is the thermal conductivity, and $T_{a,b}$ and $r_{a,b}$ are the temperatures and radii of the two RNEMD regions, respectively. |
| 264 |
> |
|
| 265 |
> |
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$: |
| 266 |
> |
|
| 267 |
> |
\begin{equation} |
| 268 |
> |
\lambda = \frac{r_a}{r_b} \frac{q_r}{4 \pi \langle \frac{dT}{dr} \rangle} |
| 269 |
> |
\label{eq:lambda} |
| 270 |
|
\end{equation} |
| 271 |
|
|
| 272 |
+ |
The rate of heat transfer between the two RNEMD regions is the amount of transferred kinetic energy over the length of the simulation, t |
| 273 |
+ |
|
| 274 |
+ |
\begin{equation} |
| 275 |
+ |
q_r = \frac{KE}{t} |
| 276 |
+ |
\label{eq:heat} |
| 277 |
+ |
\end{equation} |
| 278 |
+ |
|
| 279 |
|
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% |
| 280 |
|
% INTERFACIAL THERMAL CONDUCTANCE |
| 281 |
|
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% |
| 282 |
|
\subsection{Interfacial thermal conductance} |
| 283 |
|
|
| 284 |
< |
A thermal flux is created using VSS-RNEMD moves, and the resulting temperature |
| 284 |
> |
A thermal flux is created using VSS-RNEMD moves, and the temperature in each of the radial shells is recorded. The resulting temperature |
| 285 |
|
profiles are analyzed to yield information about the interfacial thermal |
| 286 |
|
conductance. As the simulation progresses, the VSS moves are performed on a regular basis, and |
| 287 |
|
the system develops a thermal or velocity gradient in response to the applied |
| 288 |
< |
flux. A definition of the interfacial conductance in terms of a temperature difference ($\Delta T$) measured at $r_0$, |
| 288 |
> |
flux. We can treat the temperature in each radial shell as discrete, making the thermal conductance, $G$, of each shell the inverse of its Kapitza resistance, $R_K$. To model the thermal conductance across an interface (or multiple interfaces) it is useful to consider the shells as resistors wired in series. The total resistance of the shells is then additive, and the interfacial thermal conductance is the inverse of the total Kapitza resistance. The thermal resistance of each shell is |
| 289 |
> |
|
| 290 |
|
\begin{equation} |
| 291 |
< |
G = \frac{J_r}{\Delta T_{r_0}}, \label{eq:G} |
| 291 |
> |
R_K = \frac{1}{q_r} \Delta T 4 \pi r^2 |
| 292 |
> |
\label{eq:RK} |
| 293 |
|
\end{equation} |
| 212 |
– |
is useful once the RNEMD approach has generated a |
| 213 |
– |
stable temperature gap across the interface. |
| 294 |
|
|
| 295 |
+ |
making the total resistance of two neighboring shells |
| 296 |
+ |
|
| 297 |
+ |
\begin{equation} |
| 298 |
+ |
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 ] |
| 299 |
+ |
\label{eq:Rtotal} |
| 300 |
+ |
\end{equation} |
| 301 |
+ |
|
| 302 |
|
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% |
| 303 |
|
% INTERFACIAL FRICTION |
| 304 |
|
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% |
| 305 |
|
\subsection{Interfacial friction} |
| 306 |
|
|
| 307 |
< |
The slip length, $\delta$, is defined as the ratio of the interfacial friction coefficient, $\kappa$, and dynamic solvent viscosity, $\eta$ |
| 307 |
> |
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 |
| 308 |
|
|
| 309 |
|
\begin{equation} |
| 310 |
< |
\delta = \frac{\eta}{\kappa} |
| 310 |
> |
\Xi^{rr} = 8 \pi \eta r^3 |
| 311 |
> |
\label{eq:Xistick}. |
| 312 |
|
\end{equation} |
| 313 |
|
|
| 314 |
< |
and is measured from a radial angular momentum profile generated from a nonperiodic VSS-RNEMD simulation. |
| 314 |
> |
where $\eta$ is the dynamic viscosity of the surrounding solvent. An $\eta$ value for UA hexane under these particular temperature and pressure conditions was determined by applying a traditional VSS-RNEMD linear momentum flux to a periodic box of solvent. |
| 315 |
|
|
| 228 |
– |
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 |
| 229 |
– |
|
| 230 |
– |
\begin{equation} |
| 231 |
– |
\Xi = 8 \pi \eta r^3 \label{eq:Xi}. |
| 232 |
– |
\end{equation} |
| 233 |
– |
|
| 234 |
– |
where $\eta$ is the dynamic viscosity of the surrounding solvent and was determined for TraPPE-UA hexane under these condition by applying a traditional VSS-RNEMD linear momentum flux to a periodic box of solvent. |
| 235 |
– |
|
| 316 |
|
For general ellipsoids with semiaxes $a$, $b$, and $c$, Perrin's extension of Stokes' law provides exact solutions for symmetric prolate $(a \geq b = c)$ and oblate $(a < b = c)$ ellipsoids. For simplicity, we define a Perrin Factor, $S$, |
| 317 |
|
|
| 318 |
|
\begin{equation} |
| 319 |
< |
S = \frac{2}{\sqrt{a^2 - b^2}} ln \left[ \frac{a + \sqrt{a^2 - b^2}}{b} \right]. \label{eq:S} |
| 319 |
> |
S = \frac{2}{\sqrt{a^2 - b^2}} ln \left[ \frac{a + \sqrt{a^2 - b^2}}{b} \right]. |
| 320 |
> |
\label{eq:S} |
| 321 |
|
\end{equation} |
| 322 |
|
|
| 323 |
|
For a prolate ellipsoidal rod, demonstrated here, the rotational resistance tensor $\Xi$ is a $3 \times 3$ diagonal matrix with elements |
| 324 |
|
\begin{equation} |
| 325 |
|
\Xi^{rr}_a = \frac{32 \pi}{2} \eta \frac{ \left( a^2 - b^2 \right) b^2}{2a - b^2 S} |
| 326 |
|
\label{eq:Xia} |
| 327 |
+ |
\end{equation}\vspace{-0.45in}\\ |
| 328 |
+ |
\begin{equation} |
| 329 |
+ |
\Xi^{rr}_{b,c} = \frac{32 \pi}{2} \eta \frac{ \left( a^4 - b^4 \right)}{ \left( 2a^2 - b^2 \right)S - 2a}. |
| 330 |
+ |
\label{eq:Xibc} |
| 331 |
|
\end{equation} |
| 332 |
+ |
|
| 333 |
+ |
The effective rotational friction coefficient at the interface can be extracted from nonperiodic VSS-RNEMD simulations quite easily using the applied torque ($\tau$) and the observed angular velocity of the gold structure ($\omega_{Au}$) |
| 334 |
+ |
|
| 335 |
|
\begin{equation} |
| 336 |
< |
\Xi^{rr}_{b,c} = \frac{32 \pi}{2} \eta \frac{ \left( a^4 - b^4 \right)}{ \left( 2a^2 - b^2 \right)S - 2a}. \label{eq:Xibc} |
| 336 |
> |
\Xi^{rr}_{\mathit{eff}} = \frac{\tau}{\omega_{Au}} |
| 337 |
> |
\label{eq:Xieff} |
| 338 |
|
\end{equation} |
| 339 |
|
|
| 340 |
< |
% However, the friction between hexane solvent and gold more likely falls within ``slip'' boundary conditions. Hu and Zwanzig\cite{Zwanzig} investigated the rotational friction coefficients for spheroids under slip boundary conditions and obtained numerial results for a scaling factor as a function of $\tau$, the ratio of the shorter semiaxes and the longer semiaxis of the spheroid. For the sphere and prolate ellipsoid shown here, the values of $\tau$ are $1$ and $0.3939$, respectively. According to the values tabulated by Hu and Zwanzig, the sphere friction coefficient approaches $0$, while the ellipsoidal friction coefficient must be scaled by a factor of $0.880$ to account for the reduced interfacial friction under ``slip'' boundary conditions. |
| 340 |
> |
The applied torque required to overcome the interfacial friction and maintain constant rotation of the gold is |
| 341 |
|
|
| 342 |
< |
Another useful quantity is the friction factor $f$, or the friction coefficient of a non-spherical shape divided by the friction coefficient of a sphere of equivalent volume. The nanoparticle and ellipsoidal nanorod dimensions used here were specifically chosen to create structures with equal volumes. |
| 342 |
> |
\begin{equation} |
| 343 |
> |
\tau = \frac{L}{2 t} |
| 344 |
> |
\label{eq:tau} |
| 345 |
> |
\end{equation} |
| 346 |
|
|
| 347 |
+ |
where $L$ is the total angular momentum exchanged between the two RNEMD regions and $t$ is the length of the simulation. |
| 348 |
+ |
|
| 349 |
+ |
Previous VSS-RNEMD simulations of the interfacial friction of the planar Au(111) / hexane interface have shown that the interface exists within ``slip'' boundary conditions.\cite{Kuang2012} 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 $\tau$, the ratio of the shorter semiaxes and the longer semiaxis of the spheroid. For the sphere and prolate ellipsoid shown here, the values of $\tau$ are $1$ and $0.3939$, respectively. According to the values tabulated by Hu and Zwanzig, $\Xi^{rr}_{\mathit{slip}}$ for any sphere approaches $0$, while the ellipsoidal $\Xi^{rr}_{\mathit{slip}}$ is the analytical $\Xi^{rr}_{\mathit{stick}}$ result scaled by a factor of $0.359$ to account for the reduced interfacial friction under ``slip'' boundary conditions. |
| 350 |
+ |
|
| 351 |
|
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% |
| 352 |
|
% **TESTS AND APPLICATIONS** |
| 353 |
|
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% |
| 358 |
|
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% |
| 359 |
|
\subsection{Thermal conductivities} |
| 360 |
|
|
| 361 |
< |
Calculated values for the thermal conductivity of a 40 \AA$~$ radius gold nanoparticle (15707 atoms) at different kinetic energy flux values are shown in Table \ref{table:goldconductivity}. For these calculations, the hot and cold slabs were excluded from the linear regression of the thermal gradient. The measured linear slope $\langle dT / dr \rangle$ is linearly dependent on the applied kinetic energy flux $J_r$. Calculated thermal conductivity values compare well with previous bulk QSC values of 1.08 -- 1.26 W m$^{-1}$ K$^{-1}$\cite{Kuang2010}, though still significantly lower than the experimental value of 320 W m$^{-1}$ K$^{-1}$, as the QSC force field neglects significant electronic contributions to heat conduction. The small increase relative to previous simulated bulk values is due to a slight increase in gold density -- as expected, an increase in density results in higher thermal conductivity values. The increased density is a result of nanoparticle curvature relative to an infinite bulk slab, which introduces surface tension that increases ambient density. |
| 361 |
> |
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:goldTC}. For these calculations, the hot and cold slabs were excluded from the linear regression of the thermal gradient. The measured linear slope $\langle \frac{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 {\footnotesize W / m $\cdot$ K}\cite{Kuang2010}, though still significantly lower than the experimental value of 320 {\footnotesize W / m $\cdot$ K}, as the QSC force field neglects significant electronic contributions to heat conduction. |
| 362 |
|
|
| 363 |
+ |
% 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. |
| 364 |
+ |
|
| 365 |
|
\begin{longtable}{ccc} |
| 366 |
|
\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.} |
| 367 |
|
\\ \hline \hline |
| 368 |
< |
{$J_r$} & {$\langle dT / dr \rangle$} & {$\boldsymbol \lambda$}\\ |
| 369 |
< |
{\small(kcal fs$^{-1}$ \AA$^{-2}$)} & {\small(K \AA$^{-1}$)} & {\small(W m$^{-1}$ K$^{-1}$)}\\ \hline |
| 370 |
< |
3.25$\times 10^{-6}$ & 0.11435 & 1.9753 \\ |
| 371 |
< |
6.50$\times 10^{-6}$ & 0.2324 & 1.9438 \\ |
| 372 |
< |
1.30$\times 10^{-5}$ & 0.44922 & 2.0113 \\ |
| 373 |
< |
3.25$\times 10^{-5}$ & 1.1802 & 1.9139 \\ |
| 374 |
< |
6.50$\times 10^{-5}$ & 2.339 & 1.9314 |
| 368 |
> |
{$J_r$} & {$\langle \frac{dT}{dr} \rangle$} & {$\boldsymbol \lambda$}\\ |
| 369 |
> |
{\footnotesize(kcal / fs $\cdot$ \AA$^{2}$)} & {\footnotesize(K / \AA)} & {\footnotesize(W / m $\cdot$ K)}\\ \hline |
| 370 |
> |
3.25$\times 10^{-6}$ & 0.11435 & 1.0143 \\ |
| 371 |
> |
6.50$\times 10^{-6}$ & 0.2324 & 0.9982 \\ |
| 372 |
> |
1.30$\times 10^{-5}$ & 0.44922 & 1.0328 \\ |
| 373 |
> |
3.25$\times 10^{-5}$ & 1.1802 & 0.9828 \\ |
| 374 |
> |
6.50$\times 10^{-5}$ & 2.339 & 0.9918 \\ |
| 375 |
> |
\hline |
| 376 |
> |
This work & & 1.0040 |
| 377 |
|
\\ \hline \hline |
| 378 |
< |
\label{table:goldconductivity} |
| 378 |
> |
\label{table:goldTC} |
| 379 |
|
\end{longtable} |
| 380 |
|
|
| 381 |
< |
Calculated values for the thermal conductivity of a cluster of 6912 SPC/E water molecules are shown in Table \ref{table:waterconductivity}. |
| 381 |
> |
Calculated values for the thermal conductivity of a cluster of 6912 SPC/E water molecules are shown in Table \ref{table:waterTC}. As with the gold nanoparticle thermal conductivity calculations, the RNEMD regions were excluded from the $\langle \frac{dT}{dr} \rangle$ fit. Again, the measured slope is linearly dependent on the applied kinetic energy flux $J_r$. The average calculated thermal conductivity from this work, $0.8841$ {\footnotesize W / m $\cdot$ K}, compares very well to previous nonequilibrium molecular dynamics results\cite{Romer2012, Zhang2005} and experimental values.\cite{WagnerKruse} |
| 382 |
|
|
| 383 |
|
\begin{longtable}{ccc} |
| 384 |
< |
\caption{Calculated thermal conductivity of a cluster of 6912 SPC/E water molecules. Calculations were performed at 300 K and ambient density.} |
| 384 |
> |
\caption{Calculated thermal conductivity of a cluster of 6912 SPC/E water molecules. Calculations were performed at 300 K and 5 atm.} |
| 385 |
|
\\ \hline \hline |
| 386 |
< |
{$J_r$} & {$\langle dT / dr \rangle$} & {$\boldsymbol \lambda$}\\ |
| 387 |
< |
{\small(kcal fs$^{-1}$ \AA$^{-2}$)} & {\small(K \AA$^{-1}$)} & {\small(W m$^{-1}$ K$^{-1}$)}\\ \hline |
| 388 |
< |
1.00$\times 10^{-5}$ & 0.38683 & 1.79665486\\ |
| 389 |
< |
3.00$\times 10^{-5}$ & 1.1643 & 1.79077557\\ |
| 390 |
< |
6.00$\times 10^{-5}$ & 2.2262 & 1.87314707\\ |
| 391 |
< |
\hline \hline |
| 392 |
< |
\label{table:waterconductivity} |
| 386 |
> |
{$J_r$} & {$\langle \frac{dT}{dr} \rangle$} & {$\boldsymbol \lambda$}\\ |
| 387 |
> |
{\footnotesize(kcal / fs $\cdot$ \AA$^{2}$)} & {\footnotesize(K / \AA)} & {\footnotesize(W / m $\cdot$ K)}\\ \hline |
| 388 |
> |
1$\times 10^{-5}$ & 0.38683 & 0.8698 \\ |
| 389 |
> |
3$\times 10^{-5}$ & 1.1643 & 0.9098 \\ |
| 390 |
> |
6$\times 10^{-5}$ & 2.2262 & 0.8727 \\ |
| 391 |
> |
\hline |
| 392 |
> |
This work & & 0.8841 \\ |
| 393 |
> |
Zhang, et al\cite{Zhang2005} & & 0.81 \\ |
| 394 |
> |
R$\ddot{\mathrm{o}}$mer, et al\cite{Romer2012} & & 0.87 \\ |
| 395 |
> |
Experiment\cite{WagnerKruse} & & 0.61 |
| 396 |
> |
\\ \hline \hline |
| 397 |
> |
\label{table:waterTC} |
| 398 |
|
\end{longtable} |
| 399 |
|
|
| 400 |
|
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% |
| 402 |
|
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% |
| 403 |
|
\subsection{Interfacial thermal conductance} |
| 404 |
|
|
| 405 |
+ |
Calculated interfacial thermal conductance ($G$) |
| 406 |
+ |
|
| 407 |
|
\begin{longtable}{ccc} |
| 408 |
< |
\caption{Caption.} |
| 408 |
> |
\caption{Calculated interfacial thermal conductance ($G$) values for gold nanoparticles of varying radii solvated in explicit TraPPE-UA hexane. The nanoparticle $G$ values are compared to previous results for a gold slab in TraPPE-UA hexane, revealing increased interfacial thermal conductance for non-planar interfaces.} |
| 409 |
|
\\ \hline \hline |
| 410 |
< |
{Nanoparticle Radius} & {$\boldsymbol \lambda$}\\ |
| 411 |
< |
{\small(\AA)} & {\small(W m$^{-1}$ K$^{-1}$)}\\ \hline |
| 412 |
< |
20 & 59.66\\ |
| 413 |
< |
30 & 57.88\\ |
| 414 |
< |
40 & \\ |
| 415 |
< |
$\infty$ & \\ |
| 416 |
< |
\hline \hline |
| 417 |
< |
\label{table:waterconductivity} |
| 410 |
> |
{Nanoparticle Radius} & {$G$}\\ |
| 411 |
> |
{\footnotesize(\AA)} & {\footnotesize(MW / m$^{2}$ $\cdot$ K)}\\ \hline |
| 412 |
> |
20 & {47.1} \\ |
| 413 |
> |
30 & {45.4} \\ |
| 414 |
> |
40 & {46.5} \\ |
| 415 |
> |
slab & {30.2} |
| 416 |
> |
\\ \hline \hline |
| 417 |
> |
\label{table:interfacialconductance} |
| 418 |
|
\end{longtable} |
| 419 |
|
|
| 420 |
|
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% |
| 422 |
|
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% |
| 423 |
|
\subsection{Interfacial friction} |
| 424 |
|
|
| 425 |
< |
Table \ref{table:interfacialfrictionstick} shows the calculated interfacial friction coefficients $\kappa$ for spherical gold nanoparticles and a prolate ellipsoidal gold nanorod in TraPPE-UA hexane. An angular momentum flux was applied between the A and B regions defined as the gold structure and hexane molecules beyond a certain radius, respectively. The resulting angular velocity gradient resulted in the gold structure rotating about the prescribed axis within the solvent. |
| 425 |
> |
Table \ref{table:couple} shows the calculated rotational friction coefficients $\Xi^{rr}$ for spherical gold nanoparticles and a prolate ellipsoidal gold nanorod in TraPPE-UA hexane. An angular momentum flux was applied between the A and B regions defined as the gold structure and hexane molecules beyond a certain radius, respectively. The resulting angular velocity gradient causes the gold structure to rotate about the prescribed axis. |
| 426 |
|
|
| 427 |
|
\begin{longtable}{lccccc} |
| 428 |
< |
\caption{Calculated ``stick'' interfacial friction coefficients ($\kappa$) and friction factors ($f$) of gold nanostructures solvated in TraPPE-UA hexane at 230 K. The ellipsoid is oriented with the long axis along the $z$ direction.} |
| 428 |
> |
\caption{Comparison of rotational friction coefficients under ideal ``slip'' ($\Xi^{rr}_{\mathit{slip}}$) and ``stick'' conditions ($\Xi^{rr}_{\mathit{stick}}$) and effective rotational friction coefficients ($\Xi^{rr}_{\mathit{eff}}$) of gold nanostructures solvated in TraPPE-UA hexane at 230 K. The ellipsoid is oriented with the long axis along the $z$ direction.} |
| 429 |
|
\\ \hline \hline |
| 430 |
< |
{Structure} & {Axis of Rotation} & {$\kappa_{VSS}$} & {$\kappa_{calc}$} & {$f_{VSS}$} & {$f_{calc}$}\\ |
| 431 |
< |
{} & {} & {\small($10^{-29}$ Pa s m$^{3}$)} & {\small($10^{-29}$ Pa s m$^{3}$)} & {} & {}\\ \hline |
| 432 |
< |
{Sphere (r = 20 \AA)} & {$x = y = z$} & {} & {6.13239} & {1} & {1}\\ |
| 433 |
< |
{Sphere (r = 30 \AA)} & {$x = y = z$} & {} & {20.6968} & {1} & {1}\\ |
| 434 |
< |
{Sphere (r = 40 \AA)} & {$x = y = z$} & {} & {49.0591} & {1} & {1}\\ |
| 435 |
< |
{Prolate Ellipsoid} & {$x = y$} & {} & {8.22846} & {} & {1.92477}\\ |
| 436 |
< |
{Prolate Ellipsoid} & {$z$} & {} & {3.28634} & {} & {0.768726}\\ |
| 437 |
< |
\hline \hline |
| 438 |
< |
\label{table:interfacialfrictionstick} |
| 430 |
> |
{Structure} & {Axis of Rotation} & {$\Xi^{rr}_{\mathit{slip}}$} & {$\Xi^{rr}_{\mathit{eff}}$} & {$\Xi^{rr}_{\mathit{stick}}$} & {$\Xi^{rr}_{\mathit{eff}}$ / $\Xi^{rr}_{\mathit{stick}}$}\\ |
| 431 |
> |
{} & {} & {\footnotesize(amu A$^2$ / fs)} & {\footnotesize(amu A$^2$ / fs)} & {\footnotesize(amu A$^2$ / fs)} & \\ \hline |
| 432 |
> |
Sphere (r = 20 \AA) & {$x = y = z$} & 0 & {2386} & {3314} & {0.720}\\ |
| 433 |
> |
Sphere (r = 30 \AA) & {$x = y = z$} & 0 & {8415} & {11749} & {0.716}\\ |
| 434 |
> |
Sphere (r = 40 \AA) & {$x = y = z$} & 0 & {47544} & {34464} & {1.380}\\ |
| 435 |
> |
Prolate Ellipsoid & {$x = y$} & {1792} & {3128} & {4991} & {0.627}\\ |
| 436 |
> |
Prolate Ellipsoid & {$z$} & {716} & {1590} & {1993} & {0.798} |
| 437 |
> |
\\ \hline \hline |
| 438 |
> |
\label{table:couple} |
| 439 |
|
\end{longtable} |
| 440 |
|
|
| 441 |
< |
% \begin{longtable}{lccc} |
| 335 |
< |
% \caption{Calculated ``slip'' interfacial friction coefficients ($\kappa$) and friction factors ($f$) of gold nanostructures solvated in TraPPE-UA hexane. The ellipsoid is oriented with the long axis along the $z$ direction.} |
| 336 |
< |
% \\ \hline \hline |
| 337 |
< |
% {Structure} & {Axis of rotation} & {$\kappa_{VSS}$} & {$\kappa_{calc}$}\\ |
| 338 |
< |
% {} & {} & {\small($10^{-29}$ Pa s m$^{3}$)} & {\small($10^{-29}$ Pa s m$^{3}$)}\\ \hline |
| 339 |
< |
% {Sphere} & {$x = y = z$} & {} & {0}\\ |
| 340 |
< |
% {Prolate Ellipsoid} & {$x = y$} & {} & {0}\\ |
| 341 |
< |
% {Prolate Ellipsoid} & {$z$} & {} & {7.9295392}\\ \hline \hline |
| 342 |
< |
% \label{table:interfacialfrictionslip} |
| 343 |
< |
% \end{longtable} |
| 441 |
> |
The results for $\Xi^{rr}_{\mathit{eff}}$ show that, contrary to the flat Au(111) / hexane interface, gold structures solvated by hexane do not exist in the ``slip'' boundary conditions. At this length scale, the nanostructures are not perfect spheroids due to atomic `roughening' of the surface and therefore experience increased interfacial friction which deviates 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 may at first seem counterintuitive. However, the `propellor' motion caused by rotation about 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 more than the ``stick'' limit of interfacial friction, which may be a consequence of surface features or anomalous solvent behaviors that are not fully understood at this time. |
| 442 |
|
|
| 443 |
|
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% |
| 444 |
|
% **DISCUSSION** |
