A Homogeneous Non-Equilibrium Molecular Dynamics Method for Calculating the Heat Transport Coefficient of Mixtures and Alloys Kranthi K. Mandadapu,1 Reese E. Jones,2, ∗ and Panayiotis Papadopoulos3 1

Department of Mechanical Engineering,

University of California, Berkeley, CA 94720-1740, USA 2

Sandia National Laboratories, Livermore, CA 94551-0969, USA 3

Department of Mechanical Engineering,

University of California, Berkeley, CA 94720-1740, USA

Abstract This work generalizes Evans’ homogeneous non-equilibrium method for estimating heat transport coefficient to multi-species molecular systems described by general multi-body potentials. The proposed method, in addition to being compatible with periodic boundary conditions, is shown to satisfy all the requirements of Evans’ original method, namely adiabatic incompressibility of phase space, equivalence of the dissipative and heat fluxes and momentum preservation. The difference between the new equations of motion, suitable for mixtures and alloys, and those of Evans’ original work are quantified by means of simulations for fluid Ar-Kr and solid GaN test systems.

∗

Corresponding author; Electronic address: [email protected]

1

I.

INTRODUCTION

The heat transport coefficient relating heat flux to the temperature gradient in the linear irreversible thermodynamic regime is an important phenomenological parameter for the modeling of heat conduction. For single-species systems and alloys without diffusion, this transport coefficient reduces to the usual thermal conductivity tensor κ. Molecular dynamics (MD) simulations can be employed to estimate this coefficient given the atomic structure and an interatomic potential that accurately describes the interactions among the atoms [1–7]. The two most commonly used methods for estimating the heat transport coefficient are: (i) the direct method [3, 6] and (ii) the Green-Kubo (GK) method [2, 3]. An alternative to these methods, the so-called Homogeneous Non-Equilibrium Molecular Dynamics (HNEMD) method [4, 7–12] has advantages over the GK method in that it is free of difficulties involving the calculation and integration of the heat flux autocorrelation tensor, and the direct method which has to contend with strong size effects and requires unrealistically large temperature gradients. In addition, the HNEMD method typically yields better statistical averages than both the direct and the GK method at a lower overall computational cost. In the HNEMD method, thermal transport is modeled by means of a mechanical analogue, where a fictitious external force mimics the temperature gradient. Using the linear response theory [5, 7], the long-time ensemble average of the heat flux vector can be shown to be proportional to the external force field (when the latter is sufficiently small), with the constant of proportionality being the heat transport coefficient tensor. The HNEMD method was originally developed for systems modeled by Lennard-Jones (LJ) pair potentials [4]. It was subsequently applied to molecular liquids with distance and angle constraints using an intra-molecular fourbody potential [10, 11, 13], and further extended to systems modeled by three-body potentials [7] and recently to M-body potentials [12]. The HNEMD equations of motion were developed to satisfy three basic conditions: adiabatic incompressibility of phase space (AIΓ), preservation of zero total momentum, and equivalence of dissipative and heat fluxes. These equations are suitable for single-species systems, since they yield the heat transport coefficient without having to calculate any average quantities from equilibrium simulations. If these HNEMD equations of motion are employed for multi-species systems, e.g., a semi-conductor Gallium-Nitride system 2

[14, 15], the long-time average of the heat flux vector still remains proportional to the external field. However, the constant of proportionality is not just equal to the heat transport coefficient, but involves additional correlation integrals, since the HNEMD equations do not satisfy the AIΓ condition and the dissipative flux is not equivalent to the heat flux vector as in the case of single-species systems. Hence, these correlation integrals need to be evaluated from an equilibrium simulation in addition to running the HNEMD algorithm. An earlier attempt to obtain non-equilibrium equations of motion specifically tailored to multi-species also requires equilibrium MD simulations in addition to running the non-equilibrium algorithm as the AIΓ condition is not satisfied [16]. As a result, the development of a non-equilibrium method to obtain the heat transport coefficient for multi-species systems without resorting to additional equilibrium simulations remains an open problem. This paper proposes precisely such a method, termed the Mixture-HNEMD (M-HNEMD) method. This method preserves all three conditions cited earlier and, in addition, remains compatible with periodic boundary conditions. The paper is organized as follows: in Section II, a review of basic irreversible thermodynamics is given within the context of heat conduction for binary mixtures and alloys. The linear response theory is derived in Section III for multi-species systems modeled by multi-body potentials. In Section IV, the estimate of the heat transport coefficient for multi-species systems is analyzed for the HNEMD method. Next, the M-HNEMD equations of motion are developed in Section V, followed by an application of HNEMD and MHNEMD equations of motion to a pure Argon (Ar) system, an Argon-Krypton (Ar-Kr) binary mixture and a perfect Gallium-Nitride (GaN) crystal in Section VI.

II.

LINEAR IRREVERSIBLE THERMODYNAMICS

Following the theory of irreversible processes presented in [17, Section 44], the total rate of entropy production σ for a mixture of n species capable of diffusion, heat conduction and the cross phenomena (Sor´et and Dufour effects) in the absence of external forces on any individual species is given by σ =

n X

Ji · X i + JQ · X Q .

i=1

3

(1)

µ 

∇T are the generalized forces driving the irreversible T T2 processes of diffusion and heat transfer, respectively, where µi is the concentration of species

Here, Xi = −∇

i

and XQ = −

i and T is the thermodynamic temperature. Also, JQ is the heat flux vector and Ji are the diffusive flux vectors, which satisfy

n X

Ji = 0 .

(2)

i=1

Taking into account (2), the total rate of entropy production can be equivalently expressed as σ =

n−1 X

Ji · (Xi − Xn ) + JQ · XQ .

(3)

i=1

Reducing (3) to a binary mixture (n = 2) results in σ = J1 · (X1 − X2 ) + JQ · XQ .

(4)

Assuming the irreversible processes are close to thermodynamic equilibrium, the linear phenomenological relations J1 = L11 (X1 − X2 ) + L1Q XQ

(5)

JQ = LQ1 (X1 − X2 ) + LQQ XQ

(6)

are proposed, where L11 , L1Q , LQ1 and LQQ are phenomenological coefficient matrices that satisfy Onsager’s reciprocal relations LT11 = L11 , LT1Q = LQ1 and LTQQ = LQQ , see [17–19]. Herein, the superscript T denotes the matrix transpose. ˜ 1 (t) and JQ = J ˜ Q (t), where J ˜ 1 (t) and J ˜ Q (t) represent the phase functions J1 Setting J1 = J and JQ as functions of time (see Section III), and following the derivations in [20], the GreenKubo relations for L11 , L1Q and LQQ are given in terms of the respective correlation functions by L11 L1Q LQQ

