1 |
\documentclass[11pt]{article} |
2 |
\usepackage{amsmath} |
3 |
\usepackage{amssymb} |
4 |
\usepackage{setspace} |
5 |
\usepackage{endfloat} |
6 |
\usepackage{caption} |
7 |
%\usepackage{tabularx} |
8 |
\usepackage{graphicx} |
9 |
\usepackage{multirow} |
10 |
%\usepackage{booktabs} |
11 |
%\usepackage{bibentry} |
12 |
%\usepackage{mathrsfs} |
13 |
%\usepackage[ref]{overcite} |
14 |
\usepackage[square, comma, sort&compress]{natbib} |
15 |
\usepackage{url} |
16 |
\pagestyle{plain} \pagenumbering{arabic} \oddsidemargin 0.0cm |
17 |
\evensidemargin 0.0cm \topmargin -21pt \headsep 10pt \textheight |
18 |
9.0in \textwidth 6.5in \brokenpenalty=10000 |
19 |
|
20 |
% double space list of tables and figures |
21 |
\AtBeginDelayedFloats{\renewcommand{\baselinestretch}{1.66}} |
22 |
\setlength{\abovecaptionskip}{20 pt} |
23 |
\setlength{\belowcaptionskip}{30 pt} |
24 |
|
25 |
%\renewcommand\citemid{\ } % no comma in optional reference note |
26 |
\bibpunct{[}{]}{,}{n}{}{;} |
27 |
\bibliographystyle{achemso} |
28 |
|
29 |
\begin{document} |
30 |
|
31 |
\title{Velocity Shearing and Scaling RNEMD: a minimally perturbing |
32 |
method for simulating temperature and momentum gradients} |
33 |
|
34 |
\author{Shenyu Kuang and J. Daniel |
35 |
Gezelter\footnote{Corresponding author. \ Electronic mail: gezelter@nd.edu} \\ |
36 |
Department of Chemistry and Biochemistry,\\ |
37 |
University of Notre Dame\\ |
38 |
Notre Dame, Indiana 46556} |
39 |
|
40 |
\date{\today} |
41 |
|
42 |
\maketitle |
43 |
|
44 |
\begin{doublespace} |
45 |
|
46 |
\begin{abstract} |
47 |
We present a new method for introducing stable nonequilibrium |
48 |
velocity and temperature gradients in molecular dynamics simulations |
49 |
of heterogeneous and interfacial systems. This method conserves both |
50 |
the linear momentum and total energy of the system and improves |
51 |
previous reverse non-equilibrium molecular dynamics (RNEMD) methods |
52 |
while retaining equilibrium thermal velocity distributions in each |
53 |
region of the system. The new method avoids the thermal anisotropy |
54 |
produced by previous methods by using isotropic velocity scaling and |
55 |
shearing (VSS) on all of the molecules in specific regions. To test |
56 |
the method, we have computed the thermal conductivity and shear |
57 |
viscosity of model liquid systems as well as the interfacial |
58 |
friction coeefficients of a series of solid / liquid interfaces. The |
59 |
method's ability to combine the thermal and momentum gradients |
60 |
allows us to obtain shear viscosity data for a range of temperatures |
61 |
from a single trajectory. |
62 |
\end{abstract} |
63 |
|
64 |
\newpage |
65 |
|
66 |
%\narrowtext |
67 |
|
68 |
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% |
69 |
% BODY OF TEXT |
70 |
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% |
71 |
|
72 |
\section{Introduction} |
73 |
One of the standard ways to compute transport coefficients such as the |
74 |
viscosities and thermal conductivities of liquids is to use |
75 |
imposed-flux non-equilibrium molecular dynamics |
76 |
methods.\cite{MullerPlathe:1997xw,ISI:000080382700030,kuang:164101} |
77 |
These methods establish stable non-equilibrium conditions in a |
78 |
simulation box using an applied momentum or thermal flux between |
79 |
different regions of the box. The corresponding temperature or |
80 |
velocity gradients which develop in response to the applied flux is |
81 |
related (via linear response theory) to the transport coefficient of |
82 |
interest. These methods are quite efficient, in that they do not need |
83 |
many trajectories to provide information about transport properties. To |
84 |
date, they have been utilized in computing thermal and mechanical |
85 |
transfer of both homogeneous liquids as well as heterogeneous systems |
86 |
such as solid-liquid |
87 |
interfaces.\cite{garde:nl2005,garde:PhysRevLett2009,kuang:AuThl} |
88 |
|
89 |
The reverse non-equilibrium molecular dynamics (RNEMD) methods utilize |
90 |
additional constraints that ensure conservation of linear momentum and |
91 |
total energy of the system while imposing the desired flux. The RNEMD |
92 |
methods are therefore capable of sampling various thermodynamically |
93 |
relevent ensembles, including the micro-canonical (NVE) ensemble, |
94 |
without resorting to an external thermostat. The original approaches |
95 |
proposed by M\"{u}ller-Plathe {\it et |
96 |
al.}\cite{MullerPlathe:1997xw,ISI:000080382700030} utilize simple |
97 |
momentum swapping moves for generating energy/momentum fluxes. The |
98 |
swapping moves can also be made compatible with particles of different |
99 |
identities. Although the swapping moves are simple to implement in |
100 |
molecular simulations, Tenney and Maginn have shown that they create |
101 |
nonthermal velocity distributions.\cite{Maginn:2010} Furthermore, this |
102 |
approach is not particularly efficient for kinetic energy transfer |
103 |
between particles of different identities, particularly when the mass |
104 |
difference between the particles becomes significant. This problem |
105 |
makes applying swapping-move RNEMD methods on heterogeneous |
106 |
interfacial systems somewhat difficult. |
107 |
|
108 |
Recently, we developed a somewhat different approach to applying |
109 |
thermal fluxes in RNEMD simulation using a non-isotropic velocity |
110 |
scaling (NIVS) algorithm.\cite{kuang:164101} This algorithm scales all |
111 |
atomic velocity vectors in two separate regions of a simulated system |
112 |
using two diagonal scaling matrices. The scaling matrices are |
113 |
determined by solving single quartic equation which includes linear |
114 |
momentum and kinetic energy conservation constraints as well as the |
115 |
target thermal flux between the regions. The NIVS method is able to |
116 |
effectively impose a wide range of kinetic energy fluxes without |
117 |
significant perturbation to the velocity distributions away from ideal |
118 |
Maxwell-Boltzmann distributions, even in the presence of heterogeneous |
119 |
interfaces. We successfully applied this approach in studying the |
120 |
interfacial thermal conductance at chemically-capped metal-solvent |
121 |
interfaces.\cite{kuang:AuThl} |
122 |
|
123 |
The NIVS approach works very well for preparing stable {\it thermal} |
124 |
gradients. However, as we pointed out in the original paper, it has |
125 |
limited application in imposing {\it linear} momentum fluxes (which |
126 |
are required for measuring shear viscosities). The reason for this is |
127 |
that linear momentum flux was being imposed by scaling random |
128 |
fluctuations of the center of the velocity distributions. Repeated |
129 |
application of the original NIVS approach resulted in temperature |
130 |
anisotropy, i.e. the width of the velocity distributions depended on |
131 |
coordinates perpendicular to the desired gradient direction. For this |
132 |
reason, combining thermal and momentum fluxes was difficult with the |
133 |
original NIVS algorithm. However, combinations of thermal and |
134 |
velocity gradients would provide a means to simulate thermal-linear |
135 |
coupled processes such as Soret effect in liquid flows. Therefore, |
136 |
developing improved approaches to the scaling imposed-flux methods |
137 |
would be useful. |
138 |
|
139 |
In this paper, we improve the RNEMD methods by introducing a novel |
140 |
approach to creating imposed fluxes. This approach separates the |
141 |
shearing and scaling of the velocity distributions in different |
142 |
spatial regions and can apply both transformations within a single |
143 |
time step. The approach retains desirable features of previous RNEMD |
144 |
approaches and is simpler to implement compared to the earlier NIVS |
145 |
method. In what follows, we first present the velocity shearing and |
146 |
scaling (VSS) RNEMD method and its implementation in a |
147 |
simulation. Then we compare the VSS-RNEMD method in bulk fluids to |
148 |
previous methods. We also compute interfacial frictions are computed |
149 |
for a series of heterogeneous interfaces. |
150 |
|
151 |
\section{Methodology} |
152 |
In an approach similar to the earlier NIVS method,\cite{kuang:164101} |
153 |
we consider a periodic system which has been divided into a series of |
154 |
slabs along a single axis (e.g. $z$). The unphysical thermal and |
155 |
momentum fluxes are applied from one of the end slabs to the center |
156 |
slab, and thus the thermal flux produces a higher temperature in the |
157 |
center slab than in the end slab, and the momentum flux results in a |
158 |
center slab moving along the positive $y$ axis (see Fig. \ref{cartoon}). |
159 |
The applied fluxes can be set to negative values if the reverse |
160 |
gradients are desired. For convenience the center slab is denoted as |
161 |
the {\it hot} or {\it ``H''} slab, while the end slab is denoted {\it |
162 |
``C''} (or {\it cold}). |
163 |
|
164 |
\begin{figure} |
165 |
\includegraphics[width=\linewidth]{cartoon} |
166 |
\caption{The VSS-RNEMD approach imposes unphysical transfer of both |
167 |
linear momentum and kinetic energy between a ``hot'' slab and a |
168 |
``cold'' slab in the simulation box. The system responds to this |
169 |
imposed flux by generating both momentum and temperature gradients. |
170 |
The slope of the gradients can then be used to compute transport |
171 |
properties (e.g. shear viscosity and thermal conductivity).} |
172 |
\label{cartoon} |
173 |
\end{figure} |
174 |
|
175 |
To impose these fluxes, we periodically apply a set of operations on |
176 |
the velocities of particles {$i$} within the cold slab and a separate |
177 |
operation on particles {$j$} within the hot slab. |
178 |
\begin{eqnarray} |
179 |
\vec{v}_i & \leftarrow & c\cdot\left(\vec{v}_i - \langle\vec{v}_c |
180 |
\rangle\right) + \left(\langle\vec{v}_c\rangle + |
181 |
\vec{a}_c\right) \label{eq:xformc} \\ |
182 |
\vec{v}_j & \leftarrow & h\cdot\left(\vec{v}_j - \langle\vec{v}_h |
183 |
\rangle\right) + \left(\langle\vec{v}_h\rangle + \vec{a}_h\right) \label{eq:xformh} |
184 |
\end{eqnarray} |
185 |
where $\langle\vec{v}_c\rangle$ and $\langle\vec{v}_h\rangle$ denote |
186 |
the instantaneous average velocity of the molecules within slabs $C$ |
187 |
and $H$ respectively. When a momentum flux $\vec{j}_z(\vec{p})$ is |
188 |
present, these slab-averaged velocities also get corresponding |
189 |
incremental changes ($\vec{a}_c$ and $\vec{a}_h$ respectively) that |
190 |
are applied to all particles within each slab. The incremental |
191 |
changes are obtained using Newton's second law: |
192 |
\begin{eqnarray} |
193 |
M_c \vec{a}_c & = & -\vec{j}_z(\vec{p}) \Delta t \label{eq:newton1} \\ |
194 |
M_h \vec{a}_h & = & \vec{j}_z(\vec{p}) \Delta t |
195 |
\label{eq:newton2} |
196 |
\end{eqnarray} |
197 |
where $M$ denotes total mass of particles within a slab: |
198 |
\begin{eqnarray} |
199 |
M_c & = & \sum_{i = 1}^{N_c} m_i \\ |
200 |
M_h & = & \sum_{j = 1}^{N_h} m_j |
201 |
\end{eqnarray} |
202 |
and $\Delta t$ is the interval between two separate operations. |
203 |
|
204 |
The operations in Eqs. \ref{eq:newton1} and \ref{eq:newton2} already |
205 |
conserve the linear momentum of a periodic system. To further satisfy |
206 |
total energy conservation as well as to impose the thermal flux $J_z$, |
207 |
the following constraint equations must be solved for the two scaling |
208 |
variables $c$ and $h$: |
209 |
\begin{eqnarray} |
210 |
K_c - J_z\Delta t & = & c^2 (K_c - \frac{1}{2}M_c \langle\vec{v}_c |
211 |
\rangle^2) + \frac{1}{2}M_c (\langle \vec{v}_c \rangle + \vec{a}_c)^2 \\ |
212 |
K_h + J_z\Delta t & = & h^2 (K_h - \frac{1}{2}M_h \langle\vec{v}_h |
213 |
\rangle^2) + \frac{1}{2}M_h (\langle \vec{v}_h \rangle + \vec{a}_h)^2. |
214 |
\label{constraint} |
215 |
\end{eqnarray} |
216 |
Here $K_c$ and $K_h$ denote the translational kinetic energy of slabs |
217 |
$C$ and $H$ respectively. These conservation equations are sufficient |
218 |
to ensure total energy conservation, as the operations applied in our |
219 |
method do not change the kinetic energy related to orientational |
220 |
degrees of freedom or the potential energy of a system (as long as the |
221 |
potential energy is independent of particle velocity). |
222 |
|
223 |
Equations \ref{eq:newton1}-\ref{constraint} are sufficient to |
224 |
determine the velocity scaling coefficients ($c$ and $h$) as well as |
225 |
$\vec{a}_c$ and $\vec{a}_h$. Note that there are usually two roots |
226 |
respectively for $c$ and $h$. However, the positive roots (which are |
227 |
closer to 1) are chosen so that the perturbations to the system are |
228 |
minimal. Figure \ref{method} illustrates the implementation of this |
229 |
algorithm in an individual step. |
230 |
|
231 |
\begin{figure} |
232 |
\centering |
233 |
\includegraphics[width=5in]{method} |
234 |
\caption{A single step implementation of the VSS algorithm starts with |
235 |
the velocity distributions for the two slabs in a shearing fluid |
236 |
(solid lines). Equations \ref{eq:xformc} and \ref{eq:xformh} are |
237 |
used to scale and shear velocities in both slabs to mimic the effect |
238 |
of both a thermal and a momentum flux. Gaussian distributions are |
239 |
preserved by both the scaling and shearing operations.} |
240 |
\label{method} |
241 |
\end{figure} |
242 |
|
243 |
By implementing these operations at a fixed frequency, stable thermal |
244 |
and momentum fluxes can both be applied and the corresponding |
245 |
temperature and momentum gradients can be established. |
246 |
|
247 |
Compared to the previous NIVS method, the VSS-RNEMD approach is |
248 |
computationally simpler in that only quadratic equations are involved |
249 |
to determine a set of scaling coefficients, while the NIVS method |
250 |
required solution of quartic equations. Furthermore, this method |
251 |
implements {\it isotropic} scaling of velocities in respective slabs, |
252 |
unlike NIVS, which required extra flexibility in the choice of scaling |
253 |
coefficients to allow for the energy and momentum constraints. Most |
254 |
importantly, separating the linear momentum flux from velocity scaling |
255 |
avoids the underlying cause of the thermal anisotropy in NIVS. In |
256 |
later sections, we demonstrate that this can improve the calculated |
257 |
shear viscosities from RNEMD simulations. |
258 |
|
259 |
The VSS-RNEMD approach has advantages over the original momentum |
260 |
swapping in many respects. In the swapping method, the velocity |
261 |
vectors involved are usually very different (or the generated flux |
262 |
would be quite small), thus the swapping tends to create strong |
263 |
perturbations in the neighborhood of the particles involved. In |
264 |
comparison, the VSS approach distributes the flux widely to all |
265 |
particles in a slab so that perturbations in the flux generating |
266 |
region are minimized. Additionally, momentum swapping steps tend to |
267 |
result in nonthermal velocity distributions when the imposed flux is |
268 |
large and diffusion from the neighboring slabs cannot carry momentum |
269 |
away quickly enough. In comparison, the scaling and shearing moves |
270 |
made in the VSS approach preserve the shapes of the equilibrium |
271 |
velocity disributions (e.g. Maxwell-Boltzmann). The results presented |
272 |
in later sections will illustrate that this is quite helpful in |
273 |
retaining reasonable thermal distributions in a simulation. |
274 |
|
275 |
\section{Computational Details} |
276 |
The algorithm has been implemented in our MD simulation code, |
277 |
OpenMD.\cite{Meineke:2005gd,openmd} We will first compare the method |
278 |
with previous RNEMD methods and equilibrium MD in homogeneous |
279 |
fluids (Lennard-Jones and SPC/E water). We have also used the new |
280 |
method to simulate the interfacial friction of different heterogeneous |
281 |
interfaces (Au (111) with organic solvents, Au(111) with SPC/E water, |
282 |
and the ice Ih - liquid water interface). |
283 |
|
284 |
\subsection{Simulation Protocols} |
285 |
The systems we investigated were set up in orthorhombic simulation |
286 |
cells with periodic boundary conditions in all three dimensions. The |
287 |
$z$-axes of these cells were typically quite long and served as the |
288 |
temperature and/or momentum gradient axes. The cells were evenly |
289 |
divided into $N$ slabs along this axis, with $N$ varying for the |
290 |
individual system. The $x$ and $y$ axes were of similar lengths in all |
291 |
simulations. In all cases, before introducing a nonequilibrium method |
292 |
to establish steady thermal and/or momentum gradients, equilibration |
293 |
simulations were run under the canonical ensemble with a Nos\'e-Hoover |
294 |
thermostat\cite{hoover85} followed by further equilibration using |
295 |
standard constant energy (NVE) conditions. For SPC/E water |
296 |
simulations, isobaric-isothermal equilibrations\cite{melchionna93} |
297 |
were performed before equilibration to reach standard densities at |
298 |
atmospheric pressure (1 bar); for interfacial systems, similar |
299 |
equilibrations with anisotropic box relaxations are used to relax the |
300 |
surface tension in the $xy$ plane. |
301 |
|
302 |
While homogeneous fluid systems can be set up with random |
303 |
configurations, interfacial systems are typically prepared with a |
304 |
single crystal face presented perpendicular to the $z$-axis. This |
305 |
crystal face is aligned in the x-y plane of the periodic cell, and the |
306 |
solvent occupies the region directly above and below a crystalline |
307 |
slab. The preparation and equilibration of butanethiol covered |
308 |
Au(111) surfaces, as well as the solvation and equilibration processes |
309 |
used for these interfaces are described in detail in reference |
310 |
\cite{kuang:AuThl}. |
311 |
|
312 |
For the ice / liquid water interfaces, the basal surface of ice |
313 |
lattice was first constructed. Hirsch and Ojam\"{a}e |
314 |
\cite{doi:10.1021/jp048434u} explored the energetics of ice lattices |
315 |
with all possible proton ordered configurations. We utilized Hirsch |
316 |
and Ojam\"{a}e's structure 6 ($P2_12_12_1$) which is an orthorhombic |
317 |
cell that produces a proton-ordered version of ice Ih. The basal face |
318 |
of ice in this unit cell geometry is the $\{0~0~1\}$ face. Although |
319 |
experimental solid/liquid coexistant temperature under normal pressure |
320 |
are close to 273K, Bryk and Haymet's simulations of ice/liquid water |
321 |
interfaces with different models suggest that for SPC/E, the most |
322 |
stable interface is observed at 225$\pm$5K.\cite{bryk:10258} |
323 |
Therefore, our ice/liquid water simulations were carried out at an |
324 |
average temperature of 225K. Molecular translation and orientational |
325 |
restraints were applied in the early stages of equilibration to |
326 |
prevent melting of the ice slab. These restraints were removed during |
327 |
NVT equilibration, well before data collection was carried out. |
328 |
|
329 |
\subsection{Force Field Parameters} |
330 |
For comparing the VSS-RNEMD method with previous work, we retained |
331 |
force field parameters consistent with previous simulations. Argon is |
332 |
the Lennard-Jones fluid used here, and its results are reported in |
333 |
reduced units for purposes of direct comparison with previous |
334 |
calculations. |
335 |
|
336 |
For our water simulations, we utilized the SPC/E model throughout this |
337 |
work. Previous work for transport properties of SPC/E water model is |
338 |
available\cite{Bedrov:2000,10.1063/1.3330544,Medina2011} so direct |
339 |
comparison with previous calculation methods is possible. |
340 |
|
341 |
The Au-Au interaction parameters in all simulations are described by |
342 |
the quantum Sutton-Chen (QSC) formulation.\cite{PhysRevB.59.3527} The |
343 |
QSC potentials include zero-point quantum corrections and are |
344 |
reparametrized for accurate surface energies compared to the |
345 |
Sutton-Chen potentials.\cite{Chen90} For gold/water interfaces, the |
346 |
Spohr potential was adopted\cite{ISI:000167766600035} to depict |
347 |
Au-H$_2$O interactions. |
348 |
|
349 |
For our gold/organic liquid interfaces, the small organic molecules |
350 |
included in our simulations are the Au surface capping agent |
351 |
butanethiol as well as liquid hexane and liquid toluene. The |
352 |
United-Atom |
353 |
models\cite{TraPPE-UA.thiols,TraPPE-UA.alkanes,TraPPE-UA.alkylbenzenes} |
354 |
for these components were used in this work for computational |
355 |
efficiency, while maintaining good accuracy. We refer readers to our |
356 |
previous work\cite{kuang:AuThl} for further details of these models, |
357 |
as well as the interactions between Au and the above organic molecule |
358 |
components. |
359 |
|
360 |
\subsection{Thermal conductivities} |
361 |
When the linear momentum flux $\vec{j}_z(\vec{p})$ is set to zero and |
362 |
the target $J_z$ is non-zero, VSS-RNEMD imposes kinetic energy |
363 |
transfer between the slabs, which can be used for computation of |
364 |
thermal conductivities. As with other RNEMD methods, we assume that we |
365 |
are in the linear response regime of the temperature gradient with |
366 |
respect to the thermal flux. The thermal conductivity ($\lambda$) can |
367 |
then be calculated using the imposed kinetic energy flux and the |
368 |
measured thermal gradient: |
369 |
\begin{equation} |
370 |
J_z = -\lambda \frac{\partial T}{\partial z} |
371 |
\end{equation} |
372 |
Like other imposed-flux methods, the energy flux was calculated using |
373 |
the total non-physical energy transferred (${E_{total}}$) from slab |
374 |
{\it C} to slab {\it H}, which was recorded throughout the simulation, and |
375 |
the total data collection time $t$: |
376 |
\begin{equation} |
377 |
J_z = \frac{E_{total}}{2 t L_x L_y}. |
378 |
\end{equation} |
379 |
Here, $L_x$ and $L_y$ denote the dimensions of the plane in the |
380 |
simulation cell that is perpendicular to the thermal gradient, and a |
381 |
factor of two in the denominator is necessary as the heat transport |
382 |
occurs in both the $+z$ and $-z$ directions. The average temperature |
383 |
gradient ${\langle\partial T/\partial z\rangle}$ can be obtained by a |
384 |
linear regression of the temperature profile, which is recorded during |
385 |
a simulation for each slab in a cell. For Lennard-Jones simulations, |
386 |
thermal conductivities are reported in reduced units |
387 |
(${\lambda^*=\lambda \sigma^2 m^{1/2} k_B^{-1}\varepsilon^{-1/2}}$). |
388 |
|
389 |
\subsection{Shear viscosities} |
390 |
Alternatively, when the linear momentum flux $\vec{j}_z(\vec{p})$ is |
391 |
non-zero in either the $x$ or $y$ directions, the method can be used |
392 |
to compute the shear viscosity. For isotropic systems, the direction |
393 |
of $\vec{j}_z(\vec{p})$ on the $xy$ plane should not matter, but the |
394 |
ability to arbitarily specify the vector direction in our method could |
395 |
be convenient when working with anisotropic interfaces. |
396 |
|
397 |
In a manner similar to the thermal conductivity calculations, a linear |
398 |
response of the momentum gradient with respect to the shear stress is |
399 |
assumed, and the shear viscosity ($\eta$) can be obtained with the |
400 |
imposed momentum flux (e.g. in $x$ direction) and the measured |
401 |
velocity gradient: |
402 |
\begin{equation} |
403 |
j_z(p_x) = -\eta \frac{\partial v_x}{\partial z} |
404 |
\end{equation} |
405 |
where the flux is similarly defined: |
406 |
\begin{equation} |
407 |
j_z(p_x) = \frac{P_x}{2 t L_x L_y} |
408 |
\end{equation} |
409 |
with $P_x$ being the total non-physical momentum transferred within |
410 |
the data collection time. Also, the averaged velocity gradient |
411 |
${\langle\partial v_x/\partial z\rangle}$ can be obtained using linear |
412 |
regression of the $x$ component of the mean velocity ($\langle |
413 |
v_x\rangle$) in each of the bins. For Lennard-Jones simulations, shear |
414 |
viscosities are also reported in reduced units |
415 |
(${\eta^*=\eta\sigma^2(\varepsilon m)^{-1/2}}$). |
416 |
|
417 |
Although $J_z$ may be switched off for shear viscosity simulations, |
418 |
the VSS-RNEMD method allows the user the ability to {\it |
419 |
simultaneously} impose both a thermal and a momentum flux during a |
420 |
single simulation. This creates configurations with coincident |
421 |
temperature and a velocity gradients. Since the viscosity is |
422 |
generally a function of temperature, the local viscosity depends on |
423 |
the local temperature in the fluid. Therefore, in a single simulation, |
424 |
viscosity at $z$ (corresponding to a certain $T$) can be computed with |
425 |
the applied shear flux and the local velocity gradient (which can be |
426 |
obtained by finite difference approximation). This means that the |
427 |
temperature dependence of the viscosity can be mapped out in only one |
428 |
trajectory. Results for shear viscosity computations of SPC/E water |
429 |
will demonstrate VSS-RNEMD's efficiency in this respect. |
430 |
|
431 |
\subsection{Interfacial friction and slip length} |
432 |
Shear stress creates a velocity gradient within bulk fluid phases, but |
433 |
at solid-liquid interfaces, the effects of a shear stress depend on |
434 |
the molecular details of the interface. The interfacial friction |
435 |
coefficient, $\kappa$, relates the shear stress (e.g. along the |
436 |
$x$-axis) with the relative velocity of the fluid tangent to the |
437 |
interface: |
438 |
\begin{equation} |
439 |
j_z(p_x) = \kappa \left[v_x(fluid) - v_x(solid)\right] |
440 |
\end{equation} |
441 |
where $v_x(fluid)$ and $v_x(solid)$ are the velocities measured |
442 |
directly adjacent to the interface. Under ``stick'' boundary |
443 |
condition, $\Delta v_x|_\mathrm{interface} \rightarrow 0$, which leads |
444 |
to $\kappa\rightarrow\infty$. However, for ``slip'' boundary |
445 |
conditions at a solid-liquid interface, $\kappa$ becomes finite. To |
446 |
characterize the interfacial boundary conditions, the slip length, |
447 |
$\delta$, is defined by the ratio of the fluid-phase viscosity to the |
448 |
friction coefficient of the interface: |
449 |
\begin{equation} |
450 |
\delta = \frac{\eta}{\kappa} |
451 |
\end{equation} |
452 |
In ``stick'' boundary conditions, $\delta\rightarrow 0$, so $\delta$ |
453 |
is a measure of how ``slippery'' an interface is. Figure |
454 |
\ref{slipLength} illustrates how this quantity is defined and computed |
455 |
for a solid-liquid interface. |
456 |
|
457 |
\begin{figure} |
458 |
\includegraphics[width=\linewidth]{defDelta} |
459 |
\caption{The slip length ($\delta$) can be obtained from a velocity |
460 |
profile along an axis perpendicular to the interface. The data shown |
461 |
is for an Au / hexane interface -- the crystalline region (Au) is |
462 |
moving as a block (lower dots), while the measured velocity gradient |
463 |
in the hexane phase is discontinuous at the interface.} |
464 |
\label{slipLength} |
465 |
\end{figure} |
466 |
|
467 |
Since the VSS-RNEMD method can be applied for interfaces as well as |
468 |
for bulk materials, the shear stress is applied in an identical manner |
469 |
to the shear viscosity computations, e.g. by applying an unphysical |
470 |
momentum flux, $j_z(\vec{p})$. With the correct choice of $\vec{p}$ |
471 |
in the $x-y$ plane, one can compute friction coefficients and slip |
472 |
lengths for a number of different dragging vectors on a given slab. |
473 |
The corresponding velocity profiles can be obtained as shown in Figure |
474 |
\ref{slipLength}, in which the velocity gradients within the liquid |
475 |
phase and the velocity difference at the liquid-solid interface can be |
476 |
easily measured from saved simulation data. |
477 |
|
478 |
\section{Results and Discussions} |
479 |
\subsection{Lennard-Jones fluid} |
480 |
Our orthorhombic simulation cell for the Lennard-Jones fluid has |
481 |
identical parameters to previous work to facilitate |
482 |
comparison.\cite{kuang:164101} Thermal conductivities and shear |
483 |
viscosities were computed with the new VSS algorithm, and the results |
484 |
were compared with the previous NIVS algorithm. However, since the |
485 |
NIVS algorithm produces temperature anisotropy under shear stress, |
486 |
these results are also compared to the previous momentum swapping |
487 |
approaches. Table \ref{LJ} lists these values with various fluxes in |
488 |
reduced units. |
489 |
|
490 |
\begin{table*} |
491 |
\begin{minipage}{\linewidth} |
492 |
\begin{center} |
493 |
|
494 |
\caption{Thermal conductivity ($\lambda^*$) and shear viscosity |
495 |
($\eta^*$) (in reduced units) of a Lennard-Jones fluid at |
496 |
${\langle T^* \rangle = 0.72}$ and ${\rho^* = 0.85}$ computed |
497 |
at various momentum fluxes. The new VSS method yields similar |
498 |
results to previous RNEMD methods. All results are reported in |
499 |
reduced unit. Uncertainties are indicated in parentheses.} |
500 |
|
501 |
\begin{tabular}{cc|ccc|cc} |
502 |
\hline\hline |
503 |
\multicolumn{2}{c}{Momentum Exchange} & |
504 |
\multicolumn{3}{c}{$\lambda^*$} & |
505 |
\multicolumn{2}{c}{$\eta^*$} \\ |
506 |
\hline |
507 |
Swap Interval & Equivalent $J_z^*$ or $j_z^*(p_x)$ & Swapping & |
508 |
NIVS\cite{kuang:164101} & This work & Swapping & This work \\ |
509 |
\hline |
510 |
0.116 & 0.16 & 7.03(0.34) & 7.30(0.10) & 6.25(0.23) & 3.57(0.06) & 3.53(0.16)\\ |
511 |
0.232 & 0.09 & 7.03(0.14) & 6.95(0.09) & 8.02(0.56) & 3.64(0.05) & 3.43(0.17)\\ |
512 |
0.463 & 0.047 & 6.91(0.42) & 7.19(0.07) & 6.48(0.15) & 3.52(0.16) & 3.51(0.08)\\ |
513 |
0.926 & 0.024 & 7.52(0.15) & 7.19(0.28) & 7.02(0.13) & 3.72(0.05) & 3.79(0.11)\\ |
514 |
1.158 & 0.019 & 7.41(0.29) & 7.98(0.33) & 7.37(0.23) & 3.42(0.06) & 3.53(0.06)\\ |
515 |
\hline\hline |
516 |
\end{tabular} |
517 |
\label{LJ} |
518 |
\end{center} |
519 |
\end{minipage} |
520 |
\end{table*} |
521 |
|
522 |
\subsubsection{Thermal conductivity} |
523 |
Our thermal conductivity calculations yield comparable results to the |
524 |
previous NIVS algorithm. This indicates that the thermal gradients |
525 |
produced using this method are also quite close to previous RNEMD |
526 |
methods. Simulations with moderately large thermal fluxes tend to |
527 |
yield more reliable thermal gradients than under NIVS and thus avoid |
528 |
large errors. Extreme values for the applied thermal flux can |
529 |
introduce side effects such as non-linear temperature gradients and |
530 |
inadvertent phase transitions, so these are avoided. |
531 |
|
532 |
Since the scaling operation is isotropic in VSS, the temperature |
533 |
anisotropy that was observed under the earlier NIVS approach should be |
534 |
absent. Furthermore, the VSS method avoids unintended momentum flux |
535 |
when only thermal flux is being imposed. This was not always possible |
536 |
with swapping or NIVS approaches. The thermal energy exchange in |
537 |
swapping ($\vec{p}_i$ in slab {\it C} swapped with $\vec{p}_j$ in slab |
538 |
{\it H}) or NIVS (total slab momentum components $P^\alpha$ scaled to |
539 |
$\alpha P^\alpha$) do not achieve this unless thermal flux vanishes |
540 |
(i.e. $\vec{p}_i=\vec{p}_j$ or $\alpha=1$). In this sense, the VSS |
541 |
method achieves minimal perturbation to a simulation while imposing |
542 |
thermal flux. |
543 |
|
544 |
\subsubsection{Shear viscosity} |
545 |
Table \ref{LJ} also compares our shear viscosity results with the |
546 |
momentum swapping approach. Our calculations show that the VSS method |
547 |
predicted similar values for shear viscosities to the momentum |
548 |
swapping approach, as well as providing similar velocity gradient |
549 |
profiles. Moderately large momentum fluxes are helpful in reducing the |
550 |
errors in measured velocity gradients and thus the shear viscosity |
551 |
values. However, it should be noted that the momentum swapping |
552 |
approach tends to produce non-thermal velocity distributions in the |
553 |
two slabs.\cite{Maginn:2010} |
554 |
|
555 |
To show that the temperature is isotropic within each slab under VSS, |
556 |
we measured the three one-dimensional temperatures in each of the |
557 |
slabs (Figure \ref{tempXyz}). Note that the $x$-dimensional |
558 |
temperatures were calculated after subtracting the contribution from |
559 |
bulk (shearing) velocities of the slabs. The one-dimensional |
560 |
temperature profiles show no observable differences between the three |
561 |
box dimensions. This ensures that VSS automatically preserves |
562 |
temperature isotropy during shear viscosity calculations. |
563 |
|
564 |
\begin{figure} |
565 |
\includegraphics[width=\linewidth]{tempXyz} |
566 |
\caption{Unlike the previous NIVS algorithm, the VSS method does not |
567 |
produce thermal anisotropy under shearing stress. Temperature |
568 |
differences along the different box axes were significantly smaller |
569 |
than the error bars. Note that the two regions of elevated |
570 |
temperature are caused by the shear stress. This is the same |
571 |
frictional heating reported previously by Tenney and |
572 |
Maginn.\cite{Maginn:2010}} |
573 |
\label{tempXyz} |
574 |
\end{figure} |
575 |
|
576 |
The velocity distribution profiles were further tested by imposing a |
577 |
large shear stress. Figure \ref{vDist} demonstrates how the VSS method |
578 |
is able to maintain nearly ideal Maxwell-Boltzmann velocity |
579 |
distributions even under large imposed fluxes. Previous swapping |
580 |
methods tend to deplete particles with positive velocities in the |
581 |
negative velocity slab ({\it C}) and deplete particles with negative |
582 |
velocities in the {\it H} slab, leaving significant `notches' in the |
583 |
velocity distributions. The problematic velocity distributions become |
584 |
significant when the imposed-flux becomes so large that diffusion from |
585 |
neighboring slabs cannot offset the depletion. Simultaneously, |
586 |
abnormal peaks appear in the other regions of the velocity |
587 |
distributions as high (or low) momentum particles are introduced in to |
588 |
the corresponding slab. These nonthermal distributions limit |
589 |
application of the swapping approach in shear stress simulations. The |
590 |
VSS method avoids the above problematic distributions by altering the |
591 |
means of applying momentum flux. Comparatively, velocity distributions |
592 |
recorded from simulations with the VSS method is so close to the ideal |
593 |
Maxwell-Boltzmann distributions that no obvious difference is visible |
594 |
in Figure \ref{vDist}. |
595 |
|
596 |
\begin{figure} |
597 |
\centering |
598 |
\includegraphics[width=5.5in]{velDist} |
599 |
\caption{Velocity distributions that develop under VSS more closely |
600 |
resemble ideal Maxwell-Boltzmann distributions than earlier |
601 |
momentum-swapping RNEMD approaches. This data was obtained from |
602 |
Lennard-Jones simulations with $j_z^*(p_x)\sim 0.4$ (equivalent to a |
603 |
swapping interval of 20 time steps). This is a relatively large flux |
604 |
which demonstrates some of the non-thermal velocity distributions |
605 |
that can develop under swapping RNEMD.} |
606 |
\label{vDist} |
607 |
\end{figure} |
608 |
|
609 |
\subsection{Bulk SPC/E water} |
610 |
We also computed the thermal conductivity and shear viscosity of SPC/E |
611 |
water. A simulation cell with 1000 molecules was set up in a similar |
612 |
manner to Ref. \cite{kuang:164101}. For thermal conductivity |
613 |
calculations, measurements were taken to compare with previous RNEMD |
614 |
methods; for shear viscosity computations, simulations were run under |
615 |
a series of temperatures (with corresponding pressure relaxation using |
616 |
the isobaric-isothermal ensemble\cite{melchionna93}), and results were |
617 |
compared to available data from equilibrium molecular dynamics |
618 |
calculations.\cite{10.1063/1.3330544,Medina2011} Additionally, a |
619 |
simulation with {\it both} applied thermal and momentum gradients was |
620 |
carried out to map out shear viscosity as a function of temperature. |
621 |
|
622 |
\subsubsection{Thermal conductivity} |
623 |
Table \ref{spceThermal} summarizes the thermal conductivity of SPC/E |
624 |
under different temperatures and thermal gradients compared with the |
625 |
previous NIVS results\cite{kuang:164101} and experimental |
626 |
measurements.\cite{WagnerKruse} Note that no appreciable drift of |
627 |
total system energy or temperature was observed when the VSS method |
628 |
was applied, which indicates that our algorithm conserves total energy |
629 |
well for systems involving electrostatic interactions. |
630 |
|
631 |
Measurements using the VSS method established similar temperature |
632 |
gradients to the previous NIVS method. Our simulation results are in |
633 |
fairly good agreement with those from previous simulations. Both |
634 |
methods yield values in reasonable agreement with experimental |
635 |
values. Simulations using larger thermal gradients and those with |
636 |
longer gradient axes ($z$) have less measurable noise. |
637 |
|
638 |
\begin{table*} |
639 |
\begin{minipage}{\linewidth} |
640 |
\begin{center} |
641 |
|
642 |
\caption{Thermal conductivity of SPC/E water under various |
643 |
imposed thermal gradients. Uncertainties are indicated in |
644 |
parentheses.} |
645 |
|
646 |
\begin{tabular}{cc|ccc} |
647 |
\hline\hline |
648 |
$\langle T\rangle$ & $\langle dT/dz\rangle$ & \multicolumn{3}{c} |
649 |
{$\lambda (\mathrm{W m}^{-1} \mathrm{K}^{-1})$} \\ |
650 |
(K) & (K/\AA) & This work & NIVS\cite{kuang:164101} & |
651 |
Experiment\cite{WagnerKruse} \\ |
652 |
\hline |
653 |
300 & 0.8 & 0.815(0.027) & 0.770(0.008) & 0.61 \\ \hline |
654 |
318 & 0.8 & 0.801(0.024) & 0.750(0.032) & 0.64 \\ |
655 |
& 1.6 & 0.766(0.007) & 0.778(0.019) & \\ |
656 |
& 0.8 & 0.786(0.009)\footnote{Simulation with 2000 SPC/E molecules} & & \\ |
657 |
\hline\hline |
658 |
\end{tabular} |
659 |
\label{spceThermal} |
660 |
\end{center} |
661 |
\end{minipage} |
662 |
\end{table*} |
663 |
|
664 |
\subsubsection{Shear viscosity} |
665 |
Shear viscosity was computed for SPC/E water at a range of |
666 |
temperatures that span the liquidus range for water under atmospheric |
667 |
pressure. VSS-RNEMD simulations predict a similar trend of $\eta$ |
668 |
vs. $T$ to equilibrium molecular dynamics (EMD) results and are |
669 |
presented in table \ref{spceShear}. Our results show no significant |
670 |
differences from the earlier EMD calculations. Since each value |
671 |
reported using our method takes a single trajectory of simulation, |
672 |
instead of average from many trajectories when using EMD, the VSS |
673 |
method provides an efficient means for computing the shear viscosity. |
674 |
|
675 |
\begin{table*} |
676 |
\begin{minipage}{\linewidth} |
677 |
\begin{center} |
678 |
|
679 |
\caption{Computed shear viscosity of SPC/E water under different |
680 |
temperatures. Results are compared to those obtained with |
681 |
equilibrium molecular dynamics.\cite{Medina2011} Uncertainties |
682 |
are indicated in parentheses.} |
683 |
|
684 |
\begin{tabular}{cc|cc} |
685 |
\hline\hline |
686 |
$T$ & $\partial v_x/\partial z$ & \multicolumn{2}{c} |
687 |
{$\eta (\mathrm{mPa}\cdot\mathrm{s})$} \\ |
688 |
(K) & (10$^{10}$s$^{-1}$) & This work & Previous |
689 |
simulations\cite{Medina2011} \\ |
690 |
\hline |
691 |
273 & 1.12 & 1.218(0.004) & 1.282(0.048) \\ |
692 |
& 1.79 & 1.140(0.012) & \\ |
693 |
\hline |
694 |
303 & 2.09 & 0.646(0.008) & 0.643(0.019) \\ |
695 |
\hline |
696 |
318 & 2.50 & 0.536(0.007) & \\ |
697 |
& 5.25 & 0.510(0.007) & \\ |
698 |
& 2.82 & 0.474(0.003)\footnote{Simulation with 2000 SPC/E |
699 |
molecules.} & \\ |
700 |
\hline |
701 |
333 & 3.10 & 0.428(0.002) & 0.421(0.008) \\ |
702 |
\hline |
703 |
363 & 2.34 & 0.279(0.014) & 0.291(0.005) \\ |
704 |
& 4.26 & 0.306(0.001) & \\ |
705 |
\hline\hline |
706 |
\end{tabular} |
707 |
\label{spceShear} |
708 |
\end{center} |
709 |
\end{minipage} |
710 |
\end{table*} |
711 |
|
712 |
A more effective way to map out $\eta$ vs $T$ is to combine a momentum |
713 |
flux with a thermal flux. Figure \ref{Tvxdvdz} shows the thermal and |
714 |
velocity gradient in one such simulation. At different positions with |
715 |
different temperatures, the velocity gradient is not a constant but |
716 |
can be computed locally. With the data provided in Figure |
717 |
\ref{Tvxdvdz}, a series of $\eta$ is calculated as in Figure |
718 |
\ref{etaT} and a linear fit was performed to $\partial v_x/\partial z$ |
719 |
vs. $z$ so that the resulted $\eta$ can be present as a curve as |
720 |
well. For comparison, other results are also mapped in the figure. |
721 |
|
722 |
\begin{figure} |
723 |
\centering |
724 |
\includegraphics[width=5.5in]{tvxdvdz} |
725 |
\caption{With a combination of a thermal and a momentum flux, a |
726 |
simulation can exhibit both a temperature (top) and a velocity |
727 |
(middle) gradient. $\partial v_x/\partial z$ depends on the local |
728 |
temperature, so it varies along the gradient axis. This derivative |
729 |
can be computed using finite difference approximations (lower) and |
730 |
can be used to calculate $\eta$ vs $T$ (Figure \ref{etaT}).} |
731 |
\label{Tvxdvdz} |
732 |
\end{figure} |
733 |
|
734 |
From Figure \ref{etaT}, one can see that the generated curve agrees |
735 |
well with the above RNEMD simulations at different temperatures, as |
736 |
well as results reported using EMD |
737 |
methods\cite{10.1063/1.3330544,Medina2011} in much of the temperature |
738 |
range simulated. However, this curve has relatively large error in |
739 |
lower temperature regions and has some difference in predicting $\eta$ |
740 |
near 273K. Provided that this curve only takes one trajectory to |
741 |
generate, these results are of satisfactory efficiency and |
742 |
accuracy. Since previous work already pointed out that the SPC/E model |
743 |
tends to predict lower viscosity compared to experimental |
744 |
data,\cite{Medina2011} experimental comparison are not given here. |
745 |
|
746 |
\begin{figure} |
747 |
\includegraphics[width=\linewidth]{etaT} |
748 |
\caption{The temperature dependence of the shear viscosity generated |
749 |
by a {\it single} VSS-RNEMD simulation. Applying both thermal and |
750 |
momentum gradient predicts reasonable values in much of the |
751 |
temperature range being tested.} |
752 |
\label{etaT} |
753 |
\end{figure} |
754 |
|
755 |
\subsection{Interfacial frictions and slip lengths} |
756 |
Another attractive aspect of our method is the ability to apply |
757 |
momentum and/or thermal flux in nonhomogeneous systems, where |
758 |
molecules of different identities (or phases) are segregated in |
759 |
different regions. We have previously studied the interfacial thermal |
760 |
transport of a series of metal gold-liquid |
761 |
surfaces\cite{kuang:164101,kuang:AuThl}, and would like to further |
762 |
investigate the relationship between this phenomenon and the |
763 |
interfacial frictions. |
764 |
|
765 |
Table \ref{etaKappaDelta} includes these computations and previous |
766 |
calculations of corresponding interfacial thermal conductance. For |
767 |
bare Au(111) surfaces, slip boundary conditions were observed for both |
768 |
organic and aqueous liquid phases, corresponding to previously |
769 |
computed low interfacial thermal conductance. In comparison, the |
770 |
butanethiol covered Au(111) surface appeared to be sticky to the |
771 |
organic liquid layers in our simulations. We have reported conductance |
772 |
enhancement effect for this surface capping agent,\cite{kuang:AuThl} |
773 |
and these observations have a qualitative agreement with the thermal |
774 |
conductance results. This agreement also supports discussions on the |
775 |
relationship between surface wetting and slip effect and thermal |
776 |
conductance of the |
777 |
interface.\cite{PhysRevLett.82.4671,doi:10.1080/0026897031000068578,garde:PhysRevLett2009} |
778 |
|
779 |
\begin{table*} |
780 |
\begin{minipage}{\linewidth} |
781 |
\begin{center} |
782 |
|
783 |
\caption{Computed interfacial friction coefficient values for |
784 |
interfaces with various components for liquid and solid |
785 |
phase. Error estimates are indicated in parentheses.} |
786 |
|
787 |
\begin{tabular}{llcccccc} |
788 |
\hline\hline |
789 |
Solid & Liquid & $T$ & $j_z(p_x)$ & $\eta_{liquid}$ & $\kappa$ |
790 |
& $\delta$ & $G$\footnote{References \cite{kuang:AuThl} and |
791 |
\cite{kuang:164101}.} \\ |
792 |
surface & molecules & K & MPa & mPa$\cdot$s & |
793 |
10$^4$Pa$\cdot$s/m & nm & MW/m$^2$/K \\ |
794 |
\hline |
795 |
Au(111) & hexane & 200 & 1.08 & 0.197(0.009) & 5.30(0.36) & |
796 |
3.72(0.25) & 46.5 \\ |
797 |
& & & 2.15 & 0.141(0.002) & 5.31(0.26) & |
798 |
2.76(0.14) & \\ |
799 |
Au-SC$_4$H$_9$ & hexane & 200 & 2.16 & 0.286(0.019) & $\infty$ |
800 |
& 0 & 131 \\ |
801 |
& & & 5.39 & 0.320(0.006) & $\infty$ |
802 |
& 0 & \\ |
803 |
\hline |
804 |
Au(111) & toluene & 200 & 1.08 & 0.722(0.035) & 15.7(0.7) & |
805 |
4.60(0.22) & 70.1 \\ |
806 |
& & & 2.16 & 0.544(0.030) & 11.2(0.5) & |
807 |
4.86(0.27) & \\ |
808 |
Au-SC$_4$H$_9$ & toluene & 200 & 5.39 & 0.980(0.057) & |
809 |
$\infty$ & 0 & 187 \\ |
810 |
& & & 10.8 & 0.995(0.005) & |
811 |
$\infty$ & 0 & \\ |
812 |
\hline |
813 |
Au(111) & water & 300 & 1.08 & 0.399(0.050) & 1.928(0.022) & |
814 |
20.7(2.6) & 1.65 \\ |
815 |
& & & 2.16 & 0.794(0.255) & 1.895(0.003) & |
816 |
41.9(13.5) & \\ |
817 |
\hline |
818 |
ice(basal) & water & 225 & 19.4 & 15.8(0.2) & $\infty$ & 0 & \\ |
819 |
\hline\hline |
820 |
\end{tabular} |
821 |
\label{etaKappaDelta} |
822 |
\end{center} |
823 |
\end{minipage} |
824 |
\end{table*} |
825 |
|
826 |
An interesting effect alongside the surface friction change is |
827 |
observed on the shear viscosity of liquids in the regions close to the |
828 |
solid surface. In our results, $\eta$ measured near a ``slip'' surface |
829 |
tends to be smaller than that near a ``stick'' surface. This may |
830 |
suggest the influence from an interface on the dynamic properties of |
831 |
liquid within its neighbor regions. It is known that diffusions of |
832 |
solid particles in liquid phase is affected by their surface |
833 |
conditions (stick or slip boundary).\cite{10.1063/1.1610442} Our |
834 |
observations could provide a support to this phenomenon. |
835 |
|
836 |
In addition to these previously studied interfaces, we attempt to |
837 |
construct ice-water interfaces and the basal plane of ice lattice was |
838 |
studied here. In contrast to the Au(111)/water interface, where the |
839 |
friction coefficient is substantially small and large slip effect |
840 |
presents, the ice/liquid water interface demonstrates strong |
841 |
solid-liquid interactions and appears to be sticky. The supercooled |
842 |
liquid phase is an order of magnitude more viscous than measurements |
843 |
in previous section. It would be of interst to investigate the effect |
844 |
of different ice lattice planes (such as prism and other surfaces) on |
845 |
interfacial friction and the corresponding liquid viscosity. |
846 |
|
847 |
\section{Conclusions} |
848 |
Our simulations demonstrate the validity of our method in RNEMD |
849 |
computations of thermal conductivity and shear viscosity in atomic and |
850 |
molecular liquids. Our method maintains thermal velocity distributions |
851 |
and avoids thermal anisotropy in previous NIVS shear stress |
852 |
simulations, as well as retains attractive features of previous RNEMD |
853 |
methods. There is no {\it a priori} restrictions to the method to be |
854 |
applied in various ensembles, so prospective applications to |
855 |
extended-system methods are possible. |
856 |
|
857 |
Our method is capable of effectively imposing thermal and/or momentum |
858 |
flux accross an interface. This facilitates studies that relates |
859 |
dynamic property measurements to the chemical details of an |
860 |
interface. Therefore, investigations can be carried out to |
861 |
characterize interfacial interactions using the method. |
862 |
|
863 |
Another attractive feature of our method is the ability of |
864 |
simultaneously introducing thermal and momentum gradients in a |
865 |
system. This facilitates us to effectively map out the shear viscosity |
866 |
with respect to a range of temperature in single trajectory of |
867 |
simulation with satisafactory accuracy. Complex systems that involve |
868 |
thermal and momentum gradients might potentially benefit from |
869 |
this. For example, the Soret effects under a velocity gradient might |
870 |
be models of interest to purification and separation researches. |
871 |
|
872 |
\section{Acknowledgments} |
873 |
Support for this project was provided by the National Science |
874 |
Foundation under grant CHE-0848243. Computational time was provided by |
875 |
the Center for Research Computing (CRC) at the University of Notre |
876 |
Dame. |
877 |
|
878 |
\newpage |
879 |
|
880 |
\bibliography{stokes} |
881 |
|
882 |
\end{doublespace} |
883 |
\end{document} |