Generalization of the Homogeneous Non-Equilibrium Molecular Dynamics Method for Calculating Thermal Conductivity to Multi-body Potentials Kranthi K. Mandadapu,1 Reese E. Jones,2 and Panayiotis Papadopoulos1 1

Department of Mechanical Engineering, University of California, Berkeley, CA 94720-1740, USA 2 Sandia National Laboratories, Livermore, CA 94551-0969, USA

This work provides a generalization of Evans’ homogeneous non-equilibrium method for estimating thermal conductivity to molecular systems that are described by general multi-body potentials. A perturbed form of the usual Nose-Hoover equations of motion is formally constructed and is shown to satisfy the requirements of Evans’ original method. These include adiabatic incompressibility of phase space, equivalence of the dissipative and heat fluxes and momentum preservation.

I.

INTRODUCTION

Estimating the thermal conductivity κ is essential for modeling the flow of heat using Fourier’s law. The thermal conductivity can be computed using Molecular Dynamics (MD) simulations, given sufficient information regarding the inherent molecular structure and interactions among the particles in the system [1–3]. Commonly, such interactions are modeled with empirical potentials describing the bonding energy of neighboring atoms. While many materials can be adequately described with pair potentials based simply on interatomic distances, semiconductor materials, such as Si and GaN, require additional terms that reflect the relative orientation between triplets of neighboring atoms [4–6]. Moreover, four-body terms measuring dihedral angles and torsion are also employed to describe polymers [7–9] and the surface and cluster phases of silicon [10]. In this paper, a Homogeneous Non-Equilibrium Molecular Dynamics (HNEMD) method is formulated for any system consisting of N identical particles described by potentials containing up to M -body terms. This method was originally proposed by Evans for systems modeled by pair potentials [11], applied to molecular liquids with distance and angle constraints as well as an intra-molecular four body energy in [12, 13], generalized to mixtures of arbitrarily complex non-identical molecules in [14], and was later extended to solid systems modeled by threebody potentials in [15]. The HNEMD method is synthetic, in the sense that a fictitious external force field is applied to the system in order to mimic the thermal transport process. Since the method requires only the time average of the heat flux, it is free of difficulties involving the calculation and integration of the heat flux autocorrelation tensor necessary for estimating κ using the Green-Kubo (GK) method. Also, it is free of strong size effects, like the GK method, and does not involve large temperature gradients both of which plague the direct method [2, 15]. For these reasons, the HNEMD method is a practical means for computing statistically accurate estimates of thermal conductivity. The most challenging task in the HNEMD method, which is addressed in this paper, is to formulate the field-dependent equations of motion so that the flux induced is commensurate with the microscopic heat flux vector [16, Chapter

6] The organization of the paper is as follows: The linear response theory of Nose-Hoover (NH) thermostatted systems is derived systematically in Section II. This is accomplished in a straightforward manner by appealing to a direct linearization of the Liouville equation, as in the case of adiabatic linear response theory [16, Chapter 5.1], instead of linearizing the propagator of the perturbed distribution function, as in [17]. The main result of the paper is obtained in Section III, where the equations of motion necessary for computing thermal conductivity are derived for an M -body potential consistent with the HNEMD method of Section II. II.

LINEAR RESPONSE THEORY WITH NH THERMOSTAT

The generic form of the NH-thermostatted equations of motion linearly perturbed by an external field Fe is pi + Ci (Γ)Fe (t) , r˙ i = m p˙ i = Fi + Di (Γ)Fe (t) − ζpi , (1) X  p · p 1 i i − 3N kB T . ζ˙ = Q m i Here, kB is the Boltzmann constant, T is the temperature, Γ = {(ri , pi ) , i = 1, 2, . . . , N } is the phase space with ri and pi being the position and momentum vectors of atom i respectively. Also, Fi is the interatomic force, Ci (Γ) and Di (Γ) are the tensor phase variables which describe the coupling of the system to the applied external field Fe , ζ is the thermodynamic friction coefficient associated with the NH thermostat, and Q is a parameter that affects the amplitude and period of thermal fluctuations. It is well-known that the NH thermostat preserves the canonical distribution in the absence of external field [17, 18]. When Fe = 0, the phase space distribution follows the extended canonical distribution [18]
2 1 e−β H0 (Γ)+ 2 Qζ fc (Γ, ζ) = , (2) z(β) where H0 is the total internal energy of the system and R −β H0 (Γ)+ 1 Qζ 2
2 dΓ dζ is the partition function z(β) = e