Z ∞

 V ˜ 1 (t) ⊗ J ˜ 1 (0) dt , = J c kB 0 Z ∞

 V ˜ 1 (t) ⊗ J ˜ Q (0) dt , J = c kB 0 Z ∞

 V ˜ Q (t) ⊗ J ˜ Q (0) dt , = J c kB 0 4

(7) (8) (9)

where V is the volume, kB is the Boltzmann constant and equilibrium phase-space distribution.

 · c denotes the average over

It is important to emphasize that in the case of binary liquid mixtures the evaluation of LQQ is not sufficient for estimating the thermal conductivity κ. Indeed, experiments measuring thermal conductivity of binary liquid mixtures are typically conducted using the stationary state (J1 = 0). In this case, (5) implies X1 − X2 = −L−1 11 L1Q XQ , which reduces (6) to JQ = (LQQ − LQ1 L−1 11 L1Q )XQ .

(10)

Therefore, for binary liquid mixtures, admitting Fourier’s law in the form JQ = −κ∇T and recalling the definition of XQ leads to κ =

1 (LQQ − LQ1 L−1 11 L1Q ) . T2

(11)

In contrast to binary liquids, if in a perfect crystalline solid no significant mass diffusion exists (e.g., GaN at room temperature) one may neglect diffusion and the cross phenomena. In this case, XQ is the only significant generalized force, therefore the heat flux vector in the linear irreversible regime follows the phenomenological law JQ = LQQ XQ , thus yielding the thermal 1 conductivity as κ = 2 LQQ . T III.

LINEAR RESPONSE THEORY

In linear response theory [5, Chapter 5], the general form of the equations of motion, when perturbed by an external field Fe in the presence of a Nos´e-Hoover (NH) thermostat, is r˙ i =

pi + Ci (Γ)Fe , mi

p˙ i = Fi + Di (Γ)Fe − ζpi , X  N pi · pi ˙ζ = 1 − 3NkB T , Q i=1 mi

(12)

where mi is the mass of atom i and Γ = {(ri , pi ) , i = 1, 2, . . . , N} is the phase space with ri and pi being respectively the position and momentum vectors of the i-th atom in an N-atom system. Here, Fi is the interatomic force and Ci (Γ), Di (Γ) are the second-order tensor phase 5

variables which describe the coupling of the system to the applied external field Fe [21],[5, Chapter 5]. Also, in the rate equation for the thermodynamic friction coefficient ζ associated with the Nos´e-Hoover thermostat [22], Q is a parameter chosen to yield the canonical phasespace distribution in the absence of external force fields [37]. When Fe = 0, the phase-space distribution f (Γ, ζ, t) becomes the (extended) canonical distribution fc ,
−β H0 (Γ)+ 21 Qζ 2 e fc (Γ, ζ) = , Z(β) R −β H0 (Γ)+ 1 Qζ 2
1 2 dΓ dζ and H0 is the total internal energy where β = , Z(β) = e kB T H0 =

N X pi · pi i=1

2mi

+ Φ(r) ,

(13)

(14)

see [22]. The potential Φ(r) that describes the M-body interactions among the atoms is of the general form Φ(r) =

1 X 1 1 X u2 (ri1 , ri2 ) + u3 (ri1 , ri2 , ri3 ) + . . . + 2! 3! M! (i1 ,i2 )

(i1 ,i2 ,i3 )

X

uM (ri1 , ri2 , . . . , riM ) ,

(i1 ,i2 ,...,iM )

where M ≤ N, uK (ri1 , ri2 , . . . , riK ) characterizes all K-body interactions [23], and

(15) X

(i1 ,i2 ,...,ir )

denotes summation over all

N! (N −r)!

permutations of (1, 2, . . . , r). Taking into consideration (14)

and (15), the internal energy Ei of atom i may be defined as Ei =

pi · pi 1 X 1 X 1 + u2 (ri , ri2 ) + u3 (ri , ri2 , ri3 ) + . . . + 2mi 2! 3! M! (i2 )i

where

X

(i2 ,i3 )i

denotes summation over all

(N −1)! (N −r)!

X

uM (ri , ri2 , . . . , riM ) ,

(i2 ,...,iM )i

(16)

permutations of (1, 2, . . . , r) excluding i, see

(i2 ,i3 ,...,ir )i

[12]. Here, it can be seen that the sum of Ei over i yields the total internal energy H0 given by (14). The total inter-atomic force Fi on atom i is then given by Fi = −

X 1 X 1 ∂Φ = Fii2 + Fii2 i3 + . . . + ∂ri 2! (M − 1)! (i2 )i

(i2 ,i3 )i

X

Fii2 ...iM ,

(17)

(i2 ,...,iM )i

∂ uK (ri , ri2 , . . . , riK ) is the K-body force contribution on atom i and 2 ≤ ∂ri K ≤ M. Following the Irving-Kirkwood procedure for field-free equations of motion and zero

where Fii2 ...iK = −

6

˜ Q (t) in the total momentum [24], the macroscopic instantaneous heat flux vector JQ (Γ) = J absence of external forces on any individual species is obtained in tensorial form as N 1 X T pi ν , JQ (Γ) = V i=1 i mi

(18)

where ν i = Ei I+

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

(i2 ,i3 )i

X

Fii2 ...iM ⊗(rii2 +. . .+riiM ) ,

(i2 ,...,iM )i

(19) also see [25]. The second order tensor ν i is clearly a combination of the 3 × 3 identity tensor I scaled by the energy associated with atom i and the virial associated with atom i. Assuming that the external field Fe is small and time-independent, the distribution function f (Γ, ζ, t) can be obtained as f (Γ, ζ, t) = fc (Γ, ζ) −

Z

t 0


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

(20)

where i∆L(s)fc (Γ, ζ) = iL(s)fc (Γ, ζ) − iL0 fc (Γ, ζ), with iL, iL0 being the Liouvilleans corresponding to field-dependent (Fe 6= 0) and field-independent (Fe = 0) equations of motion, respectively [5, Section 5.1]. Since fc (Γ, ζ) is the steady-state solution of the equation ∂ f ∂t

