| 53 |
|
scaling (NIVS) on the molecules in specific regions of a system, it |
| 54 |
|
is possible to impose momentum or thermal flux between regions of a |
| 55 |
|
simulation while conserving the linear momentum and total energy of |
| 56 |
< |
the system. To test the methods, we have computed the thermal |
| 56 |
> |
the system. To test the method, we have computed the thermal |
| 57 |
|
conductivity of model liquid and solid systems as well as the |
| 58 |
|
interfacial thermal conductivity of a metal-water interface. We |
| 59 |
|
find that the NIVS-RNEMD improves the problematic velocity |
| 114 |
|
relations\cite{daivis:541,mondello:9327} or the Helfand moment |
| 115 |
|
approach of Viscardy {\it et |
| 116 |
|
al.}\cite{Viscardy:2007bh,Viscardy:2007lq}) because these rely on |
| 117 |
< |
computing difficult to measure quantities. |
| 117 |
> |
computing and integrating long-time correlation functions that are |
| 118 |
> |
subject to noise issues. |
| 119 |
|
|
| 120 |
|
Another attractive feature of RNEMD is that it conserves both total |
| 121 |
|
linear momentum and total energy during the swaps (as long as the two |
| 290 |
|
isotropic fluid as possible. With this in mind, we would like the |
| 291 |
|
kinetic energies in the three different directions could become as |
| 292 |
|
close as each other as possible after each scaling. Simultaneously, |
| 293 |
< |
one would also like each scaling as gentle as possible, i.e. ${x,y,z |
| 294 |
< |
\rightarrow 1}$, in order to avoid large perturbation to the system. |
| 295 |
< |
To do this, we pick the intersection point which maintains the three |
| 296 |
< |
scaling variables ${x, y, z}$ as well as the ratio of kinetic energies |
| 297 |
< |
${K_c^z/K_c^x, K_c^z/K_c^y}$ as close as possible to 1. |
| 293 |
> |
one would also like each scaling to be as gentle as possible, i.e. |
| 294 |
> |
${x,y,z \rightarrow 1}$, in order to avoid large perturbations to the |
| 295 |
> |
system. To do this, we pick the intersection point which maintains |
| 296 |
> |
the three scaling variables ${x, y, z}$ as well as the ratio of |
| 297 |
> |
kinetic energies ${K_c^z/K_c^x, K_c^z/K_c^y}$ as close as possible to |
| 298 |
> |
1. |
| 299 |
|
|
| 300 |
|
After the valid scaling parameters are arrived at by solving geometric |
| 301 |
|
intersection problems in $x, y, z$ space in order to obtain cold slab |
| 385 |
|
NIVS-RNEMD simulation was carried out using a target momentum flux set |
| 386 |
|
to produce the same flux experienced in the swapping simulation. |
| 387 |
|
|
| 388 |
< |
To test the temperature homogeneity, directional momentum and |
| 389 |
< |
temperature distributions were accumulated for molecules in each of |
| 390 |
< |
the slabs. Transport coefficients were computed using the temperature |
| 391 |
< |
(and momentum) gradients across the $z$-axis as well as the total |
| 392 |
< |
momentum flux the system experienced during data collection portion of |
| 393 |
< |
the simulation. |
| 388 |
> |
To test the temperature homogeneity, momentum and temperature |
| 389 |
> |
distributions (for all three dimensions) were accumulated for |
| 390 |
> |
molecules in each of the slabs. Transport coefficients were computed |
| 391 |
> |
using the temperature (and momentum) gradients across the $z$-axis as |
| 392 |
> |
well as the total momentum flux the system experienced during data |
| 393 |
> |
collection portion of the simulation. |
| 394 |
|
|
| 395 |
|
\subsection{Shear viscosities} |
| 396 |
|
|
| 431 |
|
be approximated as: |
| 432 |
|
|
| 433 |
|
\begin{equation} |
| 434 |
< |
G = \frac{E_{total}}{2 t L_x L_y \left( \langle T_{gold}\rangle - |
| 435 |
< |
\langle T_{water}\rangle \right)} |
| 434 |
> |
G = \frac{E_{total}}{2 t L_x L_y \left( \langle T_\mathrm{hot}\rangle - |
| 435 |
> |
\langle T_\mathrm{cold}\rangle \right)} |
| 436 |
|
\label{interfaceCalc} |
| 437 |
|
\end{equation} |
| 438 |
|
where ${E_{total}}$ is the imposed non-physical kinetic energy |
| 439 |
< |
transfer and ${\langle T_{gold}\rangle}$ and ${\langle |
| 440 |
< |
T_{water}\rangle}$ are the average observed temperature of gold and |
| 441 |
< |
water phases respectively. If the interfacial conductance is {\it |
| 442 |
< |
not} small, it is also be possible to compute the interfacial |
| 443 |
< |
thermal conductivity using this method utilizing the change in the |
| 444 |
< |
slope of the thermal gradient ($\partial^2 \langle T \rangle / \partial |
| 445 |
< |
z^2$) at the interface. |
| 439 |
> |
transfer and ${\langle T_\mathrm{hot}\rangle}$ and ${\langle |
| 440 |
> |
T_\mathrm{cold}\rangle}$ are the average observed temperature of the |
| 441 |
> |
two separated phases. If the interfacial conductance is {\it not} |
| 442 |
> |
small, it would also be possible to compute the interfacial thermal |
| 443 |
> |
conductivity using this method by computing the change in the slope of |
| 444 |
> |
the thermal gradient ($\partial^2 \langle T \rangle / |
| 445 |
> |
\partial z^2$) at the interface. |
| 446 |
|
|
| 447 |
|
\section{Results} |
| 448 |
|
|
| 466 |
|
Our thermal conductivity calculations show that the NIVS method agrees |
| 467 |
|
well with the swapping method. Five different swap intervals were |
| 468 |
|
tested (Table \ref{LJ}). Similar thermal gradients were observed with |
| 469 |
< |
similar thermal flux under the two different methods (Figure |
| 469 |
> |
similar thermal flux under the two different methods (Fig. |
| 470 |
|
\ref{thermalGrad}). Furthermore, the 1-d temperature profiles showed |
| 471 |
< |
no observable differences between the $x$, $y$ and $z$ axes (Figure |
| 472 |
< |
\ref{thermalGrad} c), so even though we are using a non-isotropic |
| 473 |
< |
scaling method, none of the three directions are experience |
| 474 |
< |
disproportionate heating due to the velocity scaling. |
| 471 |
> |
no observable differences between the $x$, $y$ and $z$ axes (lower |
| 472 |
> |
panel of Fig. \ref{thermalGrad}), so even though we are using a |
| 473 |
> |
non-isotropic scaling method, none of the three directions are |
| 474 |
> |
experience disproportionate heating due to the velocity scaling. |
| 475 |
|
|
| 476 |
|
\begin{table*} |
| 477 |
|
\begin{minipage}{\linewidth} |
| 521 |
|
|
| 522 |
|
\subsubsection*{Velocity Distributions} |
| 523 |
|
|
| 524 |
< |
During these simulations, velocities were recorded every 1000 steps |
| 525 |
< |
and were used to produce distributions of both velocity and speed in |
| 526 |
< |
each of the slabs. From these distributions, we observed that under |
| 527 |
< |
relatively high non-physical kinetic energy flux, the speed of |
| 528 |
< |
molecules in slabs where swapping occured could deviate from the |
| 524 |
> |
To test the effects on the velocity distributions, we accumulated |
| 525 |
> |
velocities every 100 steps and produced distributions of both velocity |
| 526 |
> |
and speed in each of the slabs. From these distributions, we observed |
| 527 |
> |
that under high non-physical kinetic energy flux, the speed of |
| 528 |
> |
molecules in slabs where {\it swapping} occured could deviate from the |
| 529 |
|
Maxwell-Boltzmann distribution. This behavior was also noted by Tenney |
| 530 |
|
and Maginn.\cite{Maginn:2010} Figure \ref{thermalHist} shows how these |
| 531 |
|
distributions deviate from an ideal distribution. In the ``hot'' slab, |
| 532 |
|
the probability density is notched at low speeds and has a substantial |
| 533 |
< |
shoulder at higher speeds relative to the ideal MB distribution. In |
| 534 |
< |
the cold slab, the opposite notching and shouldering occurs. This |
| 535 |
< |
phenomenon is more obvious at higher swapping rates. |
| 533 |
> |
shoulder at higher speeds relative to the ideal distribution. In the |
| 534 |
> |
cold slab, the opposite notching and shouldering occurs. This |
| 535 |
> |
phenomenon is more obvious at high swapping rates. |
| 536 |
|
|
| 537 |
|
The peak of the velocity distribution is substantially flattened in |
| 538 |
|
the hot slab, and is overly sharp (with truncated wings) in the cold |
| 539 |
|
slab. This problem is rooted in the mechanism of the swapping method. |
| 540 |
|
Continually depleting low (high) speed particles in the high (low) |
| 541 |
< |
temperature slab is not complemented by diffusions of low (high) speed |
| 542 |
< |
particles from neighboring slabs, unless the swapping rate is |
| 541 |
> |
temperature slab is not complemented by diffusion of low (high) speed |
| 542 |
> |
particles from neighboring slabs unless the swapping rate is |
| 543 |
|
sufficiently small. Simutaneously, surplus low speed particles in the |
| 544 |
< |
low temperature slab do not have sufficient time to diffuse to |
| 545 |
< |
neighboring slabs. Since the thermal exchange rate must reach a |
| 546 |
< |
minimum level to produce an observable thermal gradient, the |
| 547 |
< |
swapping-method RNEMD has a relatively narrow choice of exchange times |
| 546 |
< |
that can be utilized. |
| 544 |
> |
cold slab do not have sufficient time to diffuse to neighboring slabs. |
| 545 |
> |
Since the thermal exchange rate must reach a minimum level to produce |
| 546 |
> |
an observable thermal gradient, the swapping-method RNEMD has a |
| 547 |
> |
relatively narrow choice of exchange times that can be utilized. |
| 548 |
|
|
| 549 |
|
For comparison, NIVS-RNEMD produces a speed distribution closer to the |
| 550 |
< |
Maxwell-Boltzmann curve (Figure \ref{thermalHist}). The reason for |
| 551 |
< |
this is simple; upon velocity scaling, a Gaussian distribution remains |
| 550 |
> |
Maxwell-Boltzmann curve (Fig. \ref{thermalHist}). The reason for this |
| 551 |
> |
is simple; upon velocity scaling, a Gaussian distribution remains |
| 552 |
|
Gaussian. Although a single scaling move is non-isotropic in three |
| 553 |
|
dimensions, our criteria for choosing a set of scaling coefficients |
| 554 |
|
helps maintain the distributions as close to isotropic as possible. |
| 555 |
< |
Consequently, NIVS-RNEMD is able to impose an unphysical thermal flux |
| 556 |
< |
as the previous RNEMD methods but without large perturbations to the |
| 557 |
< |
velocity distributions in the two slabs. |
| 555 |
> |
Consequently, NIVS-RNEMD is able to impose a non-physical thermal flux |
| 556 |
> |
without large perturbations to the velocity distributions in the two |
| 557 |
> |
slabs. |
| 558 |
|
|
| 559 |
|
\begin{figure} |
| 560 |
|
\includegraphics[width=\linewidth]{thermalHist} |
| 573 |
|
|
| 574 |
|
\subsubsection*{Shear Viscosity} |
| 575 |
|
Our calculations (Table \ref{LJ}) show that velocity-scaling RNEMD |
| 576 |
< |
predicted comparable shear viscosities to swap RNEMD method. The |
| 577 |
< |
average molecular momentum gradients of these samples are shown in |
| 578 |
< |
Figure \ref{shear} (a) and (b). |
| 576 |
> |
predicted similar values for shear viscosities to the swapping RNEMD |
| 577 |
> |
method. The average molecular momentum gradients of these samples are |
| 578 |
> |
shown in the upper two panels of Fig. \ref{shear}. |
| 579 |
|
|
| 580 |
|
\begin{figure} |
| 581 |
|
\includegraphics[width=\linewidth]{shear} |
| 589 |
|
|
| 590 |
|
Observations of the three one-dimensional temperatures in each of the |
| 591 |
|
slabs shows that NIVS-RNEMD does produce some thermal anisotropy, |
| 592 |
< |
particularly in the hot and cold slabs. Figure \ref{shear} (c) |
| 593 |
< |
indicates that with a relatively large imposed momentum flux, |
| 594 |
< |
$j_z(p_x)$, the $x$ direction approaches a different temperature from |
| 595 |
< |
the $y$ and $z$ directions in both the hot and cold bins. This is an |
| 596 |
< |
artifact of the scaling constraints. After the momentum gradient has |
| 597 |
< |
been established, $P_c^x < 0$. Consequently, the scaling factor $x$ |
| 598 |
< |
is nearly always greater than one in order to satisfy the constraints. |
| 599 |
< |
This will continually increase the kinetic energy in $x$-dimension, |
| 600 |
< |
$K_c^x$. If there is not enough time for the kinetic energy to |
| 601 |
< |
exchange among different directions and different slabs, the system |
| 602 |
< |
will exhibit the observed thermal anisotropy in the hot and cold bins. |
| 592 |
> |
particularly in the hot and cold slabs. The lower panel of Fig. |
| 593 |
> |
\ref{shear} indicates that with a relatively large imposed momentum |
| 594 |
> |
flux, $j_z(p_x)$, the $x$ direction approaches a different temperature |
| 595 |
> |
from the $y$ and $z$ directions in both the hot and cold bins. This |
| 596 |
> |
is an artifact of the scaling constraints. After a momentum gradient |
| 597 |
> |
has been established, $P_c^x$ is always less than zero. Consequently, |
| 598 |
> |
the scaling factor $x$ is always greater than one in order to satisfy |
| 599 |
> |
the constraints. This will continually increase the kinetic energy in |
| 600 |
> |
$x$-dimension, $K_c^x$. If there is not enough time for the kinetic |
| 601 |
> |
energy to exchange among different directions and different slabs, the |
| 602 |
> |
system will exhibit the observed thermal anisotropy in the hot and |
| 603 |
> |
cold bins. |
| 604 |
|
|
| 605 |
|
Although results between scaling and swapping methods are comparable, |
| 606 |
|
the inherent temperature anisotropy does make NIVS-RNEMD method less |
| 607 |
|
attractive than swapping RNEMD for shear viscosity calculations. We |
| 608 |
< |
note that this problem appears only when momentum flux is applied, and |
| 609 |
< |
does not appear in thermal flux calculations. |
| 608 |
> |
note that this problem appears only when a large {\it linear} momentum |
| 609 |
> |
flux is applied, and does not appear in {\it thermal} flux |
| 610 |
> |
calculations. |
| 611 |
|
|
| 612 |
|
\subsection{Bulk SPC/E water} |
| 613 |
|
|
| 614 |
|
We compared the thermal conductivity of SPC/E water using NIVS-RNEMD |
| 615 |
|
to the work of Bedrov {\it et al.}\cite{Bedrov:2000}, which employed |
| 616 |
|
the original swapping RNEMD method. Bedrov {\it et |
| 617 |
< |
al.}\cite{Bedrov:2000} argued that exchange of the molecule |
| 618 |
< |
center-of-mass velocities instead of single atom velocities in a |
| 619 |
< |
molecule conserves the total kinetic energy and linear momentum. This |
| 620 |
< |
principle is also adopted Fin our simulations. Scaling was applied to |
| 621 |
< |
the center-of-mass velocities of the rigid bodies of SPC/E model water |
| 619 |
< |
molecules. |
| 617 |
> |
al.}\cite{Bedrov:2000} argued that exchange of the molecular |
| 618 |
> |
center-of-mass velocities instead of single atom velocities conserves |
| 619 |
> |
the total kinetic energy and linear momentum. This principle is also |
| 620 |
> |
adopted in our simulations. Scaling was applied to the center-of-mass |
| 621 |
> |
velocities of SPC/E water molecules. |
| 622 |
|
|
| 623 |
|
To construct the simulations, a simulation box consisting of 1000 |
| 624 |
|
molecules were first equilibrated under ambient pressure and |
| 719 |
|
in volume. Following adjustment of the box volume, equilibrations in |
| 720 |
|
both the canonical and microcanonical ensembles were carried out. With |
| 721 |
|
the simulation cell divided evenly into 10 slabs, different thermal |
| 722 |
< |
gradients were established by applying a set of target thermal |
| 721 |
< |
transfer fluxes. |
| 722 |
> |
gradients were established by applying a set of target thermal fluxes. |
| 723 |
|
|
| 724 |
|
The results for the thermal conductivity of gold are shown in Table |
| 725 |
|
\ref{AuThermal}. In these calculations, the end and middle slabs were |
| 726 |
< |
excluded in thermal gradient linear regession. EAM predicts slightly |
| 727 |
< |
larger thermal conductivities than QSC. However, both values are |
| 728 |
< |
smaller than experimental value by a factor of more than 200. This |
| 726 |
> |
excluded from the thermal gradient linear regession. EAM predicts |
| 727 |
> |
slightly larger thermal conductivities than QSC. However, both values |
| 728 |
> |
are smaller than experimental value by a factor of more than 200. This |
| 729 |
|
behavior has been observed previously by Richardson and Clancy, and |
| 730 |
|
has been attributed to the lack of electronic contribution in these |
| 731 |
|
force fields.\cite{Clancy:1992} It should be noted that the density of |
| 786 |
|
relax under ambient temperature and pressure. A separate (but |
| 787 |
|
identically sized) box of SPC/E water was also equilibrated at ambient |
| 788 |
|
conditions. The two boxes were combined by removing all water |
| 789 |
< |
molecules within 3 \AA radius of any gold atom. The final |
| 789 |
> |
molecules within 3 \AA~ radius of any gold atom. The final |
| 790 |
|
configuration contained 1862 SPC/E water molecules. |
| 791 |
|
|
| 792 |
|
The Spohr potential was adopted in depicting the interaction between |
| 800 |
|
phase.\cite{Cahill:793,RevModPhys.61.605} Experiments on the thermal |
| 801 |
|
conductivity of gold nanoparticles and nanorods in solvent and in |
| 802 |
|
glass cages have predicted values for $G$ between 100 and 350 |
| 803 |
< |
(MW/m$^2$/K). The experiments typically have multiple gold surfaces |
| 804 |
< |
that have been protected by a capping agent (citrate or CTAB) or which |
| 805 |
< |
are in direct contact with various glassy solids. In these cases, the |
| 806 |
< |
acoustic impedance mismatch would be substantially reduced, leading to |
| 807 |
< |
much higher interfacial conductances. It is also possible, however, |
| 808 |
< |
that the lack of electronic effects that gives rise to the low thermal |
| 809 |
< |
conductivity of EAM gold is also causing a low reading for this |
| 810 |
< |
particular interface. |
| 803 |
> |
(MW/m$^2$/K), two orders of magnitude larger than the value reported |
| 804 |
> |
here. The experiments typically have multiple surfaces that have been |
| 805 |
> |
protected by ionic surfactants, either |
| 806 |
> |
citrate\cite{Wilson:2002uq,plech:195423} or cetyltrimethylammonium |
| 807 |
> |
bromide (CTAB), or which are in direct contact with various glassy |
| 808 |
> |
solids. In these cases, the acoustic impedance mismatch would be |
| 809 |
> |
substantially reduced, leading to much higher interfacial |
| 810 |
> |
conductances. It is also possible, however, that the lack of |
| 811 |
> |
electronic effects that gives rise to the low thermal conductivity of |
| 812 |
> |
EAM gold is also causing a low reading for this particular interface. |
| 813 |
|
|
| 814 |
< |
Under this low thermal conductance, both gold and water phase have |
| 815 |
< |
sufficient time to eliminate temperature difference inside |
| 816 |
< |
respectively (Figure \ref{interface} b). With indistinguishable |
| 817 |
< |
temperature difference within respective phase, it is valid to assume |
| 818 |
< |
that the temperature difference between gold and water on surface |
| 819 |
< |
would be approximately the same as the difference between the gold and |
| 820 |
< |
water phase. This assumption enables convenient calculation of $G$ |
| 818 |
< |
using Eq. \ref{interfaceCalc} instead of measuring temperatures of |
| 819 |
< |
thin layer of water and gold close enough to surface, which would have |
| 820 |
< |
greater fluctuation and lower accuracy. Reported results (Table |
| 821 |
< |
\ref{interfaceRes}) are of two orders of magnitude smaller than our |
| 822 |
< |
calculations on homogeneous systems, and thus have larger relative |
| 823 |
< |
errors than our calculation results on homogeneous systems. |
| 814 |
> |
Under this low thermal conductance, both gold and water phases have |
| 815 |
> |
sufficient time to eliminate local temperature differences (Fig. |
| 816 |
> |
\ref{interface}). With flat thermal profiles within each phase, it is |
| 817 |
> |
valid to assume that the temperature difference between gold and water |
| 818 |
> |
surfaces would be approximately the same as the difference between the |
| 819 |
> |
gold and water bulk regions. This assumption enables convenient |
| 820 |
> |
calculation of $G$ using Eq. \ref{interfaceCalc}. |
| 821 |
|
|
| 822 |
|
\begin{figure} |
| 823 |
|
\includegraphics[width=\linewidth]{interface} |
| 857 |
|
|
| 858 |
|
|
| 859 |
|
\section{Conclusions} |
| 860 |
< |
NIVS-RNEMD simulation method is developed and tested on various |
| 861 |
< |
systems. Simulation results demonstrate its validity in thermal |
| 862 |
< |
conductivity calculations, from Lennard-Jones fluid to multi-atom |
| 863 |
< |
molecule like water and metal crystals. NIVS-RNEMD improves |
| 864 |
< |
non-Boltzmann-Maxwell distributions, which exist inb previous RNEMD |
| 865 |
< |
methods. Furthermore, it develops a valid means for unphysical thermal |
| 860 |
> |
|
| 861 |
> |
Our simulations demonstrate that validity of non-isotropic velocity |
| 862 |
> |
scaling (NIVS) in RNEMD calculations of thermal conductivity in atomic |
| 863 |
> |
and molecular liquids and solids. NIVS-RNEMD improves the problematic |
| 864 |
> |
velocity distributions which can develop in other RNEMD methods. |
| 865 |
> |
Furthermore, it provides a means for carrying out non-physical thermal |
| 866 |
|
transfer between different species of molecules, and thus extends its |
| 867 |
< |
applicability to interfacial systems. Our calculation of gold/water |
| 868 |
< |
interfacial thermal conductivity demonstrates this advantage over |
| 869 |
< |
previous RNEMD methods. NIVS-RNEMD has also limited application on |
| 870 |
< |
shear viscosity calculations, but could cause temperature difference |
| 871 |
< |
among different dimensions under high momentum flux. Modification is |
| 872 |
< |
necessary to extend the applicability of NIVS-RNEMD in shear viscosity |
| 873 |
< |
calculations. |
| 867 |
> |
applicability to interfacial systems. Our calculation of the gold / |
| 868 |
> |
water interfacial thermal conductivity demonstrates this advantage |
| 869 |
> |
over previous RNEMD methods. NIVS-RNEMD also has limited applications |
| 870 |
> |
for shear viscosity calculations, but has the potential to cause |
| 871 |
> |
temperature anisotropy under high momentum fluxes. Further work will |
| 872 |
> |
be necessary to eliminate the one-dimensional heating if shear |
| 873 |
> |
viscosities are required. |
| 874 |
|
|
| 875 |
|
\section{Acknowledgments} |
| 876 |
|
The authors would like to thank Craig Tenney and Ed Maginn for many |