2 with β =

1 . The total internal energy H0 is given by kB T X pi · pi H0 = + Φ(r) , (3) 2m i

where Φ(r) is the potential that describes interactions among the particles such that the interatomic force is ∂ Φ. Assuming the external field given by Fi = − ∂x i Fe is applied at time t = 0, the perturbed distribution f (Γ, ζ, t) is obtained by solving the Liouville equation   ∂f (Γ, ζ, t) ∂ ∂ ∂ ˙ ∂ ˙ ˙ ˙ ζ f (Γ, ζ, t) = − ·Γ+Γ· +ζ + ∂t ∂Γ ∂Γ ∂ζ ∂ζ = −iLf (Γ, ζ, t) , (4)

where i∆L(t) and ∆f (Γ, ζ, t) are linear perturbations to the field-free Liouvillean iL0 and the extended canonical distribution fc (Γ, ζ), respectively, due to the presence of the external field Fe . Here, iL0 corresponds to the fieldfree equations of motion, i.e., equation (1) with Fe = 0, ∂ f = −iL0 f . Using the approximation while fc satisfies ∂t (5), the Liouville equation (4) can be linearized as ∂ ∆f (Γ, ζ, t) + iL0 ∆f (Γ, ζ, t) = −i∆L(t)fc (Γ, ζ) . (6) ∂t The solution ∆f (Γ, ζ, t) of the linearized Liouville equation (6) can be obtained formally as ∆f (Γ, ζ, t) = −

Z

t

0

with initial condition f (Γ, ζ, 0) = fc (Γ, ζ). Here, iL represents the Liouvillean. Assuming the external field Fe to be small enough, the Liouvillean iL and the perturbed distribution function f (Γ, ζ, t) may be approximated by iL = iL0 + i∆L(t) , f (Γ, ζ, t) = fc (Γ, ζ) + ∆f (Γ, ζ, t)

(5)


exp −iL0 (t − s) i∆L(s)fc (Γ, ζ) ds ,

(7) where i∆L(s)fc (Γ, ζ) = iL(s)fc (Γ, ζ) − iL0 fc (Γ, ζ) [16, Chapter 5.1]. Since fc (Γ, ζ) is the steady-state solution of ∂ f = −iL0 f , it follows that iL0 fc (Γ, ζ) = the equation ∂t 0 [18], hence i∆L(s)fc (Γ, ζ) = iL(s)fc (Γ, ζ). Using (4) and (2), it can be shown that for the field-dependent equations of motion (1)



 ∂ ˙ ∂ ∂ ∂ ˙ ˙ ˙ i∆L(s)fc (Γ, ζ) = ζ fc (Γ, ζ) ·Γ+Γ· +ζ + ∂Γ ∂Γ ∂ζ ∂ζ X    X ∂ ∂ T pi T = · Ci Fe + · Di Fe fc (Γ, ζ) + β − Di + Ci Fi · Fe (t)fc (Γ, ζ) . ∂qi ∂pi m i i

At this stage, the equations of motion are assumed to satisfy the condition of adiabatic incompressibility of phase space (AIΓ) [16, Chapter 5], i.e.,  X ∂   X ∂  · Ci Fe + · Di Fe = 0 . ∂qi ∂pi i i

(9)

(8)

f (Γ, ζ, t) is obtained as f (Γ, ζ, t) = fc (Γ, ζ) Z t
  exp −iL0 (t − s) J(Γ) · Fe (s) fc (Γ, ζ) ds . +β 0

(11)