= −iL0 f , it follows that iL0 fc (Γ, ζ) = 0 [22]. Hence,

i∆L(s)fc (Γ, ζ) = iL(s)fc (Γ, ζ)   ∂ ˙ ∂ ∂ ∂ ˙ ˙ ˙ = ·Γ+Γ· +ζ + ζ fc (Γ, ζ) ∂Γ ∂Γ ∂ζ ∂ζ X    X N N N N X X ∂ ∂ T pi T = Di Ci Fi · Fe fc (Γ, ζ) · Ci Fe + · Di Fe fc (Γ, ζ) + β − + ∂r ∂p m i i i=1 i=1 i=1   i=1 = βV B(Γ) − J(Γ) · Fe fc (Γ, ζ) , (21)

where 1 B(Γ) = βV

N N X X ∂ ∂ T · Ci + · DTi ∂r ∂p i i i=1 i=1

7

!

(22)

is the phase space compression factor [5, Section 3.3], and N X

1 J(Γ) = V

i=1

N

pi X T Di T − Ci Fi mi i=1

!

(23)

is defined as the “dissipative” flux in [5, Section 5.1]. In earlier works [4, 7, 12, 21], it was assumed that the equations of motion (12) satisfy the condition of adiabatic incompressibility of phase space (AIΓ), i.e., B = 0. However, in this paper, B 6= 0 is allowed in the derivation of the linear response theory. Using (20) and (21), the long-time ensemble average of the heat ˜ Q (t) evolving with the field-dependent equations of motion (12) is given by flux vector J

 ˜ Q (∞) J T

=



V kB T 2

Z

0

∞

  h i ˜ Q ((s)0 ) ⊗ J(0) ˜ ˜ J − B(0) ds Fe ,

(24)

c

˜ Q ((s)0 ) is obtained by solving the field-independent where the subscript “0” emphasizes that J equations of motion [7, 12]. Two conditions are now imposed on the equations of motion (12), and consequently on Ci and Di , for the determination of LQQ using (9). These are: (i) preservation of zero total momentum, N X

pi = 0 ,

(25)

i=1

and (ii) equivalency of fluxes, J − B = JQ .

(26)

Taking into account these two conditions, (24) becomes

   Z ∞ ˜ Q (∞)

J V ˜ Q ((s)0 ) ⊗ J ˜ Q (0)ic ds Fe . = J (27) T kB T 2 0

 ˜ Q (∞) J is directly proportional to As noted in [7, 12], assuming that Fe is small, the quantity T the external field Fe with constant of proportionality equal to the phenomenological coefficient LQQ matrix 2 . It should be mentioned here that the condition given by (26) can also be imposed T by constructing Ci and Di so that B = 0 and J = JQ . As a consequence of this alternative

construction, the equations of motion (12) also satisfy adiabatic incompressibility of phase space AIΓ. 8

IV.

THERMAL CONDUCTIVITY ESTIMATION BY THE HNEMD METHOD

As mentioned in Section I, the HNEMD equations of motion have been employed in the estimation of GaN thermal conductivity in both bulk and nanowire systems [14, 15]. However, it is shown here that the original HNEMD method [4, 12] is not consistent for mixtures and alloys due to the non-equivalency of the fluxes J − B and JQ (i.e., J − B 6= JQ ). To appreciate this point, recall that Ci and Di in the original HNEMD equations of motion [12] are given in the form

N 1 X Di = ν i − νj . N j=1

Ci = 0 ,

These choices are compatible with preservation of zero total momentum, since

(28) P

i

Di = 0. On

the other hand, in this case the fluxes B and J defined respectively in (22) and (23) are given by 1 B = βV

J =

1 V

where

N X pi mi i=1

!

N −1 , N ! N N N X X X pi 1 pi − ν Tj = JQ − A , ν Ti mi N j=1 mi i=1 i=1 N N 1 X T X pi ˜ A(Γ) = A(t) . νj V N j=1 m i i=1

(29)

(30)

(31)

is the difference between the heat flux JQ and the dissipative flux J. It can be seen from (29), (30) and (31) that the difference in the fluxes J and B is not equal to the heat flux JQ , which violates the second condition (26) in Section III. Therefore, in the HNMED method, the longtime average of the macroscopic heat flux vector JQ in (24) can be expressed with the aid (29) and (30) as D E ˜ Q (∞) J T

or

=



V kB T 2

Z

0

∞

  ˜ ˜ JQ ((s)0 ) ⊗ JQ (0) ds − c

D

˜ Q (∞) J T

E

=

V kB T 2

Z

∞ 0



 h i ˜ ˜ ˜ JQ ((s)0 ) ⊗ A(0) + B(0) ds Fe

1 (LQQ − LQA − LQB )Fe , T2 9

c

(32)

(33)

where LQA

V = kB

Z

V = kB

Z

and LQB

∞ 0

0

∞



˜ Q ((s)0 ) ⊗ A(0) ˜ J



ds ,

(34)

  ˜ ˜ JQ ((s)0 ) ⊗ B(0) ds ,

(35)

c

c

Hence, estimation of LQQ requires the calculation of LQA + LQB in addition to running the hJ˜Q (∞)i standard HNEMD algorithm to obtain . These extra correlation terms need to be T evaluated by means of a Green-Kubo method, which naturally involves the errors related to the integration and calculation of correlation function. If upon employing the Green-Kubo method LQA + LQB is found to be small compared to LQQ , then it is sufficient to use the classical HNEMD method given by (28) for estimating LQQ .

V.

M-HNEMD ALGORITHM

In this section, equations of motion are derived for a system of n different species which satisfy the two basic conditions (25) and (26) discussed in Section III. These equations are of the form (12), with particular choices for the tensor variables Ci and Di . In the following, a sequence of proposed forms for Ci and Di is suggested and refined to satisfy (25) and (26), hence showing the development of the final expressions given at the end of this section. As observed in Section IV, using Di in the form (28) violates the equivalence of fluxes (26). To circumvent this problem, Di may be initially modified as Di = ν i −

where m ¯ =

N X

N mi X νj , m ¯ j=1

(36)

mi . Using (36) with Ci = 0, the flux vector J in (23) becomes

i=1

 N  X  N N 1 X T pi 1 X T J = ν − ν pi V i=1 i mi m ¯ j=1 j i=1 and

N X

(37)

pi = 0, provided that the total initial momentum is zero, hence J = JQ . However, the

i=1

10

flux vector B in (22) is now equal to 1 B = βV

