| 83 |
|
\section{Introduction} |
| 84 |
|
|
| 85 |
|
Non-equilibrium Molecular Dynamics (NEMD) methods impose a temperature or velocity {\it gradient} on a |
| 86 |
< |
system,\cite{ASHURST:1975tg,Evans:1982zk,ERPENBECK:1984sp,MAGINN:1993hc,Berthier:2002ij,Evans:2002ai,Schelling:2 |
| 87 |
< |
002dp,PhysRevA.34.1449,JiangHao_jp802942v} and use linear response theory to connect the resulting thermal or |
| 88 |
< |
momentum flux to transport coefficients of bulk materials. However, for heterogeneous systems, such as phase |
| 86 |
> |
system,\cite{ASHURST:1975tg,Evans:1982zk,ERPENBECK:1984sp,MAGINN:1993hc,Berthier:2002ij,Evans:2002ai,Schelling:2002dp,PhysRevA.34.1449,JiangHao_jp802942v} and use linear response theory to connect the resulting thermal or |
| 87 |
> |
momentum flux to transport coefficients of bulk materials. However, for heterogeneous systems, such as phase |
| 88 |
|
boundaries or interfaces, it is often unclear what shape of gradient should be imposed at the boundary between |
| 89 |
|
materials. |
| 90 |
|
|
| 99 |
|
properties, and they have become widely used to compute thermal and mechanical transport in both homogeneous |
| 100 |
|
liquids and solids~\cite{MullerPlathe:1997xw,ISI:000080382700030,Maginn:2010} as well as heterogeneous |
| 101 |
|
interfaces.\cite{garde:nl2005,garde:PhysRevLett2009,kuang:AuThl} |
| 102 |
+ |
|
| 103 |
+ |
The strengths of specific algorithms for imposing the flux between two |
| 104 |
+ |
different slabs of the simulation cell has been the subject of some |
| 105 |
+ |
renewed interest. The original RNEMD approach used kinetic energy or |
| 106 |
+ |
momentum exchange between particles in the two slabs, either through |
| 107 |
+ |
direct swapping of momentum vectors or via virtual elastic collisions |
| 108 |
+ |
between atoms in the two regions. There have been recent |
| 109 |
+ |
methodological advances which involve scaling all particle velocities |
| 110 |
+ |
in both slabs. Constraint equations are simultaneously imposed to |
| 111 |
+ |
require the simulation to conserve both total energy and total linear |
| 112 |
+ |
momentum. The most recent and simplest of the velocity scaling |
| 113 |
+ |
approaches allows for simultaneous shearing (to provide viscosity |
| 114 |
+ |
estimates) as well as scaling (to provide information about thermal |
| 115 |
+ |
conductivity). |
| 116 |
+ |
|
| 117 |
+ |
To date, however, the RNEMD methods have only been usable in periodic |
| 118 |
+ |
simulation cells where the exchange regions are physically separated |
| 119 |
+ |
along one of the axes of the simulation cell. This limits the |
| 120 |
+ |
applicability to infinite planar interfaces. |
| 121 |
|
|
| 122 |
+ |
In order to model steady-state non-equilibrium distributions for |
| 123 |
+ |
curved surfaces (e.g. hot nanoparticles in contact with colder |
| 124 |
+ |
solvent), or for regions that are not planar slabs, the method |
| 125 |
+ |
requires some generalization for non-parallel exchange regions. In |
| 126 |
+ |
the following sections, we present the Velocity Shearing and Scaling |
| 127 |
+ |
(VSS) RNEMD algorithm which has been explicitly designed for |
| 128 |
+ |
non-periodic simulations, and use the method to compute some thermal |
| 129 |
+ |
transport and solid-liquid friction at the surfaces of spherical and |
| 130 |
+ |
ellipsoidal nanoparticles, and discuss how the method can be extended |
| 131 |
+ |
to provide other kinds of non-equilibrium fluxes. |
| 132 |
|
|
| 133 |
|
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% |
| 134 |
|
% **METHODOLOGY** |
| 135 |
|
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% |
| 136 |
|
\section{Velocity Shearing and Scaling (VSS) for non-periodic systems} |
| 137 |
|
|
| 138 |
< |
The VSS-RNEMD approach uses a series of simultaneous velocity shearing and scaling exchanges between the two |
| 139 |
< |
slabs.\cite{Kuang2012} This method imposes energy and momentum conservation constraints while simultaneously |
| 140 |
< |
creating a desired flux between the two slabs. These constraints ensure that all configurations are sampled |
| 141 |
< |
from the same microcanonical (NVE) ensemble. |
| 138 |
> |
The periodic VSS-RNEMD approach uses a series of simultaneous velocity |
| 139 |
> |
shearing and scaling exchanges between the two slabs.\cite{Kuang2012} |
| 140 |
> |
This method imposes energy and momentum conservation constraints while |
| 141 |
> |
simultaneously creating a desired flux between the two slabs. These |
| 142 |
> |
constraints ensure that all configurations are sampled from the same |
| 143 |
> |
microcanonical (NVE) ensemble. |
| 144 |
|
|
| 145 |
|
\begin{figure} |
| 146 |
|
\includegraphics[width=\linewidth]{figures/npVSS} |
| 276 |
|
the Langevin Hull methodology\cite{Vardeman2011} -- before any non-equilibrium methods were introduced. For |
| 277 |
|
heterogeneous systems, the gold nanoparticles and ellipsoid were first created from a bulk lattice and |
| 278 |
|
thermally equilibrated before being solvated in hexane. Packmol\cite{packmol} was used to solvate previously |
| 279 |
< |
equilibrated gold nanostructures within a droplet of hexane. |
| 279 |
> |
equilibrated gold nanostructures within a spherical droplet of hexane. |
| 280 |
|
|
| 281 |
|
Once fully equilibrated, a thermal or angular momentum flux was applied for 1 - 2 |
| 282 |
|
ns, until a stable temperature or angular velocity gradient had developed. Systems containing liquids were run |
| 283 |
< |
under moderate pressure (5 atm) and temperatures (230 K) to avoid the formation of a substantial vapor phase. 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 |
| 283 |
> |
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 |
| 284 |
|
imposed flux. |
| 285 |
|
|
| 286 |
|
To stabilize the gold nanoparticle under the imposed angular momentum flux we altered the gold atom at the |
| 287 |
|
designated coordinate origin to have $10,000$ times its original mass. The nonbonded interactions remain |
| 288 |
< |
unchanged. The heavy atom is excluded from the angular momentum exchange. For rotation of the ellipsoid about |
| 288 |
> |
unchanged and the heavy atom is excluded from the angular momentum exchange. For rotation of the ellipsoid about |
| 289 |
|
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. |
| 290 |
|
|
| 291 |
|
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% |
| 356 |
|
\end{equation} |
| 357 |
|
|
| 358 |
|
This series can be expanded for any number of adjacent shells, allowing for the calculation of the interfacial |
| 359 |
< |
thermal conductance for interfaces of significant thickness, such as self-assembled ligand monolayers on a |
| 359 |
> |
thermal conductance for interfaces of considerable thickness, such as self-assembled ligand monolayers on a |
| 360 |
|
metal surface. |
| 361 |
|
|
| 362 |
|
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% |
| 387 |
|
momentum flux to a periodic box of solvent. |
| 388 |
|
|
| 389 |
|
For general ellipsoids with semiaxes $a$, $b$, and $c$, Perrin's extension of Stokes' law provides exact |
| 390 |
< |
solutions for symmetric prolate $(a \geq b = c)$ and oblate $(a < b = c)$ ellipsoids. For simplicity, we define |
| 390 |
> |
solutions for symmetric prolate $(a \geq b = c)$ and oblate $(a < b = c)$ ellipsoids under ideal ``stick'' conditions. For simplicity, we define |
| 391 |
|
a Perrin Factor, $S$, |
| 392 |
|
|
| 393 |
|
\begin{equation} |
| 405 |
|
\label{eq:Xibc} |
| 406 |
|
\end{equation} |
| 407 |
|
|
| 408 |
< |
The effective rotational friction coefficient at the interface can be extracted from non-periodic VSS-RNEMD simulations quite easily using the applied torque ($\tau$) and the observed angular velocity of the gold structure ($\omega_{Au}$) |
| 408 |
> |
corresponding to rotation about the long axis ($a$), and each of the equivalent short axes ($b$ and $c$), respectively. |
| 409 |
|
|
| 410 |
+ |
Previous VSS-RNEMD simulations of the interfacial friction of the planar Au(111) / hexane interface have shown |
| 411 |
+ |
that the interface exists within ``slip'' boundary conditions.\cite{Kuang2012} Hu and Zwanzig\cite{Zwanzig} |
| 412 |
+ |
investigated the rotational friction coefficients for spheroids under slip boundary conditions and obtained |
| 413 |
+ |
numerial results for a scaling factor to be applied to $\Xi^{rr}_{\mathit{stick}}$ as a function of $\tau$, the |
| 414 |
+ |
ratio of the shorter semiaxes and the longer semiaxis of the spheroid. For the sphere and prolate ellipsoid |
| 415 |
+ |
shown here, the values of $\tau$ are $1$ and $0.3939$, respectively. Under ``slip'' conditions, |
| 416 |
+ |
$\Xi^{rr}_{\mathit{slip}}$ for any sphere and rotation of the prolate ellipsoid about its long axis approaches |
| 417 |
+ |
$0$, as no solvent is displaced by either of these rotations. $\Xi^{rr}_{\mathit{slip}}$ for rotation of the |
| 418 |
+ |
prolate ellipsoid about its short axis is $35.9\%$ of the analytical $\Xi^{rr}_{\mathit{stick}}$ result, |
| 419 |
+ |
accounting for the reduced interfacial friction under ``slip'' boundary conditions. |
| 420 |
+ |
|
| 421 |
+ |
The effective rotational friction coefficient, $\Xi^{rr}_{\mathit{eff}}$ at the interface can be extracted from non-periodic VSS-RNEMD simulations quite easily using the applied torque ($\tau$) and the observed angular velocity of the gold structure ($\omega_{Au}$) |
| 422 |
+ |
|
| 423 |
|
\begin{equation} |
| 424 |
|
\Xi^{rr}_{\mathit{eff}} = \frac{\tau}{\omega_{Au}} |
| 425 |
|
\label{eq:Xieff} |
| 434 |
|
|
| 435 |
|
where $L$ is the total angular momentum exchanged between the two RNEMD regions and $t$ is the length of the simulation. |
| 436 |
|
|
| 394 |
– |
Previous VSS-RNEMD simulations of the interfacial friction of the planar Au(111) / hexane interface have shown |
| 395 |
– |
that the interface exists within ``slip'' boundary conditions.\cite{Kuang2012} Hu and Zwanzig\cite{Zwanzig} |
| 396 |
– |
investigated the rotational friction coefficients for spheroids under slip boundary conditions and obtained |
| 397 |
– |
numerial results for a scaling factor to be applied to $\Xi^{rr}_{\mathit{stick}}$ as a function of $\tau$, the |
| 398 |
– |
ratio of the shorter semiaxes and the longer semiaxis of the spheroid. For the sphere and prolate ellipsoid |
| 399 |
– |
shown here, the values of $\tau$ are $1$ and $0.3939$, respectively. Under ``slip'' conditions, |
| 400 |
– |
$\Xi^{rr}_{\mathit{slip}}$ for any sphere and rotation of the prolate ellipsoid about its long axis approaches |
| 401 |
– |
$0$, as no solvent is displaced by either of these rotations. $\Xi^{rr}_{\mathit{slip}}$ for rotation of the |
| 402 |
– |
prolate ellipsoid about its short axis is $35.9\%$ of the analytical $\Xi^{rr}_{\mathit{stick}}$ result, |
| 403 |
– |
accounting for the reduced interfacial friction under ``slip'' boundary conditions. |
| 404 |
– |
|
| 437 |
|
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% |
| 438 |
|
% **TESTS AND APPLICATIONS** |
| 439 |
|
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% |
| 465 |
|
\end{longtable} |
| 466 |
|
|
| 467 |
|
The measured linear slope $\langle \frac{dT}{dr} \rangle$ is linearly dependent on the applied kinetic energy |
| 468 |
< |
flux $J_r$. Calculated thermal conductivity values compare well with previous bulk QSC values of 1.08 -- 1.26 |
| 469 |
< |
{\footnotesize W / m $\cdot$ K}\cite{Kuang2010}, though still significantly lower than the experimental value |
| 438 |
< |
of 320 {\footnotesize W / m $\cdot$ K}, as the QSC force field neglects significant electronic contributions to |
| 468 |
> |
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 |
| 469 |
> |
of 320 W / m $\cdot$ K, as the QSC force field neglects significant electronic contributions to |
| 470 |
|
heat conduction. |
| 471 |
|
|
| 472 |
|
Calculated values for the thermal conductivity of a cluster of 6912 SPC/E water molecules are shown in Table |
| 491 |
|
\end{longtable} |
| 492 |
|
|
| 493 |
|
Again, the measured slope is linearly dependent on the applied kinetic energy flux $J_r$. The average |
| 494 |
< |
calculated thermal conductivity from this work, $0.8841$ {\footnotesize W / m $\cdot$ K}, compares very well to |
| 494 |
> |
calculated thermal conductivity from this work, $0.8841$ W / m $\cdot$ K, compares very well to |
| 495 |
|
previous non-equilibrium molecular dynamics results\cite{Romer2012, Zhang2005} and experimental |
| 496 |
|
values.\cite{WagnerKruse} |
| 497 |
|
|
| 504 |
|
surface solvated in TraPPE-UA hexane are shown in Table \ref{table:G}. |
| 505 |
|
|
| 506 |
|
\begin{longtable}{ccc} |
| 507 |
< |
\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 Au(111) interface in TraPPE-UA hexane.} |
| 507 |
> |
\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.} |
| 508 |
|
\\ \hline \hline |
| 509 |
|
{Nanoparticle Radius} & {$G$}\\ |
| 510 |
|
{\footnotesize(\AA)} & {\footnotesize(MW / m$^{2}$ $\cdot$ K)}\\ \hline |
| 520 |
|
The introduction of surface curvature increases the interfacial thermal conductance by a factor of |
| 521 |
|
approximately $1.5$ relative to the flat interface. There are no significant differences in the $G$ values for |
| 522 |
|
the varying nanoparticle sizes. It seems likely that for the range of nanoparticle sizes represented here, any |
| 523 |
< |
particle size effects are not evident. |
| 523 |
> |
particle size effects are not evident. The simulation of larger nanoparticles may demonstrate an approach to the $G$ value of a flat Au(111) slab but would require prohibitively costly numbers of atoms. |
| 524 |
|
|
| 525 |
|
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% |
| 526 |
|
% INTERFACIAL FRICTION |
| 547 |
|
\end{longtable} |
| 548 |
|
|
| 549 |
|
The results for $\Xi^{rr}_{\mathit{eff}}$ show that, contrary to the flat Au(111) / hexane interface, gold |
| 550 |
< |
structures solvated by hexane do not exist in the ``slip'' boundary conditions. At this length scale, the |
| 550 |
> |
structures solvated by hexane do not exist in the ``slip'' boundary condition. At this length scale, the |
| 551 |
|
nanostructures are not perfect spheroids due to atomic `roughening' of the surface and therefore experience |
| 552 |
|
increased interfacial friction which deviates from the ideal ``slip'' case. The 20 and 30 \AA$\,$ radius |
| 553 |
|
nanoparticles experience approximately 70\% of the ideal ``stick'' boundary interfacial friction. Rotation of |
| 554 |
|
the ellipsoid about its long axis more closely approaches the ``stick'' limit than rotation about the short |
| 555 |
< |
axis, which may at first seem counterintuitive. However, the `propellor' motion caused by rotation about the |
| 555 |
> |
axis, which at first seems counterintuitive. However, the `propellor' motion caused by rotation about the |
| 556 |
|
short axis may exclude solvent from the rotation cavity or move a sufficient amount of solvent along with the |
| 557 |
|
gold that a smaller interfacial friction is actually experienced. The largest nanoparticle (40 \AA$\,$ radius) |
| 558 |
|
appears to experience more than the ``stick'' limit of interfacial friction, which may be a consequence of |
