--- trunk/nivsRnemd/nivsRnemd.tex 2010/03/10 02:00:06 3570 +++ trunk/nivsRnemd/nivsRnemd.tex 2010/07/14 15:20:31 3611 @@ -6,6 +6,7 @@ \usepackage{caption} %\usepackage{tabularx} \usepackage{graphicx} +\usepackage{multirow} %\usepackage{booktabs} %\usepackage{bibentry} %\usepackage{mathrsfs} @@ -38,7 +39,22 @@ Notre Dame, Indiana 46556} \begin{doublespace} \begin{abstract} - + We present a new method for introducing stable non-equilibrium + velocity and temperature distributions in molecular dynamics + simulations of heterogeneous systems. This method extends earlier + Reverse Non-Equilibrium Molecular Dynamics (RNEMD) methods which use + momentum exchange swapping moves that 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 and stable + thermal and momentum gradients can then be established. The scaling + method we have developed conserves the total 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 +65,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 +119,29 @@ 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 +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 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 +149,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 +157,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,7 +170,7 @@ 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 @@ -170,13 +199,15 @@ where the kinetic energies, $K_h^\alpha$ and $K_c^\alp \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}: +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,230 +219,359 @@ 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. \begin{figure} \includegraphics[width=\linewidth]{ellipsoids} -\caption{Scaling points 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 hot bin are scaled by only those +\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} -One may also define momentum flux (say along the x-direction) as: +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, [CITATIONS NEEDED] but numerically finding the roots +of high-degree polynomials is generally an ill-conditioned +problem.[CITATION NEEDED] One simplification is to maintain velocity +scalings that are {\it as isotropic as possible}. To do this, we +impose $x=y$, and to 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. +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. -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 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 +after each MD step. We have tested the method in a variety of +different systems, including homogeneous fluids (Lennard-Jones and +SPC/E water), crystalline solids ({\sc eam}~\cite{PhysRevB.33.7983} and +quantum Sutton-Chen ({\sc q-sc})~\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 using velocity scaling moves, we can +even obtain estimates of the interfacial thermal conductivities ($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 ``cold'' slab, while the +central slab ($n= N/2 + 1$) was designated the ``hot'' 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. -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$: +\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 a the same momentum or energy flux exhibited in the +swapping simulation. + +To test the temperature homogeneity (and to compute transport +coefficients), directional momentum and temperature distributions were +accumulated for molecules in each of the slabs. + +\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 occurs in two directions due to our periodic +boundary conditions. The velocity gradient ${\langle \partial v_x + /\partial z \rangle}$ was obtained using linear regression of the +velocity profiles in the bins. For Lennard-Jones simulations, shear +viscosities are reporte 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 similarly 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 materials with a relatively low interfacial conductance, and in +cases where the flux between the materials is small, 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. +In calculating the interfacial thermal conductivity $G$, this +assumption was made, and the conductance can be approximated as: +\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. + +\section{Results} + +\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 first run under reduced +temperature ${\langle T^*\rangle = 0.72}$ in the microcanonical +ensemble, but other temperatures ([XXX, YYY, and ZZZ]) were also +sampled. 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. + +\subsubsection*{Thermal Conductivity} + +Our thermal conductivity calculations show that the NIVS method agrees +well with the swapping method. Four different swap intervals were +tested (Table \ref{thermalLJRes}). With a fixed 10 fs [WHY NOT REDUCED +UNITS???] scaling interval, the target exchange kinetic energy +produced equivalent kinetic energy flux as in the swapping method. +Similar thermal gradients were observed with similar thermal flux +under the two different methods (Figure \ref{thermalGrad}). + \begin{table*} -\begin{minipage}{\linewidth} -\begin{center} + \begin{minipage}{\linewidth} + \begin{center} -\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. } + \caption{Thermal conductivity (in reduced units) of a + Lennard-Jones fluid at ${\langle T^* \rangle = 0.72}$ and + ${\rho^* = 0.85}$ for the swapping and scaling methods at + various kinetic energy exchange rates. Uncertainties are + indicated in parentheses.} + + \begin{tabular}{|cc|cc|} + \hline + \multicolumn{2}{|c|}{Swapping RNEMD} & + \multicolumn{2}{|c|}{NIVS-RNEMD} \\ + \hline + Swap Interval (fs) & $\lambda^*_{swap}$ & Equilvalent $J_z^*$ & $\lambda^*_{scale}$\\ + \hline + 250 & 7.03(0.34) & 0.16 & 7.30(0.10)\\ + 500 & 7.03(0.14) & 0.09 & 6.95(0.09)\\ + 1000 & 6.91(0.42) & 0.047 & 7.19(0.07)\\ + 2000 & 7.52(0.15) & 0.024 & 7.19(0.28)\\ + \hline + \end{tabular} + \label{thermalLJRes} + \end{center} + \end{minipage} +\end{table*} -\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)\\ -\hline -\end{tabular} -\label{shearRate} -\end{center} -\end{minipage} -\end{table*} +\begin{figure} +\includegraphics[width=\linewidth]{thermalGrad} +\caption{NIVS-RNEMD method creates similar temperature gradients + compared with the swapping method under a variety of imposed kinetic + energy flux values.} +\label{thermalGrad} +\end{figure} +During these simulations, velocities were recorded every 1000 steps +and was 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 spee 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. + +In the velocity distributions, the ideal Gaussian peak 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]{shearGrad} -\caption{Average momentum gradients of shear viscosity simulations} -\label{shearGrad} +\includegraphics[width=\linewidth]{thermalHist} +\caption{Speed distribution for thermal conductivity using a) + ``swapping'' and b) NIVS- RNEMD methods. Shown is from the + simulations with an exchange or equilvalent exchange interval of 250 + fs. In circled areas, distributions from ``swapping'' RNEMD + simulation have deviation from ideal Maxwell-Boltzmann distribution + (curves fit for each distribution).} +\label{thermalHist} \end{figure} + +\subsubsection*{Shear Viscosity} +Our calculations (Table \ref{shearRate}) show that velocity-scaling +RNEMD predicted comparable shear viscosities to swap RNEMD 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{shear} (a) and (b). + +\begin{table*} + \begin{minipage}{\linewidth} + \begin{center} + + \caption{Shear viscosities of Lennard-Jones fluid at ${T^* = + 0.72}$ and ${\rho^* = 0.85}$ using swapping and NIVS methods + at various momentum exchange rates. Uncertainties are + indicated in parentheses. } + + \begin{tabular}{ccccc} + Swapping method & & & NIVS-RNEMD & \\ + \hline + Swap Interval (fs) & $\eta^*_{swap}$ & & Equilvalent $j_p^*(v_x)$ & + $\eta^*_{scale}$\\ + \hline + 500 & 3.64(0.05) & & 0.09 & 3.76(0.09)\\ + 1000 & 3.52(0.16) & & 0.046 & 3.66(0.06)\\ + 2000 & 3.72(0.05) & & 0.024 & 3.32(0.18)\\ + 2500 & 3.42(0.06) & & 0.019 & 3.43(0.08)\\ + \hline + \end{tabular} + \label{shearRate} + \end{center} + \end{minipage} +\end{table*} + \begin{figure} -\includegraphics[width=\linewidth]{shearTempScale} -\caption{Temperature profile for scaling RNEMD simulation.} -\label{shearTempScale} + \includegraphics[width=\linewidth]{shear} + \caption{Average momentum gradients in shear viscosity simulations, + using (a) ``swapping'' method and (b) NIVS-RNEMD method + respectively. (c) Temperature difference among x and y, z dimensions + observed when using NIVS-RNEMD with equivalent exchange interval of + 500 fs.} + \label{shear} \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 +two slabs which were scaled. Figure \ref{shear} (c) 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 @@ -431,146 +591,285 @@ attractive than swapping RNEMD in shear viscosity calc 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 results -agree 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{thermalGrad} exhibits similar thermal -gradients of respective similar kinetic energy flux. +\subsection{Bulk SPC/E water} -\begin{table*} -\begin{minipage}{\linewidth} -\begin{center} +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 in our simulations. Scaling was applied to +the center-of-mass velocities of the rigid bodies of SPC/E model water +molecules. -\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.} +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 +[XXX ps] period under constant-NVE conditions without any momentum +flux. [YYY ps] was allowed to stabilize the system with an imposed +momentum transfer to create a gradient, and [ZZZ ps] was alotted for +data collection under RNEMD. -\begin{tabular}{ccc} -\hline -Series & $\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)\\ -\hline -\end{tabular} -\label{thermal} -\end{center} -\end{minipage} -\end{table*} +As shown in Figure \ref{spceGrad}, temperature gradients were +established similar to the previous work. However, the average +temperature of our system is 300K, while that in Bedrov {\it et al.} +is 318K, which would be attributed for part of the difference between +the final calculation results (Table \ref{spceThermal}). [WHY DIDN'T +WE DO 318 K?] Both methods yield values in reasonable agreement with +experiment [DONE AT WHAT TEMPERATURE?] \begin{figure} -\includegraphics[width=\linewidth]{thermalGrad} -\caption{Temperature gradients of thermal conductivity simulations} -\label{thermalGrad} + \includegraphics[width=\linewidth]{spceGrad} + \caption{Temperature gradients in SPC/E water thermal conductivity + simulations.} + \label{spceGrad} \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 -swapping 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. +\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|ccc|} + \hline + $\langle dT/dz\rangle$(K/\AA) & \multicolumn{3}{|c|}{$\lambda + (\mathrm{W m}^{-1} \mathrm{K}^{-1})$} \\ + & This work (300K) & Previous simulations (318K)\cite{Bedrov:2000} & + Experiment\cite{WagnerKruse}\\ + \hline + 0.38 & 0.816(0.044) & & 0.64\\ + 0.81 & 0.770(0.008) & 0.784 & \\ + 1.54 & 0.813(0.007) & 0.730 & \\ + \hline + \end{tabular} + \label{spceThermal} + \end{center} + \end{minipage} +\end{table*} -\begin{figure} -\includegraphics[width=\linewidth]{histSwap} -\caption{Speed distribution for thermal conductivity using swapping RNEMD.} -\label{histSwap} -\end{figure} +\subsection{Crystalline Gold} -\begin{figure} -\includegraphics[width=\linewidth]{histScale} -\caption{Speed distribution for thermal conductivity using scaling RNEMD.} -\label{histScale} -\end{figure} +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 ({\sc eam})~\cite{PhysRevB.33.7983} and Sutton-Chen ({\sc + 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 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 ({\sc q-sc}) formulation matches these properties while +including zero-point quantum corrections for different transition +metals.\cite{PhysRevB.59.3527} The {\sc eam} functional forms differ +slightly from {\sc sc} but the overall method is very similar. -\subsection{Interfaciel Thermal Conductivity} +In this work, we have utilized both the {\sc eam} and the {\sc q-sc} +potentials to test the behavior of scaling RNEMD. -\begin{figure} -\includegraphics[width=\linewidth]{spceGrad} -\caption{Temperature gradients for SPC/E water thermal conductivity.} -\label{spceGrad} -\end{figure} +A face-centered-cubic (FCC) lattice was prepared containing 2880 Au +atoms. [LxMxN 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 OR LINEAR DIMENSIONS ?]. 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. {\sc eam} predicts +slightly larger thermal conductivities than {\sc q-sc}. 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 effects in +these force fields.\cite{Clancy:1992} The non-equilibrium MD method +employed in their simulations produced comparable results to ours. It +should be noted that the density of the metal being simulated also +greatly affects the thermal conductivity. 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 [PAGE +NUMBERS?].\cite{AshcroftMermin} + \begin{table*} -\begin{minipage}{\linewidth} -\begin{center} + \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. [CLANCY COMPARISON? SWAPPING + COMPARISON?]} + + \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}{*}{\sc q-sc} & \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}{*}{\sc 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*} -\caption{Calculation results for thermal conductivity of SPC/E water - at ${\langle T\rangle}$ = 300K at various thermal exchange rates. Errors of - calculations in parentheses. } +\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 withing 3 \AA radius of any gold atom. The final +configuration contained 1862 SPC/E water molecules. -\begin{tabular}{cccc} -\hline -$\langle dT/dz\rangle$(K/\AA) & & $\lambda$(W/m/K) & \\ -& This work & Previous simulations$^a$ & Experiment$^b$\\ -\hline -0.3 & 0.82() & 0.784 & 0.64\\ -0.8 & 0.770() & 0.730\\ -1.5 & 0.813() & \\ -\hline -\end{tabular} -\label{spceThermal} -\end{center} -\end{minipage} -\end{table*} +After simulations of bulk water and crystal gold, a mixture system was +constructed, consisting of 1188 Au atoms and 1862 H$_2$O +molecules. Spohr potential was adopted in depicting the interaction +between metal atom and water molecule.\cite{ISI:000167766600035} A +similar protocol of equilibration was followed. Several thermal +gradients was built under different target thermal flux. It was found +out that compared to our previous simulation systems, the two phases +could have large temperature difference even under a relatively low +thermal flux. +After simulations of homogeneous water and gold systems using +NIVS-RNEMD method were proved valid, calculation of gold/water +interfacial thermal conductivity was followed. It is found out that +the low interfacial conductance is probably due to the hydrophobic +surface in our system. Figure \ref{interface} (a) demonstrates mass +density change along $z$-axis, which is perpendicular to the +gold/water interface. It is observed that water density significantly +decreases when approaching the surface. 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]{AuGrad} -\caption{Temperature gradients for crystal gold thermal conductivity.} -\label{AuGrad} +\includegraphics[width=\linewidth]{interface} +\caption{Simulation results for Gold/Water interfacial thermal + conductivity: (a) Significant water density decrease is observed on + crystalline gold surface, which indicates low surface contact and + leads to low thermal conductance. (b) Temperature profiles for a + series of simulations. Temperatures of different slabs in the same + phase show no significant differences.} +\label{interface} \end{figure} \begin{table*} \begin{minipage}{\linewidth} \begin{center} -\caption{Calculation results for thermal conductivity of crystal gold - at ${\langle T\rangle}$ = 300K at various thermal exchange rates. Errors of - calculations in parentheses. } +\caption{Calculation results for interfacial thermal conductivity + at ${\langle T\rangle \sim}$ 300K at various thermal exchange + rates. Errors of calculations in parentheses. } -\begin{tabular}{ccc} +\begin{tabular}{cccc} \hline -$\langle dT/dz\rangle$(K/\AA) & $\lambda$(W/m/K)\\ -& This work & Previous simulations$^a$ \\ +$J_z$ (MW/m$^2$) & $T_{gold}$ (K) & $T_{water}$ (K) & $G$ +(MW/m$^2$/K)\\ \hline -1.4 & 1.10() & \\ -2.8 & 1.08() & \\ -5.1 & 1.15() & \\ +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{AuThermal} +\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 in 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} 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 -\bibliographystyle{jcp2} +\bibliographystyle{aip} \bibliography{nivsRnemd} + \end{doublespace} \end{document}