N N X pi 1 X pi − m m ¯ i i=1 i=1

!

=

N 1 X pi 6= 0 , βV i=1 mi

(38)

resulting, again, in J − B 6= JQ . One way to achieve B = 0 is to choose Ci = −ri ⊗

pi ; mi

(39)

however, this introduces additional terms in J, such that J = JQ +

N  X pi i=1

mi

 ⊗ ri Fi 6= JQ .

(40)

To simultaneously eliminate B and set J = JQ , one should also change Di to N N mi X mi X ν j − (Fi · ri )I + (Fj · rj )I . Di = ν i − m ¯ j=1 m ¯ j=1

(41)

While (39) and (41) enforce the conditions (25) and (26), they do not satisfy frame invariance, i.e., they do not yield the same behavior under arbitrary translations and rotations of the system. To impose frame invariance, one should subtract from the position vector ri of every atom a suitably defined vector ¯r, which convects with the body. In this case, (39) and (41) become Ci = − (ri − ¯r) ⊗

pi , mi

N mi X νj − Di = ν i − m ¯ j=1

! N mi X Fi · (ri − ¯r) − Fj · (rj − ¯r) I m ¯ j=1

(42)

There are various possible choices for ¯r which enforce frame invariance without violating N 1 X (25) and (26). One such choice is ¯r = rCM = mi ri , namely the position vector of the m ¯ i=1 mass center of the system. Other choices include the position vector of any fixed point in the system, provided that this point convects with any arbitrary translation and rotation of the whole system. One such choice would be, for example, the left-bottom corner ra of the system. However, all the aforementioned choices are incompatible with periodic boundary conditions (PBCs). Indeed, in the case of PBCs, all atoms that exit the system from one side enter it again 11

from the opposite side with the same momentum, see, e.g., [26, Section 1.5.2]. In this so-called wrapping process, the terms (42) with ¯r = rCM or ¯r = ra create large jumps in velocities and forces, as ri − ¯r is discontinuous in time as an atom reenters the periodic cell with different values of Ci and Di than those with which it exited, see Section VI for numerical evidence of this problem. In this work, the effect of any discontinuities due to incompatibility of Ci and Di with PBCs is minimized by choosing for each atom an ¯r that depends only on atoms in its neighborhood. Specifically, the neighborhood of an atom is now comprised only of atoms within the cut-off radius rc > 0 of the potential. This treatment is still time-discontinuous due to neighboring atoms entering or exiting the cutoff region of a given atom. Clearly, the magnitude of this time-discontinuity depends crucially on the thermally-driven motion in the material [38]. For instance, in solid alloys transitions of atoms to/from the cut-off region can be quite rare and become virtually non-existent for low enough temperatures. On the other hand, in dense fluids an atom may change neighbors much more frequently, see Section VI for further discussion. For this choice of ¯r, equation (42) can be written as pi , Ci = − (ri − ¯ri ) ⊗ mi ! N N   m X   mi X i Di = ν i − ν j − Fi · ri − ¯ri − Fj · rj − ¯rj I m ¯ j=1 m ¯ j=1

(43)

Ni 1 X rj and Ni is the number of atoms in the neighborhood of atom i within the where ¯ri = Ni j=1 j6=i

cut-off radius. Despite the time-discontinuity, the choice of (43) preserves total momentum and satisfies (26) with B = 0 since the partial derivatives of Ci and Di with respect to the phase variables are unaffected. Therefore, the equations of motion (12) with (43) can be used to evaluate LQQ and are referred to as the M-HNEMD equations of motion under periodic boundary conditions. Also, Ci in (43) remains bounded with increasing system size, unlike the choice of ¯r = rCM . Since B = 0, the equations of motion also satisfy the adiabatic incompressibility of phase space, as in the case of HNEMD method for single-species systems. A casual review reveals that, when reduced to a single-species system (mi = m), the MHNEMD algorithm is not identical to the HNEMD algorithm. This demonstrates the non12

uniqueness of the NEMD algorithms for evaluating thermal conductivity.

VI.

RESULTS

The HNEMD and M-HNEMD algorithms are applied to estimate the thermal conductivity LQQ κ of Ar and GaN systems and the transport coefficient 2 for an Ar-Kr system. A sequence of T the heat flux averages hJ˜Qx,k (∞)i corresponding to a decreasing sequence of Fex,k is computed to determine the linear regime of the system under the action of the external field Fe = (Fex,k , 0, 0). LQQ,xx In this case, κxx is estimated for Ar and GaN and is estimated for Ar-Kr systems using T2 the slope method described in [7]. Before the external field Fe is switched on, the system is equilibriated using the NH thermostat which generates the (extended) canonical ensemble given in (13). When the external field is applied, the duration of each of the simulations in the sequence is chosen to make the variance of the estimated hJ˜Qx,k (∞)i negligible and hence the reported error estimates are based on the regression errors. The time integration scheme used for the HNEMD method is described in [7] and is modified for the M-HNEMD method, as in Appendix A.

A.

Argon

As an initial test, the M-HNEMD algorithm with (43) is applied to a single-species (n = 1) ˚ The Ar system at its triple point (T = 86.5 K, face-centered cubic lattice of length 5.719 A). system is modeled using a Lennard-Jones potential (M = 2 in (15)) with N = 256 atoms, ˚ and periodic boundary conditions on the 4 × 4 × 4 unit cells. In cut-off radius of rc = 13A the potential, the characteristic energy is ǫ = 1.6539 × 10−21 J and the characteristic length is ˚ [27]. Each run is carried out for 106 time steps with a step-size of ∆t = 4.0f s. σ = 3.405A ˚−1 < The computationally tractable linear regime of Fex is found to be between 0.001 A ˜ ˚−1 , as shown in Figure 1(a), where κxx,k = hJQx,k (∞)i is independent of Fex,k . It is Fex < 0.1 A T Fex,k concluded from Figure 1(a) that the results from the M-HNEMD and HNEMD algorithms are practically identical, thus validating the M-HNEMD algorithm for a single-species system. The value of thermal conductivity using the M-HNEMD algorithm with the slope method is found 13