Note that, in the absence of a thermostat, the AIΓ condition is identical to the incompressibility of the full phase ∂ · Γ˙ = 0. Also, the AIΓ condition is always space, i.e., ∂Γ satisfied if the equations of motion are derivable from a Hamiltonian. However, here the AIΓ is assumed to hold even if, as in the present case, the equations of motion (1) are not derivable from a general Hamiltonian. Using (9), equation (8) may be further reduced to   T pi T i∆L(s)fc (Γ, ζ) = β − Di + Ci Fi · Fe (t)fc (Γ, ζ) m = −βJ(Γ) · Fe (t)fc (Γ, ζ) , (10)

The ensemble average of the microscopic virial heat flux ˜ Q (t) derived from the field-free equavector JQ (Γ(t)) = J tions of motion using the Irving-Kirkwood procedure [19], but evolving with the field-dependent equations of motion (1), is given by Z ˜ Q (t)i = hJ JQ (Γ)f (Γ, ζ, t) dΓ dζ Z tZ ˜ JQ (Γ) = hJQ (0)ic + β 0  i h 
 exp − iL0 (t − s) J(Γ) · Fe fc (Γ, ζ) dΓ dζ ds ,

where J(Γ) is defined as the dissipative flux. Substituting (10) into (7), the perturbed distribution function

˜ Q (0)ic = JQ (Γ)fc (Γ, ζ) dΓ dζ. Assuming the where hJ external field Fe to be independent of time, the inner

(12)

R

3 integral on the right-hand side of (12) can be written as

Z

i h 
 JQ (Γ) exp − iL0 (t − s) J(Γ) · Fe fc (Γ, ζ) dΓ dζ

 ˜ Q ((t − s)0 ) ⊗ J(0) ˜ = J F , c e (13)

III.

The interaction potential of a system of N identical particles is of the general form Φ(r) =

1 X 1 X u2 (ri1 , ri2 ) + u3 (ri1 , ri2 , ri3 ) + . . . 2! i ,i 3! i ,i ,i 1

