| 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/VSS} |
| 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,Kuang:2010uq} 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 |
|
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% |
| 123 |
|
% **METHODOLOGY** |
| 124 |
|
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% |
| 125 |
< |
\section{Methodology} |
| 125 |
> |
\section{Velocity Shearing and Scaling (VSS) 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{2012MolPh.110..691K} |
| 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 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 |
> |
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 |
> |
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}} \\ \\ |
| 155 |
> |
\mathbf{v}_i $~~~$\leftarrow & |
| 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 |
> |
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 |
> |
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{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 |
+ |
To simultaneously impose a thermal flux ($J_r$) between the shells we |
| 181 |
+ |
use energy conservation constraints, |
| 182 |
+ |
\begin{eqnarray} |
| 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 |
+ |
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$) and shear viscosity ($\eta$), |
| 208 |
+ |
\begin{equation} |
| 209 |
+ |
J_r = -\lambda \frac{\partial T}{\partial |
| 210 |
+ |
r} \hspace{2in} j_r(\mathbf{L}_z) = -\eta \frac{\partial |
| 211 |
+ |
\omega_z}{\partial r} |
| 212 |
+ |
\end{equation} |
| 213 |
+ |
of a liquid cluster. |
| 214 |
+ |
|
| 215 |
|
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% |
| 216 |
|
% NON-PERIODIC DYNAMICS |
| 217 |
|
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% |
| 237 |
|
of rectangular slabs. A temperature profile along the $r$ coordinate is created by recording the average temperature of concentric spherical shells. |
| 238 |
|
|
| 239 |
|
\begin{figure} |
| 240 |
< |
\center{\includegraphics[width=7in]{figures/VSS}} |
| 240 |
> |
\center{\includegraphics[width=7in]{figures/npVSS2}} |
| 241 |
|
\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.} |
| 242 |
|
\label{fig:VSS} |
| 243 |
|
\end{figure} |
| 412 |
|
\\ \hline \hline |
| 413 |
|
{$J_r$} & {$\langle dT / dr \rangle$} & {$\boldsymbol \lambda$}\\ |
| 414 |
|
{\small(kcal fs$^{-1}$ \AA$^{-2}$)} & {\small(K \AA$^{-1}$)} & {\small(W m$^{-1}$ K$^{-1}$)}\\ \hline |
| 415 |
< |
\\ \hline \hline |
| 415 |
> |
1.00$\times 10^{-5}$ & 0.38683 & 1.79665486\\ |
| 416 |
> |
3.00$\times 10^{-5}$ & 1.1643 & 1.79077557\\ |
| 417 |
> |
6.00$\times 10^{-5}$ & 2.2262 & 1.87314707\\ |
| 418 |
> |
\hline \hline |
| 419 |
|
\label{table:waterconductivity} |
| 420 |
|
\end{longtable} |
| 421 |
|
|
| 424 |
|
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% |
| 425 |
|
\subsection{Interfacial thermal conductance} |
| 426 |
|
|
| 427 |
+ |
\begin{longtable}{ccc} |
| 428 |
+ |
\caption{Caption.} |
| 429 |
+ |
\\ \hline \hline |
| 430 |
+ |
{Nanoparticle Radius} & $J_r$ & {G}\\ |
| 431 |
+ |
{\small(\AA)} & {\small(kcal fs$^{-1}$ \AA$^{-2}$)} & {\small(W m$^{-2}$ K$^{-1}$)}\\ \hline |
| 432 |
+ |
20 & & 59.66 \\ |
| 433 |
+ |
30 & & 57.88 \\ |
| 434 |
+ |
40 & & 37.48 \\ |
| 435 |
+ |
$\infty$ & & 30.2 \\ |
| 436 |
+ |
\hline \hline |
| 437 |
+ |
\label{table:interfacialconductance} |
| 438 |
+ |
\end{longtable} |
| 439 |
+ |
|
| 440 |
|
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% |
| 441 |
|
% INTERFACIAL FRICTION |
| 442 |
|
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% |
| 443 |
|
\subsection{Interfacial friction} |
| 444 |
|
|
| 445 |
< |
Table \ref{table:interfacialfrictionstick} shows the calculated interfacial friction coefficients $\kappa$ for a spherical gold nanoparticle and prolate ellipsoidal gold nanorod in TraPPE-UA hexane. An angular momentum flux was applied between the A and B regions defined as the gold structure and hexane molecules beyond a certain radius, respectively. The resulting angular velocity gradient resulted in the gold structure rotating about the prescribed axis within the solvent. |
| 445 |
> |
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. |
| 446 |
|
|
| 447 |
|
\begin{longtable}{lccccc} |
| 448 |
< |
\caption{Calculated ``stick'' 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.} |
| 448 |
> |
\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.} |
| 449 |
|
\\ \hline \hline |
| 450 |
< |
{Structure} & {Axis of rotation} & {$\kappa_{VSS}$} & {$\kappa_{calc}$} & {$f_{VSS}$} & {$f_{calc}$}\\ |
| 450 |
> |
{Structure} & {Axis of Rotation} & {$\kappa_{VSS}$} & {$\kappa_{calc}$} & {$f_{VSS}$} & {$f_{calc}$}\\ |
| 451 |
|
{} & {} & {\small($10^{-29}$ Pa s m$^{3}$)} & {\small($10^{-29}$ Pa s m$^{3}$)} & {} & {}\\ \hline |
| 452 |
< |
{Sphere} & {$x = y = z$} & {} & {5.37237} & {1} & {1}\\ |
| 453 |
< |
{Prolate Ellipsoid} & {$x = y$} & {} & {3.59881} & {} & {0.768726}\\ |
| 454 |
< |
{Prolate Ellipsoid} & {$z$} & {} & {9.01084} & {} & {1.92477}\\ \hline \hline |
| 452 |
> |
{Sphere (r = 20 \AA)} & {$x = y = z$} & {} & {6.13239} & {1} & {1}\\ |
| 453 |
> |
{Sphere (r = 30 \AA)} & {$x = y = z$} & {} & {20.6968} & {1} & {1}\\ |
| 454 |
> |
{Sphere (r = 40 \AA)} & {$x = y = z$} & {} & {49.0591} & {1} & {1}\\ |
| 455 |
> |
{Prolate Ellipsoid} & {$x = y$} & {} & {8.22846} & {} & {1.92477}\\ |
| 456 |
> |
{Prolate Ellipsoid} & {$z$} & {} & {3.28634} & {} & {0.768726}\\ |
| 457 |
> |
\hline \hline |
| 458 |
|
\label{table:interfacialfrictionstick} |
| 459 |
|
\end{longtable} |
| 460 |
|
|