to be 0.1305 ± 0.0002 W/mK, see Figure 1(b). This compares favorably to the HNEMD result 0.1302 ± 0.0002 W/mK and the Green-Kubo result 0.129 W/mK, see Figure 2. As an aside, it is also observed that when using ¯r = rCM in (42) the thermal conductivity κxx shows a decreasing trend for increasing values of Fex until it reaches a regime of severe discontinuity, see Figure 3. As argued in Section V, the discontinuity of κxx (Fex ) can be attributed to the non-smooth transition during the re-mapping due to the PBCs. This problem is eliminated with the use of eq. (43). Interestingly, the use of equations (43) also eliminates the decreasing trend of κxx , which facilitates the identification of a linear regime.

B.

Argon-Krypton

The M-HNEMD algorithm (43) is applied to the binary Ar-Kr mixture system at T = 115.6 K. The system is comprised of a total of N = 256 atoms (128 Ar and 128 Kr) in a cube of volume 24.1923 ˚ A3 with PBCs. It is modeled using a Lennard-Jones potential, where the ˚ [28]. The two parameters for Ar-Kr parameters for Kr are ǫ = 2.3056 × 10−21J and σ = 3.670A pairs are determined from the Lorentz-Berthelot mixing rules [29]. The equations of motion are integrated with step-size of ∆t = 4.0 f s for 106 time steps. ˚−1 < Fex < Figure 4(a) shows that the linear regime is found approximately when 0.01 A ˚−1 for both the M-HNEMD and HNEMD algorithms. The M-HNEMD estimate for the 0.085 A transport coefficient LQQ,xx /T 2 is found to be 0.08743 ± 0.00011 W/mK, see Figure 4(b). This is in good agreement with the Green-Kubo estimate of 0.08720 W/mK at a correlation time of t = 1ps, see Figure 5. The HNEMD estimate is found to be 0.08618 ± 0.00013 W/mK, which is very close to the Green-Kubo and M-HNEMD estimates, see Figure 5. As mentioned in Section LQA,xx + LQB,xx IV, since = 0.00043 W/mK is small compared to the Green-Kubo estimate T2 from Figure 5, either the HNEMD or the M-HNEMD algorithm may be employed to obtain an accurate estimate of LQQ,xx for Ar-Kr in the given state. Figure 6 shows the magnitude of the discontinuities in the trajectory of a typical ¯ri relative to the trajectory of the associated atom ri . These trajectories were captured before the field Fex was applied to allow for comparison of ¯ri to an unperturbed ri . Even in a dense fluid, it is apparent that the discontinuities are quite small relative to the normal fluctuations in position 14

0.16 0.14

κ [W/m−K]

0.12 HNEMD: κxx HNEMD: κyx HNEMD: κzx M−HNEMD: κxx M−HNEMD: κyx M−HNEMD: κzx

0.1 0.08 0.06 0.04 0.02 0 −0.02 0

0.02 0.04 0.06 0.08 0.1 Fex [1/Å]

0.12 0.14 0.16

(a) Variation of thermal conductivities κxx , κyx and κzx with Fex . 0.025

HNEMD: JQx HNEMD: JQy HNEMD: JQz M−HNEMD: JQx M−HNEMD: JQy M−HNEMD: JQz

〈JQ〉/T [W/m−K−Å]

0.02 0.015 0.01 0.005 0 −0.005 0

0.02 0.04 0.06 0.08 0.1 Fex [1/Å]

(b) Variation of the components of

0.12 0.14 0.16

˜ Q (∞)i hJ with Fex . T

FIG. 1: Comparison of the M-HNEMD and HNEMD algorithms for Ar at the triple point. The dashed and solid lines correspond respectively to the M-HNEMD and HNEMD estimates of κxx .

for the chosen time step.

C.

Gallium-Nitride

The HNEMD and M-HNEMD methods are employed to estimate the thermal conductivity of a perfect GaN system (wurtzite crystal structure) using an orthogonal cell with lattice parame-

15

NORMALIZED CORRELATION

1 0.8 0.6 0.4 0.2 0 −0.2 0

2

4

6

8

10

t [ps]

(a) Normalized Ar auto-correlation function hJ˜Qx (0)J˜Qx (t)i . hJ˜Qx (0)J˜Qx (0)i 0.14 0.12

κ [W/m−K]

0.1 0.08 0.06 0.04 0.02 0 0

2

4

6

8

10

t [ps]

(b) Integral of the auto-correlation function leading to thermal conductivity estimate κxx .

FIG. 2: Green-Kubo estimate of the thermal conductivity of Ar based on the direct integration of the heat flux auto-correlation function.

˚ and c = 5.20A ˚ at T=500K, see [6] for details on the lattice structure. The system ters a = 3.19A is modeled using 512 atoms with 4 unit cells per direction (4 × 4 × 4) and a Stillinger-Weber potential (M = 3 in eq. (15)), see [30, 31], under periodic boundary conditions. It is observed that 4 × 4 × 4 system is sufficient for obtaining the value of thermal conductivity using both HNEMD and M-HNEMD algorithms. A time-step of 0.5 fs is used for 8 × 106 time steps. The 16

2

κxx κyx κzx

κ [W/m−K]

1.5 1

0.2 0.15 0.1 0.05 0

0.5

0

0.04

0.08

0.04

0.06 Fex [1/Å]

0 −0.5 −1 0

0.02

0.08

0.1

0.12

FIG. 3: Variation of thermal conductivities κxx , κyx and κzx with Fex for Ar at the triple point for ˚−1 and a significant decreasing ¯r = rCM , demonstrating an instability at approximately Fex = 0.089 A trend (see inset).

˚−1 < Fez < 1.3 × 10−4 A ˚−1 linear regime in the HNEMD method is found in the range 7 × 10−5 A and the thermal conductivity is estimated via the slope method as 133.65 ± 2.50 W/mK, see Figure 7(a) and Figure 7(b). Using the M-HNEMD method, the thermal conductivity is found to be 132.02 ±5.85 W/mK, see Figure 7(a) and Figure 7(b), which is comparable to the HNEMD result. Both the HNEMD and M-HNEMD estimates are comparable to the estimate of the direct method estimate, which ˜ Q (t)i is 120.59 W/mK [6]. It can also be seen in Figure 8 that the average of the heat flux hJ ˚−1 thus demonstrating that the duconverges after 4 × 106 time steps for Fex = 1.0 × 10−4A ration of the simulations is sufficient. Finally, the Green-Kubo approach yields an estimate of

17

LQQ/T2 [W/m−K]

0.1 0.09 0.08 0.07 0.06 0.05 0.04 0.03 0.02 0.01 0 −0.01

2