˜ Q ((t − s)0 ) ⊗ where ⊗ denotes the tensor product, J  ˜ J(0) is the correlation function of the heat flux vector c JQ and the dissipative flux J with respect to the extended ˜ = J(Γ(t)) [15, Section 2]. canonical distribution, and J(t) ˜ − s)0 ) in the The subscript ‘0’ in (13) indicates that J((t correlation function is obtained by solving the field-free equations of motion, see [16, Chapter 5] and [15, Section 2]. Using (13), equation (12) reduces to

  Z t

 ˜ Q ((s)0 ) ⊗ J(0) ˜ ˜ Q (t)i = β ds Fe , J hJ c

1

1

2

3

M

2

(16) where M ≤ N and uM (ri1 , ri2 , . . . , riM ) describes the Mbody interactions [20]. The total inter-atomic force Fi on atom i is given by 1 X 1 X ∂Φ = Fii2 + Fii2 i3 ∂ri 1! i 2! i ,i 2 2 3 X 1 + ... + Fii2 ...iM . (M − 1)! i ,...,i

Fi = −

(14)

The HNEMD method requires that the equations of motion (1) are constructed so that the dissipative flux is equivalent to the heat flux vector obtained by using the field-free equations of motion, i.e., J = JQ , and equation (14) can be written as

2

X 1 + uM (ri1 , ri2 , . . . , riM ) , M ! i ,i ,···i

2

0



˜ Q (0) = 0 since there is no heat flux when the where J c system is in canonical equilibrium.

EQUATIONS OF MOTION

(17)

M

where Fii2 ...iM = − ∂r∂ i uM (ri , ri2 , . . . , riM ) is the M-body force contribution on atom i. The energy Ei of atom i is taken to be of the form Ei =

1 X 1 X pi · pi u2 (ri , ri2 ) + u3 (ri , ri2 , ri3 ) + 2m 2! i 3! i ,i 2

2

3

1 X + ...+ uM (ri , ri2 , . . . , riM ) . M ! i ,...,i 2

M

(18)

 ˜ Q (∞) J VT

=



1 V kB T 2

Z

0

∞



 ˜ Q ((s)0 ) ⊗ J ˜ Q (0)ic ds Fe , J

(15) where V is the volume of the system. The dynamics ˜ Q ((s)0 ) in (15) generates the ensemble of used to obtain J starting states following fc (Γ, ζ), thus making the theory ergodically consistent. In the range of Fe where 

the lin˜ Q (∞) J ear response theory is valid, (15) indicates that VT is linearly proportional to Fe once the steady state is reached with the constant of proportionality being the thermal conductivity tensor κ as in the GK method. Since the system has been  assumed to be ergodic, the ˜ Q (∞) is equivalent to the time averensemble average J ˜ Q (t) with respect to the field dependent dynamics age of J once the linear non-equilibrium steady state is reached. In the following section, the equations of motion, consistent with the preceding HNEMD method, to compute the thermal conductivity of a system consisting of N identical particles modeled by a potential with up to M-body terms are proposed.

based on the assumption that the energy from any Kbody term uK (ri1 , ri2 , . . . , riK ) is divided equally among the atoms i1 , i2 , . . . , iK . For this choice of the energy, the heat flux vector with the field-free equations of motion is given by JQ =

X pi Ei m

i

+

Fi i · pi1 1 X ri1 i2 ( 1 2 ) 2! i ,i m 1

2

Fi i i · pi1 1 X (ri1 i2 + ri1 i3 )( 1 2 3 ) + ... + 6! i ,i ,i m 1

2

3

X Fi i ...i · pi1 1 (ri1 i2 + ri1 i3 + . . . + ri1 iM )( 1 2 M ), + M ! i ,i ,...i m 1

2

M

(19) where ri1 i2 = ri1 − ri2 , provided the system momentum is zero [19]. (Note that other partitions of energy are possible but the choice does not appear to significantly affect the thermal conductivity estimate in practice [21].) The heat flux vector can be rewritten in tensorial form

4 X ¯ = 1 where E Ei is the average instantaneous inN i ternal energy. This form of Di in (23) induces no net momentum and also satisfies the AIΓ condition

as JQ =

X

Ei1 I +

i1

1 X (ri1 i2 ⊗ Fi1 i2 ) 2! i 2

 1 X + (ri1 i2 + ri1 i3 ) ⊗ Fi1 i2 i3 + . . . 6! i ,i 2

+

3

1 X  (ri1 i2 + ri1 i3 M ! i ,...i 2

M

  pi1 + . . . + ri1 iM ) ⊗ Fi1 i2 ...iM m (20)

where I is the identity tensor. The criteria for generating the field-dependent equations of motion (1) are: (a) adiabatic incompressibility (AIΓ), (b) equivalence of virial heat flux and dissipative heat flux (J = JQ ), and (c) momentum preservation. Recalling P P theT definition for dissipative flux J = T pi i Di m − i Ci Fi in equation (10) and comparing it with the tensorial form of the heat flux vector in (20), suggests the form 1 X Fii2 ⊗ rii2 2! i

Di = Ei I +

2

1 X Fii2 i3 ⊗ (rii2 + rii3 ) + . . . + 3! i ,i 2

+

(21)

3

1 X Fii2 ...iM ⊗ (rii2 + . . . + riiM ) M ! i ,...,i M

2

with Ci = 0. However, this form violates momentum preservation X X X Di Fe = 0 . (22) Fi + p˙ i = i

i

1 X Fii2 ⊗ rii2 2! i 2 1 X 1 X − Fi1 i2 ⊗ ri1 i2 + Fii2 i3 ⊗ (rii2 + rii3 ) 2!N i ,i 3! i ,i

¯ + Di = (Ei − E)I

2

(24)

since Ci = 0. Using (23), the dissipative heat flux may be now written as pi X T Ci Fi − m i i X 1 X 1 X Ei I + = rii2 ⊗ Fii2 + (rii2 + rii3 ) ⊗ Fii2 i3 2! i 3! i ,i i 2 2 3  pi 1 X (rii2 + . . . + riiM ) ⊗ Fii2 ...iM + ...+ M ! i ,...,i m 2 M  1 X − E¯ + ri i ⊗ Fi1 i2 2!N i ,i 1 2 1 2 1 X + (ri1 i2 + ri1 i3 ) ⊗ Fi1 i2 i3 3!N i ,i ,i 1 2 3 X X 1 pi + (ri1 i2 + . . . + ri1 iM ) ⊗ Fi1 i2 ...iM . M !N i ,i ...,i m i

J=

X

Di T

1

2

M

(25)

i

P which requires i Di = 0. This can be remedied by subtracting the mean of each term, as in Evans’ original formula [11], leading to

1

 X ∂   X ∂  · Ci Fe + · Di Fe ∂qi ∂pi i i X ∂ ¯ · Fe = 0 , (Ei − E) = ∂p i i

2

X pi

= 0, (25) reduces to (19), hence J = JQ . m Thus, given (24), (25) and (22), the form (23) with Ci = 0 is a viable HNEMD method for calculating the thermal conductivity of a material described with multibody potentials.

Since

i

3

1 X − Fi1 i2 i3 ⊗ (ri1 i2 + ri1 i3 ) 3!N i ,i ,i 1

2

3

1 X Fii2 ...iM ⊗ (rii2 + . . . + riiM ) + ... + M ! i ,...,i 2 M X 1 − Fi1 i2 ...iM ⊗ (ri1 i2 + . . . + ri1 iM ) , M !N i ,i ...,i 1

2

M

(23)

[1] W. G. Hoover, Computational Statistical Mechanics (Elsevier, Amsterdam, 1991).

Acknowledgement

Sandia is a multiprogram laboratory operated by Sandia Corporation, a Lockheed Martin Company, for the United States Department of Energy under contract DEACO4-94AL85000.

[2] P. K. Schelling, S. R. Phillpot, and P. Keblinski, Phys.

5 Rev. B 65, 144306 (2002). [3] X. W. Zhou, S. Aubry, R. E. Jones, A. Greenstein, and P. K. Schelling, Phys. Rev. B 79, 115201 (2009). [4] F. H. Stillinger and T. A. Weber, Phys. Rev. B 31, 5262 (1985). [5] J. Tersoff, Phys. Rev. B 37, 6991 (1988). [6] J. Tersoff, Phys. Rev. B 39, 5566 (1989). [7] S. Weiner, P. Kollman, D. Case, U. C. Singh, C. Ghio, G. Alagona, S. P. . Jr., and P. Weiner, J. Am. Chem. Soc. 106, 765 (1984). [8] A. MacKerell, D. Bashford, M. Bellott, R. Dunbrack, J. Evanseck, M. Field, S. Fischer, J. Gao, H. Guo, S. Ha, et al., J. Phys. Chem. B 102, 3586 (1998). [9] S. Stuart, A. Tutein, and J. Harrison, J. Chem Phys 112, 6472 (2000). [10] A. D. Mistriotis, G. E. Froudakis, P. Vendras, and N. Flytzanis, Phys. Rev. B 47, 10648 (1993). [11] D. J. Evans, Phys. Let. A 91, 457 (1982). [12] P. Daivis and D. Evans, Molecular Physics 81, 1289 (1994), ISSN 0026-8976, URL http: //permalink.lanl.gov/object/view?what=info: lanl-repo/isi%/A1994NK82900003.

[13] P. Daivis and D. Evans, Chemical Physics 198, 25 (1995). [14] A. Perronace, J.-M. Simon, B. Rousseau, and G. Ciccotti, Molecular Physics 99, 1139 (2001), ISSN 00268976, URL http://permalink.lanl.gov/object/view? what=info:lanl-repo/ins%pec/7003368. [15] K. K. Mandadapu, R. E. Jones, and P. Papadopoulos, J. Chem. Phys. 130, 204106 (2009). [16] D. J. Evans and G. P. Morris, Statistical Mechanics of Non-equilibrium Liquids (Academic Press, 1990). [17] D. J. Evans and B. L. Holian, J. Chem. Phys. 83, 4069 (1985). [18] W. G. Hoover, Phys. Rev. A 31, 1695 (1985). [19] J. H. Irving and J. G. Kirkwood, J. Chem. Phys. 18, 817 (1950). [20] J. W. Martin, J. Phys. C: Solid State Phys. 8, 2837 (1975). [21] G. Marechal and J.-P. Ryckaert, Chemical Physics Letters 101, 548 (1983), ISSN 0009-2614, URL http://permalink.lanl.gov/object/view?what=info: lanl-repo/ins%pec/2165965.