| 128 |
|
interfaces. |
| 129 |
|
|
| 130 |
|
\section{Methodology} |
| 131 |
< |
Similar to the NIVS methodology,\cite{kuang:164101} we consider a |
| 131 |
> |
Similar to the NIVS method,\cite{kuang:164101} we consider a |
| 132 |
|
periodic system divided into a series of slabs along a certain axis |
| 133 |
|
(e.g. $z$). The unphysical thermal and/or momentum flux is designated |
| 134 |
< |
from the center slab to one of the end slabs, and thus the center slab |
| 135 |
< |
would have a lower temperature than the end slab (unless the thermal |
| 136 |
< |
flux is negative). Therefore, the center slab is denoted as ``$c$'' |
| 137 |
< |
while the end slab as ``$h$''. |
| 134 |
> |
from the center slab to one of the end slabs, and thus the thermal |
| 135 |
> |
flux results in a lower temperature of the center slab than the end |
| 136 |
> |
slab, and the momentum flux results in negative center slab momentum |
| 137 |
> |
with positive end slab momentum (unless these fluxes are set |
| 138 |
> |
negative). Therefore, the center slab is denoted as ``$c$'', while the |
| 139 |
> |
end slab as ``$h$''. |
| 140 |
|
|
| 141 |
< |
To impose these fluxes, we periodically apply separate operations to |
| 142 |
< |
velocities of particles {$i$} within the center slab and of particles |
| 143 |
< |
{$j$} within the end slab: |
| 141 |
> |
To impose these fluxes, we periodically apply different set of |
| 142 |
> |
operations on velocities of particles {$i$} within the center slab and |
| 143 |
> |
those of particles {$j$} within the end slab: |
| 144 |
|
\begin{eqnarray} |
| 145 |
|
\vec{v}_i & \leftarrow & c\cdot\left(\vec{v}_i - \langle\vec{v}_c |
| 146 |
|
\rangle\right) + \left(\langle\vec{v}_c\rangle + \vec{a}_c\right) \\ |
| 149 |
|
\end{eqnarray} |
| 150 |
|
where $\langle\vec{v}_c\rangle$ and $\langle\vec{v}_h\rangle$ denotes |
| 151 |
|
the instantaneous bulk velocity of slabs $c$ and $h$ respectively |
| 152 |
< |
before an operation occurs. When a momentum flux $\vec{j}_z(\vec{p})$ |
| 152 |
> |
before an operation is applied. When a momentum flux $\vec{j}_z(\vec{p})$ |
| 153 |
|
presents, these bulk velocities would have a corresponding change |
| 154 |
|
($\vec{a}_c$ and $\vec{a}_h$ respectively) according to Newton's |
| 155 |
|
second law: |
| 157 |
|
M_c \vec{a}_c & = & -\vec{j}_z(\vec{p}) \Delta t \\ |
| 158 |
|
M_h \vec{a}_h & = & \vec{j}_z(\vec{p}) \Delta t |
| 159 |
|
\end{eqnarray} |
| 160 |
< |
where |
| 160 |
> |
where $M$ denotes total mass of particles within a slab: |
| 161 |
|
\begin{eqnarray} |
| 162 |
|
M_c & = & \sum_{i = 1}^{N_c} m_i \\ |
| 163 |
|
M_h & = & \sum_{j = 1}^{N_h} m_j |
| 164 |
|
\end{eqnarray} |
| 165 |
< |
and $\Delta t$ is the interval between two operations. |
| 165 |
> |
and $\Delta t$ is the interval between two separate operations. |
| 166 |
|
|
| 167 |
< |
The above operations conserve the linear momentum of a periodic |
| 168 |
< |
system. To satisfy total energy conservation as well as to impose a |
| 169 |
< |
thermal flux $J_z$, one would have |
| 170 |
< |
[SUPPORT INFO MIGHT BE NECESSARY TO PUT EXTRA MATH IN] |
| 167 |
> |
The above operations already conserve the linear momentum of a |
| 168 |
> |
periodic system. To further satisfy total energy conservation as well |
| 169 |
> |
as to impose the thermal flux $J_z$, the following equations are |
| 170 |
> |
included as well: |
| 171 |
> |
[MAY PUT EXTRA MATH IN SUPPORT INFO OR APPENDIX] |
| 172 |
|
\begin{eqnarray} |
| 173 |
|
K_c - J_z\Delta t & = & c^2 (K_c - \frac{1}{2}M_c \langle\vec{v}_c |
| 174 |
|
\rangle^2) + \frac{1}{2}M_c (\langle \vec{v}_c \rangle + \vec{a}_c)^2 \\ |
| 176 |
|
\rangle^2) + \frac{1}{2}M_h (\langle \vec{v}_h \rangle + \vec{a}_h)^2 |
| 177 |
|
\end{eqnarray} |
| 178 |
|
where $K_c$ and $K_h$ denotes translational kinetic energy of slabs |
| 179 |
< |
$c$ and $h$ respectively before an operation occurs. These |
| 179 |
> |
$c$ and $h$ respectively before an operation is applied. These |
| 180 |
|
translational kinetic energy conservation equations are sufficient to |
| 181 |
< |
ensure total energy conservation, as the operations applied do not |
| 182 |
< |
change the potential energy of a system, given that the potential |
| 181 |
> |
ensure total energy conservation, as the operations applied in our |
| 182 |
> |
method do not change the kinetic energy related to other degrees of |
| 183 |
> |
freedom or the potential energy of a system, given that its potential |
| 184 |
|
energy does not depend on particle velocity. |
| 185 |
|
|
| 186 |
|
The above sets of equations are sufficient to determine the velocity |
| 187 |
|
scaling coefficients ($c$ and $h$) as well as $\vec{a}_c$ and |
| 188 |
< |
$\vec{a}_h$. Note that two roots of $c$ and $h$ exist |
| 189 |
< |
respectively. However, to avoid dramatic perturbations to a system, |
| 190 |
< |
the positive roots (which are closer to 1) are chosen. Figure |
| 191 |
< |
\ref{method} illustrates the implementation of this algorithm in an |
| 192 |
< |
individual step. |
| 188 |
> |
$\vec{a}_h$. Note that there are two roots respectively for $c$ and |
| 189 |
> |
$h$. However, the positive roots (which are closer to 1) are chosen so |
| 190 |
> |
that the perturbations to a system can be reduced to a minimum. Figure |
| 191 |
> |
\ref{method} illustrates the implementation sketch of this algorithm |
| 192 |
> |
in an individual step. |
| 193 |
|
|
| 194 |
|
\begin{figure} |
| 195 |
|
\includegraphics[width=\linewidth]{method} |
| 196 |
|
\caption{Illustration of the implementation of the algorithm in a |
| 197 |
|
single step. Starting from an ideal velocity distribution, the |
| 198 |
< |
transformation is used to apply both thermal and momentum flux from |
| 199 |
< |
the ``c'' slab to the ``h'' slab. As the figure shows, the thermal |
| 200 |
< |
distributions preserve after this operation.} |
| 198 |
> |
transformation is used to apply the effect of both a thermal and a |
| 199 |
> |
momentum flux from the ``c'' slab to the ``h'' slab. As the figure |
| 200 |
> |
shows, thermal distributions can preserve after this operation.} |
| 201 |
|
\label{method} |
| 202 |
|
\end{figure} |
| 203 |
|
|
| 205 |
|
thermal and/or momentum flux can be applied and the corresponding |
| 206 |
|
temperature and/or momentum gradients can be established. |
| 207 |
|
|
| 208 |
< |
This approach is more computationaly efficient compared to the |
| 209 |
< |
previous NIVS method, in that only quadratic equations are involved, |
| 210 |
< |
while the NIVS method needs to solve a quartic equations. Furthermore, |
| 211 |
< |
the method implements isotropic scaling of velocities in respective |
| 212 |
< |
slabs, unlike the NIVS, where an extra criteria function is necessary |
| 213 |
< |
to choose a set of coefficients that performs the most isotropic |
| 214 |
< |
scaling. More importantly, separating the momentum flux imposing from |
| 215 |
< |
velocity scaling avoids the underlying cause that NIVS produced |
| 216 |
< |
thermal anisotropy when applying a momentum flux. |
| 208 |
> |
Compared to the previous NIVS method, this approach is computationally |
| 209 |
> |
more efficient in that only quadratic equations are involved to |
| 210 |
> |
determine a set of scaling coefficients, while the NIVS method needs |
| 211 |
> |
to solve quartic equations. Furthermore, this method implements |
| 212 |
> |
isotropic scaling of velocities in respective slabs, unlike the NIVS, |
| 213 |
> |
where an extra criteria function is necessary to choose a set of |
| 214 |
> |
coefficients that performs a scaling as isotropic as possible. More |
| 215 |
> |
importantly, separating the means of momentum flux imposing from |
| 216 |
> |
velocity scaling avoids the underlying cause to thermal anisotropy in |
| 217 |
> |
NIVS when applying a momentum flux. And later sections will |
| 218 |
> |
demonstrate that this can improve the performance in shear viscosity |
| 219 |
> |
simulations. |
| 220 |
|
|
| 221 |
< |
The advantages of the approach over the original momentum swapping |
| 222 |
< |
approach lies in its nature to preserve a Gaussian |
| 223 |
< |
distribution. Because the momentum swapping tends to render a |
| 224 |
< |
nonthermal distribution, when the imposed flux is relatively large, |
| 225 |
< |
diffusion of the neighboring slabs could no longer remedy this effect, |
| 226 |
< |
and nonthermal distributions would be observed. Results in later |
| 227 |
< |
section will illustrate this effect. |
| 221 |
> |
The advantages of this approach over the original momentum swapping |
| 222 |
> |
one lies in its nature of preserve a Maxwell-Boltzmann distribution |
| 223 |
> |
(mathimatically a Gaussian function). However, because the momentum |
| 224 |
> |
swapping steps tend to result in a nonthermal distribution, when an |
| 225 |
> |
imposed flux is relatively large, diffusions from the neighboring |
| 226 |
> |
slabs could no longer remedy this effect, problematic distributions |
| 227 |
> |
would be observed. Results in later section will illustrate this |
| 228 |
> |
effect in more detail. |
| 229 |
|
|
| 230 |
|
\section{Computational Details} |
| 231 |
|
The algorithm has been implemented in our MD simulation code, |
| 232 |
|
OpenMD\cite{Meineke:2005gd,openmd}. We compare the method with |
| 233 |
< |
previous RNEMD methods or equilibrium MD methods in homogeneous fluids |
| 233 |
> |
previous RNEMD methods or equilibrium MD (EMD) methods in homogeneous fluids |
| 234 |
|
(Lennard-Jones and SPC/E water). And taking advantage of the method, |
| 235 |
|
we simulate the interfacial friction of different heterogeneous |
| 236 |
|
interfaces (gold-organic solvent and gold-SPC/E water and ice-liquid |
| 237 |
|
water). |
| 238 |
|
|
| 239 |
|
\subsection{Simulation Protocols} |
| 240 |
< |
The systems to be investigated are set up in a orthorhombic simulation |
| 241 |
< |
cell with periodic boundary conditions in all three dimensions. The |
| 242 |
< |
$z$ axis of these cells were longer and was set as the gradient axis |
| 243 |
< |
of temperature and/or momentum. Thus the cells were divided into $N$ |
| 244 |
< |
slabs along this axis, with various $N$ depending on individual |
| 245 |
< |
system. The $x$ and $y$ axis were usually of the same length in |
| 246 |
< |
homogeneous systems or close to each other where interfaces |
| 247 |
< |
presents. In all cases, before introducing a nonequilibrium method to |
| 248 |
< |
establish steady thermal and/or momentum gradients for further |
| 249 |
< |
measurements and calculations, canonical ensemble with a Nos\'e-Hoover |
| 250 |
< |
thermostat\cite{hoover85} and microcanonical ensemble equilibrations |
| 251 |
< |
were used to prepare systems ready for data |
| 252 |
< |
collections. Isobaric-isothermal equilibrations are performed before |
| 253 |
< |
this for SPC/E water systems to reach normal pressure (1 bar), while |
| 254 |
< |
similar equilibrations are used for interfacial systems to relax the |
| 255 |
< |
surface tensions. |
| 240 |
> |
The systems to be investigated are set up in orthorhombic simulation |
| 241 |
> |
cells with periodic boundary conditions in all three dimensions. The |
| 242 |
> |
$z$ axis of these cells were longer and set as the temperature and/or |
| 243 |
> |
momentum gradient axis. And the cells were evenly divided into $N$ |
| 244 |
> |
slabs along this axis, with various $N$ depending on individual |
| 245 |
> |
system. The $x$ and $y$ axis were of the same length in homogeneous |
| 246 |
> |
systems or had length scale close to each other where heterogeneous |
| 247 |
> |
interfaces presents. In all cases, before introducing a nonequilibrium |
| 248 |
> |
method to establish steady thermal and/or momentum gradients for |
| 249 |
> |
further measurements and calculations, canonical ensemble with a |
| 250 |
> |
Nos\'e-Hoover thermostat\cite{hoover85} and microcanonical ensemble |
| 251 |
> |
equilibrations were used before data collections. For SPC/E water |
| 252 |
> |
simulations, isobaric-isothermal equilibrations\cite{melchionna93} are |
| 253 |
> |
performed before the above to reach normal pressure (1 bar); for |
| 254 |
> |
interfacial systems, similar equilibrations are used to relax the |
| 255 |
> |
surface tensions of the $xy$ plane. |
| 256 |
|
|
| 257 |
< |
While homogeneous fluid systems can be set up with random |
| 258 |
< |
configurations, our interfacial systems needs extra steps to ensure |
| 259 |
< |
the interfaces be established properly for computations. The |
| 257 |
> |
While homogeneous fluid systems can be set up with rather random |
| 258 |
> |
configurations, our interfacial systems needs a series of steps to |
| 259 |
> |
ensure the interfaces be established properly for computations. The |
| 260 |
|
preparation and equilibration of butanethiol covered gold (111) |
| 261 |
|
surface and further solvation and equilibration process is described |
| 262 |
< |
as in reference \cite{kuang:AuThl}. |
| 262 |
> |
in details as in reference \cite{kuang:AuThl}. |
| 263 |
|
|
| 264 |
|
As for the ice/liquid water interfaces, the basal surface of ice |
| 265 |
|
lattice was first constructed. Hirsch {\it et |
| 266 |
|
al.}\cite{doi:10.1021/jp048434u} explored the energetics of ice |
| 267 |
< |
lattices with different proton orders. We refer to their results and |
| 268 |
< |
choose the configuration of the lowest energy after geometry |
| 269 |
< |
optimization as the unit cells of our ice lattices. Although |
| 270 |
< |
experimental solid/liquid coexistant temperature near normal pressure |
| 271 |
< |
is 273K, Bryk and Haymet's simulations of ice/liquid water interfaces |
| 272 |
< |
with different models suggest that for SPC/E, the most stable |
| 273 |
< |
interface is observed at 225$\pm$5K. Therefore, all our ice/liquid |
| 274 |
< |
water simulations were carried out under 225K. To have extra |
| 275 |
< |
protection of the ice lattice during initial equilibration (when the |
| 276 |
< |
randomly generated liquid phase configuration could release large |
| 277 |
< |
amount of energy in relaxation), a constraint method (REF?) was |
| 278 |
< |
adopted until the high energy configuration was relaxed. |
| 279 |
< |
[MAY ADD A FIGURE HERE FOR BASAL PLANE, MAY INCLUDE PRISM IF POSSIBLE] |
| 267 |
> |
lattices with all possible proton order configurations. We refer to |
| 268 |
> |
their results and choose the configuration of the lowest energy after |
| 269 |
> |
geometry optimization as the unit cell for our ice lattices. Although |
| 270 |
> |
experimental solid/liquid coexistant temperature under normal pressure |
| 271 |
> |
should be close to 273K, Bryk and Haymet's simulations of ice/liquid |
| 272 |
> |
water interfaces with different models suggest that for SPC/E, the |
| 273 |
> |
most stable interface is observed at 225$\pm$5K.\cite{bryk:10258} |
| 274 |
> |
Therefore, our ice/liquid water simulations were carried out at |
| 275 |
> |
225K. To have extra protection of the ice lattice during initial |
| 276 |
> |
equilibration (when the randomly generated liquid phase configuration |
| 277 |
> |
could release large amount of energy in relaxation), restraints were |
| 278 |
> |
applied to the ice lattice to avoid inadvertent melting by the heat |
| 279 |
> |
dissipated from the high enery configurations. |
| 280 |
> |
[MAY ADD A SNAPSHOT FOR BASAL PLANE] |
| 281 |
|
|
| 282 |
|
\subsection{Force Field Parameters} |
| 283 |
|
For comparison of our new method with previous work, we retain our |
| 284 |
< |
force field parameters consistent with the results we will compare |
| 285 |
< |
with. The Lennard-Jones fluid used here for argon , and reduced unit |
| 286 |
< |
results are reported for direct comparison purpose. |
| 284 |
> |
force field parameters consistent with previous simulations. Argon is |
| 285 |
> |
the Lennard-Jones fluid used here, and its results are reported in |
| 286 |
> |
reduced unit for direct comparison purpose. |
| 287 |
|
|
| 288 |
|
As for our water simulations, SPC/E model is used throughout this work |
| 289 |
|
for consistency. Previous work for transport properties of SPC/E water |
| 298 |
|
Spohr potential was adopted\cite{ISI:000167766600035} to depict |
| 299 |
|
Au-H$_2$O interactions. |
| 300 |
|
|
| 301 |
< |
The small organic molecules included in our simulations are the Au |
| 302 |
< |
surface capping agent butanethiol and liquid hexane and toluene. The |
| 303 |
< |
United-Atom |
| 301 |
> |
For our gold/organic liquid interfaces, the small organic molecules |
| 302 |
> |
included in our simulations are the Au surface capping agent |
| 303 |
> |
butanethiol and liquid hexane and toluene. The United-Atom |
| 304 |
|
models\cite{TraPPE-UA.thiols,TraPPE-UA.alkanes,TraPPE-UA.alkylbenzenes} |
| 305 |
|
for these components were used in this work for better computational |
| 306 |
|
efficiency, while maintaining good accuracy. We refer readers to our |
| 328 |
|
\end{equation} |
| 329 |
|
where $L_x$ and $L_y$ denotes the dimensions of the plane in a |
| 330 |
|
simulation cell perpendicular to the thermal gradient, and a factor of |
| 331 |
< |
two in the denominator is present for the heat transport occurs in |
| 332 |
< |
both $+z$ and $-z$ directions. The temperature gradient |
| 331 |
> |
two in the denominator is necessary for the heat transport occurs in |
| 332 |
> |
both $+z$ and $-z$ directions. The average temperature gradient |
| 333 |
|
${\langle\partial T/\partial z\rangle}$ can be obtained by a linear |
| 334 |
|
regression of the temperature profile, which is recorded during a |
| 335 |
|
simulation for each slab in a cell. For Lennard-Jones simulations, |
| 338 |
|
|
| 339 |
|
\subsection{Shear viscosities} |
| 340 |
|
Alternatively, the method can carry out shear viscosity calculations |
| 341 |
< |
by switching off $J_z$. One can specify the vector |
| 342 |
< |
$\vec{j}_z(\vec{p})$ by choosing the three components |
| 341 |
> |
by specify a momentum flux. In our algorithm, one can specify the |
| 342 |
> |
three components of the flux vector $\vec{j}_z(\vec{p})$ |
| 343 |
|
respectively. For shear viscosity simulations, $j_z(p_z)$ is usually |
| 344 |
< |
set to zero. Although for isotropic systems, the direction of |
| 345 |
< |
$\vec{j}_z(\vec{p})$ on the $xy$ plane should not matter, the ability |
| 346 |
< |
of arbitarily specifying the vector direction in our method provides |
| 347 |
< |
convenience in anisotropic simulations. |
| 344 |
> |
set to zero. For isotropic systems, the direction of |
| 345 |
> |
$\vec{j}_z(\vec{p})$ on the $xy$ plane should not matter, but the |
| 346 |
> |
ability of arbitarily specifying the vector direction in our method |
| 347 |
> |
could provide convenience in anisotropic simulations. |
| 348 |
|
|
| 349 |
< |
Similar to thermal conductivity computations, linear response of the |
| 350 |
< |
momentum gradient with respect to the shear stress is assumed, and the |
| 351 |
< |
shear viscosity ($\eta$) can be obtained with the imposed momentum |
| 352 |
< |
flux (e.g. in $x$ direction) and the measured gradient: |
| 349 |
> |
Similar to thermal conductivity computations, for a homogeneous |
| 350 |
> |
system, linear response of the momentum gradient with respect to the |
| 351 |
> |
shear stress is assumed, and the shear viscosity ($\eta$) can be |
| 352 |
> |
obtained with the imposed momentum flux (e.g. in $x$ direction) and |
| 353 |
> |
the measured gradient: |
| 354 |
|
\begin{equation} |
| 355 |
|
j_z(p_x) = -\eta \frac{\partial v_x}{\partial z} |
| 356 |
|
\end{equation} |
| 359 |
|
j_z(p_x) = \frac{P_x}{2 t L_x L_y} |
| 360 |
|
\end{equation} |
| 361 |
|
with $P_x$ being the total non-physical momentum transferred within |
| 362 |
< |
the data collection time. Also, the velocity gradient |
| 362 |
> |
the data collection time. Also, the averaged velocity gradient |
| 363 |
|
${\langle\partial v_x/\partial z\rangle}$ can be obtained using linear |
| 364 |
< |
regression of the $x$ component of the mean velocity, $\langle |
| 365 |
< |
v_x\rangle$, in each of the bins. For Lennard-Jones simulations, shear |
| 366 |
< |
viscosities are reported in reduced units |
| 364 |
> |
regression of the $x$ component of the mean velocity ($\langle |
| 365 |
> |
v_x\rangle$) in each of the bins. For Lennard-Jones simulations, shear |
| 366 |
> |
viscosities are also reported in reduced units |
| 367 |
|
(${\eta^*=\eta\sigma^2(\varepsilon m)^{-1/2}}$). |
| 368 |
+ |
|
| 369 |
+ |
[COMBINE JZKE W JZPX] |
| 370 |
|
|
| 371 |
|
\subsection{Interfacial friction and Slip length} |
| 372 |
|
While the shear stress results in a velocity gradient within bulk |
| 440 |
|
\multicolumn{2}{c}{$\eta^*$} \\ |
| 441 |
|
\hline |
| 442 |
|
Swap Interval $t^*$ & Equivalent $J_z^*$ or $j_z^*(p_x)$ & |
| 443 |
< |
NIVS & This work & Swapping & This work \\ |
| 443 |
> |
NIVS\footnote{Reference \cite{kuang:164101}.} & This work & |
| 444 |
> |
Swapping & This work \\ |
| 445 |
|
\hline |
| 446 |
|
0.116 & 0.16 & 7.30(0.10) & 6.25(0.23) & 3.57(0.06) & 3.53(0.16)\\ |
| 447 |
|
0.232 & 0.09 & 6.95(0.09) & 8.02(0.56) & 3.64(0.05) & 3.43(0.17)\\ |