HNEMD: LQQ,xx/T 2 HNEMD: LQQ,yx/T 2 HNEMD: LQQ,zx/T 2 M−HNEMD: LQQ,xx/T 2 M−HNEMD: LQQ,yx/T 2 M−HNEMD: LQQ,zx/T

0

0.02 0.04 0.06 0.08 0.1 Fex [1/Å]

0.12 0.14 0.16

(a) Variation of the transport coefficients LQQ,xx , LQQ,yx and LQQ,zx with Fex . 0.016

HNEMD: JQx HNEMD: JQy HNEMD: JQz M−HNEMD: JQx M−HNEMD: JQy M−HNEMD: JQz

〈JQ〉/T [W/m−K−Å]

0.014 0.012 0.01 0.008 0.006 0.004 0.002 0 −0.002 0

0.02 0.04 0.06 0.08 0.1 Fex [1/Å]

0.12 0.14 0.16

(b) Variation of the components of the heat flux vector ˜ Q (∞)i hJ with Fex . T

FIG. 4: Comparison of the M-HNEMD and HNEMD algorithms for Ar-Kr. The dashed and solid lines correspond respectively to the M-HNEMD method and HNEMD estimates of LQQ,xx /T 2 .

128.09 W/mK using 8×106 time steps and 1000 atoms (5×5×5). In the Green-Kubo approach, it is observed that a system at least as large as 5 × 5 × 5 is needed to obtain a limiting thermal conductivity value. In all cases, the Green-Kubo estimate is based on the criterion of the first ˜ Q (t)J ˜ Q (0)i plateau attained by the integral of the correlation, as in [3]. Figure 9 shows that hJ decays in a complicated manner with correlation time, unlike the Ar-Kr system, which makes 18

NORMALIZED CORRELATION

1 0.8 0.6

〈JQx(0).JQx(t)〉

0.4

〈JQx(0).Ax(t)〉 〈JQx(0).Bx(t)〉

0.2 0 −0.2 0

2

4

6

8

10

t [ps]

CORRELATION INTEGRAL [W/m−K]

(a) Correlations normalized by hJ˜Qx (0)J˜Qx (0)i. 0.1 0.09 0.08 0.07 0.06

LQQ,xx/T2

0.05

LQA,xx/T2

0.04

LQB,xx/T2

0.03 0.02 0.01 0 −0.01 0

2

4

6

8

10

t [ps]

(b) Integrals of the correlation functions.

˜x (t)i, as well as FIG. 5: Ar-Kr correlation functions hJ˜Qx (0)J˜Qx (t)i, hJ˜Qx (0)A˜x (t)i, and hJ˜Qx (0)B LQQ,xx LQA,xx LQB,xx integration of the correlation functions leading to , and estimates using the 2 2 T T T2 Green-Kubo method.

this criterion somewhat ambiguous. As seen in Figure 9, the noise past a correlation time of approximately 200 ps is quite significant, thus requiring considerably more simulation time than the M-HNEMD method to obtain a similar quality thermal conductivity estimate. However, (LQA,xx + LQB,xx ) /T 2 = 0.051 W/mK is extremely small compared to the Green-Kubo estimate leading to the conclusion that the HNEMD and M-HNEMD methods yield comparable results for GaN.

19

15

X POSITION [A]

14 13 12 11 10

− ri ri

9 0

20

40

60

80 100 120 140 160 180 200 TIME [ps]

FIG. 6: Trajectories ri and ¯ri of a typical atom i in an Ar-Kr system. VII.

DISCUSSION

As observed in Section VI, the Green-Kubo method may not be viable for the estimation of thermal conductivity in systems like GaN. In such systems, it is difficult to establish that the complex autocorrelation function has decayed sufficiently to permit accurate integration. In addition, given that the decay of the autocorrelation function is fairly long (at least 150-200 ps in the case of GaN), it takes considerable simulation time (of the order of 8 ns in GaN) to obtain the thermal conductivity within reasonable statistical certainty. These problems may be avoided by using the HNEMD and M-HNEMD methods, which involve the average of the heat flux vector alone as opposed to the integration of the noisy autocorrelation of the heat flux vector. Also, these non-equilibrium methods yield results of reasonable accuracy for relatively small total simulation times, as the signal-to-noise ratio is high due to the action of the finite non-zero external field Fe . For single-component systems, the HNEMD and M-HNEMD algorithms, while not mathematically identical, yield statistically identical results, which demonstrates the non-uniqueness in the construction of non-equilibrium molecular dynamics methods. Interestingly, even when applied to multi-component Ar-Kr and GaN systems, the HNEMD and M-HNEMD algorithms yield very similar results, as LQA + LQB is very small compared to the heat transport coefficient tensor LQQ . However, the M-HNEMD is theoretically consistent with the linear response 20

800

2

HNEMD: LQQ,xx/T 2 HNEMD: LQQ,yx/T 2 HNEMD: LQQ,zx/T 2 M−HNEMD: LQQ,xx/T 2 M−HNEMD: LQQ,yx/T 2 M−HNEMD: LQQ,zx/T

700

κ [W/m−K]

600 500 400 300 200 100 0 −100 0

5e−05

0.0001 Fex [1/Å]

0.00015

0.0002

(a) Variation of the thermal conductivities κxx , κyx and κzx with Fex . 0.05

HNEMD: JQx HNEMD: JQy HNEMD: JQz M−HNEMD: JQx M−HNEMD: JQy M−HNEMD: JQz

〈JQ〉/T [W/m−K−Å]

0.04 0.03 0.02 0.01 0 −0.01 0

5e−05

0.0001 Fex [1/Å]

0.00015

0.0002

(b) Variation of the components of the heat flux vector ˜ Q (∞)i hJ with Fex . T

FIG. 7: Comparison of the M-HNEMD and HNEMD algorithms for GaN. The dashed and solid lines correspond respectively to the M-HNEMD method and HNEMD estimates of κxx .

theory for multi-component system, in the sense that, in addition to satisfying preservation of total momentum (see eq. (25)) and being compatible with periodic boundary conditions, it also satisfies equivalency of fluxes (see eq. (26)). In contrast, for multi-component systems the HNEMD method violates eq. (26), thereby requiring additional equilibrium simulations to establish that LQA + LQB is small compared to LQQ . Although these results have not an-

21

0.007 0.006 〈 JQ 〉 [eV/Å2−ps]

0.005 0.004 0.003

〈JQx〉 〈JQy〉 〈JQz〉

