A Homogeneous Non-Equilibrium Molecular Dynamics Method for Calculating Thermal Conductivity with a Three-Body Potential Kranthi K. Mandadapu Department of Mechanical Engineering, University of California, Berkeley, CA 94720-1740, USA
Reese E. Jones
∗
Sandia National Laboratories, Livermore, CA 94551-0969, USA
Panayiotis Papadopoulos Department of Mechanical Engineering, University of California, Berkeley, CA 94720-1740, USA
August 14, 2009
Abstract In this work, Evans’ homogeneous non-equilibrium molecular dynamics method for estimating thermal conductivity is extended to systems employing three-body potentials. This extension is put on a firm theoretical basis and applied to a silicon lattice modeled by the Stillinger-Weber potential. Two new methods are suggested for estimating the thermal conductivity based on a range of values of the fictitious force. Also, kinetic theory is used to estimate the linear range of the fictitious force necessary to bias the heat flow, thereby potentially reducing the number of simulations needed to estimate thermal conductivity.
1
Introduction
Thermal conductivity estimates are crucial in understanding and simulating the process of heat transfer by conduction. In problems related to mechanics, the thermal conductivity tensor κ is often strongly dependent on the deformation gradient F and the temperature T . The task of designing experiments that measure κ(F, T ) is generally complex and expensive, owing to the size of the parametric space of F and T . Moreover, it is often difficult to identify from the experimental results the contribution of, say, a distribution of defects or impurities. Molecular dynamics (MD) simulations provide an efficient supplement to such experiments. There are two well-known methods for calculating the thermal conductivity through MD simulations, namely the Green-Kubo (GK) method and the direct method. In the GK method, fluctuations of the heat flux from equilibrium MD simulations are employed to calculate the thermal conductivity via the fluctuation-dissipation theorem [1, 2, 3, 4]. Typically, the GK method requires very long simulation times to yield values of thermal conductivity that are converged to within statistical uncertainties. On the other hand, the direct method is a non-equilibrium molecular dynamics (NEMD) method, which attempts ∗
Corresponding author:
[email protected]
1
to mimic an experiment by imposing a temperature gradient across the system [4, 5, 6]. This method also has difficulties with finite-size effects and unrealistically large temperature gradients. A review of the application of the GK and direct methods to silicon is given in [4]. In this paper, the homogeneous non-equilibrium molecular dynamics (HNEMD) method proposed by Evans [7, 8, 3] for calculating the thermal transport coefficients is extended to three-body interatomic potentials. Evans’ method generates smaller statistical errors than the GK and requires minimum system sizes of the order of the GK method (which are considerably smaller than those required for the direct method). The HNEMD method was initially developed for thermal conductivity of systems with pair-wise interactions and attempts have been made to apply it to systems with higher-order interactions among the particles, [9, 10]. However, in both [9] and [10], it is unclear how the thermal conductivity is deduced when it varies nonlinearly within the considered range of the external field. Here, a theory similar to Evans’ is developed and applied to a silicon system governed by the three-body Stillinger-Weber (SW) potential [11]. The form of the heat flux vector, including the necessary three body terms, has already been developed for the related GK method [12, 4]. To this end, it is shown that the equations of motion require additional terms arising from the three-body interactions in order to be consistent with the theory. Moreover, it is demonstrated that these additional terms significantly affect the resulting conductivity estimates. The HNEMD method of Evans defines a mechanical analogue to the thermal transport process and uses the linear response theory to calculate the transport coefficients. This method is synthetic, in the sense that a fictitious force field is used to mimic the effect of a thermal gradient, thereby reducing the thermal transport problem to a mechanical problem. If this force field is sufficiently small, the response of a phase variable corresponding to the resultant non-equilibrium system involves only equilibrium time correlation functions. Finally, the thermal conductivity at a desired temperature and deformation gradient can be computed from the response of the virial heat flux vector to the fictitious force field. Here, two methods are proposed for estimating the thermal conductivity based on a range of values of the fictitious force, as alternatives to the extrapolation to zero fictitious force originally suggested by Evans [7]. Further, it is shown that kinetic theory can be used to estimate the linear range of the fictitious force needed for the linear response assumption. Both of these developments add significant practical value to the HNEMD method.
2
Theory
In linear response theory [3], the equations of motion for a system of N identical atoms are taken to be of the form pi + Ci (Γ)Fe (t) , m p˙ i = Fi + Di (Γ)Fe (t) , r˙ i =
(1)
where ri and pi are the position and momentum vectors of atom i respectively, Γ = {(ri , pi ) , i = 1, 2, . . . , N } is the phase space, F i is the interatomic force on atom i, and Ci (Γ), Di (Γ) are tensor phase variables which describe the coupling of the system to the fictitious external force field Fe . It is assumed that in the absence of the external field F e , 2
the phase-space of the system follows the canonical distribution, that is fc (Γ) =
e−βH0 (Γ) , z(β)
(2)
1 , with kB being kB T the Boltzmann constant. The equipartition theorem [13] furnishes a kinetic temperature relation + * X pi · p i 3 , (3) N kB T = 2 2m where H0 is the total energy, z(β) is the partition function, β =
i
c
where m is the mass of each atom and h·i c denotes the canonical phase space average. For the classical system under consideration, the total internal energy H 0 (which coincides with the total energy of the unperturbed system in (2)) takes the form H0 =
X pi · p i i
2m
+ Φ(r) ,
(4)
where Φ(r) is the potential that describes particle interactions. In this case, the time derivative of the total internal energy with respect to the new field-dependent equations of motion (1) is X pi X pi · p˙ i ˙ − Fi · r˙ i = · Di Fe (t) − Fi · Ci Fe (t) H0 (Γ, t) = m m i i (5) X T T pi − Ci Fi · Fe (t) = J(Γ) · Fe (t) , = Di m i
∂Φ are the interatomic forces and J(Γ) is termed the dissipative flux. It ∂ri is remarked here that the equations of motion (1) need not be derivable from a general Hamiltonian. ˜ Let B(Γ(t)) = B(t) be a general phase variable. Then, the ensemble average of B evolving with (1) in the Schr¨odinger representation is given by Z ˜ hB(t)i = B(Γ)f (Γ, t) dΓ , (6) where Fi = −
where f (Γ, t) is the perturbed phase-space distribution function obtained by solving the Liouville equation ∂ ˙ ∂ ∂f (Γ, t) ˙ = − ·Γ+Γ· f (Γ, t) = −iLf (Γ, t) , (7) ∂t ∂Γ ∂Γ with iL being the Liouvillean and the initial condition given by f (Γ, 0) = f c (Γ). At this stage, it is assumed that the equations of motion (1) satisfy the condition of Adiabatic ∂ Incompressibility of phase space (AIΓ), i.e., · Γ˙ = 0. Considering the force field Fe ∂Γ
3
to be sufficiently small, the Liouville equation (7) can be linearized and solved for the perturbed distribution function Z t exp − iL0 (t − s) β J(Γ) · Fe (s) fc (Γ) ds , (8) f (Γ, t) = fc (Γ) + 0
where iL0 is the Liouvillean for the field-free equations of motion, i.e., (1) with F e = 0, see [3, Chapter 5] for a detailed derivation. Substituting (8) into (6), and assuming that the external field Fe is time-independent, the ensemble average of B is obtained as Z t Z i h ˜ ˜ B(Γ) exp − iL0 (t − s) J(Γ) · Fe fc (Γ) dΓ ds , (9) hB(t)i = hB(0)ic + β 0
˜ where hB(0)i B(Γ)fc (Γ) dΓ. The inner integral on the right-hand side of (9) can be c = written as Z Z h i i h J(Γ)·Fe fc (Γ) exp iL0 (t−s) B(Γ) dΓ , B(Γ) exp −iL0 (t−s) J(Γ)·Fe fc (Γ) dΓ = R
(10)
using the identity Z Z h i h i B(Γ) exp − iL0 t g(Γ, 0) dΓ = g(Γ, 0) exp iL0 t B(Γ) dΓ ,
(11)
where iL0 is the phase variable p-Liouvillean of the field-free equations of motion and g(Γ, 0) is chosen to be equal to J(Γ)·Fe fc (Γ) [3, Section 3.3]. In equation (10), exp iL0 (t−s) B(Γ) is the value of B at time t − s starting with an initial value Γ chosen from J(Γ) · F e fc (Γ). ˜ Introducing the notation exp iL0 (t − s) B(Γ) = B((t − s)0 ; Γ), where the subscript ‘0’ in (t − s)0 indicates that the propagation is by field-free equations of motion, (10) can be rewritten as Z Z ˜ J(Γ) · Fe fc (Γ) exp iL0 (t − s) B(Γ) dΓ = J(Γ) · Fe fc (Γ)B((t − s)0 ; Γ) dΓ . (12) It is clear that J(Γ) in (12) is the value of J at the starting time and hence can be ˜ Γ). Hence, (12) reduces to written as J(Γ) = J(0; Z Z ˜ ˜ Γ) · Fe fc (Γ)B((t − s)0 ; Γ) dΓ J(Γ) · Fe fc (Γ) exp iL0 (t − s) B(Γ) dΓ = J(0; Z ˜ ˜ Γ) Fe fc (Γ) dΓ = B((t − s)0 ; Γ) ⊗ J(0;
˜ ˜ = B((t − s)0 ) ⊗ J(0) Fe , c
(13)
˜ ˜ where ‘⊗’ denotes tensor product and B((t − s)0 ) ⊗ J(0) is the correlation function of c B and J with respect to the canonical distribution. Substituting (10) and (13) in (9), the ensemble average of B becomes Z t
˜ ˜ ˜ ˜ hB(t)i = hB(0)ic + β B((t − s)0 ) ⊗ J(0) F ds c e 0 (14) Z t
˜ ˜ ˜ B((s) = hB(0)i 0 ) ⊗ J(0) c ds Fe . c+ β 0
4
˜ ˜ It is noted that, since B((s) 0 ) ⊗ J(0) c in (14) is independent of F e , the ensemble ˜ average hB(t)i is linearly dependent on the external field F e once the non-equilibrium steady state is reached. To deduce the thermal conductivity from (14), one needs to follow two steps. First, the phase variable B(Γ) is chosen as the virial heat flux vector B(Γ) = J Q (Γ). The latter is defined in Section 3 for silicon using a three-body potential. Second, the equations of motion are written such that the dissipative flux J(Γ) in (5) is equivalent to the virial heat flux vector J(Γ) = JQ (Γ) and further satisfy the AIΓ condition. Using these two steps, (14) reduces to Z t
˜ ˜ ˜ ˜ JQ (t) = JQ (0) c + β JQ ((s)0 ) ⊗ JQ (0) c ds Fe , (15) 0
˜ Q (0)ic = 0, as there is no heat flux when the system is in canonical equilibrium. where hJ Hence, (15) can be rewritten as
Z ∞ ˜ Q (∞)
J 1 ˜ ˜ = JQ ((s)0 ) ⊗ JQ (0)ic ds Fe , (16) VT V kB T 2 0 where V is the volume of the system. Thus, in the linear non-equilibrium steady state, ˜ Q (∞)i hJ is proportional to the external field F e with the constant of proportionality being VT the classical Green-Kubo expression for the thermal conductivity tensor. In other words, the thermal conductivity tensor κ can be calculated by the HNEMD method by tracking ˜ Q (∞)i hJ the response of to Fe rather than calculating the autocorrelation function directly VT from an equilibrium simulation as in the GK method. Assuming the system to be ergodic, ˜ Q (∞)i is equivalent to the time average with respect to the fieldthe ensemble average hJ dependent dynamics, i.e., Z 1 t0 +Y ˜ ˜ hJQ (∞)i = lim JQ (t) dt, (17) Y →∞ Y t0
where t0 denotes the time to attain the steady state [14]. Observing the similarity between the GK formula and (16), Evans [3, Section 6.5] proposed an expression for the (zz-component of the) thermal conductivity by setting F e = (0, 0, Fez ) so that
Z ∞ ˜Qz (∞)
J 1 κ = . (18) J˜Qz ((s)0 )J˜Qz (0)ic ds = lim Fez →0 V T Fez V kB T 2 0
This method of obtaining the limit and alternatives for estimating κ will be discussed in Section 4. It is noted here that, as it stands, the theory is ergodically inconsistent, in the ˜ Q ((s)0 ) in the equilibrium heat flux autocorrelation sense that the dynamics used to obtain J tensor of (16) does not generate the ensemble of starting states following f c (Γ). This inconsistency is removed in Section 4 by application of a thermostat.
3
Thermal conductivity for a three-body potential
In this section, the HNEMD method is developed for estimating the thermal conductivity of materials, such as bulk silicon, which are modeled using a three-body potential to account 5
for non-pairwise interactions. Three-body potentials for systems consisting of identical particles are of the general form 1 X 1 X u3 (ri , rj , rk ), (19) Φ(r) = u2 (ri , rj ) + 2! 3! i,j
i,j,k
where u2 (ri , rj ) describes the pair-wise interaction between atoms i and j and u 3 (ri , rj , rk ) describes the interaction of a triplet of atoms i, j and k. The total internal energy of this system can be written as X pi · p i 1 X 1 X H0 = + u2 (ri , rj ) + u3 (ri , rj , rk ) . (20) 2m 2! 3! i
i,j
i,j,k
The energy of each individual atom i is taken to be of the form 1 X pi · p i 1 X + Ei = u2 (ri , rj ) + u3 (ri , rj , rk ) , 2m 2! 3! j
(21)
j,k
and it can easily be seen that the total internal energy H 0 is recovered by summing the individual atomic energies. In this case, the heat flux is given by X pi Fijk · pi Fij · pi 1 X 1 X d X Ei + ri Ei = rij + , (rij + rik ) JQ (Γ) = dt m 2! m 3! m i
i
i,j
i,j,k
(22) ∂ ∂ u2 (ri , rj ) is the pair force and Fijk = − u3 (ri , rj , rk ) = where rij = ri − rj , Fij = − ∂ri ∂ri Fikj is the three body force [15, 12, 4]. It is noted here that there exist alternative choices for defining the energy Ei , see Appendix A for the corresponding heat flux expressions. The field-dependent equations of motion (1) for the three-body potential are chosen such that Ci Fe = 0 and X 1 X ¯ e+1 Di Fe = (Ei − E)F Fij (rij · Fe ) − Fjk (rjk · Fe ) 2 2N j j,k (23) 1X 1 X + Fijk [(rij + rik ) · Fe ] − Fljk [(rlj + rlk ) · Fe ] , 6 6N j,k
where Fi = −
l,j,k
X X 1X ∂Φ ¯ = 1 = Fij + Fijk and E Ei is the average instantaneous ∂ri 2 N j
i
j,k
internal energy. This choice of Di Fe extends Evans’ in [3, Section 6.5] by the addition of the last two terms involving three-body interactions. It is not necessary for C i (Γ) ¯ e in (23) implies that a heat current is to be zero, see, e.g., [16]. The term (E i − E)F induced by driving the particles with energy E i greater (respectively, smaller) than the ¯ in the direction (respectively, against the direction) of the external field average energy E X 1 X 1 X Fe . The terms 2!1 Fij (rij · Fe ) = (Fij ⊗ rij )Fe and Fijk [(rij + rik ) · Fe ] = 2! 3! j j j,k 1 X [Fijk ⊗ (rij + rik )]Fe imply that the heat current is also driven based on the deviation 3! j,k
of the virial part of the pressure tensor of atom i from the average pressure. 6
The equations of motion (1) and (23) are tailored to match the dissipative flux J(Γ) from the time derivative of total energy in (5) with the actual heat flux vector J Q (Γ) in (22). Further, it can be shown that they satisfy the AIΓ condition, which allows the application of the linear response theory, X as described in Section 2. In addition, they are both momentum-preserving, i.e., p˙ i = 0, and compatible with periodic boundary i
conditions. Using the momentum-preservation property, it is shown in Appendix B that X pi 1 X pi 1 X pi ˙ rij Fij · + (rij + rik ) Fijk · Ei + · Fe = JQ (Γ) · Fe . H0 = m 2! m 3! m i,j
i
i,j,k
(24) Therefore, the linear response theory of Section 2 implies that the average of heat flux vector JQ under the application of the external field F e is given by (16). In the remainder of this paper, the SW potential is used to describe the interactions among Si atoms, since it is widely accepted that elastic constants, phonon dispersion curves and thermal expansion coefficients can be accurately modeled using this potential, see, e.g., [17, 18, 19, 20].
4
Results
A scalar heat conduction coefficient κ is estimated by setting F e = (0, 0, Fez ). In this case, (16) reduces to Z ∞ hJ˜Qz (∞)i 1 ˜Qz ((s)0 )J˜Qz (0)ic ds Fez = κFez . = h J (25) VT V kB T 2 0 The analysis in Section 2 corresponds to an unthermostatted MD system and the resulting linear response theory is referred to as adiabatic. In practice, it is important for the MD system to be thermostatted, in order to remove the heat due to the work done by fictitious external field Fe . This ensures that the system reaches a steady non-equilibrium state and ˜ Q (∞)i in (16) remains finite. Here, the Nose-Hoover (NH) thermostat the long time-limit hJ is employed first to generate a canonical system before the external field F e is applied and then to remove the excess heat once F e is in place. Moreover, the NH thermostat ensures that the dynamics used to obtain J˜Q ((s)0 ) in the equilibrium heat flux autocorrelation tensor of (25) generates a canonical ensemble of starting states, hence enforcing ergodic consistency. The linear response theory for thermostatted systems is discussed in [21] and [3, Chapter 5]. Here, the equations of motion (1) and (23) with the NH thermostat are integrated by the operator split method, as described in Appendix C. The following three methods are considered to deduce the thermal conductivity from estimates of hJ˜Qz (∞)i corresponding to a set of force fields {F ez,j , j = 1, . . . , n}: 1. Extrapolation: This is the approach used
in [7], where a least squares fit of the ˜ JQz,j (∞) vs. Fez,j is extrapolated to Fez = 0 to function κ + mFez to κj = V T Fez,j determine a value for the thermal conductivity κ. The variance in the intercept n n X σ2 X 1 2 2 σκ = Fez,j is used to estimate the error, with σ = n−2 (κj − κ − mFez,j ) ∆ j=1
j=1
7
being the fit error and ∆ = n
n X
Fez,j
j=1
2
−
n X j=1
2 Fez,j , see [22].
2. Gradient: In this method, the slope of a least squares fit of the function κF ez to hJ˜Qz,j (∞)i vs. Fez,j in (25) is identified as the thermal conductivity. A zero intercept, TV σ2 n i.e. lim JQz = 0, is assumed. The error in the slope is estimated via σ κ2 = Fez →0 ∆ where, again, σ 2 is the fit error and ∆ is the variance of F ez,j . 3. Mean: Here, the conductivity is taken to be the average over κ j , namely κ =
1 n
n X
κj ,
j=1
and the error is estimated by σκ2 =
1 1 n n−1
n X
(κj − κ).
j=1
4.1
Argon using the Lennard-Jones potential
As an initial test, the algorithm is applied to an argon system at its triple point (86.5 K, 5.719 ˚ A) using a Lennard-Jones (LJ) pair potential with 256 atoms in an FCC lattice and a time step-size of 4.0 f s. Figure 1 shows the variation of κ j with the external field ˚−1 < Fez < 0.1 A ˚−1 , in which Fez,j . It is observed that there exists a finite range 0.04 A κj is essentially independent of Fez . Outside this range, there is strong dependence of κj on Fez and the linear response theory is not applicable. The thermal conductivity obtained by the extrapolation method using the estimates in the range of validity of the linear response theory is found to be 0.126 ± 0.001 W/mK, where the experimental result is 0.126 W/mK. The results in Table 1 show that this method can also be applied to calculate the thermal conductivity at different temperatures and volumetric deformations yielding results comparable to those obtained using the Green-Kubo method, see, e.g., [23].
4.2
Silicon using the Stillinger-Weber potential
Applying the HNEMD method to bulk silicon (1000 K, 5.43 ˚ A), a sequence of simulations is performed using the Stillinger-Weber potential with 1728 atoms in a diamond cubic lattice for about 4 million time steps, each of size 0.55 f s. The system is chosen large enough to be free of the size effects noted for the GK method in [4]. It is concluded from Figure 2 ˚−1 < Fez < 2 × 10−4 A ˚−1 . Since that the linear response theory is valid for 0.3 × 10 −4 A the external field is applied only along the z-direction, it is expected for the diamond cubic system that the average heat fluxes h J˜Qx (∞)i and hJ˜Qy (∞)i in the x and y directions approach zero respectively, as shown in Figure 3. The value of thermal conductivity by linear extrapolation to Fez → 0 is 50 ± 4 W/mK, which is lower than but consistent with the Green-Kubo result 66 ± 16 W/mK and the direct method result 65 ± 16 W/mK [4]. The gradient method gives a conductivity estimate of 55 ± 2 W/mK, see Figure 4. On the ˚−1 < Fez < 1 × 10−4 A ˚−1 gives other hand, the mean of values in the range 0.3 × 10 −4 A a conductivity estimate in the range 53 ± 2 W/mK, see Figure 2. Hence, all the three approaches for evaluating κ from the HNEMD method yield consistent results. To quantify the significance of three-body terms on deducing thermal conductivity, the simulations were repeated with only two-body terms included in the equations of motion 8
(1) and (23), as well as in the heat flux vector J Q (Γ) of (22). A regime of Fez was found where the linear response theory is valid, as shown in Figure 5. The value of thermal conductivity from both the gradient and mean methods is found to be 26 ± 1 W/mK which is approximately half of the value obtained from simulations including three-body terms, see Figure 6. This clearly demonstrates the significance of the three-body terms in the equations of motion and the heat flux vector when estimating the thermal conductivity of silicon.
5
Discussion
The HNEMD method for silicon requires 1-2 million time steps to obtain the thermal conductivity from a heat flux average that is converged to within statistical uncertainties. This is small relative to the GK method, which requires on the order of 6 million time steps to reach a comparable confidence interval [4, 6]. This is because the GK method relies on the integration of heat flux autocorrelation functions with long decay times. As argued in [4], bypassing the long-time correlation noise in the GK method by using an exponential fit to the short time largely underestimates the thermal conductivity in the GK method. However, the alternative approach of direct quadrature requires longer simulation times to obtain an accurate rate of decay. An additional complication in the GK method stems from the effect of round-off, which is more prominent when evaluating equilibrium fluctuations and increases with the number of time steps. Since the HNEMD method estimates the average hJ˜Qz (∞)i for a range of Fez , it eschews errors related to autocorrelation calculations. Each of these estimates may be obtained via simple time averaging, assuming the system to be ergodic, thereby avoiding the ensemble averaging needed in GK methods using NVE dynamics. In addition, time averaging in the linear non-equilibrium steady-state is less prone to round-off errors, due to the nature of the averaging and the smaller number of time steps, when compared to the GK method. The HNEMD method has two potential shortcomings. First, it requires simulations for multiple values of Fez in order to identify the range of the linear response theory. This is alleviated by the fact that such simulations are independent of each other and can be performed in parallel. Second, the method is inefficient for extremely small values of Fez . In this range, the system is very close to the equilibrium state and the effect of the external force field on the system becomes comparable to the system at equilibrium. hJ˜ (∞)i as Therefore, hJ˜Qz (∞)i approaches zero making it difficult to estimate the ratio Qz T Fez Fez approaches zero, see Figure 7. Hence, it is important to determine a range of F ez that hJ˜
(∞)i
is large enough to obtain reasonable values of Qz and small enough for the system to T Fez be in the linear non-equilibrium range. The kinetic theory of gases can be employed to estimate the linear response range of a system to the external perturbation F ez . It is known for the case of gases that Fourier’s Law (JQ = κ∇T ) is valid when the relative variation of temperature is small within a length equal to the mean free path ϕ [24, 25], namely ϕ
|∇T | 1. T
(26)
Assuming that silicon at 1000K is a gas of phonons, (26) can be employed with ϕ being the phonon mean free path. In the linear range, it can be seen from (25) and Fourier’s 9
law that the effect of external field F ez is identical to that of a logarithmic temperature ∂ ln T 1 ∂T gradient, i.e., Fez = = . Therefore, appealing to (26), it is concluded that ∂z T ∂z 1 ˚−1 , if the system is in the linear range, then F ez . Hence, for silicon Fez 0.0015 A ϕ since the phonon mean free path of Si at 1000K is estimated from the kinetic theory to be 600 ˚ A, see [4]. In fact, it can be observed from Figure 2 that the linear range of F ez ˚−1 . It follows that the number of simulations is approximately 1/100 to 1/10 of 0.0015 A that must be performed for different F ez to obtain a reliable conductivity estimate with the HNEMD method can be significantly reduced. At a minimum, the proposed extrapolation method requires only two conductivity estimates, whereas the gradient method requires a single flux measurement in the linear range since zero flux at zero F ez is assumed at the outset. The estimation of conductivity by the mean method can be also potentially accomplished with a single measurement, but it is more sensitive than the gradient method to the accurate knowledge of the lower limit of the linear range.
Acknowledgements REJ would like to acknowledge the support of the Laboratory Directed Research and Development program at Sandia National Laboratories. Sandia is a multiprogram laboratory operated by Sandia Corporation, a Lockheed Martin Company, for the United States Department of Energy under contract DE-ACO4-94AL85000. The authors would like to thank Professor Denis J. Evans of the Australian National University for insightful correspondence concerning the HNEMD method and Dr. Philip C. Howell of Siemens for discussions concerning Appendix A.
References [1] M. S. Green. Markoff random processes and the statistical mechanics of timedependent phenomena. 2. Irreversible processes in fluids. J. Chem. Phys., 22:398–413, 1954. [2] R. Kubo. Statistical-mechanical theory of irreversible processes. 1. General theory and simple applications to magnetic and conduction problems. J. Phys. Soc. Jpn., 12:570–586, 1957. [3] D. J. Evans and G. P. Morris. Statistical Mechanics of Non-equilibrium Liquids. Academic Press, 1990. [4] P. K. Schelling, S. R. Phillpot, and P. Keblinski. Comparison of atomic-level simulation methods for computing thermal conductivity. Phys. Rev. B, 65(14):144306–12, 2002. [5] W. G. Hoover. Computational Statistical Mechanics. Elsevier, Amsterdam, 1991. [6] X. W. Zhou, S. Aubry, R. E. Jones, A. Greenstein, and P. K. Schelling. Towards more accurate molecular dynamics calculation of thermal conductivity: Case study of GaN bulk crystals. Phys. Rev. B, 79(11):115201, 2009.
10
[7] D. J. Evans. Homogeneous NEMD algorithm for thermal conductivity–application of non-canonical linear response theory. Phys. Let. A, 91(9):457 – 460, 1982. [8] D. P. Hansen and D. J. Evans. A generalized heat flow algorithm. Mol. Phys., 81(4):767 – 779, 1994. [9] S. Berber, Y. K. Kwon, and D. Tomanek. Unusually high thermal conductivity of carbon nanotubes. Phys. Rev. Lett, 84:4613–5616, 2000. [10] J. R. Lukes and H. Zhong. Thermal conductivity of individual single-wall carbon nanotubes. J. Heat Transfer, 129(6):705–716, 2007. [11] F. H. Stillinger and T. A. Weber. Computer-simulation of local order in condensed phases of silicon. Phys. Rev. B, 31:5262–5271, 1985. [12] S. G. Volz and G. Chen. Molecular-dynamics simulation of thermal conductivity of Silicon crystals. Phys. Rev. B, 61(4):2651–2656, 2000. [13] J.H. Weiner. Statistical Mechanics of Elasticity. John Wiley, New York, 1983. [14] A.I. Khinchin. Mathematical foundations of statistical mechanics. Dover, New York, 1949. [15] J. H. Irving and J. G. Kirkwood. The statistical mechanical theory of transport processes. IV. The equations of hydrodynamics. J. Chem. Phys., 18:817–829, 1950. [16] M. J. Gillan and M. Dixon. The calculation of thermal conductivities by perturbed molecular dynamics simulation. J. Phys. C, 16(5):869–878, 1983. [17] M. D. Kluge, J. R. Ray, and A. Rahman. Molecular dynamic calculation of elastic constants of silicon. J. Chem. Phys., 85:4028–4031, 1986. [18] J. Q. Broughton and X. P. Li. Phase diagram of silicon by molecular dynamics. Phys. Rev. B, 35:9120–9127, 1987. [19] S. J. Cook and P. Clancy. Comparison of semiempirical potential functions for silicon and germanium. Phys. Rev. B, 47:7686–7699, 1993. [20] J. S. Kallman, W. G. Hoover, C. G. Hoover, A. J. DeGroot, S. M. Lee, and F. Wooten. Molecular-dynamics of silicon indentation. Phys. Rev. B., 47:7705–7709, 1993. [21] D. J. Evans and B. L. Holian. The Nose-Hoover thermostat. J. Chem. Phys., 83:4069– 4074, 1985. [22] P. R. Bevington and D. K. Robinson. Data Reduction and Error Analysis for the physical sciences. McGraw-Hill, Boston, 3rd edition, 2003. [23] H. Kaburaki, J. Li, S. Yip, and H. Kimizuka. Dynamical thermal conductivity of argon crystal. J. Appl. Phys., 102(4):043514, 2007. [24] S. Chapman and T. G. Cowling. The Mathematical Theory of Non-uniform Gases. Cambridge University Press, New York, 1939.
11
[25] I. Prigogine. Introduction to Thermodynamics of Irreversible Processes. John Wiley and Sons, New York, 1967. [26] G. J. Martyna, M. E. Tuckerman, D. J. Tobias, and M. L. Klein. Explicit reversible integrators for extended systems dynamics. Mol. Phys., 87:1117–1157(41), 1996. [27] S. J. Plimpton. Fast parallel algorithms for short-range molecular dynamics. J. Comp. Phys., 117:1–19, 1995. See also http://lammps.sandia.gov/.
12
A
Partition of energy to atoms and the heat flux vector for the SW potential
The total potential energy from pair and triplet bonds is given in (19). For this appendix, a shorthand notation is adopted by letting φ ij = u2 (ri , rj ) and φijk = u3 (ri , rj , rk ). For the SW potential, the three-body term φ ijk is taken to be the sum of three triplet bonds, ij jk ki φijk = φjk i + φj + φk , where φi = h(rij , rik , θjik ), with θjik representing the angle between rij and rik . Here, φjk i is the potential with i as the vertex, and j, k as the legs. In this case, the partition of total energy to individual atoms may be described using a one-parameter family 1 X pi · p i 1 X + αφ (A.1) φij + (1 − 3α)φjk + Ei = ijk , i 2m 2 2 j
j,k
where α ∈ [0, 31 ] determines the weights of φjk i and φijk . For example, α = 0 implies that jk 1 φi is assigned entirely to atom i and α = 3 asserts that φjk i is equally divided among the atoms i, j and k, thus leading to (21). Taking into account (1) for F 3 = 0 and (A.1), the heat flux vector JQ is given by d X ri Ei dt i X d pi · p i 1 X X pi 1X jk φij + (1 − 3α)φi + αφijk Ei + + ri = m dt 2m 2 2 j i i j,k X pi · p˙ i X 1 X ∂φij pi ∂φij pj X pi · · Ei + ri + = + ri m m 2 ∂ri m ∂rj m i i i j jk jk ∂φi ∂φjk pj 1X pi ∂φi pk i (1 − 3α) · · + + + · 2 ∂ri m ∂rj m ∂rk m j,k ijk ∂φi pj ∂φijk pi ∂φijk pk i i · + · + · +α ∂ri m ∂rj m ∂rk m X X pi ∂φij pi pi · p˙ i 1X Ei + ri ri + r j · = + m m 2 ∂ri m i i i,j X ∂φik ∂φji ∂φjk 1 pi pi pi j i k + (1 − 3α) ri ( · ) + rj ( · ) + rk ( · ) 2 ∂ri m ∂ri m ∂ri m i,j,k ∂φijk pk i · . (A.2) + α ri + r j + r k ∂rk m
JQ =
Letting
ji
Fij = −
∂φ ∂φijk ∂φij , Fji , = − k , Fijk = − ki ∂ri ∂ri ∂ri
13
(A.3)
substituting (A.3) into (A.2) and using (1) with F e = 0, the heat flux vector can be written as X pi X Fij · pi 1 X Fijk · pi 1 X Fij · pi JQ = + ri − ri + r j Ei + ri m m 2 m 2 m i i,j i,j i,j,k pi pi pi 1X · + rj Fik + rk Fji (1 − 3α) ri Fjk − ji · ii · ki 2 m m m i,j,k pi + α ri + rj + rk Fijk · . (A.4) m The three-body force Fijk may also be written as ji ik Fijk = (1 − 3α)Fijk + 3αFijk = (1 − 3α) Fjk + F + F ji ii ki + 3αFijk ,
(A.5)
where the right-hand side is obtained by using (A.3) and taking the gradient of φ ijk = ij ik φjk i + φj + φk with respect to ri . Using (A.5) and recalling the definition of r ij , the heat flux vector of (A.4) can be simplified to X X pi 1X Fij · pi ki pi Ei + rij + (1 − 3α)rij Fji · JQ = m 2 m m i,j i i,j,k 1X pi + α rij + rik Fijk · . (A.6) 2 m i,j,k
It is observed in [4] that the choices α = 0 and α = 13 yield similar results for the thermal conductivity of silicon systems. In this paper, all results are obtained for α = 13 , which corresponds to (22).
B
Correspondence of virial and dissipative heat fluxes
In this appendix, it is proved that the equations of motion (1) and (23) are written in such a way that the dissipative heat flux vector J coincides with the virial heat flux vector J Q of (22). To this end, it is noted that the rate of change of the total internal energy H 0 is X pi 1 X ∂φij pi ∂φij pj ˙ · · H0 = · p˙ i + + m 2 ∂ri m ∂rj m i i,j ∂φijk pk 1 X ∂φijk pi ∂φijk pj + · · + + · 6 ∂ri m ∂rj m ∂rk m i,j,k X pi 1 X ∂φij pi 1 X ∂φijk pi · p˙ i + · · = + . (B.1) m 2 ∂ri m 2 ∂ri m i
i,j
i,j,k
14
Using the equations of motion (1) and (23), and the definitions (A.3), it follows from (B.1) that X X pi pi ¯ 1 X pi ˙ H0 = Ei · F e − E · Fe + · Fij rij · Fe m m 2 m i i,j i X pi 1 X 1 X pi − Fjk rjk · Fe + · Fijk (rij + rik ) · Fe m 2N 6 m i j,k i,j,k X pi 1 X − Fljk rlj + rlk · Fe . (B.2) m 6N i
l,j,k
Using standard properties of tensors, it is readily concluded that p i · Fij (rij · F) = pi ·[(Fij ⊗rij )F] = [(rij ⊗Fij )pi ]·F = rij (Fij ·pi )·F and, similarly, pi ·Fijk [(rij +rik )·F] = (rij + rik )(Fijk · pi ) · F. Since the equations of motion are momentum-preserving and the X pi initial momentum is zero, = 0 at all times. Using these deductions, (B.2) simplifies m i to X pi pi pi 1X 1X ˙ rij (Fij · ) + (rij + rik )(Fijk · ) · Fe = JQ · Fe . (B.3) H0 = Ei + m 2 m 6 m i,j
i,j,k
Is it noted that for (B.3) to be used, the momentum-preservation property should be inherited by the discrete counterparts of the equations of motion (1) and (23).
C
Numerical integration algorithm
The governing system of equations (1) and (23) with the addition of a Nose-Hoover (NH) thermostat can be decomposed into 1 r 0 0 0 0 m pi d i + −ζpi + Di (r, p, F(r))Fe + Fi (r) + 0 (C.1) pi = 0 dt 1 0 ζ (T (p)/T0 − 1) 0 0 0 τ2 where T (p) is the temperature (3), T 0 is the expected temperature, ζ is the NH control variable, Di Fe is the Evans bias force (23), and Fi is the interatomic force. To integrate this system of ordinary differential equations the operator split method [26] is employed. This method is based on the decomposition of the p-Liouvillean operator that propagates the initial state Γ(t) = exp(iLt)Γ(0) (C.2)
into L = Lζ +Lp1 +Lp2 +Lp3 +Lx which correspond to sequence of vectors on the right-hand side of (C.1), so that ∆t ∆t ∆t ∆t exp(iL∆ t) = exp iLζ exp iLp3 exp iLp2 exp iLp1 exp(iLx ∆ t) 2 2 2 2 ∆t ∆t ∆t ∆t exp iLp1 exp iLp2 exp iLp3 exp iLζ + h.o.t. 2 2 2 2 (C.3) 15
In applying the operator split methodology the action of each of the evolution operators Lζ , Lp1 . . . on {ri , pi , ζ} is integrated in turn, preferrably exactly; however, in the case of Evans force Di Fe , which depends on p only through the energy E i , it is possible but not feasible to integrate p˙ = Di (r, p, F(r))Fe exactly given the quadratic dependence of Di Fe on p. Therefore, given the relative expense of the operations, where the interatomic force evaluation is the most expensive followed by the Evans force, the integration of D i Fe is lumped with F, see (C.8) and (C.9) below. With this choice, the resulting update is essentially the same as NVT algorithm in [26, Equation 26] ∆ t 2 T (pn ) ζn+1/2 = ζn − τ −1 (C.4) 2 T0 ∆t ˜ n = exp ζ p pn (C.5) 2 n+1/2 ∆t ˜n + Fn (C.6) pn+1/2 = p 2m ∆t p (C.7) rn+1 = rn + m n+1/2 ˜ n+1 = F(rn+1 ) F (C.8) ˜ n+1 + D(rn+1 , F ˜ n+1 , pn+1/2 )Fe Fn+1 = F (C.9) ∆t ˜ n+1 + Fn+1 pn+1 = p 2m ∆t ˜ n+1 = exp p ζ pn+1/2 2 n+1/2 ∆ t 2 T (pn+1 ) τ −1 ζn+1 = ζn+1/2 − 2 T0
(C.10) (C.11) (C.12)
where (C.4) and (C.12) are the exact integration of L ζ over ∆2 t . Likewise, (C.5) and (C.11) are the exact integration of Lp3 over ∆2 t , (C.6) and (C.10) are the integration of L p2 +p1 over ∆2 t and, finally, (C.7) is the exact integration of L x over ∆ t. Here, the subscripts enumerate time-steps. It also should be noted that this operator split method is not second order accurate, unlike the similar schemes in [26], due to its inability to integrate the Evans force Di Fe exactly. However, the method does preserve momentum exactly, assuming the total initial momentum is zero. By summing either (C.5) or (C.11) over the atoms, it is clear that the exponential factor merely scales the original momentum, i.e. zero, and therefore does not change it. Likewise, (C.6) and (C.10) do not change the total momentum since both the interatomic and field-dependent forces sum to zero. The sum of the interatomic forces F(r) is zero since the system is periodic and the sum of the field-dependent forces D(r, F, p)Fe is zero since these forces are formulated as a deviation from an average. This algorithm was incorporated as part of the MD code LAMMPS [27].
16
Interatomic distance (˚ A) 5.376 5.376 5.376 5.376 5.719
Temperature (K) 10 20 30 50 86.5
HNEMD κ (W/mK) 1.97 0.86 0.60 0.37 0.13
GK κ (W/mK) 1.91 0.95 0.63 0.39 0.13
Table 1: Comparison of thermal conductivity (κ) obtained by the HNEMD and GK methods at different temperatures and volumetric deformations
17
List of Figures: 1. Variation of thermal conductivity κ(F ez ) with Fez for argon at its triple point. The intercept at Fez = 0 is the estimate of the true thermal conductivity κ. Note that the errors in the estimates for smaller F ez are small enough to make the error bars indistinguishable. 2. Variation of thermal conductivity κ(F ez ) with Fez for Si system at T = 1000K. The intercept at Fez = 0 is the estimate of the true thermal conductivity κ. Note that the errors in the estimates for smaller F ez are small enough to make the error bars indistinguishable. 3. The running average of components of the heat flux J Q in the x and z directions for ˚−1 . Note that the Si system in the non-equilibrium steady state with F ez = 8 × 10−5 A the y direction is omitted for clarity. 4. Heat flux Jz as a function of Fez used to obtain the thermal conductivity via the gradient method. It is assumed that the linear fit passes through the origin. 5. Comparison of thermal conductivity κ(F ez ) variation with Fez for silicon at T=1000K. The intercept at Fez = 0 is the estimate of the true thermal conductivity κ. The method using the full heat flux and fictitious force is compared with the same method using only the two body terms. 6. Comparison of heat flux Jz estimates for silicon at T=1000K using the gradient method. The method using the full heat flux and fictitious force is compared with the same method using only the two body terms. 7. The running average of components of the heat flux J Q for the silicon system in the ˚−1 . non-equilibrium steady state with F ez = 1 × 10−5 A
18
0.2 0.18 0.16 κ [W/mK]
0.14 0.12 0.1 0.08 0.06 0.04 0.02 0 0
0.5
1
1.5
2
2.5
3
3.5
4
9
Fez [10 /m]
Figure 1: Variation of thermal conductivity κ(F ez ) with Fez for argon at its triple point. The intercept at Fez = 0 is the estimate of the true thermal conductivity κ. Note that the errors in the estimates for smaller F ez are small enough to make the error bars indistinguishable.
19
600 500
κ [W/mK]
400 300 200 100 0 0
0.5
1
1.5
2
2.5
3
3.5
4
6
Fez [10 /m]
Figure 2: Variation of thermal conductivity κ(F ez ) with Fez for Si system at T = 1000K. The intercept at Fez = 0 is the estimate of the true thermal conductivity κ. Note that the errors in the estimates for smaller F ez are small enough to make the error bars indistinguishable.
20
200
〈Jz〉 〈Jx〉
JQ [W/mK]
150 100 50 0 -50 0
0.5
1
1.5
2
2.5
6
t [10 fs]
Figure 3: The running average of components of the heat flux J Q in the x and z directions ˚−1 . Note that for the Si system in the non-equilibrium steady state with F ez = 8 × 10−5 A the y direction is omitted for clarity.
21
1800
〈JQz〉/VT [W/m2K]
1600 1400 1200 1000 800 600 400 200 0 0
0.5
1
1.5 2 2.5 Fez [106/m]
3
3.5
4
Figure 4: Heat flux Jz as a function of Fez used to obtain the thermal conductivity via the gradient method. It is assumed that the linear fit passes through the origin.
22
600
3-body 2-body only
500
κ [W/mK]
400 300 200 100 0 0
0.5
1
1.5
2
2.5
3
3.5
4
6
Fez [10 /m]
Figure 5: Comparison of thermal conductivity κ(F ez ) variation with Fez for silicon at T=1000K. The intercept at Fez = 0 is the estimate of the true thermal conductivity κ. The method using the full heat flux and fictitious force is compared with the same method using only the two body terms.
23
1800
3-body 2-body only
〈JQz〉/VT [W/m2K]
1600 1400 1200 1000 800 600 400 200 0 0
0.5
1
1.5 2 2.5 6 Fez [10 /m]
3
3.5
4
Figure 6: Comparison of heat flux Jz estimates for silicon at T=1000K using the gradient method. The method using the full heat flux and fictitious force is compared with the same method using only the two body terms.
24
200
〈Jz〉 〈Jx〉
JQ [W/mK]
150 100 50 0 -50 0
0.5
1
1.5
2
2.5
t [106 fs]
Figure 7: The running average of components of the heat flux J Q for the silicon system in ˚−1 . the non-equilibrium steady state with F ez = 1 × 10−5 A
25
