--- trunk/nivsRnemd/nivsRnemd.tex 2010/03/09 20:37:52 3565 +++ trunk/nivsRnemd/nivsRnemd.tex 2010/08/12 14:48:03 3630 @@ -6,10 +6,13 @@ \usepackage{caption} %\usepackage{tabularx} \usepackage{graphicx} +\usepackage{multirow} %\usepackage{booktabs} %\usepackage{bibentry} %\usepackage{mathrsfs} -\usepackage[ref]{overcite} +%\usepackage[ref]{overcite} +\usepackage[square, comma, sort&compress]{natbib} +\usepackage{url} \pagestyle{plain} \pagenumbering{arabic} \oddsidemargin 0.0cm \evensidemargin 0.0cm \topmargin -21pt \headsep 10pt \textheight 9.0in \textwidth 6.5in \brokenpenalty=10000 @@ -19,7 +22,9 @@ \setlength{\abovecaptionskip}{20 pt} \setlength{\belowcaptionskip}{30 pt} -\renewcommand\citemid{\ } % no comma in optional referenc note +%\renewcommand\citemid{\ } % no comma in optional referenc note +\bibpunct{[}{]}{,}{s}{}{;} +\bibliographystyle{aip} \begin{document} @@ -38,7 +43,21 @@ Notre Dame, Indiana 46556} \begin{doublespace} \begin{abstract} - + We present a new method for introducing stable non-equilibrium + velocity and temperature gradients in molecular dynamics simulations + of heterogeneous systems. This method extends earlier Reverse + Non-Equilibrium Molecular Dynamics (RNEMD) methods which use + momentum exchange swapping moves. The standard swapping moves can + create non-thermal velocity distributions and are difficult to use + for interfacial calculations. By using non-isotropic velocity + scaling (NIVS) on the molecules in specific regions of a system, it + is possible to impose momentum or thermal flux between regions of a + simulation while conserving the linear momentum and total energy of + the system. To test the methods, we have computed the thermal + conductivity of model liquid and solid systems as well as the + interfacial thermal conductivity of a metal-water interface. We + find that the NIVS-RNEMD improves the problematic velocity + distributions that develop in other RNEMD methods. \end{abstract} \newpage @@ -49,39 +68,51 @@ Notre Dame, Indiana 46556} % BODY OF TEXT %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% - - \section{Introduction} The original formulation of Reverse Non-equilibrium Molecular Dynamics (RNEMD) obtains transport coefficients (thermal conductivity and shear viscosity) in a fluid by imposing an artificial momentum flux between two thin parallel slabs of material that are spatially separated in the simulation cell.\cite{MullerPlathe:1997xw,ISI:000080382700030} The -artificial flux is typically created by periodically ``swapping'' either -the entire momentum vector $\vec{p}$ or single components of this -vector ($p_x$) between molecules in each of the two slabs. If the two -slabs are separated along the z coordinate, the imposed flux is either -directional ($j_z(p_x)$) or isotropic ($J_z$), and the response of a -simulated system to the imposed momentum flux will typically be a -velocity or thermal gradient. The transport coefficients (shear -viscosity and thermal conductivity) are easily obtained by assuming -linear response of the system, +artificial flux is typically created by periodically ``swapping'' +either the entire momentum vector $\vec{p}$ or single components of +this vector ($p_x$) between molecules in each of the two slabs. If +the two slabs are separated along the $z$ coordinate, the imposed flux +is either directional ($j_z(p_x)$) or isotropic ($J_z$), and the +response of a simulated system to the imposed momentum flux will +typically be a velocity or thermal gradient (Fig. \ref{thermalDemo}). +The shear viscosity ($\eta$) and thermal conductivity ($\lambda$) are +easily obtained by assuming linear response of the system, \begin{eqnarray} j_z(p_x) & = & -\eta \frac{\partial v_x}{\partial z}\\ -J & = & \lambda \frac{\partial T}{\partial z} +J_z & = & \lambda \frac{\partial T}{\partial z} \end{eqnarray} -RNEMD has been widely used to provide computational estimates of thermal -conductivities and shear viscosities in a wide range of materials, -from liquid copper to monatomic liquids to molecular fluids -(e.g. ionic liquids).\cite{ISI:000246190100032} +RNEMD has been widely used to provide computational estimates of +thermal conductivities and shear viscosities in a wide range of +materials, from liquid copper to both monatomic and molecular fluids +(e.g. ionic +liquids).\cite{Bedrov:2000-1,Bedrov:2000,Muller-Plathe:2002,ISI:000184808400018,ISI:000231042800044,Maginn:2007,Muller-Plathe:2008,ISI:000258460400020,ISI:000258840700015,ISI:000261835100054} -RNEMD is preferable in many ways to the forward NEMD methods because -it imposes what is typically difficult to measure (a flux or stress) -and it is typically much easier to compute momentum gradients or -strains (the response). For similar reasons, RNEMD is also preferable -to slowly-converging equilibrium methods for measuring thermal -conductivity and shear viscosity (using Green-Kubo relations or the -Helfand moment approach of Viscardy {\it et +\begin{figure} +\includegraphics[width=\linewidth]{thermalDemo} +\caption{RNEMD methods impose an unphysical transfer of momentum or + kinetic energy between a ``hot'' slab and a ``cold'' slab in the + simulation box. The molecular system responds to this imposed flux + by generating a momentum or temperature gradient. The slope of the + gradient can then be used to compute transport properties (e.g. + shear viscosity and thermal conductivity).} +\label{thermalDemo} +\end{figure} + +RNEMD is preferable in many ways to the forward NEMD +methods\cite{ISI:A1988Q205300014,hess:209,Vasquez:2004fk,backer:154503,ISI:000266247600008} +because it imposes what is typically difficult to measure (a flux or +stress) and it is typically much easier to compute the response +(momentum gradients or strains). For similar reasons, RNEMD is also +preferable to slowly-converging equilibrium methods for measuring +thermal conductivity and shear viscosity (using Green-Kubo +relations\cite{daivis:541,mondello:9327} or the Helfand moment +approach of Viscardy {\it et al.}\cite{Viscardy:2007bh,Viscardy:2007lq}) because these rely on computing difficult to measure quantities. @@ -91,29 +122,30 @@ Recently, Tenney and Maginn\cite{ISI:000273472300004} typically samples from the same manifold of states in the microcanonical ensemble. -Recently, Tenney and Maginn\cite{ISI:000273472300004} have discovered -some problems with the original RNEMD swap technique. Notably, large +Recently, Tenney and Maginn\cite{Maginn:2010} have discovered some +problems with the original RNEMD swap technique. Notably, large momentum fluxes (equivalent to frequent momentum swaps between the -slabs) can result in ``notched'', ``peaked'' and generally non-thermal momentum -distributions in the two slabs, as well as non-linear thermal and -velocity distributions along the direction of the imposed flux ($z$). -Tenney and Maginn obtained reasonable limits on imposed flux and -self-adjusting metrics for retaining the usability of the method. +slabs) can result in ``notched'', ``peaked'' and generally non-thermal +momentum distributions in the two slabs, as well as non-linear thermal +and velocity distributions along the direction of the imposed flux +($z$). Tenney and Maginn obtained reasonable limits on imposed flux +and proposed self-adjusting metrics for retaining the usability of the +method. In this paper, we develop and test a method for non-isotropic velocity -scaling (NIVS-RNEMD) which retains the desirable features of RNEMD +scaling (NIVS) which retains the desirable features of RNEMD (conservation of linear momentum and total energy, compatibility with periodic boundary conditions) while establishing true thermal -distributions in each of the two slabs. In the next section, we -develop the method for determining the scaling constraints. We then -test the method on both single component, multi-component, and -non-isotropic mixtures and show that it is capable of providing +distributions in each of the two slabs. In the next section, we +present the method for determining the scaling constraints. We then +test the method on both liquids and solids as well as a non-isotropic +liquid-solid interface and show that it is capable of providing reasonable estimates of the thermal conductivity and shear viscosity -in these cases. +in all of these cases. \section{Methodology} -We retain the basic idea of Muller-Plathe's RNEMD method; the periodic -system is partitioned into a series of thin slabs along a particular +We retain the basic idea of M\"{u}ller-Plathe's RNEMD method; the +periodic system is partitioned into a series of thin slabs along one axis ($z$). One of the slabs at the end of the periodic box is designated the ``hot'' slab, while the slab in the center of the box is designated the ``cold'' slab. The artificial momentum flux will be @@ -121,7 +153,7 @@ moves. For molecules $\{i\}$ located within the cold hot slab. Rather than using momentum swaps, we use a series of velocity scaling -moves. For molecules $\{i\}$ located within the cold slab, +moves. For molecules $\{i\}$ located within the cold slab, \begin{equation} \vec{v}_i \leftarrow \left( \begin{array}{ccc} x & 0 & 0 \\ @@ -129,9 +161,10 @@ where ${x, y, z}$ are a set of 3 scaling variables for 0 & 0 & z \\ \end{array} \right) \cdot \vec{v}_i \end{equation} -where ${x, y, z}$ are a set of 3 scaling variables for each of the -three directions in the system. Likewise, the molecules $\{j\}$ -located in the hot slab will see a concomitant scaling of velocities, +where ${x, y, z}$ are a set of 3 velocity-scaling variables for each +of the three directions in the system. Likewise, the molecules +$\{j\}$ located in the hot slab will see a concomitant scaling of +velocities, \begin{equation} \vec{v}_j \leftarrow \left( \begin{array}{ccc} x^\prime & 0 & 0 \\ @@ -141,20 +174,20 @@ Conservation of linear momentum in each of the three d \end{equation} Conservation of linear momentum in each of the three directions -($\alpha = x,y,z$) ties the values of the hot and cold bin scaling +($\alpha = x,y,z$) ties the values of the hot and cold scaling parameters together: \begin{equation} P_h^\alpha + P_c^\alpha = \alpha^\prime P_h^\alpha + \alpha P_c^\alpha \end{equation} where \begin{eqnarray} -P_c^\alpha & = & \sum_{i = 1}^{N_c} m_i \left[\vec{v}_i\right]_\alpha \\ -P_h^\alpha & = & \sum_{j = 1}^{N_h} m_j \left[\vec{v}_j\right]_\alpha +P_c^\alpha & = & \sum_{i = 1}^{N_c} m_i v_{i\alpha} \\ +P_h^\alpha & = & \sum_{j = 1}^{N_h} m_j v_{j\alpha} \label{eq:momentumdef} \end{eqnarray} Therefore, for each of the three directions, the hot scaling parameters are a simple function of the cold scaling parameters and -the instantaneous linear momentum in each of the two slabs. +the instantaneous linear momenta in each of the two slabs. \begin{equation} \alpha^\prime = 1 + (1 - \alpha) p_\alpha \label{eq:hotcoldscaling} @@ -167,16 +200,18 @@ Conservation of total energy also places constraints o Conservation of total energy also places constraints on the scaling: \begin{equation} -\sum_{\alpha = x,y,z} K_h^\alpha + K_c^\alpha = \sum_{\alpha = x,y,z} -\left(\alpha^\prime\right)^2 K_h^\alpha + \alpha^2 K_c^\alpha -\end{equation} -where the kinetic energies, $K_h^\alpha$ and $K_c^\alpha$, are computed -for each of the three directions in a similar manner to the linear momenta -(Eq. \ref{eq:momentumdef}). Substituting in the expressions for the -hot scaling parameters ($\alpha^\prime$) from Eq. (\ref{eq:hotcoldscaling}), -we obtain the {\it constraint ellipsoid equation}: +\sum_{\alpha = x,y,z} \left\{ K_h^\alpha + K_c^\alpha \right\} = \sum_{\alpha = x,y,z} +\left\{ \left(\alpha^\prime\right)^2 K_h^\alpha + \alpha^2 K_c^\alpha \right\} +\end{equation} +where the translational kinetic energies, $K_h^\alpha$ and +$K_c^\alpha$, are computed for each of the three directions in a +similar manner to the linear momenta (Eq. \ref{eq:momentumdef}). +Substituting in the expressions for the hot scaling parameters +($\alpha^\prime$) from Eq. (\ref{eq:hotcoldscaling}), we obtain the +{\it constraint ellipsoid}: \begin{equation} -\sum_{\alpha = x,y,z} a_\alpha \alpha^2 + b_\alpha \alpha + c_\alpha = 0 +\sum_{\alpha = x,y,z} \left( a_\alpha \alpha^2 + b_\alpha \alpha + + c_\alpha \right) = 0 \label{eq:constraintEllipsoid} \end{equation} where the constants are obtained from the instantaneous values of the @@ -188,317 +223,668 @@ This ellipsoid equation defines the set of cold slab s c_\alpha & = & K_h^\alpha p_\alpha^2 + 2 K_h^\alpha p_\alpha - K_c^\alpha \label{eq:constraintEllipsoidConsts} \end{eqnarray} -This ellipsoid equation defines the set of cold slab scaling -parameters which can be applied while preserving both linear momentum -in all three directions as well as kinetic energy. +This ellipsoid (Eq. \ref{eq:constraintEllipsoid}) defines the set of +cold slab scaling parameters which, when applied, preserve the linear +momentum of the system in all three directions as well as total +kinetic energy. -The goal of using velocity scaling variables is to transfer linear -momentum or kinetic energy from the cold slab to the hot slab. If the -hot and cold slabs are separated along the z-axis, the energy flux is -given simply by the decrease in kinetic energy of the cold bin: +The goal of using these velocity scaling variables is to transfer +kinetic energy from the cold slab to the hot slab. If the hot and +cold slabs are separated along the z-axis, the energy flux is given +simply by the decrease in kinetic energy of the cold bin: \begin{equation} (1-x^2) K_c^x + (1-y^2) K_c^y + (1-z^2) K_c^z = J_z\Delta t \end{equation} The expression for the energy flux can be re-written as another ellipsoid centered on $(x,y,z) = 0$: \begin{equation} -x^2 K_c^x + y^2 K_c^y + z^2 K_c^z = K_c^x + K_c^y + K_c^z - J_z\Delta t +\sum_{\alpha = x,y,z} K_c^\alpha \alpha^2 = \sum_{\alpha = x,y,z} +K_c^\alpha -J_z \Delta t \label{eq:fluxEllipsoid} \end{equation} -The spatial extent of the {\it flux ellipsoid equation} is governed -both by a targetted value, $J_z$ as well as the instantaneous values of the -kinetic energy components in the cold bin. +The spatial extent of the {\it thermal flux ellipsoid} is governed +both by the target flux, $J_z$ as well as the instantaneous values of +the kinetic energy components in the cold bin. To satisfy an energetic flux as well as the conservation constraints, -it is sufficient to determine the points ${x,y,z}$ which lie on both -the constraint ellipsoid (Eq. \ref{eq:constraintEllipsoid}) and the -flux ellipsoid (Eq. \ref{eq:fluxEllipsoid}), i.e. the intersection of -the two ellipsoids in 3-dimensional space. +we must determine the points ${x,y,z}$ that lie on both the constraint +ellipsoid (Eq. \ref{eq:constraintEllipsoid}) and the flux ellipsoid +(Eq. \ref{eq:fluxEllipsoid}), i.e. the intersection of the two +ellipsoids in 3-dimensional space. -One may also define momentum flux (say along the x-direction) as: +\begin{figure} +\includegraphics[width=\linewidth]{ellipsoids} +\caption{Velocity scaling coefficients which maintain both constant + energy and constant linear momentum of the system lie on the surface + of the {\it constraint ellipsoid} while points which generate the + target momentum flux lie on the surface of the {\it flux ellipsoid}. + The velocity distributions in the cold bin are scaled by only those + points which lie on both ellipsoids.} +\label{ellipsoids} +\end{figure} + +Since ellipsoids can be expressed as polynomials up to second order in +each of the three coordinates, finding the the intersection points of +two ellipsoids is isomorphic to finding the roots a polynomial of +degree 16. There are a number of polynomial root-finding methods in +the literature,\cite{Hoffman:2001sf,384119} but numerically finding +the roots of high-degree polynomials is generally an ill-conditioned +problem.\cite{Hoffman:2001sf} One simplification is to maintain +velocity scalings that are {\it as isotropic as possible}. To do +this, we impose $x=y$, and treat both the constraint and flux +ellipsoids as 2-dimensional ellipses. In reduced dimensionality, the +intersecting-ellipse problem reduces to finding the roots of +polynomials of degree 4. + +Depending on the target flux and current velocity distributions, the +ellipsoids can have between 0 and 4 intersection points. If there are +no intersection points, it is not possible to satisfy the constraints +while performing a non-equilibrium scaling move, and no change is made +to the dynamics. + +With multiple intersection points, any of the scaling points will +conserve the linear momentum and kinetic energy of the system and will +generate the correct target flux. Although this method is inherently +non-isotropic, the goal is still to maintain the system as close to an +isotropic fluid as possible. With this in mind, we would like the +kinetic energies in the three different directions could become as +close as each other as possible after each scaling. Simultaneously, +one would also like each scaling as gentle as possible, i.e. ${x,y,z + \rightarrow 1}$, in order to avoid large perturbation to the system. +To do this, we pick the intersection point which maintains the three +scaling variables ${x, y, z}$ as well as the ratio of kinetic energies +${K_c^z/K_c^x, K_c^z/K_c^y}$ as close as possible to 1. + +After the valid scaling parameters are arrived at by solving geometric +intersection problems in $x, y, z$ space in order to obtain cold slab +scaling parameters, Eq. (\ref{eq:hotcoldscaling}) is used to +determine the conjugate hot slab scaling variables. + +\subsection{Introducing shear stress via velocity scaling} +It is also possible to use this method to magnify the random +fluctuations of the average momentum in each of the bins to induce a +momentum flux. Doing this repeatedly will create a shear stress on +the system which will respond with an easily-measured strain. The +momentum flux (say along the $x$-direction) may be defined as: \begin{equation} (1-x) P_c^x = j_z(p_x)\Delta t \label{eq:fluxPlane} \end{equation} -The above {\it flux equation} is essentially a plane which is -perpendicular to the x-axis, with its position governed both by a -targetted value, $j_z(p_x)$ as well as the instantaneous value of the -momentum along the x-direction. +This {\it momentum flux plane} is perpendicular to the $x$-axis, with +its position governed both by a target value, $j_z(p_x)$ as well as +the instantaneous value of the momentum along the $x$-direction. -Similarly, to satisfy a momentum flux as well as the conservation -constraints, it is sufficient to determine the points ${x,y,z}$ which -lie on both the constraint ellipsoid (Eq. \ref{eq:constraintEllipsoid}) -and the flux plane (Eq. \ref{eq:fluxPlane}), i.e. the intersection of -an ellipsoid and a plane in 3-dimensional space. +In order to satisfy a momentum flux as well as the conservation +constraints, we must determine the points ${x,y,z}$ which lie on both +the constraint ellipsoid (Eq. \ref{eq:constraintEllipsoid}) and the +flux plane (Eq. \ref{eq:fluxPlane}), i.e. the intersection of an +ellipsoid and a plane in 3-dimensional space. -To summarize, by solving respective equation sets, one can determine -possible sets of scaling variables for cold slab. And corresponding -sets of scaling variables for hot slab can be determine as well. - -The following problem will be choosing an optimal set of scaling -variables among the possible sets. Although this method is inherently -non-isotropic, the goal is still to maintain the system as isotropic -as possible. Under this consideration, one would like the kinetic -energies in different directions could become as close as each other -after each scaling. Simultaneously, one would also like each scaling -as gentle as possible, i.e. ${x,y,z \rightarrow 1}$, in order to avoid -large perturbation to the system. Therefore, one approach to obtain the -scaling variables would be constructing an criteria function, with -constraints as above equation sets, and solving the function's minimum -by method like Lagrange multipliers. +In the case of momentum flux transfer, we also impose another +constraint to set the kinetic energy transfer as zero. In other +words, we apply Eq. \ref{eq:fluxEllipsoid} and let ${J_z = 0}$. With +one variable fixed by Eq. \ref{eq:fluxPlane}, one now have a similar +set of quartic equations to the above kinetic energy transfer problem. -In order to save computation time, we have a different approach to a -relatively good set of scaling variables with much less calculation -than above. Here is the detail of our simplification of the problem. +\section{Computational Details} -In the case of kinetic energy transfer, we impose another constraint -${x = y}$, into the equation sets. Consequently, there are two -variables left. And now one only needs to solve a set of two {\it - ellipses equations}. This problem would be transformed into solving -one quartic equation for one of the two variables. There are known -generic methods that solve real roots of quartic equations. Then one -can determine the other variable and obtain sets of scaling -variables. Among these sets, one can apply the above criteria to -choose the best set, while much faster with only a few sets to choose. +We have implemented this methodology in our molecular dynamics code, +OpenMD,\cite{Meineke:2005gd,openmd} performing the NIVS scaling moves +with a variable frequency after the molecular dynamics (MD) steps. We +have tested the method in a variety of different systems, including: +homogeneous fluids (Lennard-Jones and SPC/E water), crystalline +solids, using both the embedded atom method +(EAM)~\cite{PhysRevB.33.7983} and quantum Sutton-Chen +(QSC)~\cite{PhysRevB.59.3527} models for Gold, and heterogeneous +interfaces (QSC gold - SPC/E water). The last of these systems would +have been difficult to study using previous RNEMD methods, but the +current method can easily provide estimates of the interfacial thermal +conductivity ($G$). -In the case of momentum flux transfer, we impose another constraint to -set the kinetic energy transfer as zero. In another word, we apply -Eq. \ref{eq:fluxEllipsoid} and let ${J_z = 0}$. After that, with one -variable fixed by Eq. \ref{eq:fluxPlane}, one now have a similar set -of equations on the above kinetic energy transfer problem. Therefore, -an approach similar to the above would be sufficient for this as well. +\subsection{Simulation Cells} -\section{Computational Details} -Our simulation consists of a series of systems. All of these -simulations were run with the OpenMD simulation software -package\cite{Meineke:2005gd} integrated with RNEMD methods. +In each of the systems studied, the dynamics was carried out in a +rectangular simulation cell using periodic boundary conditions in all +three dimensions. The cells were longer along the $z$ axis and the +space was divided into $N$ slabs along this axis (typically $N=20$). +The top slab ($n=1$) was designated the ``hot'' slab, while the +central slab ($n= N/2 + 1$) was designated the ``cold'' slab. In all +cases, simulations were first thermalized in canonical ensemble (NVT) +using a Nos\'{e}-Hoover thermostat.\cite{Hoover85} then equilibrated in +microcanonical ensemble (NVE) before introducing any non-equilibrium +method. -A Lennard-Jones fluid system was built and tested first. In order to -compare our method with swapping RNEMD, a series of simulations were -performed to calculate the shear viscosity and thermal conductivity of -argon. 2592 atoms were in a orthorhombic cell, which was ${10.06\sigma - \times 10.06\sigma \times 30.18\sigma}$ by size. The reduced density -${\rho^* = \rho\sigma^3}$ was thus 0.85, which enabled direct -comparison between our results and others. These simulations used -velocity Verlet algorithm with reduced timestep ${\tau^* = - 4.6\times10^{-4}}$. +\subsection{RNEMD with M\"{u}ller-Plathe swaps} -For shear viscosity calculation, the reduced temperature was ${T^* = - k_B T/\varepsilon = 0.72}$. Simulations were first equilibrated in canonical -ensemble (NVT), then equilibrated in microcanonical ensemble -(NVE). Establishing and stablizing momentum gradient were followed -also in NVE ensemble. For the swapping part, Muller-Plathe's algorithm was -adopted.\cite{ISI:000080382700030} The simulation box was under -periodic boundary condition, and devided into ${N = 20}$ slabs. In each swap, -the top slab ${(n = 1)}$ exchange the most negative $x$ momentum with the -most positive $x$ momentum in the center slab ${(n = N/2 + 1)}$. Referred -to Tenney {\it et al.}\cite{ISI:000273472300004}, a series of swapping -frequency were chosen. According to each result from swapping -RNEMD, scaling RNEMD simulations were run with the target momentum -flux set to produce a similar momentum flux and shear -rate. Furthermore, various scaling frequencies can be tested for one -single swapping rate. To compare the performance between swapping and -scaling method, temperatures of different dimensions in all the slabs -were observed. Most of the simulations include $10^5$ steps of -equilibration without imposing momentum flux, $10^5$ steps of -stablization with imposing momentum transfer, and $10^6$ steps of data -collection under RNEMD. For relatively high momentum flux simulations, -${5\times10^5}$ step data collection is sufficient. For some low momentum -flux simulations, ${2\times10^6}$ steps were necessary. +In order to compare our new methodology with the original +M\"{u}ller-Plathe swapping algorithm,\cite{ISI:000080382700030} we +first performed simulations using the original technique. At fixed +intervals, kinetic energy or momentum exchange moves were performed +between the hot and the cold slabs. The interval between exchange +moves governs the effective momentum flux ($j_z(p_x)$) or energy flux +($J_z$) between the two slabs so to vary this quantity, we performed +simulations with a variety of delay intervals between the swapping moves. -After each simulation, the shear viscosity was calculated in reduced -unit. The momentum flux was calculated with total unphysical -transferred momentum ${P_x}$ and data collection time $t$: +For thermal conductivity measurements, the particle with smallest +speed in the hot slab and the one with largest speed in the cold slab +had their entire momentum vectors swapped. In the test cases run +here, all particles had the same chemical identity and mass, so this +move preserves both total linear momentum and total energy. It is +also possible to exchange energy by assuming an elastic collision +between the two particles which are exchanging energy. + +For shear stress simulations, the particle with the most negative +$p_x$ in the hot slab and the one with the most positive $p_x$ in the +cold slab exchanged only this component of their momentum vectors. + +\subsection{RNEMD with NIVS scaling} + +For each simulation utilizing the swapping method, a corresponding +NIVS-RNEMD simulation was carried out using a target momentum flux set +to produce the same flux experienced in the swapping simulation. + +To test the temperature homogeneity, directional momentum and +temperature distributions were accumulated for molecules in each of +the slabs. Transport coefficients were computed using the temperature +(and momentum) gradients across the $z$-axis as well as the total +momentum flux the system experienced during data collection portion of +the simulation. + +\subsection{Shear viscosities} + +The momentum flux was calculated using the total non-physical momentum +transferred (${P_x}$) and the data collection time ($t$): \begin{equation} j_z(p_x) = \frac{P_x}{2 t L_x L_y} \end{equation} -And the velocity gradient ${\langle \partial v_x /\partial z \rangle}$ -can be obtained by a linear regression of the velocity profile. From -the shear viscosity $\eta$ calculated with the above parameters, one -can further convert it into reduced unit ${\eta^* = \eta \sigma^2 - (\varepsilon m)^{-1/2}}$. +where $L_x$ and $L_y$ denote the $x$ and $y$ lengths of the simulation +box. The factor of two in the denominator is present because physical +momentum transfer between the slabs occurs in two directions ($+z$ and +$-z$). The velocity gradient ${\langle \partial v_x /\partial z + \rangle}$ was obtained using linear regression of the mean $x$ +component of the velocity, $\langle v_x \rangle$, in each of the bins. +For Lennard-Jones simulations, shear viscosities are reported in +reduced units (${\eta^* = \eta \sigma^2 (\varepsilon m)^{-1/2}}$). -For thermal conductivity calculation, simulations were first run under -reduced temperature ${T^* = 0.72}$ in NVE ensemble. Muller-Plathe's -algorithm was adopted in the swapping method. Under identical -simulation box parameters, in each swap, the top slab exchange the -molecule with least kinetic energy with the molecule in the center -slab with most kinetic energy, unless this ``coldest'' molecule in the -``hot'' slab is still ``hotter'' than the ``hottest'' molecule in the ``cold'' -slab. According to swapping RNEMD results, target energy flux for -scaling RNEMD simulations can be set. Also, various scaling -frequencies can be tested for one target energy flux. To compare the -performance between swapping and scaling method, distributions of -velocity and speed in different slabs were observed. +\subsection{Thermal Conductivities} -For each swapping rate, thermal conductivity was calculated in reduced -unit. The energy flux was calculated similarly to the momentum flux, -with total unphysical transferred energy ${E_{total}}$ and data collection -time $t$: +The energy flux was calculated in a similar manner to the momentum +flux, using the total non-physical energy transferred (${E_{total}}$) +and the data collection time $t$: \begin{equation} J_z = \frac{E_{total}}{2 t L_x L_y} \end{equation} -And the temperature gradient ${\langle\partial T/\partial z\rangle}$ -can be obtained by a linear regression of the temperature -profile. From the thermal conductivity $\lambda$ calculated, one can -further convert it into reduced unit ${\lambda^*=\lambda \sigma^2 - m^{1/2} k_B^{-1}\varepsilon^{-1/2}}$. +The temperature gradient ${\langle\partial T/\partial z\rangle}$ was +obtained by a linear regression of the temperature profile. For +Lennard-Jones simulations, thermal conductivities are reported in +reduced units (${\lambda^*=\lambda \sigma^2 m^{1/2} + k_B^{-1}\varepsilon^{-1/2}}$). -Another series of our simulation is to calculate the interfacial -thermal conductivity of a Au/H${_2}$O system. Respective calculations of -water (SPC/E) and gold (QSC) thermal conductivity were performed and -compared with current results to ensure the validity of -NIVS-RNEMD. After that, the mixture system was simulated. +\subsection{Interfacial Thermal Conductivities} -\section{Results And Discussion} -\subsection{Shear Viscosity} -Our calculations (Table \ref{shearRate}) shows that scale RNEMD method -produced comparable shear viscosity to swap RNEMD method. In Table -\ref{shearRate}, the names of the calculated samples are devided into -two parts. The first number refers to total slabs in one simulation -box. The second number refers to the swapping interval in swap method, or -in scale method the equilvalent swapping interval that the same -momentum flux would theoretically result in swap method. All the scale -method results were from simulations that had a scaling interval of 10 -time steps. The average molecular momentum gradients of these samples -are shown in Figure \ref{shearGrad}. +For interfaces with a relatively low interfacial conductance, the bulk +regions on either side of an interface rapidly come to a state in +which the two phases have relatively homogeneous (but distinct) +temperatures. The interfacial thermal conductivity $G$ can therefore +be approximated as: -\begin{table*} -\begin{minipage}{\linewidth} -\begin{center} +\begin{equation} +G = \frac{E_{total}}{2 t L_x L_y \left( \langle T_{gold}\rangle - + \langle T_{water}\rangle \right)} +\label{interfaceCalc} +\end{equation} +where ${E_{total}}$ is the imposed non-physical kinetic energy +transfer and ${\langle T_{gold}\rangle}$ and ${\langle + T_{water}\rangle}$ are the average observed temperature of gold and +water phases respectively. If the interfacial conductance is {\it + not} small, it is also be possible to compute the interfacial +thermal conductivity using this method utilizing the change in the +slope of the thermal gradient ($\partial^2 \langle T \rangle / \partial +z^2$) at the interface. -\caption{Calculation results for shear viscosity of Lennard-Jones - fluid at ${T^* = 0.72}$ and ${\rho^* = 0.85}$, with swap and scale - methods at various momentum exchange rates. Results in reduced - unit. Errors of calculations in parentheses. } +\section{Results} -\begin{tabular}{ccc} -\hline -Series & $\eta^*_{swap}$ & $\eta^*_{scale}$\\ -\hline -20-500 & 3.64(0.05) & 3.76(0.09)\\ -20-1000 & 3.52(0.16) & 3.66(0.06)\\ -20-2000 & 3.72(0.05) & 3.32(0.18)\\ -20-2500 & 3.42(0.06) & 3.43(0.08) -\end{tabular} -\label{shearRate} -\end{center} -\end{minipage} -\end{table*} +\subsection{Lennard-Jones Fluid} +2592 Lennard-Jones atoms were placed in an orthorhombic cell +${10.06\sigma \times 10.06\sigma \times 30.18\sigma}$ on a side. The +reduced density ${\rho^* = \rho\sigma^3}$ was thus 0.85, which enabled +direct comparison between our results and previous methods. These +simulations were carried out with a reduced timestep ${\tau^* = + 4.6\times10^{-4}}$. For the shear viscosity calculations, the mean +temperature was ${T^* = k_B T/\varepsilon = 0.72}$. For thermal +conductivity calculations, simulations were run under reduced +temperature ${\langle T^*\rangle = 0.72}$ in the microcanonical +ensemble. The simulations included $10^5$ steps of equilibration +without any momentum flux, $10^5$ steps of stablization with an +imposed momentum transfer to create a gradient, and $10^6$ steps of +data collection under RNEMD. -\begin{figure} -\includegraphics[width=\linewidth]{shearGrad} -\caption{Average momentum gradients of shear viscosity simulations} -\label{shearGrad} -\end{figure} +\subsubsection*{Thermal Conductivity} -\begin{figure} -\includegraphics[width=\linewidth]{shearTempScale} -\caption{Temperature profile for scaling RNEMD simulation.} -\label{shearTempScale} -\end{figure} -However, observations of temperatures along three dimensions show that -inhomogeneity occurs in scaling RNEMD simulations, particularly in the -two slabs which were scaled. Figure \ref{shearTempScale} indicate that with -relatively large imposed momentum flux, the temperature difference among $x$ -and the other two dimensions was significant. This would result from the -algorithm of scaling method. From Eq. \ref{eq:fluxPlane}, after -momentum gradient is set up, $P_c^x$ would be roughly stable -($<0$). Consequently, scaling factor $x$ would most probably larger -than 1. Therefore, the kinetic energy in $x$-dimension $K_c^x$ would -keep increase after most scaling steps. And if there is not enough time -for the kinetic energy to exchange among different dimensions and -different slabs, the system would finally build up temperature -(kinetic energy) difference among the three dimensions. Also, between -$y$ and $z$ dimensions in the scaled slabs, temperatures of $z$-axis -are closer to neighbor slabs. This is due to momentum transfer along -$z$ dimension between slabs. +Our thermal conductivity calculations show that the NIVS method agrees +well with the swapping method. Five different swap intervals were +tested (Table \ref{LJ}). Similar thermal gradients were observed with +similar thermal flux under the two different methods (Figure +\ref{thermalGrad}). Furthermore, the 1-d temperature profiles showed +no observable differences between the $x$, $y$ and $z$ axes (Figure +\ref{thermalGrad} c), so even though we are using a non-isotropic +scaling method, none of the three directions are experience +disproportionate heating due to the velocity scaling. -Although results between scaling and swapping methods are comparable, -the inherent temperature inhomogeneity even in relatively low imposed -exchange momentum flux simulations makes scaling RNEMD method less -attractive than swapping RNEMD in shear viscosity calculation. - -\subsection{Thermal Conductivity} - -Our thermal conductivity calculations also show that scaling method -agrees with swapping method. Table \ref{thermal} lists our simulation -results with similar manner we used in shear viscosity -calculation. All the data reported from scaling method were obtained -by simulations of 10-step exchange frequency, and the target exchange -kinetic energy were set to produce equivalent kinetic energy flux as -in swapping method. Figure \ref{thermalGradSwap} and -\ref{thermalGradScale} exhibit similar thermal gradients of respective -similar kinetic energy flux. - \begin{table*} -\begin{minipage}{\linewidth} -\begin{center} + \begin{minipage}{\linewidth} + \begin{center} -\caption{Calculation results for thermal conductivity of Lennard-Jones - fluid at ${\langle T^* \rangle = 0.72}$ and ${\rho^* = 0.85}$, with - swap and scale methods at various kinetic energy exchange rates. Results - in reduced unit. Errors of calculations in parentheses.} - -\begin{tabular}{ccc} -\hline -Name & $\lambda^*_{swap}$ & $\lambda^*_{scale}$\\ -\hline -20-250 & 7.03(0.34) & 7.30(0.10)\\ -20-500 & 7.03(0.14) & 6.95(0.09)\\ -20-1000 & 6.91(0.42) & 7.19(0.07)\\ -20-2000 & 7.52(0.15) & 7.19(0.28) -\end{tabular} -\label{thermal} -\end{center} -\end{minipage} + \caption{Thermal conductivity ($\lambda^*$) and shear viscosity + ($\eta^*$) (in reduced units) of a Lennard-Jones fluid at + ${\langle T^* \rangle = 0.72}$ and ${\rho^* = 0.85}$ computed + at various momentum fluxes. The original swapping method and + the velocity scaling method give similar results. + Uncertainties are indicated in parentheses.} + + \begin{tabular}{|cc|cc|cc|} + \hline + \multicolumn{2}{|c}{Momentum Exchange} & + \multicolumn{2}{|c}{Swapping RNEMD} & + \multicolumn{2}{|c|}{NIVS-RNEMD} \\ + \hline + \multirow{2}{*}{Swap Interval (fs)} & Equivalent $J_z^*$ or & + \multirow{2}{*}{$\lambda^*_{swap}$} & + \multirow{2}{*}{$\eta^*_{swap}$} & + \multirow{2}{*}{$\lambda^*_{scale}$} & + \multirow{2}{*}{$\eta^*_{scale}$} \\ + & $j_z^*(p_x)$ (reduced units) & & & & \\ + \hline + 250 & 0.16 & 7.03(0.34) & 3.57(0.06) & 7.30(0.10) & 3.54(0.04)\\ + 500 & 0.09 & 7.03(0.14) & 3.64(0.05) & 6.95(0.09) & 3.76(0.09)\\ + 1000 & 0.047 & 6.91(0.42) & 3.52(0.16) & 7.19(0.07) & 3.66(0.06)\\ + 2000 & 0.024 & 7.52(0.15) & 3.72(0.05) & 7.19(0.28) & 3.32(0.18)\\ + 2500 & 0.019 & 7.41(0.29) & 3.42(0.06) & 7.98(0.33) & 3.43(0.08)\\ + \hline + \end{tabular} + \label{LJ} + \end{center} + \end{minipage} \end{table*} \begin{figure} -\includegraphics[width=\linewidth]{thermalGradSwap} -\caption{Temperature gradients of simulations using swap method.} -\label{thermalGradSwap} + \includegraphics[width=\linewidth]{thermalGrad} + \caption{The NIVS-RNEMD method creates similar temperature gradients + compared with the swapping method under a variety of imposed + kinetic energy flux values. Furthermore, the implementation of + Non-Isotropic Velocity Scaling does not cause temperature + anisotropy to develop in thermal conductivity calculations.} + \label{thermalGrad} \end{figure} +\subsubsection*{Velocity Distributions} + +During these simulations, velocities were recorded every 1000 steps +and were used to produce distributions of both velocity and speed in +each of the slabs. From these distributions, we observed that under +relatively high non-physical kinetic energy flux, the speed of +molecules in slabs where swapping occured could deviate from the +Maxwell-Boltzmann distribution. This behavior was also noted by Tenney +and Maginn.\cite{Maginn:2010} Figure \ref{thermalHist} shows how these +distributions deviate from an ideal distribution. In the ``hot'' slab, +the probability density is notched at low speeds and has a substantial +shoulder at higher speeds relative to the ideal MB distribution. In +the cold slab, the opposite notching and shouldering occurs. This +phenomenon is more obvious at higher swapping rates. + +The peak of the velocity distribution is substantially flattened in +the hot slab, and is overly sharp (with truncated wings) in the cold +slab. This problem is rooted in the mechanism of the swapping method. +Continually depleting low (high) speed particles in the high (low) +temperature slab is not complemented by diffusions of low (high) speed +particles from neighboring slabs, unless the swapping rate is +sufficiently small. Simutaneously, surplus low speed particles in the +low temperature slab do not have sufficient time to diffuse to +neighboring slabs. Since the thermal exchange rate must reach a +minimum level to produce an observable thermal gradient, the +swapping-method RNEMD has a relatively narrow choice of exchange times +that can be utilized. + +For comparison, NIVS-RNEMD produces a speed distribution closer to the +Maxwell-Boltzmann curve (Figure \ref{thermalHist}). The reason for +this is simple; upon velocity scaling, a Gaussian distribution remains +Gaussian. Although a single scaling move is non-isotropic in three +dimensions, our criteria for choosing a set of scaling coefficients +helps maintain the distributions as close to isotropic as possible. +Consequently, NIVS-RNEMD is able to impose an unphysical thermal flux +as the previous RNEMD methods but without large perturbations to the +velocity distributions in the two slabs. + \begin{figure} -\includegraphics[width=\linewidth]{thermalGradScale} -\caption{Temperature gradients of simulations using scale method.} -\label{thermalGradScale} +\includegraphics[width=\linewidth]{thermalHist} +\caption{Velocity and speed distributions that develop under the + swapping and NIVS-RNEMD methods at high flux. The distributions for + the hot bins (upper panels) and cold bins (lower panels) were + obtained from Lennard-Jones simulations with $\langle T^* \rangle = + 4.5$ with a flux of $J_z^* \sim 5$ (equivalent to a swapping interval + of 10 time steps). This is a relatively large flux which shows the + non-thermal distributions that develop under the swapping method. + NIVS does a better job of producing near-thermal distributions in + the bins.} +\label{thermalHist} \end{figure} -During these simulations, molecule velocities were recorded in 1000 of -all the snapshots. These velocity data were used to produce histograms -of velocity and speed distribution in different slabs. From these -histograms, it is observed that with increasing unphysical kinetic -energy flux, speed and velocity distribution of molecules in slabs -where swapping occured could deviate from Maxwell-Boltzmann -distribution. Figure \ref{histSwap} indicates how these distributions -deviate from ideal condition. In high temperature slabs, probability -density in low speed is confidently smaller than ideal distribution; -in low temperature slabs, probability density in high speed is smaller -than ideal. This phenomenon is observable even in our relatively low -swpping rate simulations. And this deviation could also leads to -deviation of distribution of velocity in various dimensions. One -feature of these deviated distribution is that in high temperature -slab, the ideal Gaussian peak was changed into a relatively flat -plateau; while in low temperature slab, that peak appears sharper. +\subsubsection*{Shear Viscosity} +Our calculations (Table \ref{LJ}) show that velocity-scaling RNEMD +predicted comparable shear viscosities to swap RNEMD method. The +average molecular momentum gradients of these samples are shown in +Figure \ref{shear} (a) and (b). + \begin{figure} -\includegraphics[width=\linewidth]{histSwap} -\caption{Speed distribution for thermal conductivity using swapping RNEMD.} -\label{histSwap} + \includegraphics[width=\linewidth]{shear} + \caption{Average momentum gradients in shear viscosity simulations, + using ``swapping'' method (top panel) and NIVS-RNEMD method + (middle panel). NIVS-RNEMD produces a thermal anisotropy artifact + in the hot and cold bins when used for shear viscosity. This + artifact does not appear in thermal conductivity calculations.} + \label{shear} \end{figure} -\subsection{Interfaciel Thermal Conductivity} +Observations of the three one-dimensional temperatures in each of the +slabs shows that NIVS-RNEMD does produce some thermal anisotropy, +particularly in the hot and cold slabs. Figure \ref{shear} (c) +indicates that with a relatively large imposed momentum flux, +$j_z(p_x)$, the $x$ direction approaches a different temperature from +the $y$ and $z$ directions in both the hot and cold bins. This is an +artifact of the scaling constraints. After the momentum gradient has +been established, $P_c^x < 0$. Consequently, the scaling factor $x$ +is nearly always greater than one in order to satisfy the constraints. +This will continually increase the kinetic energy in $x$-dimension, +$K_c^x$. If there is not enough time for the kinetic energy to +exchange among different directions and different slabs, the system +will exhibit the observed thermal anisotropy in the hot and cold bins. -\section{Acknowledgments} -Support for this project was provided by the National Science -Foundation under grant CHE-0848243. Computational time was provided by -the Center for Research Computing (CRC) at the University of Notre -Dame. \newpage +Although results between scaling and swapping methods are comparable, +the inherent temperature anisotropy does make NIVS-RNEMD method less +attractive than swapping RNEMD for shear viscosity calculations. We +note that this problem appears only when momentum flux is applied, and +does not appear in thermal flux calculations. -\bibliographystyle{jcp2} +\subsection{Bulk SPC/E water} + +We compared the thermal conductivity of SPC/E water using NIVS-RNEMD +to the work of Bedrov {\it et al.}\cite{Bedrov:2000}, which employed +the original swapping RNEMD method. Bedrov {\it et + al.}\cite{Bedrov:2000} argued that exchange of the molecule +center-of-mass velocities instead of single atom velocities in a +molecule conserves the total kinetic energy and linear momentum. This +principle is also adopted Fin our simulations. Scaling was applied to +the center-of-mass velocities of the rigid bodies of SPC/E model water +molecules. + +To construct the simulations, a simulation box consisting of 1000 +molecules were first equilibrated under ambient pressure and +temperature conditions using the isobaric-isothermal (NPT) +ensemble.\cite{melchionna93} A fixed volume was chosen to match the +average volume observed in the NPT simulations, and this was followed +by equilibration, first in the canonical (NVT) ensemble, followed by a +100~ps period under constant-NVE conditions without any momentum flux. +Another 100~ps was allowed to stabilize the system with an imposed +momentum transfer to create a gradient, and 1~ns was allotted for data +collection under RNEMD. + +In our simulations, the established temperature gradients were similar +to the previous work. Our simulation results at 318K are in good +agreement with those from Bedrov {\it et al.} (Table +\ref{spceThermal}). And both methods yield values in reasonable +agreement with experimental values. + +\begin{table*} + \begin{minipage}{\linewidth} + \begin{center} + + \caption{Thermal conductivity of SPC/E water under various + imposed thermal gradients. Uncertainties are indicated in + parentheses.} + + \begin{tabular}{|c|c|ccc|} + \hline + \multirow{2}{*}{$\langle T\rangle$(K)} & + \multirow{2}{*}{$\langle dT/dz\rangle$(K/\AA)} & + \multicolumn{3}{|c|}{$\lambda (\mathrm{W m}^{-1} + \mathrm{K}^{-1})$} \\ + & & This work & Previous simulations\cite{Bedrov:2000} & + Experiment\cite{WagnerKruse}\\ + \hline + \multirow{3}{*}{300} & 0.38 & 0.816(0.044) & & + \multirow{3}{*}{0.61}\\ + & 0.81 & 0.770(0.008) & & \\ + & 1.54 & 0.813(0.007) & & \\ + \hline + \multirow{2}{*}{318} & 0.75 & 0.750(0.032) & 0.784 & + \multirow{2}{*}{0.64}\\ + & 1.59 & 0.778(0.019) & 0.730 & \\ + \hline + \end{tabular} + \label{spceThermal} + \end{center} + \end{minipage} +\end{table*} + +\subsection{Crystalline Gold} + +To see how the method performed in a solid, we calculated thermal +conductivities using two atomistic models for gold. Several different +potential models have been developed that reasonably describe +interactions in transition metals. In particular, the Embedded Atom +Model (EAM)~\cite{PhysRevB.33.7983} and Sutton-Chen (SC)~\cite{Chen90} +potential have been used to study a wide range of phenomena in both +bulk materials and +nanoparticles.\cite{ISI:000207079300006,Clancy:1992,Vardeman:2008fk,Vardeman-II:2001jn,ShibataT._ja026764r,Sankaranarayanan:2005lr,Chui:2003fk,Wang:2005qy,Medasani:2007uq} +Both potentials are based on a model of a metal which treats the +nuclei and core electrons as pseudo-atoms embedded in the electron +density due to the valence electrons on all of the other atoms in the +system. The SC potential has a simple form that closely resembles the +Lennard Jones potential, +\begin{equation} +\label{eq:SCP1} +U_{tot}=\sum _{i}\left[ \frac{1}{2}\sum _{j\neq i}D_{ij}V^{pair}_{ij}(r_{ij})-c_{i}D_{ii}\sqrt{\rho_{i}}\right] , +\end{equation} +where $V^{pair}_{ij}$ and $\rho_{i}$ are given by +\begin{equation} +\label{eq:SCP2} +V^{pair}_{ij}(r)=\left( \frac{\alpha_{ij}}{r_{ij}}\right)^{n_{ij}}, \rho_{i}=\sum_{j\neq i}\left( \frac{\alpha_{ij}}{r_{ij}}\right) ^{m_{ij}}. +\end{equation} +$V^{pair}_{ij}$ is a repulsive pairwise potential that accounts for +interactions between the pseudoatom cores. The $\sqrt{\rho_i}$ term in +Eq. (\ref{eq:SCP1}) is an attractive many-body potential that models +the interactions between the valence electrons and the cores of the +pseudo-atoms. $D_{ij}$, $D_{ii}$ set the appropriate overall energy +scale, $c_i$ scales the attractive portion of the potential relative +to the repulsive interaction and $\alpha_{ij}$ is a length parameter +that assures a dimensionless form for $\rho$. These parameters are +tuned to various experimental properties such as the density, cohesive +energy, and elastic moduli for FCC transition metals. The quantum +Sutton-Chen (QSC) formulation matches these properties while including +zero-point quantum corrections for different transition +metals.\cite{PhysRevB.59.3527} The EAM functional forms differ +slightly from SC but the overall method is very similar. + +In this work, we have utilized both the EAM and the QSC potentials to +test the behavior of scaling RNEMD. + +A face-centered-cubic (FCC) lattice was prepared containing 2880 Au +atoms (i.e. ${6\times 6\times 20}$ unit cells). Simulations were run +both with and without isobaric-isothermal (NPT)~\cite{melchionna93} +pre-equilibration at a target pressure of 1 atm. When equilibrated +under NPT conditions, our simulation box expanded by approximately 1\% +in volume. Following adjustment of the box volume, equilibrations in +both the canonical and microcanonical ensembles were carried out. With +the simulation cell divided evenly into 10 slabs, different thermal +gradients were established by applying a set of target thermal +transfer fluxes. + +The results for the thermal conductivity of gold are shown in Table +\ref{AuThermal}. In these calculations, the end and middle slabs were +excluded in thermal gradient linear regession. EAM predicts slightly +larger thermal conductivities than QSC. However, both values are +smaller than experimental value by a factor of more than 200. This +behavior has been observed previously by Richardson and Clancy, and +has been attributed to the lack of electronic contribution in these +force fields.\cite{Clancy:1992} It should be noted that the density of +the metal being simulated has an effect on thermal conductance. With +an expanded lattice, lower thermal conductance is expected (and +observed). We also observed a decrease in thermal conductance at +higher temperatures, a trend that agrees with experimental +measurements.\cite{AshcroftMermin} + +\begin{table*} + \begin{minipage}{\linewidth} + \begin{center} + + \caption{Calculated thermal conductivity of crystalline gold + using two related force fields. Calculations were done at both + experimental and equilibrated densities and at a range of + temperatures and thermal flux rates. Uncertainties are + indicated in parentheses. Richardson {\it et + al.}\cite{Clancy:1992} give an estimate of 1.74 $\mathrm{W + m}^{-1}\mathrm{K}^{-1}$ for EAM gold + at a density of 19.263 g / cm$^3$.} + + \begin{tabular}{|c|c|c|cc|} + \hline + Force Field & $\rho$ (g/cm$^3$) & ${\langle T\rangle}$ (K) & + $\langle dT/dz\rangle$ (K/\AA) & $\lambda (\mathrm{W m}^{-1} \mathrm{K}^{-1})$\\ + \hline + \multirow{7}{*}{QSC} & \multirow{3}{*}{19.188} & \multirow{3}{*}{300} & 1.44 & 1.10(0.06)\\ + & & & 2.86 & 1.08(0.05)\\ + & & & 5.14 & 1.15(0.07)\\\cline{2-5} + & \multirow{4}{*}{19.263} & \multirow{2}{*}{300} & 2.31 & 1.25(0.06)\\ + & & & 3.02 & 1.26(0.05)\\\cline{3-5} + & & \multirow{2}{*}{575} & 3.02 & 1.02(0.07)\\ + & & & 4.84 & 0.92(0.05)\\ + \hline + \multirow{8}{*}{EAM} & \multirow{3}{*}{19.045} & \multirow{3}{*}{300} & 1.24 & 1.24(0.16)\\ + & & & 2.06 & 1.37(0.04)\\ + & & & 2.55 & 1.41(0.07)\\\cline{2-5} + & \multirow{5}{*}{19.263} & \multirow{3}{*}{300} & 1.06 & 1.45(0.13)\\ + & & & 2.04 & 1.41(0.07)\\ + & & & 2.41 & 1.53(0.10)\\\cline{3-5} + & & \multirow{2}{*}{575} & 2.82 & 1.08(0.03)\\ + & & & 4.14 & 1.08(0.05)\\ + \hline + \end{tabular} + \label{AuThermal} + \end{center} + \end{minipage} +\end{table*} + +\subsection{Thermal Conductance at the Au/H$_2$O interface} +The most attractive aspect of the scaling approach for RNEMD is the +ability to use the method in non-homogeneous systems, where molecules +of different identities are segregated in different slabs. To test +this application, we simulated a Gold (111) / water interface. To +construct the interface, a box containing a lattice of 1188 Au atoms +(with the 111 surface in the $+z$ and $-z$ directions) was allowed to +relax under ambient temperature and pressure. A separate (but +identically sized) box of SPC/E water was also equilibrated at ambient +conditions. The two boxes were combined by removing all water +molecules within 3 \AA radius of any gold atom. The final +configuration contained 1862 SPC/E water molecules. + +The Spohr potential was adopted in depicting the interaction between +metal atoms and water molecules.\cite{ISI:000167766600035} A similar +protocol of equilibration to our water simulations was followed. We +observed that the two phases developed large temperature differences +even under a relatively low thermal flux. + +The low interfacial conductance is probably due to an acoustic +impedance mismatch between the solid and the liquid +phase.\cite{Cahill:793,RevModPhys.61.605} Experiments on the thermal +conductivity of gold nanoparticles and nanorods in solvent and in +glass cages have predicted values for $G$ between 100 and 350 +(MW/m$^2$/K). The experiments typically have multiple gold surfaces +that have been protected by a capping agent (citrate or CTAB) or which +are in direct contact with various glassy solids. In these cases, the +acoustic impedance mismatch would be substantially reduced, leading to +much higher interfacial conductances. It is also possible, however, +that the lack of electronic effects that gives rise to the low thermal +conductivity of EAM gold is also causing a low reading for this +particular interface. + +Under this low thermal conductance, both gold and water phase have +sufficient time to eliminate temperature difference inside +respectively (Figure \ref{interface} b). With indistinguishable +temperature difference within respective phase, it is valid to assume +that the temperature difference between gold and water on surface +would be approximately the same as the difference between the gold and +water phase. This assumption enables convenient calculation of $G$ +using Eq. \ref{interfaceCalc} instead of measuring temperatures of +thin layer of water and gold close enough to surface, which would have +greater fluctuation and lower accuracy. Reported results (Table +\ref{interfaceRes}) are of two orders of magnitude smaller than our +calculations on homogeneous systems, and thus have larger relative +errors than our calculation results on homogeneous systems. + +\begin{figure} +\includegraphics[width=\linewidth]{interface} +\caption{Temperature profiles of the Gold / Water interface at four + different values for the thermal flux. Temperatures for slabs + either in the gold or in the water show no significant differences, + although there is a large discontinuity between the materials + because of the relatively low interfacial thermal conductivity.} +\label{interface} +\end{figure} + +\begin{table*} + \begin{minipage}{\linewidth} + \begin{center} + + \caption{Computed interfacial thermal conductivity ($G$) values + for the Au(111) / water interface at ${\langle T\rangle \sim}$ + 300K using a range of energy fluxes. Uncertainties are + indicated in parentheses. } + + \begin{tabular}{|cccc| } + \hline + $J_z$ (MW/m$^2$) & $\langle T_{gold} \rangle$ (K) & $\langle + T_{water} \rangle$ (K) & $G$ + (MW/m$^2$/K)\\ + \hline + 98.0 & 355.2 & 295.8 & 1.65(0.21) \\ + 78.8 & 343.8 & 298.0 & 1.72(0.32) \\ + 73.6 & 344.3 & 298.0 & 1.59(0.24) \\ + 49.2 & 330.1 & 300.4 & 1.65(0.35) \\ + \hline + \end{tabular} + \label{interfaceRes} + \end{center} + \end{minipage} +\end{table*} + + +\section{Conclusions} +NIVS-RNEMD simulation method is developed and tested on various +systems. Simulation results demonstrate its validity in thermal +conductivity calculations, from Lennard-Jones fluid to multi-atom +molecule like water and metal crystals. NIVS-RNEMD improves +non-Boltzmann-Maxwell distributions, which exist inb previous RNEMD +methods. Furthermore, it develops a valid means for unphysical thermal +transfer between different species of molecules, and thus extends its +applicability to interfacial systems. Our calculation of gold/water +interfacial thermal conductivity demonstrates this advantage over +previous RNEMD methods. NIVS-RNEMD has also limited application on +shear viscosity calculations, but could cause temperature difference +among different dimensions under high momentum flux. Modification is +necessary to extend the applicability of NIVS-RNEMD in shear viscosity +calculations. + +\section{Acknowledgments} +The authors would like to thank Craig Tenney and Ed Maginn for many +helpful discussions. Support for this project was provided by the +National Science Foundation under grant CHE-0848243. Computational +time was provided by the Center for Research Computing (CRC) at the +University of Notre Dame. +\newpage + \bibliography{nivsRnemd} + \end{doublespace} \end{document}