0.002 0.001 0 −0.001 −0.002 −0.003 0

0.5

1

1.5

2 t [ns]

2.5

3

3.5

4

FIG. 8: The running average of the components of the heat flux vector JQ for GaN system for ˚−1 using the M-HNEMD method. Fex = 1.3 × 10−4 A

swered the question whether these quantities are insignificant for all material systems under all conditions conclusively, it does suggest that the original HNEMD method can provide accurate thermal conductivity estimates without strictly satisfying the equivalence of fluxes and the incompressibility of phase space. Nevertheless, in all cases, without resorting to additional simulations to calculate expensive time-correlations, the marginally more complex M-HNEMD algorithm will provide accurate results if used in the linear regime. As expected, the M-HNEMD method shares some of the same challenges that apply to other non-equilibrium methods that rely on the linear response theory. First, it requires simulations corresponding to a decreasing sequence of Fe to establish the linear regime. This can be alleviated by performing these simulations in parallel, as they are completely independent of each other. Second, it becomes extremely inefficient for very small Fe , where the signal-to-noise ratio is low. In such a case, it takes simulation time on the order of a Green-Kubo estimate to hJ˜Qx (∞)i obtain a reasonable estimate of as Fex tends to zero. Other NEMD methods may T Fex become applicable in this case, such as the one by Ciccotti et al. [32, 33], which employs a direct difference of heat flux between unperturbed and periodically perturbed trajectories to obtain the linear response of the system. In the NEMD under consideration, it is important to set Fex far enough from zero to avoid this problem, while still staying within the linear regime.

22

NORMALIZED CORRELATION

1 0.8

〈JQx(0).JQx(t)〉

0.6

〈JQx(0).Ax(t)〉

0.4

〈JQx(0).Bx(t)〉

0.2 0 1

−0.2 −0.4

0

−0.6 −1

−0.8 0

50

100

150 t [ps]

200

250

300

CORRELATION INTEGRAL [W/m−K]

(a) Correlations normalized by hJ˜Qx (0)J˜Qx (0)i. 140 120 100 80

LQQ,xx/T2

60

LQA,xx/T2

40

LQB,xx/T2

20 0 −20 0

50

100

150

200

250

t [ps]

(b) Integrals of the correlation functions

˜x (t)i, as well as inteFIG. 9: GaN correlation functions hJ˜Qx (0)J˜Qx (t)i, hJ˜Qx (0)A˜x (t)i and hJ˜Qx (0)B LQQ,xx LQA,xx LQB,xx gration of the correlations functions leading to κxx = , and estimates using the 2 2 T T T2 Green-Kubo method. Note that in (a) the data has been decimated to show each correlation clearly; the inset shows the full correlations which densely fill the envelopes described by the decimated data.

An estimate of the linear regime can be obtained from prior knowledge of the phonon mean-free path [7]. Third, and most important, the width of the linear regime decreases with an increasing size of the system, thereby making it difficult to obtain a reliable thermal conductivity estimate [34]. This problem may be rectified by adopting the methodology in [34], which has been shown to maintain a constant width of the linear regime.

23

Acknowledgments

This work was supported by 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.

[1] W. G. Hoover. Computational Statistical Mechanics. Elsevier, Amsterdam, 1991. [2] S. G. Volz and G. Chen. Molecular-dynamics simulation of thermal conductivity of Silicon crystals. Phys. Rev. B, 61(4):2651–2656, 2000. [3] 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. [4] D. J. Evans.

Homogeneous NEMD algorithm for thermal conductivity–application of non-

canonical linear response theory. Phys. Let. A, 91(9):457 – 460, 1982. [5] D. J. Evans and G. P. Morris. Statistical Mechanics of Non-equilibrium Liquids. Academic Press, 1990. [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. [7] K. K. Mandadapu, R. E. Jones, and P. Papadopoulos. A homogeneous non-equilibrium molecular dynamics method for calculating thermal conductivity with the three-body potential. J. Chem. Phys., 130:204106–1–9, 2009. [8] J. M. Luttinger. Theory of thermal transport coefficients. Phys. Rev., 135(6A):A1505–A1514, 1964. [9] M. J. Gillan and M. Dixon. The calculation of thermal conductivities by perturbed molecular dynamics simulation. J. Phys. C, 16(5):869–878, 1983. [10] P. J. Daivis and D. J. Evans.

Nonequilibrium molecular-dynamics calculation of thermal-

conductivity of flexible molecules - Butane. Molecular Physics, 81(6):1289–1295, 1994.

24

[11] A. Perronace, J.-M. Simon, B. Rousseau, and G. Ciccotti. Flux expressions and NEMD perturbations for models of semi-flexible molecules. Molecular Physics, 99(13):1139–49, 2001. [12] K. K. Mandadapu, R. E. Jones, and P. Papadopoulos. Generalization of the homogeneous nonequilibrium molecular dynamics method for calculating thermal conductivity to multi-body potentials. Phys. Rev. E, 2009. [13] P. J. Daivis and D. J. Evans. Temperature dependence of the thermal conductivity for two models of liquid Butane. Chemical Physics, 198(1-2):25–34, 1995. [14] Z. Wang, X. Zu, F. Gao, W. J. Weber, and J.-P. Crocombette. Atomistic simulation of the size and orientation dependences of thermal conductivity in GaN nanowires. Appl. Phys. Lett., 90(16):161923, 2007. [15] Z. Wang, F. Gao, J. P. Crocombette, X. T. Zu, L. Yang, and W. J. Weber. Thermal conductivity of GaN nanotubes simulated by nonequilibrium molecular dynamics. Physical Review B, 75(15):153303–1–153303–4, 2007. [16] D. J. Evans and P. T. Cummings. Non-equilibrium molecular dynamics algorithm for the calculation of thermal diffusion in simple fluid mixtures. Mol. Phys., 72(4):893–898, 1991. [17] S. R. de Groot. Thermodynamics of Irreversible Processes. Interscience Publishers Inc., New York, 1951. [18] L. Onsager. Reciprocal relations in irreversible processes. I. Phys. Rev., 37(4):405–426, 1931. [19] L. Onsager. Reciprocal relations in irreversible processes. II. Phys. Rev., 38(12):2265–2279, 1931. [20] R. Zwanzig. Elementary derivations of time-correlation formulas for transport coefficients. J. Chem. Phys., 40(9):2527–2533, 1964. [21] D. J. Evans and B. L. Holian. The Nose-Hoover thermostat. J. Chem. Phys., 83:4069–4074, 1985. [22] W. G. Hoover. Canonical dynamics: Equilibrium phase-space distributions. Phys. Rev. A, 31:1695–1697, 1985. [23] J. W. Martin. Many-body forces in metals and the Brugger elastic constants. J. Phys. C: Solid State Phys., 8:2837–2857, 1975. [24] 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.

25

[25] R. J. Hardy. Formulas for determining local properties in molecular dynamic simulations: Shock waves. J. Chem. Phys, 76(1):622–628, 1982. [26] M. P. Allen and D. J. Tildesley. Computer Simulation of Liquids. Oxford University Press, New York, 1987. [27] L. A. Rowley, D. Nicholson, and N. G. Parsonage. Monte Carlo grand canonical ensemble calculation in a gas-liquid transition for 12-6 Argon. J. Chem. Phys, 17:401 – 414, 1975. [28] M. Moosavi and E. K. Gaharshadi. Molecular dynamics simulations of some thermodynamic properties of mixtures of Argon with Neon, Krypton and Xenon using two-body and three-body interaction potentials. Fluid Phase Equilibria, 274:51–58, 2008. [29] R. C. Reid, J. M. Prausnitz, and B. E. Poling. The properties of gases and liquids. McGraw-Hill, New York, 4th edition, 1987. [30] A. Bere and A. Serra. Atomic structure of dislocation cores in GaN. Phys. Rev. B, 65(20):205323, 2002. [31] A. Bere and A. Serra. On the atomic structures, mobility and interactions of extended defects in GaN: dislocations, tilt and twin boundaries. Philos. Mag, 86:2152–2192, 2006. [32] G. Ciccotti, G. Jacucci, and I. R. McDonald. Thought-Experiments by Molecular-Dynamics. J. Stat. Phys., 21(1):1–22, 1979. [33] C. Massobrio and G. Ciccotti. Lennard-Jones Triple-Point Conductivity via Weak External Fields. Phys. Rev. A, 30(6):3191–3197, 1984. [34] D. P. Hansen and D. J. Evans. A generalized heat flow algorithm. Mol. Phys., 81(4):767 – 779, 1994. [35] G. J. Martyna, M. L. Klein, and M. E. Tuckerman. Nos´e-Hoover Chains - The Canonical Ensemble via Continuous Dynamics. J. Chem. Phys., 97(4):2635–2643, 1992. [36] 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. [37] It should be noted that the property of the Nos´e-Hoover thermostat to sample the canonical phase-space distribution relies on the system being ergodic under Nos´e-Hoover dynamics. For systems that are not ergodic in this manner, e.g., certain small or stiff systems, the technique of

26

Nos´e-Hoover chains [35] provides a remedy. It is apparent that the systems analyzed in this work are well-behaved in this respect and hence Nos´e-Hoover chains are not applied. [38] The jumps in ¯ri are proportional to

n Ni rc ,

where n is the net efflux of atoms from the cutoff region

over a time-step. Given a well-defined density, the long-time expectation of

n Ni

is zero; but on

shorter time-scales, n will be proportional to the surface area of the cutoff region and be affected by the diffusion and density fluctuation characteristics of the material at the given thermodynamic state. Specifically, n will be related to the spatio-temporal two-point density correlation function of the system integrated over the cutoff region.

Appendix A: Numerical integration algorithm

Following the work in [7], the governing system of equations (12) can be decomposed into               1 0 ri p 0 0 0 C (r, p)Fe           m i  i  d                0 pi  =   + −ζpi  + Di (r, p)Fe  + Fi (r) +  0  +   0   dt     ¯          T (p) 1 − 1 ζ 0 0 0 0 0 Q T | {z } | {z } | {z } | {z } | {z } | {z } L1

L2

L3

L4

L5

L6

(A.1) PN pi · pi 1 . The operator-split method in [36] is employed to integrate where T¯(p) = 3NkB i=1 mi the preceding system of ordinary differential equations. In this method, the propagation Γ(t) = exp(iLt)Γ(0) of the initial state is approximated according to         ∆t ∆t ∆t ∆t exp(iL∆t) = exp iL1 exp iL2 exp iL3+4 exp iL5 exp(iL6 ∆t) 2 2 2 2         ∆t ∆t ∆t ∆t exp iL3+4 exp iL2 exp iL1 , exp iL5 2 2 2 2

27

(A.2)

(A.3)

where higher-order terms have been omitted and iL is the p-Liouvillean [5, Section 3.3] . The update formulae for the individual steps are defined as   ∆t 1 T¯ (pn ) L1 : ζn+1/2 = ζn + −1 2 Q T   ∆t ˜ n+1/2 = exp − ζn+1/2 pn L2 : p 2 ∆t ˜ n+1/2 + Fn L3+4 : pn+1/2 = p 2m ∆t L5 : rn+1/2 = rn + pn+1/2 2m L6 : ˜rn+1 = rn+1/2 + ∆t C(rn+1/2 , pn+1/2 ) Fe ∆t pn+1/2 L5 : rn+1 = ˜rn+1 + 2m ˜ n+1 = F(rn+1 ) F ˜ n+1 + D(rn+1 , pn+1/2 , F ˜ n+1 )Fe Fn+1 = F ∆t ˜ n+1 = pn+1/2 + Fn+1 L3+4 : p 2m   ∆t ˜ n+1 L2 : pn+1 = exp − ζn+1/2 p 2   ∆t 1 T¯(pn+1 ) L1 : ζn+1 = ζn+1/2 + −1 2 Q T

(A.4) (A.5) (A.6) (A.7) (A.8) (A.9) (A.10) (A.11) (A.12) (A.13) (A.14)

where, in the interest of brevity, all the i subscripts referring to atoms have been omitted. Here the subscript n refers to the time-step. Also note that F0 is defined as F0 = F(r0 ), given initial conditions r0 and p0 . All the steps of this algorithm involve exact integration of the respective propagators, with the exception of (A.8) and (A.12) which involve integration of the forces Ci Fe and Di Fe . Indeed, (A.12) is exact except for the part of Di that depends on kinetic energy and, consequently on terms involving squares of the components of pi . Moreover, Di involves the momenta of each atom through the average kinetic energy embedded in Di . Likewise, (A.8) is inexact due to the dependence of Ci on a set of rj , which couples the integration for atom i to the motions of its neighboring atoms.

28