Adaptive Green-Kubo estimates of transport coefficients from molecular dynamics based on robust error analysis Reese E. Jones1, ∗ and Kranthi K. Mandadapu2, † 1

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

Department of Mechanical Engineering,

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

Abstract We present a rigorous Green-Kubo methodology for calculating transport coefficients based on on-the-fly estimates of : (a) statistical stationarity of the relevant process, and (b) error in the resulting coefficient. The methodology uses time samples efficiently across an ensemble of parallel replicas to yield accurate estimates, which is particularly useful for estimating the thermal conductivity of semi-conductors near their Debye temperatures where the characteristic decay times of the heat flux correlation functions are large. Employing and extending the error analysis of Zwanzig and Ailawadi (1969) and Frenkel (1980) to the integral of correlation, we are able to provide tight theoretical bounds for the error in the estimate of the transport coefficient. To demonstrate the performance of the method, four test cases of increasing computational cost and complexity are presented : the viscosity of Ar and water, and the thermal conductivity of Si and GaN. In addition to producing accurate estimates of the transport coefficients for these materials, this work demonstrates precise agreement of the computed variances in the estimates of the correlation and the transport coefficient with the extended theory based on the assumption that fluctuations follow a Gaussian process. The proposed algorithm in conjunction with the extended theory enable the calculation of transport coefficients with the Green-Kubo method accurately and efficiently.

∗

Corresponding author; Electronic address: [email protected]

†

currently at Sandia National Laboratories, Livermore, CA

1

I.

INTRODUCTION

The Green-Kubo (GK) method [1–3] is a well-established technique in the physics [4–6], chemistry [7, 8] and mathematics communities [9, 10] for determining transport coefficients. In essence, a GK estimate of a transport coefficient relies on the integral of an accurate time-correlation of the equilibrium fluctuations of the corresponding flux. When employed with molecular dynamics (MD), it is found empirically that Green-Kubo methods yield bulk transport coefficients with relatively little system size dependence [4, 6, 11] compared to “direct” simulation methods [12]. However, the GK method can suffer from issues related to convergence of the correlation function that results in imprecise heuristics and ad hoc estimates. In fact, reported transport coefficients estimated with GK techniques can vary considerably, e.g. for GaN [12, 13]. There are a multitude of reasons for discrepancies in the estimates. For example, mundane issues like different potentials and parameterizations may cause discrepancies in the estimates of transport coefficients [14] due to differences in the underlying properties, such as the effects of elastic wave speed on thermal conductivity. Another clearly addressable issue is poor quality estimates due to insufficient averaging time or premature truncation of the correlation integral. The intent of this work is to define and employ quantifiable convergence metrics to address exactly this last source of error. To estimate the statistical error in the transport coefficient, typically the variance of a small ensemble of long-time averaged estimates is used [4, 6, 12]. Unfortunately, the size of the ensemble, in many instances, is barely sufficient to generate viable statistics due to the expense of each of the long-time estimates. Moreover, the error estimates based on this ensemble alone ignore the variance intrinsic in each long-time estimate. Therefore, the error calculated in this fashion is an over-estimation of the actual error leading unnecessary additional computation. These coarse error estimates can also obscure the determination of the convergence of the mean integral of correlation to a constant value over the range of correlation times computed. Clearly, there is a need for accurate estimation of the error in the transport coefficient calculated by the Green-Kubo procedure. In this work, we measure and control the error due to: (a) use of data before the system has reached equilibrium from its initial conditions, (b) truncation of the integral of the correlation, and (c) finite averaging time. We also discuss various aspects of the time-discretization of the relevant process. To this end, we establish

2

a rigorous on-the-fly technique for determining a transport coefficient to within a pre-set error tolerance using a mixture of parallel replicas and ergodic time samples. We calculate the error in the mean of this mixture using a block averaging technique [15, 16] and [17, Section 6.4.1] that accounts for interdependence of the samples. This idea of mixing time and ensemble averages dates back at least to the late 1960’s [18] for calculating accurate correlations; however, it has not been put into practice to calculate transport coefficients in conjunction with rigorous error analysis. In the next section, we introduce the theory necessary for calculating estimates of the correlations, transport coefficients and their theoretical errors. Specifically, in Appendix A, we develop an extension of the Zwanzig-Aliwadi/Frenkel Gaussian process (GP) based error analysis for correlation function to the integral of correlation and hence the transport coefficient. In Section III and Appendix B, we develop the corresponding adaptive and efficient algorithms to compute the estimates using MD. In Section IV, we apply the method to calculate, directly and without assumptions, the transport coefficients and the associated errors directly for a range of materials and compare these results to the newly developed theoretical estimates.

II.

THEORY

Molecular dynamics (MD) provides the trajectories {xα (t)} of a set of atoms {α = ¨ α = fα , given initial condi1, . . . , N } by integrating Newton’s equations of motion, mα x tions, masses mα , and forces fα . The total force fα on an atom α is the sum of interatomic forces derived from an empirical potential Φ({xα }) and the Nos´e-Hoover thermostat [38],[26, Chapter 5]. The initial conditions consist of the initial positions {xα (0)} and momenta {pα (0)} = {mα vα (0)}, where vα = x˙ α are the velocities. The total interatomic potential Φ is typically sub-structured with “bonds” based on pairs (αβ), triplets (αβγ), etc., of distinct atoms so that it can be represented in the form Φ=

1 X 1 X φ2 (xα , xβ ) + φ3 (xα , xβ , xγ ) + . . . 2! αβ 3! αβγ

and may include long-range electrostatic interactions.

3

(1)

The expression for the stress σ [19–22] is given by:   1 X 1 X 1 X 1 X σα = − mα vα ⊗vα + σ=− fαβ ⊗ xαβ + fαβγ ⊗ (xαβ + xαγ ) + . . . , V α V α 2! β 3! βγ | {z } να

(2) where V is the volume of the system, xαβ ≡ xα − xβ , fαβ ≡ − ∂x∂ α φ2 (xα , xβ ), fαβγ ≡ − ∂x∂ α φ3 (xα , xβ , xγ ) and ν α is the virial for atom α. A corresponding formula [21, 22] for the heat flux J is J=


1 X εα I + ν Tα vα , V α

(3)

where the per-atom energy α is formed from the kinetic energy of the atom α and a reasonable partition of the total potential energy [4, 63]: 1 1 X 1 X φ2 (xα , xβ ) + φ3 (xα , xβ , xγ ) + . . . . εα = mα vα · vα + 2 2! β 3! βγ

(4)

With a definition of the appropriate flux, the non-equilibrium transport coefficients shear viscosity ν and thermal conductivity κij can be obtained from the Green-Kubo formulae [1–3, 23–25] ν =

Z

V kB T

κij



E ς(0) − hςi ς(t) − hςi dt

(5)

0

and Z

V = kB T 2

∞D

∞D



E Ji (0) − hJi i Jj (t) − hJj i dt ,

(6)

0

respectively [64]. Here, ς is any of the off-diagonal components of the momentum-flux/stress tensor {σ12 , σ13 , σ23 }, Ji is the heat flux in i-direction, T is the temperature, kB is the Boltzmann constant and h·i denotes the average over the equilibrium ensemble. It is important to note that the averages of stress hςi and heat flux hJi corresponding to an equilibrium system are zero by definition. Generically, the Green-Kubo formulae (5)-(6) can be written as Z

∞

Z hχ(0)χ(t)i dt = λ

ξ = λ 0

∞

C(t) dt ,

(7)

0

where ξ is the transport coefficient, χ is the fluctuation in flux, e.g. , χ(t) = ς(t) − hςi, and λ is a constant containing kB , T and V . To obtain the transport coefficient ξ, we need to calculate the auto-correlation function of χ denoted by C(t), and its integral. Clearly, the accuracy of the transport coefficient ξ depends directly on the accuracy of the calculated correlation and its integral. With MD, usually the correlation function C(t) is 4

calculated as a time-average from the data of a single simulation using the ergodic hypothesis [26, Section 3.4], [27, Section 4.2] : 1 C(t) = lim τ →∞ τ

Z

τ

χ(s) χ(s + t) ds = lim Cτ (t) . τ →∞

0

(8)

Here, Cτ (t) is an estimate of the correlation function C(t) for a finite simulation time τ . Truncation of the limit in Eq. (8) leads to what we call averaging error. Likewise, the estimate of the transport coefficient ξ is given by Z τc Cτ (t) dt , ξτ (τc ) = λ

(9)

0

which also requires truncation of the upper limit τc → ∞ in Eq. (7) for a finite time simulation. This truncation leads to what we call truncation error, which compounds the averaging error already present in Eq. (9). Clearly, τc should be much greater than the characteristic decay time τC of the correlation in order to minimize truncation error; therefore, we call τc the maximum significant correlation time. The error in the estimate of ξ depends on the variance in the estimate of the correlation, which is given by
2 var Cτ (t) ≡ h Cτ (t) − hCτ i i = hCτ (t)2 i − hCτ (t)i2 ,

(10)

where recall h·i denotes the average over the equilibrium ensemble and hCτ (t)i ≡ C(t). For a finite simulation time τ and assuming that the fluctuations χ follow a Gaussian process (GP), the variance in Cτ (t) can be calculated theoretically for the following two cases: (a) short correlation times, t → 0, var Cτ (0) = 4

τC 2 hCτ (0)i2 , τ

(11)

τC 2 hCτ (0)i2 , τ

(12)

and (b) long correlation times, t → ∞, var Cτ (∞) = 2 where τC 2 is given by

Z τC 2 =

∞

hCτ (s)i2 ds

0

hCτ (0)i2

,

(13)

refer to [18, 28], [17, Section 6.4], and [29, Appendix D]. From Eqs. (11) and (12), it can be seen that the variance in Cτ (t) changes only by a factor of 2 over the entire range of 5

correlation times, given the assumption that χ follows a Gaussian process. Note that (a) Eqs. (11) and (12) are obtained by assuming that the total simulation time τ is much longer than the characteristic decay time τC of the correlation, and (b) equal number of samples are used to obtain Cτ (t) at both shorter and longer correlation times. Also, τc , τC , and τC 2 denote different but related characteristic times. By applying GP theory to the integral of correlation in Eq. (9), we obtain the derivative of the variance of ξτ (τc ) with respect to τc for the following two cases: (a) τc → 0,

and (b) τc → ∞,

d var ξτ = 0 , dτc

(14)

ξ2 d var ξτ = 4 , dτc τ

(15)

see Appendix A for the derivation. In deriving Eq. (15), we assumed that τc and τ are much greater than the characteristic decay time of the correlation τC . Eqs. (14) and (15), and the fact that var ξτ (τc ) is a monotonically increasing function of τc , allow us to conclude that var ξτ (τc ) < 4ξ 2

τc , τ

(16)

where ξ is the true transport coefficient. Therefore, we can apply Eq. (16) to obtain the variance in estimate of the transport coefficient by replacing ξ with a sufficiently accurate ξτ . Given this approximation, the relative error in the estimate of ξ should be bounded by p 2 ττc . This extension of the Zwanzig and Ailawadi theory, which was limited to the error in the correlation function, to the integral of correlation gives us a direct estimate of the error in the calculated transport coefficient ξτ . The remainder of this work is devoted to obtaining a sufficiently accurate ξτ empirically, i.e. without the assumption of an underlying Gaussian process.

III.

METHODS

In this section, we will describe, in detail, the methodology to evaluate the two types of error in the estimate of ξ already mentioned, averaging and truncation, and discuss a third due to discretization. In Section II, we examined the GK theory in the setting of a single trajectory. It is intuitively clear, however, that averaging over parallel replica trajectories will lead to faster 6

convergence of the estimate ξτ . Hence, we expand Eq. (9) as an average over a number of replicas Nr , N r Z τc Z τ 1 X χr (s) ⊗ χr (s + t) dt ds , ξτ (τc ) = λ Nr τ r=1 0 0

(17)

each employing a time-average based on the ergodic hypothesis and stationarity of the underlying statistics. These replica samples can come from symmetry considerations, e.g. hσ12 (0)σ12 (t)i ∼ hσ23 (0)σ23 (t)i ∼ hσ13 (0)σ13 (t)i for a fluid, or parallel simulations generating trajectories from different initial conditions sampled from the appropriate equilibrium distribution. Since each MD simulation generates samples of χ only at discrete times, we now discretize Eq. (17) to approximate ξτ (τc ) as Nr X Nτ λ X χr (i∆t) ⊗ χr ((i + j)∆t) ∆t wj ξτ (Nc ∆t) ≈ Nr Nτ r=1 i=0 j=0 | {z } Nc X

(18)

Cτ (j∆t)

using a sampling interval of ∆t, quadrature weights wj , and an assumed maximum significant correlation time τc = Nc ∆t [65]. Herein, Nτ is the number of time-samples, which is assumed to be the same for each replica for simplicity. (Refer to Table I for a summary of nomenclature used in this section.) With the introduction of discrete time samples, we recognise that a third source of error is due to discretization. These errors can be due to inaccuracies in computing the trajectories of atoms, which are typically controlled by the size of the time-step and the order of accuracy of the integration method. A full treatment of errors due to numerical time integration and their effects on the flux correlations is beyond the scope of the present work. Herein, we employ the velocity Verlet integrator [36] and select a time-step sufficiently small so as not to affect the computed transport coefficient ξτ . Likewise, we employ the trapezoidal rule (w0 = 1/2 and wj = 1 ∀ j > 0) for the quadrature in Eq. (17) and select a sampling interval ∆t such that the resulting ξτ is insensitive to this choice. To obtain an accurate Green-Kubo estimate we need to know: (a) when the process has reached a steady/equilibrium state [66] after initialization, in effect when to start taking data for ξτ , (b) the maximum significant correlation time, τc , to sufficiently minimize the truncation error in ξτ (τc ), and (c) when the process has accumulated enough samples so that variance in the estimate of the transport coefficient, var ξτ (τc ), is less than a pre-selected 7

∆t sample interval (a multiple of the time-step) Na number of samples in a (block) average Na  Nτ Nb number of blocks Nb = Nτ /Na on a replica Nc number of samples stored, upper limit of correlation, τc = ∆t Nc  τC Nr number of replicas Ns total number of samples of block size Na is Ns = Nr Nb = Nr Nτ /Na Nτ total number of time samples, max time τ = ∆t Nτ TABLE I: Parameters and notation.

error bound.

A.

Detection of a steady state process

In practice, there will be a non-equilibrium transient in the distribution of a phase variable, such as χ, generated by an MD simulation. This is due to the fact that the initial positions, unlike the velocities, are typically not sampled from an equilibrium distribution at temperature T because of the non-linearity in the inter-atomic potential. The resulting practical issues for the GK method is the trade-off between the time it takes for the system to relax to equilibrium versus the extra time it takes for a time average to converge to the ensemble average due to pollution by non-equilibrium values. A more subtle issue is that the Green-Kubo formula (7) is predicated on the fact that the average flux hχi is zero. Hence, it is important to ensure that the system has reached equilibrium in order to avoid errors due to non-equilibrium fluxes. To this end, we take the steady state of the flux distribution, which we define as the state where the probability density function (PDF) of the flux ρ(χ) is unchanging in time, as a convenient and relevant indication of equilibrium [66]. To obtain a representation of the PDF, we use a (normalized) histogram which is a sum of the frequencies of χ(t) over a range of bins. Sampling the process χ(t) gives us a sequence of histograms, each of which contains all the samples of the previous histograms, from which we form a measure of (Cauchy) convergence by comparing consecutive PDFs using a metric [30, 31]. In particular, we use

8

the Hellinger metric kρ(m) (χ) − ρ(m−1) (χ)kH =

"Z q

ρ(m) (χ) −

q

#1/2 2 ρ(m−1) (χ) dχ

(19)

Ω

where Ω is the union of the support of the PDFs ρ(m−1) (χ) and ρ(m) (χ), and m is an index representing the number of samples used to generate the representation. When this metric falls below a pre-set threshold, we assume a steady-state has been reached and the system is in equilibrium [66].

B.

Estimation of the maximum significant correlation time

It is clear from Eq. (9) that the evaluation of ξτ (τc ) requires an a priori bound of the maximum significant correlation time τc . In order to minimize the truncation error in the transport coefficient ξτ (τc ), τc should be considerably greater than the characteristic decay time of the correlation τC , which can be defined as [67] Z ∞ C(t) dt ξ ξ˜ 0 = ≈ . τC = C(0) λ C(0) λ Cτ (0)

(20)

To obtain an approximate value of the characteristic decay time τC , we need a reasonable guess ξ˜ of the transport coefficient ξ and an estimate of the variance of χ, Cτ (0), which can be obtained with a short time average of χ2 once the steady state has been reached. With this in hand, we choose τc = Nc ∆t  τC . We use the value of Nc obtained in this manner to determine an optimal amount of storage for the flux samples needed to calculate an on-the-fly version of the correlation in Eq. (18). More details of our on-the-fly procedure will be given in the subsequent subsections.

C.

Convergence of the estimate with respect to averaging error

In this subsection, we adapt the block averaging technique in [15, 16] and [17, Section 6.4.1] to obtain the estimates of correlation Cτ and the transport coefficient ξτ and their respective variances, and thus quantify the averaging error. In the block averaging technique, an entire time sequence of data is divided into blocks comprised of an equal number of samples and the block averages are calculated. These block 9

averages are then used to calculate the estimate of the mean and the variance in the estimate of the mean. To develop an on-the-fly estimate of averaging error, we calculate ξτ (τc ) in Eq. (18) as a running average of block averages (r,b) ξNa (Nc ∆t)

=λ

Nc X

(r,b)

CNa (j∆t)∆t wj ,

(21)

j=0 (r,b)

where the block-averaged correlations CNa are given by (r,b) CNa (j

(b+1)Na 1 X χr ((i − j)∆t) ⊗ χr (i∆t) . ∆t) = Na i=bN

(22)

a

Here, r refers to a particular replica, b denotes the block index, and Na is the size of the (r,b)

block. To calculate the block-averaged correlation CNa (j ∆t) with Eq. (22), we use the Na samples from the present block b as well as j samples from the preceding blocks. In this way, an equal number of samples are used to obtain the averages for all correlation times j ∆t in order to connect to the GP theory given in Section II. Also, as emphasized in the previous subsection, only a history of Nc flux values are stored for each replica at all times in order (r,b)

to calculate CNa on-the-fly. Figure 1 shows an illustration of the procedure to generate (r,b)

CNa . Consequently, the estimate of the correlation and transport coefficient are simply the averages of the block averages, which are given by Ns Ns 1 X 1 X (r,b) (r,b) Cτ (j∆t) = CNa (j∆t) and ξτ (nc ∆t) = ξNa (nc ∆t) , Ns Ns (r,b)=1

(23)

(r,b)=1

respectively. Here, Ns = Nr Nb is the number of block samples over all replicas, and Nb = Nτ /Na is the number of blocks produced by a single replica (which is assumed to be the same over all replicas for simplicity). As the estimates of the correlation and the transport coefficient are averages of the block averages, we expect them to follow a Gaussian distribution due to the central limit theorem. Thus, the variance in the estimate of ξ can be obtained from the variance of its block averages: (r,b) var ξNa

Ns  2 1 X (r,b) = ξNa − ξτ , Ns

(24)

(r,b)=1

(r,b)

and likewise for var CNa . However, in calculating the block averages we see from Figure 1 (r,b)

(r,b)

that there is a potential for significant interdependence of samples of CNa and ξNa for any 10

χi

-1 N a

0N a

1N a

2N a

3N a

... NbNa

χj N c

+ N a

∆t (r,1)

C

Na

(r,2)

C

Na

(r,b)

FIG. 1: Procedure to calculate the average CNa of a block of size Na for a single replica r = 1. For each sample of the flux χ(i ∆t), a set {χ((i − j) ∆t)χ(i ∆t) : j = 0 . . . Nc } is calculated. This (r,b)

process continues until the sum of Na of these sets has been accumulated and the average CNa

is calculated. The dependence between consecutive block averages, due to the overlap in the flux data used to create them, decreases as Na becomes much larger than Nc .

given block size Na due to the overlap in the underlying flux data. Therefore, the variances in the estimates of C and ξ can not be obtained from the results of Eq. (24) simply by division by number of blocks, unlike the averages, Eq. (23) . To elaborate on this point, given a finite correlation time for the flux τC , the significant dependence between the block averages results from the flux samples near the edges of (r,i) (r,j)

the blocks. However, as the block size increases, the cross-correlations hξNa ξNa i between (r,b)

(r,b)

blocks become less significant and the block averages CNa and ξNa become statistically uncorrelated [17, Section 6.4]. Therefore, for large block sizes, the variance in the estimate of the mean ξτ may be given by [16] (r,b)

(r,b)

var ξNa var ξτ ≈ lim Na →∞ Nr Nb

Na var ξNa = lim , Na →∞ Nr Nτ

(25)

and likewise for var Cτ . Clearly, the variances decrease inversely with the number of replicas (r,b)

(r,b)

Nr . Also, for a valid limit to exist Na var CNa and Na var ξNa must become constant. The existence of such limits is supported by Eqs. (11), (12) and (15), which assume that the fluctuations χ follow a Gaussian process. Ultimately, we estimate the averaging error in the coefficient ξ using the standard deviation, i.e. the square root of the large block limit variance ∆ξτ =

p var ξτ

which represents a confidence interval of approximately 68%. 11

(26)

D.

On-the-fly implementation

The primary tasks of the proposed on-the-fly algorithm are to: (a) determine that the MD trajectories generate flux samples corresponding to a steady state, (b) estimate the maximum significant correlation time, (c) calculate block averages and the variance of the averages, (d) determine when the variances of the block averages of various block sizes are converged in time, (e) estimate the large block limit of the variances, and (f) determine that the integral of correlation is sufficiently converged with respect to averaging and truncation errors. The full algorithm is summarized together with parallel communication details in Appendix B and demonstrated in the next section. The important algorithmic parameters, including the tolerances used for determining convergence, are summarized in Table II.

IV.

RESULTS

To verify that the Green-Kubo algorithm described in Section III and Appendix B is robust and efficient and also to measure the accuracy of the GP theory developed in Section II, we apply the algorithm to compute the transport coefficients of well-known materials with well-tested empirical potentials that include electrostatic and multibody effects. Specifically, we calculate the viscosity of Lennard-Jones (LJ) liquid argon [17] and TIP4P water [32, 33] and the thermal conductivity of Stillinger-Weber parameterizations of solid silicon [34] and gallium-nitride [35]. In each simulation, the dynamics are generated with the velocity Verlet integrator [36] and the Nos´e-Hoover thermostat [37–39] [68] at constant volume (using LAMMPS [40]). The momenta for each replica simulation were initialized with different samples of the Boltzmann distribution. As described in the Appendix B and summarized in Table II, a number of tolerances are needed by the algorithm. For all cases, the tolerance to determine a steady state pdf was set to 0.001, the base multiplier S for the expansion of the smallest block size Na was set to 2, the window size to assess constancy of running averages was set to Nw = 100, and the constant integral fraction ` was set to 50%. The steady variance tolerance var was set to (1)

10% for all cases, leaving just the smallest block size Na

and ξmax to be chosen based on (1)

the particular material. The results of the algorithm are not sensitive to S or Na (m)

the sequence of block sizes Na

since

continually expands with increased running time (refer to

12

pdf tolerance on change of flux PDF (1)

Na

S

initial block size block size multiplier

Nw length of windowed averages to assess steadiness of the block variance var tolerance on steady block variance and equivalent block variances ξmax reasonable guess of the coefficient times a large multiplier ave tolerance on averaging error in the coefficient `

percentage of integral within tolerance ave TABLE II: Tolerances and algorithmic parameters.

Eq. (B2)). Also, the estimate ξτ does not have a significant truncation error as long as ξmax is large enough. Therefore, given a sufficiently large ξmax , only the performance of the (1)

algorithm, and not its accuracy, is affected by specific choices of S, Na

and ξmax .

Finally, we estimate the averaging error on-the-fly using Eq. (26). Additionally, we calculate an a posteriori estimate of the truncation error using exponential fits of the correlation Cτ (t). Specifically, we estimate the truncation error as the integral of the fitted function from τc to ∞. We also performed studies of the sensitivity of the estimates to: (a) system size, (b) time-step, and (c) sampling interval. In the following, the results we report are independent of the values of these parameters but we do not report the details of these studies for brevity.

A.

Viscosity

In this section we calculate the viscosity of two molecular fluids: argon and water. For these isotropic fluids, the viscosity ν computed from Eq. (5) using each of the three shear stresses σ12 , σ13 , σ23 is expected to be equivalent and hence we use each as a replica. Since argon has relatively short characteristic correlation decay time τC , we exploit it as an exhaustive and illustrative test case.

13

1.

Lennard-Jones fluid argon

The widely-used LJ parameterization for Ar is

 kB

= 119.8 K, σ = 3.405 ˚ A, and its mass

is 39.948 amu. With this interaction potential, we simulated the dynamics of 256 atoms at temperature T = 86.4956 K and density ρ = 0.854258 amu/˚ A3 (near the triple point) using a time-step of 1.0 f s for long time stability and accuracy. The pressure in this state was approximately 11 M P a, which, in experiments [41, 42], corresponds to a viscosity of ν ≈ 0.307 mP a-s. Here, we use argon to examine: (a) the convergence of the system to a steady state from initial conditions, (b) the behavior of the variances of the block averages as a function of number of replicas and running time, (c) the convergence of averaging error in the estimate of the transport coefficient with running time, and (d) the correspondence of the simulation results with Gaussian process theory. With a sufficiently small sampling interval ∆t = 8 f s, we examine the transient phase after initializing the velocities (for each replica) from the Boltzmann distribution by monitoring the Hellinger metric given in Eq. (19). Figure 2 shows that the Hellinger metric converges essentially in a uniform manner and appears to follow an inverse square root of the number of samples trend from the outset. Likewise, no indication of a significant nonequilibrium transient in the system was observed in subsequent simulations. During this initial phase the algorithm also estimated the characteristic correlation time τC ≈ 0.3 ps using Eq. (20) and ν¯ ≈ 0.3 mP a-s. Consequently, the maximum correlation length was set to Nc = 1900 via Eq. (B1), where we chose νmax = 50 ν¯. After a statistically steady state has been reached, the sequences of block-averaged correlations for various block sizes necessary to estimate the averaging error ∆ντ from Eq. (26) (r,b)

are computed. Figure 3 shows convergence of Na var νNa with block size Na at Nc for fixed (r,b)

number of replicas Nr = 3. Figure 4 shows the convergence of Na var νNa at Nc for a range of replicas Nr for a single block size. To facilitate comparison, both of these graphs have been normalized by the long-time value of the variance of the Na = 3200 block for Nr = 3 replicas. It can be observed from Figures 3 and 4 that there are significant initial transients before the variance of any given block average can be considered steady. Furthermore, for a few replicas (Nr < 12), there seems to be a systematic trend where the running average takes a long time to become constant. This observation certainly recommends multiple replicas

14

for attaining a steady block variance. Lastly, the inset in Figure 4 shows the error ∆ντ as a function of number of blocks Nb = Nτ /Na and number of replicas Nr . It illustrates that time samples and replica samples can be used interchangeably at a sufficiently large block size Na , so that the error is simply a function of the number of samples Ns = Nr Nb . Figure 5 shows the convergence of the estimate of the viscosity using three (symmetry) replicas with running time τ exceeding 90 ns. The inset demonstrates that errors follow the −1/2

trend line ∼ Nτ

, as expected. Here, the system was run well past the point at which

convergence was detected by the algorithm in order to explore the long-time behavior of the block-averaged samples. Now that the long-time behavior of the error has been established, we use the on-the-fly algorithm to obtain an estimate of the viscosity with a pre-selected error tolerance. With reference to Figure 5, the algorithm converged to ave = 10% in 4.23 ns and 5 % in 10.2 ns, which are comparable with the run time estimates of 2-20 ns in [5]. All subsequent results are reported at a run time of 10.2 ns. The variances of the estimate of the correlation, shown in Figure 6a, display a near perfect agreement with GP theory. Specifically, the data follows the theoretical estimates of the variance at short and long correlation times, Eqs. (11) and (12). These expressions were evaluated using τC 2 estimated from both: (a) direct integration of Eq. (13) and (b) integration of a fit of the correlation function to a series of two exponentials, independently [69]. Now, we examine the resulting variance in the estimate of viscosity as a function of correlation time. Figure 6b shows that the large block limit of the error exists and is achieved with uniform convergence for all calculated correlation times. The trend seen in Figure 6b is analogous to that observed in [16]: (a) for small block sizes the error is underestimated, and (b) for large block sizes, Na → Nτ , the error estimate is inaccurate due to insufficient sampling. Here, and in the subsequent examples, we have verified that the algorithm has calculated a valid large block limit, and, on the other hand, we have omitted plotting block sizes with insufficient sampling. For Ar, the smallest block that converged to within 10% of the larger block estimates was Na = 1600. Additionally, we calculated the average block-to-block correlation Z τc  (r,0) (r,b) C(b) = (CNa (s) − C(s))(CNa (s) − C(s)) ds 0

15

directly as a function of block size Na [70]. We found that the nearest block correlations were

C(1) C(0)

= 0.048, 0.022, 0.013 for Na = 100, 200, 400, respectively, and C(b) for b > 1 was

insignificant. This is direct empirical evidence that the block samples become uncorrelated as the block size increases. Figure 6b also shows that the variance of the integral ντ (τc ) with respect to its upper limit τc is linear after the initial significant correlation. The measured slope of large-block limit of this variance is 0.964 × 10−5 (mP a-s)2 /ns, which compares well to our theoretical estimate 1.01 × 10−5 (mP a-s)2 /ns calculated using Eq. (15). As a further investigation of the linear trend of variance in the estimate of viscosity shown in Figure 6b, we calculated the full covariance matrix embedded in Eq. (A1) : C(s, t) =

D

(r,b)

CNa (s) − C(s)

E  (r,b) . CNa (t) − C(t)

(27)

Clearly, the covariance matrix shown in Figure 7 has a finite bandwidth that becomes constant for large enough |s − t|. Hence, applying Eqs. (24) and (25) results in a linear relationship between the variance in the transport coefficient and the upper limit of the integration τc . Finally, Figure 8 shows the estimates of the correlation and viscosity with error bars corresponding to ave = 5%. It is clear from Figure 8a that the error in the estimate of the correlation is negligible when the error in the integral of correlation is less than 5 %. The error bars shown in Figure 8b demonstrate that the error in the integral of correlation is indeed less than 5% over ` = 50% of its range. The truncation error resulting from the computed Nc was estimated to be negligible at 0.00011%, based on the time-scales resolved by the a posteriori double exponential fit of the correlation. Also, Figure 8b illustrates that there is an optimal range for the upper limit of the integration of correlation Nc given a fixed amount of run time. When Nc is too small, truncation error is large and the transport coefficient is typically underestimated. In this case, however, the averaging error is relatively small. Conversely, when Nc is large, the truncation error is small; however, the averaging error is large. It is also important to notice that the averaging error ∆ντ ∼ Nc /Ns (refer to Eq. (16)); therefore, once a sufficiently large Nc is chosen to control the truncation error, the averaging error can then be reduced by a combination of longer run times and a larger set of replicas.

16

0.1

0.01

0.001

STRESS [MPa]

HELLINGER METRIC

1

30 15 0 -15 -30

σ12 σ13 σ23 0

0.0001

0.05

0.0001

0.1

0.001

0.01

0.1

RUNNING TIME [ns]

FIG. 2: LJ Ar: convergence of ρ(χ) to a steady state as measured by the Hellinger metric. The trend line ∼ 1/t. The inset shows shear stresses as function of time t which at equilibrium should have zero mean. NORMALIZED INTEGRAL VARIANCE

1.5

1

Na=100 Na=200 Na=400 Na=800 Na=1600 Na=3200 Na=6400

0.5

0 0

10

20

30

40

50

60

70

80

90

100

RUNNING TIME [ns]

(r,b)

FIG. 3: LJ Ar: convergence of block variance Na var ξNa (at Nc = 4000 normalized by final value of Na = 3200) in time for Nr = 3. 2.

TIP4P water

Here, we used the on-the-fly algorithm to compute the viscosity of water at T = 298 K and density 0.605 amu/˚ A3 with the TIP4P potential. At this state, the viscosity is expected to be ν ≈ 0.89313 mP a-s from experiments [42]. The TIP4P parameterization for water [32, 33] is essentially a 12-6 LJ pair potential between O atoms plus a long-range Coulomb interaction. The long-range Coulomb forces create significant spatial correlations and are computed using a particle-particle particle-mesh (PPPM) solver [43], with an inner cutoff of 9.0 ˚ A, and an outer cutoff of 9.5 ˚ A [71]. 17

Nr=3 Nr=6 Nr=12 Nr=24 Nr=48

1 100

Nb

NORMALIZED INTEGRAL VARIANCE

1.5

10

0.5 1 1

10 Nr

0 0

10

20

30

40

50

100 60

70

80

90

100

RUNNING TIME [ns]

(r,b)

FIG. 4: LJ Ar: convergence of block variance Na var ξNa for block Na = 3200 ( at Nc = 4000 normalized by final value of Nr = 3) over time for increasing number of replicas Nr Note that (r,b)

the number of block-averaged samples used to calculate Na var νNa increases with Nr . The inset shows the error as a function of number of replicas Nr and number of block samples Nb = Nτ /Na showing approximate ∼ 1/Ns trend for Na = Nc = 2000 samples. In this contour plot, red is assigned to the highest error and pink to the lowest. 0.5

0.3 ERROR [mPa-s]

VISCOSITY [mPa-s]

0.4

0.2

0.1

0.01 0.01 0.1 1

10 100

0 0

10

20

30

40

50

60

70

80

90

100

RUNNING TIME [ns]

FIG. 5: LJ Ar: running average ντ and error for the converged block size Na = 1600. The blue lines are ντ ± 5%. The inset shows the error ∆ντ (red) and a trend line ∼ τ −1/2 (black).

Using this widely-employed potential, we simulated a system of 512 H2 O molecules in a 24.83 ˚ A3 periodic box. We used Nr = 12 replicas each with of a time-step 1.0 f s and a sampling interval of ∆t = 4.0 f s. During the pre-steady state stage, the algorithm calculated an estimate of the characteristic correlation time τC = 0.1 ps, based on ν˜ ≈ 1 mP a-s. This estimate was then used to calculate the maximum significant correlation length Nc = 2500

18

VISCOSITY ERROR SQUARED [mPa-s]2

× 10-6 0.25

16

FREQUENCY

NORMALIZED CORRELATION VARIANCE

18

14 12

Nc = 10 Nc = 20 Nc = 50 Nc = 100

0.2 0.15 0.1 0.05 0

-3 -2 -1 0 1 2 3 4 5 6 NORMALIZED CORRELATION

10 8 6 0

2

4

6

8

10

12

14

0.00016

Na=100 Na=200 Na=400 Na=800 Na=1600 Na=3200 Na=6400 theory

0.00014 0.00012 0.00010 0.00008 0.00006 0.00004 0.00002 0.00000 0

CORRELATION TIME [ps]

2

4

6

8

10

12

14

CORRELATION TIME [ps]

(a)variance of correlation

(b)variance of integral of correlation

(r,b)

(r,b)

Na var νNa Na var CNa normalized by hCτ (0)i2 and (b) both as a function FIG. 6: LJ Ar: (a) Nr Nτ Nr Nτ of correlation time and for a range of block sizes. The large block limits of (a) and (b) yield the normalized variance of the estimate of correlation var Cτ and the variance in the estimate of the viscosity var ντ , respectively. In (a), a posteriori theoretical estimates of the short and long time variances in Cτ are denoted by horizontal black lines. The inset in (a) shows histograms of the (r,b)

steady state distribution of CNa

for selected correlation times. Note that the PDFs for short

correlation times where the average correlation is significant display distinct skew. In (b), the a posteriori theoretical estimate of the error is shown with a diagonal black line.

via Eq. (B1). The algorithm converged to 5% error in 1.55 ns of simulation time (and to 10% error in 0.46 ns) for each of the 12 replicas. The large block limit was detected for Na ≥ 1600. As with Ar, the correspondence with the GP theory is remarkable as shown in Figure 9. Not only do the short and long-time variances of the correlation match the theoretical estimates, but also the measured slope of the variance of the integral, 0.000206 (mP a-s)2 /ns compares well to the theoretical estimate 0.000251 (mP a-s)2 /ns. Interestingly, Figure 10a shows a kink in the correlation function that was not seen in Ar. As can be seen in Figure 10b, the viscosity of TIP4P water is obtained with the error less than 5% over 50% of the evaluated range of the correlation time. The estimate of the viscosity is ν = 1.08 ± 0.045 mP a-s for 5 % averaging error, (and ν = 1.11 ± 0.077 mP a-s for 10% error). Finally, using the same techniques as described in Section IV A 1, the truncation error was found to be negligible at 0.0015%. 19

3

CORRELATION TIME [ps]

2.5 2 1.5 1 0.5 0 0

FIG. 7: LJ Ar: covariance structure

0.5 1 1.5 2 2.5 CORRELATION TIME [ps]

D

3

 E (r,b) (r,b) CNa (i ∆t) − C(i ∆t) CNa (j ∆t) − C(j ∆t) for Na =

100. The contours are evenly spaced with red being the highest value (in the neighborhood of

1

0.3

0.8

0.25 VISCOSITY [mPa-s]

NORMALIZED CORRELATION

i = j = 0) and purple the lowest (which is essentially zero away from the diagonal i = j).

0.6 0.4 0.2 0

0.2 0.15 0.1 0.05

-0.2

0 0

2

4

6

8

10

12

14

0

CORRELATION TIME [ps]

2

4

6

8

10

12

14

CORRELATION TIME [ps]

(a)correlation

(b)integral of correlation

FIG. 8: LJ Ar: estimates of the (a) correlation and (b) viscosity with error bars. B.

Thermal conductivity

In this section, we calculated the thermal conductivity of two crystalline systems: (a) Si and (b) GaN, modelled by the Stillinger-Weber (SW) potential. The SW potential complements the usual pair-wise interaction terms with three-body angularly-dependent terms that stabilize these crystal structures. Silicon has (diamond) cubic symmetry; therefore, 20

VISCOSITY ERROR SQUARED [mPa-s]2

NORMALIZED CORRELATION VARIANCE

4.5

× 10-6

4 Na=250 Na=500 Na=1000 Na=2000 Na=4000 Na=8000 theory

3.5 3 2.5 2 1.5 0

1

2

3

4

5

6

7

8

9

10

0.003

Na=250 Na=500 Na=1000 Na=2000 Na=4000 Na=8000 theory

0.0025 0.002 0.0015 0.001 0.0005 0 0

1

2

CORRELATION TIME [ps]

3

4

5

6

7

8

9

10

CORRELATION TIME [ps]

(a)variance of correlation

(b)variance of integral of correlation

(r,b)

(r,b)

Na var νNa Na var CNa normalized by hCτ (0)i2 and (b) both as a FIG. 9: TIP4P water: (a) Nr Nτ Nr Nτ function of correlation time and for a range of block sizes. The large block limits of (a) and (b) yield normalized var Cτ and var ντ , respectively. The black horizontal lines in (a) are the short, Eq. (11), and long time, Eq. (12), estimates. The black diagonal line in (b) is our theoretical estimate of the slope.

we assume that the auto-correlation of the x, y, and z heat flux components (aligned with the 100, 010, 001 crystal directions) are equivalent. The GaN crystal we studied, however, has wurtzite structure, and so the flux auto-correlations in the x, y, and z directions are not equivalent. Consequently, in this case, we calculate three components of the anisotropic conductivity tensor.

1.

Stillinger-Weber silicon

The Stillinger-Weber parameterization of Si [34] is well-known and commonly used to calculate thermal conductivity. The thermal conductivity of isotopically enriched Si at our selected temperature T = 1000K is κ ≈ 35-50 W/mK which is obtained by extrapolation of experimental data [45, 46]. Given that the chosen temperature is well above the Debye temperature ( ≈ 645 K), no quantum correction is necessary. Given the size-dependence studies in [4, 7, 11], we chose a simulation of 6 × 6 × 6 unit cells (1728 atoms) with a lattice parameter 5.43095 ˚ A. Furthermore, we chose a time-step of 0.5 f s for long-time stability and a sample interval of ∆t = 2.5 f s for accuracy of the

21

1.2 1

1

0.8 0.6

0.6

0.4 0.4

0.2 0

0.2

-0.2 0

1

2

0

VISCOSITY [mPa-s]

0.8

VISCOSITY [mPa-s]

NORMALIZED CORRELATION

1

0.8 0.6 0.4

1.6 1.4 1.2 1 0.8 0.6 0.4 0

0.2

-0.2

0.5 1 1.5 RUNNING TIME [ns]

0 0

1

2

3

4

5

6

7

8

9

10

CORRELATION TIME [ps]

0

2

4

6

8

10

12

14

16

CORRELATION TIME [ps]

(a)correlation

(b)average integral of correlation

FIG. 10: TIP4P water: estimates of the (a) correlation and (b) viscosity with error bars. Inset in (a) is a close-up of the short time correlation. Inset in (b) is the evolution of the estimate with running time.

correlation integral. Using data collected across Nr = 24 parallel replicas, we applied the on-the-fly algorithm (1)

with an averaging error tolerance of ave = 10% and smallest block size of Na

= 10, 000.

In the pre-steady phase, the algorithm estimated the maximum significant correlation time τc = 6.6 ps ( Nc = 132500 ) based on an assumption of κmax ≈ 50 × 65 W/m-K. The algorithm converged to 10 % error at 6.0 ns of simulation time with an estimate of κ = 53.3 ± 5.2 W/m-K. This result agrees well with our previous estimate using Evans’ method [47]. The short and long-time variances in the estimate of the correlation, shown in Figure 11a, demonstrate that the GP theory is accurate in this case. Likewise, Figure 11b shows the same trends and convergence with block size as in the fluid viscosity examples albeit for larger block sizes. Specifically, the characteristic correlation time here is much larger than for the fluid cases, thereby larger block sizes. Here, the large block limit was obtained for Na ≥ 40, 000. With reference to Figure 11b, the slope of the converged block integral variance 0.0623 (W/m-K)2 /ps compares well to our GP estimate of 0.0726 (W/m-K)2 /ps. The inset of Figure 12a shows a slow decay of the correlation after a very fast decay (comparable to the results in [4, Figure 7]). Lastly, Figure 12b shows that the estimated Nc was sufficient to achieve a constant value over greater than 50% of the calculated range of the integral.

22

CONDUCTIVITY ERROR SQUARED [W/m-K]2

NORMALIZED CORRELATION VARIANCE

2.8

× 10-6

2.6 2.4 2.2

Na=10000 Na=20000 Na=40000 Na=80000 Na=160000 theory

2 1.8 1.6 1.4 1.2 1 0.8 0

50

100

150

200

250

300

25

Na=10000 Na=20000 Na=40000 Na=80000 Na=160000 theory

20

15

10

5

0 0

50

CORRELATION TIME [ps]

100

150

200

250

300

CORRELATION TIME [ps]

(a)variance of correlation

(b)variance of integral of correlation

(r,b)

(r,b)

Na var νNa Na var CNa normalized by hCτ (0)i2 and (b) both as a function FIG. 11: SW Si: (a) Nr Nτ Nr Nτ of correlation time and for a range of block sizes. The large block limits of (a) and (b) yield normalized var Cτ and var ντ , respectively. The horizontal black lines in (a) and the diagonal black line in (b) are estimates from Gaussian process theory. 60

0.8 0.6

0.6

0.4 0.4

0.2 0

0.2

-0.2 0

1

2

0

50 40

CONDUCTIVITY

1

0.8

CONDUCTIVITY [W/m-K]

NORMALIZED CORRELATION

1

30 20

70 60 50 40 30 20 10 0

10

1

2

3

4

5

6

RUNNING TIME [ns] -0.2

0 0

50

100

150

200

250

300

0

CORRELATION TIME [ps]

50

100

150

200

250

300

CORRELATION TIME [ps]

(a)correlation

(b)integral of correlation

FIG. 12: SW Si: estimates of (a) correlation and (b) thermal conductivity with error bars. The inset in (a) is a close-up of the short-time correlation decay. The inset in (b) is the evolution of the estimate and error with simulation time. 2.

Stillinger-Weber gallium-nitride

We use the Stillinger-Weber parameterization of GaN in [35], which has been employed in other studies of thermal conductivity, e.g. [12, 49, 57]. The selected temperature T =

23

800 K [72] is well above the Debye temperature ≈ 590 K and below the melt temperature. We expect the thermal conductivity to be κ ≈ 67-85 W/mK based on thin-film and bulk experimental data [50–55] extrapolated to T = 800K using the fits in [56]. The x, y, and z axes of the system are aligned with the c-axis [0001], [¯1100], and [11¯20] directions of the crystal, respectively. In addition, we used lattice constants a = 3.19 ˚ A and c = 5.20 ˚ A, and an internal displacement u = 0.375. As in the previous example, we chose a system of 6 × 6 × 6 cells (1728 atoms), a time-step of 0.5 f s, and a sample interval of ∆t = 2.5 f s. Using Nr = 36 replicas with the on-the-fly algorithm, the characteristic decay time was estimated as τC = 0.16 ps in the initial phase based on 100 W/mK. However, given the results in [12, 49], this value is a clear underestimation of the mean decay of the correlation. In this case, estimation of a reasonable length of correlation Nc essentially failed. Instead, we use an (arbitrary) large multiplier to obtain τc = 400 ps (Nc = 160000 samples) and minimize truncation error. With this choice, the simulation converged to a 10 % accurate estimate at 7.2 ns of simulation time: κxx = 103.3 ± 7.9 W/m-K, κyy = 95.8 ± 7.7 W/m-K, and κzz = 76.9 ± 5.9 W/m-K. Now examining the details of the converged results, Figure 13a shows a highly oscillatory variance in the estimate of the correlation. Nevertheless, the correspondence with GP theory is good. Here, we calculated τC 2 using a direct quadrature of Eq. (13). Despite the oscillations in the variance of the estimate of the correlation, the trends for the variance in the estimate of the thermal conductivity, shown in Figure 13b, are similar to all the other cases. Again, the measured slope 0.149 (W/m-K)2 /ps compares well to the theoretical estimate 0.165 (W/m-K)2 /ps. The estimate of correlation in Figure 14a displays high frequency oscillations, with both positive and negative values, similar to [7, Fig. 7a] for sphalerite, as well as the results in [12, 49]. Despite this feature, the integral shown in Figure 14b is well-behaved and is constant over more than half of the calculated range, based on our choice of Nc .

V.

DISCUSSION

The procedure of mixing time and replica samples was first proposed by Zwanzig and Ailwadi [18] to obtain an accurate correlation function. In [18], they anticipated that the √ limiting convergence rate for the error in estimate of the correlation function is ∼ 1/ Ns . In 24

× 10-6

18

70

Na=5000 Na=10000 Na=20000 Na=40000 Na=80000 Na=160000 Na=320000 theory

16 14 12

Na=5000 Na=10000 Na=20000 Na=40000 Na=80000 Na=160000 Na=320000 theory

60 50 S [W/m-K]2

NORMALIZED CORRELATION VARIANCE

20

10 8 6

40 30 20

4 10

2 0

0 0

50

100

150

200

250

300

350

400

0

50

CORRELATION TIME [ps]

100

150

200

250

300

350

400

CORRELATION TIME [ps]

(a)variance of correlation

(b)variance of integral of correlation

(r,b)

(r,b)

Na var νNa Na var CNa normalized by hCτ (0)i2 and (b) both as a function FIG. 13: SW GaN: (a) Nr Nτ Nr Nτ of correlation time and for a range of block sizes. The large block limits of (a) and (b) yield normalized var Cτ and var ντ , respectively. The horizontal black lines in (a) and the diagonal black line in (b) are estimates from Gaussian process theory. 120

0.4 0.2

0

1

2

0 -0.2 -0.4

100 80 CONDUCTIVITY

1 0.5 0 -0.5 -1

0.6

CONDUCTIVITY [W/m-K]

NORMALIZED CORRELATION

1 0.8

60 40 20

160 140 120 100 80 60 0

-0.6

1

2

3

4

5

6

RUNNING TIME [ns] -0.8

0 0

50

100

150

200

250

300

350

400

CORRELATION TIME [ps]

0

50

100

150

200

250

300

350

400

CORRELATION TIME [ps]

(a)correlation

(b)integral of correlation

FIG. 14: SW GaN: estimates of (a) correlation and (b) thermal conductivity with error bars for the c-axis direction. The inset in (a) is a close-up of the short-time correlation oscillations and decay envelope. The inset in (b) is the evolution of the estimate and error with simulation time.

addition to demonstrating this behavior, our algorithm shows how to handle the ensemble of replicas simultaneously and efficiently in parallel. Our results clearly imply that employing more replica samples relative to time samples is advantageous. Assuming the flux fluctuations follow a Gaussian process, Zwanzig and Aliwadi [18] also 25

established that the variance in the estimate of the correlation, at short and long correlation times, is proportional to the ratio of the correlation time τC 2 to the averaging time τ . In this work, we extended the theory to obtain the variance in the estimate of the transport coefficient and found that it is bounded by the ratio of upper limit of the integral of correlation τc to the total run time τ . When the GP theory is accurate, it gives nearly all the characteristics of the error in the estimates of the correlation and its integral. In all the four test cases we examined, the empirical results adhere well to the assumptions of an underlying Gaussian process. In fact, our estimate of the error in the integral of correlation proved to be a tight and conservative bound, slightly over-estimating the empirically calculated error in all cases. Consequently, the Gaussian process theory is likely to be sufficient to estimate the error in the transport coefficient without the block-averaging technique; however, the only route we see to prove this assumption for a particular case is via the algorithm. Also, in applying the GP theory to estimate the error in the transport coefficient a sufficiently accurate mean must be calculated, which provides additional motivation for the empirical algorithm. The proposed algorithm computed the averaging error directly and in parallel using the block averaging technique [15, 16] applied to the integral of correlation. Also, it computed the maximum significant correlation time to set an appropriate upper limit for the integral of correlation and thus control the truncation error. Our results demonstrate that taking this limit to be many multiples of the estimated characteristic correlation decay time is necessary to make the truncation error negligible. However, increasing the upper limit leads to the integration of the noise in the correlation (where the mean is close to zero) and therefore the averaging error increases proportionally. Hence, there is an optimal upper limit of integration to obtain the estimate of the transport coefficient ξτ (τc ) when both averaging and truncation errors are taken into consideration. The algorithm can certainly be improved and this is left for future work; however, the need for a rigorous and efficient procedure of the type proposed is clear. For example, the characteristic correlation time of the equilibrium heat flux in the insulators we studied will increase as the temperature approaches room temperature. Hence, accurate GK estimates of thermal conductivity of semiconductors and insulators near room temperature (with the usual quantum correction) are particularly expensive, especially if calculated na¨ıvely. The proposed algorithm with its error-based termination enables efficient parallel calculation of 26

these important material properties.

Acknowledgments

This work was supported by the U.S. Department of Energy Office of Science through the Applied Mathematics program in the Office of Advanced Scientific Computing Research (ASCR) and the Laboratory Directed Research and Development (LDRD) program at Sandia National Laboratories. Sandia is a multiprogram laboratory managed and operated by Sandia Corporation, a wholly owned subsidiary of Lockheed Martin Corporation, for the U.S. Department of Energy’s National Nuclear Security Administration under contract No. DE-AC04-94AL85000. We would also like to acknowledge the help of: P. Saxe, D. Rigby and others at Materials Design (who provided the original on-the-fly correlation code for LAMMPS); S. Plimpton at Sandia (for insight into LAMMPS); M. Salloum and K. Sargsyan (for suggestions of metrics to measure convergence of PDFs); H. Adalsteinsson, R. Berry, B. Debusschere, H. Najm, J. Templeton at Sandia (for feedback on the preliminary results and theory); and S. Sankararaman at Harvard (for discussions of statistics) , P. Howell at Siemens (for insights into GK practice).

[1] M. S. Green, J. Chem. Phys. 22, 398 (1954). [2] R. Kubo, J. Phys. Soc. Jpn. 12, 570 (1957). [3] R. Kubo, M. Yokota, and S. Nakajima, J. Phys. Soc. Jpn. 12, 1203 (1957). [4] P. K. Schelling, S. R. Phillpot, and P. Keblinski, Phys. Rev. B 65, 144306 (2002). [5] H. Kaburaki, J. Li, S. Yip, and H. Kimizuka, J. Appl. Phys. 102, 043514 (2007). [6] A. S. Henry and G. Chen, Journal of Computational and Theoretical Nanoscience 5, 141 (2008). [7] J. Che, T. C ¸ a˘ gin, W. Deng, and W. A. Goddard III, J. Chem. Phys. 113, 6888 (2000). [8] N. Galamba, C. A. N. de Castro, and J. F. Ely, J. Chem. Phys. 126, 204511 (2007). [9] N. V. Brilliantov and T. P¨ oschel, Chaos 15, 026108 (2005). [10] G. A. Pavliotis, IMA J. Appl. Math. 75, 951 (2010).

27

[11] L. Sun and J. Y. Murthy, Appl. Phys. Lett. 89, 171919 (2006). [12] X. W. Zhou, S. Aubry, R. E. Jones, A. Greenstein, and P. K. Schelling, Phys. Rev. B 79, 115201 (2009). [13] T. Kawamura, Y. Kangawa, and K. Kakimoto, J. Cryst. Growth 284, 197 (2005). [14] X. W. Zhou and R. E. Jones, Modelling and Simulation in Materials Science and Engineering 19, 025004 (2011). [15] R. Friedberg and J. E. Cameron, J. Chem. Phys. 52, 6049 (1970). [16] H. Flyvbjerg and H. G. Petersen, 91, 461 (1989). [17] M. P. Allen and D. J. Tildesley, Computer Simulation of Liquids (Oxford University Press, New York, 1987). [18] R. Zwanzig and N. K. Ailawadi, Phys. Rev. 182, 280 (1969). [19] R. J. E. Clausius, Philosophical Magazine 40, 122 (1870). [20] J. C. Maxwell, Transactions of the Royal Society Edinborough XXVI, 1 (1870). [21] J. H. Irving and J. G. Kirkwood, J. Chem. Phys. 18, 817 (1950). [22] K. K. Mandadapu, R. E. Jones, and P. Papadopoulos, Phys. Rev. E 80, 047702 (2009). [23] L. Onsager, Phys. Rev. 37, 405 (1931). [24] L. Onsager, Phys. Rev. 38, 2265 (1931). [25] R. Zwanzig, J. Chem. Phys. 40, 2527 (1964). [26] D. J. Evans and G. P. Morriss, Statistical Mechanics of Non-equilibrium Liquids (ANU E Press, 2007), 2nd ed. [27] J. P. Sethna, Entropy, Order Parameters, and Complexity (Oxford Univ. Press, 2010). [28] D. Frenkel, in Proceedings of the International School of Physics “Enrico Fermi”, Course LXXV (1980). [29] D. Frenkel and B. Smit, Understanding Molecular Simulation (Academic Press, San Diego, 2001). [30] P. Diaconis and S. L. Zabell, Journal of the American Statistical Association 77, 822 (1982). [31] A. L. Gibbs and F. E. Su, International Statistical Review 70, 419 (2002). [32] W. L. Jorgensen, J. Chandrasekhar, J. D. Madura, R. W. Impey, and M. L. Klein, J. Chem. Phys. 79, 926 (1983). [33] H. W. Horn, W. C. Swope, J. W. Pitera, J. D. Madura, T. J. Dick, G. L. Hura, and T. HeadGordon, J. Chem. Phys. 120, 9665 (2004).

28

[34] F. H. Stillinger and T. A. Weber, Phys. Rev. B 31, 5262 (1985). [35] A. Bere and A. Serra, Philos. Mag 86, 2152 (2006). [36] L. Verlet, Phys. Rev. 159, 98 (1967). [37] S. Nos´e, J. Chem. Phys. 81, 511 (1984). [38] W. G. Hoover, Phys. Rev. A 31, 1695 (1985). [39] G. Martyna, M. Klein, and M. Tuckerman, J. Chem. Phys. 97, 2635 (1992). [40] Sandia National Laboratories, Large-scale Atomic/Molecular Massively Parallel Simulator (LAMMPS), URL http://lammps.sandia.gov. [41] P. Malbrunot, A. Boyer, E. Charles, and H. Abachi, Phys. Rev. A 27, 1523 (1983). [42] National Institute of Testing and Standards, Thermophysical Properties of Fluid Systems, URL http://webbook.nist.gov/chemistry/fluid/. [43] R. Hockney and J. Eastwood, Computer Simulation using Particles (Taylor & Francis Group, 1988). [44] A. P. Thompson, S. J. Plimpton, and W. Mattson, J. Chem. Phys. 131, 154107 (2009). [45] W. S. Capinski, H. J. Maris, E. Bauser, I. Silier, M. Asen-Palmer, T. Ruf, M. Cardona, and E. Gmelin, Applied Physics Letters 71, 2109 (1997). [46] T. Ruf, R. Henn, M. Asen-Palmer, E. Gmelin, M. Cardona, H.-J. Pohl, G. Devyatych, and P. Sennikov, Solid State Communications 115, 243 (2000). [47] K. K. Mandadapu, R. E. Jones, and P. Papadopoulos, J. Chem. Phys. 130, 204106 (2009). [48] D. MacGowan and D. J. Evans, Phys. Rev. A 34, 2133 (1986). [49] K. K. Mandadapu, R. E. Jones, and P. Papadopoulos, J. Chem. Phys. 133, 034122 (2010). [50] G.A. and Slack, Journal of Physics and Chemistry of Solids 34, 321 (1973). [51] E. Sichel and J. Pankove, J. Phys. Chem. Solids 38, 330 (1977). [52] C. Luo, D. Clarke, and J. Dryden, Journal of Electronic Materials 30, 138 (2001). [53] A. Je˙zowski, P. Stachowiak, T. Plackowski, T. Suski, S. Krukowski, M. Bo´ckowski, I. Grzegory, B. Danilchenko, and T. Paszkiewicz, Solid State Communications 128, 69 (2003). [54] W. Liu and A. A. Balandin, Appl. Phys. Lett. 85, 5230 (2004). [55] C. Mion, J. F. Muth, E. A. Preble, and D. Hanser, 89, 092123 (2006). [56] S. Vitanov, Ph.D. thesis, TU-Wien (2010). [57] X. W. Zhou, R. E. Jones, and S. Aubry, Phys. Rev. B 81, 155321 (2010). [58] A. N. Kolmogorov and S. V. Fomin, Introductory Real Analysis (Dover, New York, 1975).

29

[59] D. E. Knuth, The Art of Computer Programming, vol. 2: Seminumerical Algorithms (AddisonWesley, Boston, 1998), 3rd ed. [60] B. P. Welford, Technometrics 4, 419 (1962). [61] S. G. Volz and G. Chen, Phys. Rev. B 61, 2651 (2000). [62] D. J. Evans and B. L. Holian, J. Chem. Phys. 83, 4069 (1985). [63] Herein, we use an equally-weighted split of the bond energy to all the atoms involved in any particular bond. Also, no example we consider in the Results section requires the Coulomb interaction energy as part of the heat flux and so we omit it here. [64] Since we consider only solid mixtures, the thermal conductivity formula omits mass diffusion terms. [65] It should be noted that there are at least three approaches to obtain the value of the integral of correlation related to the transport coefficient: direct integration [4], fit frequency dependent thermal conductivity to single relaxation time approximation [61], fit auto-correlation function directly to a small set of relaxation times [5, 7]. We chose direct quadrature method since it is amenable to on-the-fly estimation of estimate and error (unlike Fourier transform methods) and does not require fitting and extrapolation. [66] In general, a steady state is not an equilibrium state. In the present context, however, the steady state generated by the Nos´e-Hoover dynamics corresponds to the canonical phase space distribution [38], [26, Section 5.2]. Hence, once the transient relaxation from initial conditions is over, the steady state is equivalent to equilibrium state. [67] Note that hχ2 i is clearly dependent on system size, the product V hχ2 i should converge with increasing system size. [68] Note the use of a thermostat differs from the common practice of using constant energy dynamics. As proved in [62] the differences between correlations taken with Hamiltonian and Nos´e-Hoover dynamics are of the order of 1/N where N is the number of particles. Moreover, the use of thermostatted dynamics avoids the difficulty in generating initial conditions corresponding to a given temperature. [69] Although not shown here, the variance in the estimate of the auto-correlation hχi (0)χi (t)i at long correlation times also agreed well with the variance in the estimate of the cross correlation of flux hχi (0)χj (t)i, i 6= j for all t > 0. [70] Since the average block-to-block correlation displays a decreasing trend, the individual compo-

30

(r,0)

(r,b)

nents of this average h(CNa (nc ∆t) − C(nc ∆t)) · (CNa (nc ∆t) − C(nc ∆t))i must also decrease with increased block size. We expect the trend for the individual components of the block-toblock correlation to have a maximum at a correlation time (a) larger than zero, where there is no overlap of block data used to calculate the estimate, but (b) lesser than the correlation time where the estimate is very close to zero. [71] Given the long-range electrostatic interactions and periodic boundary conditions, we used the formulation for the system-averaged stress found in [44] instead of Eq. (2). [72] The temperature T = 800 K was chosen to correspond to the inconclusive results given in [12] using the direct method.

Appendix A: Gaussian process theory for the variance in the estimate of the transport coefficient

Following the derivations in [18, 28] and [29, Appendix D] , the estimate of the transport Rτ Rτ coefficient ξ(τc ), Eq.(9), is given by ξτ (τc ) = λ 0 c Cτ (s)ds where Cτ (s) = τ1 0 χ(t)χ(s + t)dt, Eq.(8), is the estimate of the correlation of the flux χ from a time-average over an interval of length τ . The variance in the estimate of the transport coefficient is then given by

 var ξτ (τc ) = [ξτ (τc ) − hξτ (τc )i]2 Z τc  Z τc  2 0 0 0 00 00 00 =λ Cτ (s ) − hCτ (s )i ds Cτ (s ) − hCτ (s )i ds 0 0 Z τc Z τc D E 0 2 ds00 [Cτ (s0 ) − hCτ (s0 )i] [Cτ (s00 ) − hCτ (s00 )i] . ds =λ 0

(A1)

0

We can see from Section IV that var ξτ (τc ) varies linearly with the upper limit of integration τc  τC . In this Appendix, we obtain a theoretical value of the slope of this linear trend by assuming that the fluctuations in χ can be modeled by a Gaussian Process. With Eq. (A1), the variance in the transport coefficient coefficient can be written as Z τc Z τc Z τc Z τc 0 00 0 00 2 0 2 var ξτ (τc ) =λ ds ds hCτ (s )Cτ (s )i − λ ds ds00 hCτ (s0 )i hCτ (s00 )i 0 0 0 Z τ Z τ Z τc 0 Z τc 1 =λ2 2 dt0 dt00 hχ(s0 + t0 )χ(t0 )χ(s00 + t00 )χ(t00 )i ds0 ds00 (A2) τ 0 0 0 0 Z τc Z τc Z τ Z τ 2 1 0 00 0 −λ 2 ds ds dt dt00 hχ(s0 + t0 )χ(t0 )i hχ(s00 + t00 )χ(t00 )i . τ 0 0 0 0

31

Now, assuming that the fluctuations χ follow a Gaussian process, Eq.(A2) reduces to Z Z Z Z λ2 τc 0 τc 00 τ 0 τ 00 dt hχ(s0 + t0 )χ(s00 + t00 )i hχ(t0 )χ(t00 )i dt ds var ξτ (τc ) = 2 ds τ 0 0 (A3) Z τc 0 Z τ 0 Z τ 2 Z τc λ 0 00 0 0 00 00 00 0 00 0 + 2 dt dt hχ(s + t )χ(t )i hχ(s + t )χ(t )i , ds ds τ 0 0 0 0 see [18] and [29, Eq. D.2.2]. Using change of variables x = t0 −t00 , y = s0 −s00 and stationarity of the flux statistics, C(t) ≡ hχ(t + t0 )χ(t0 )i = hχ(t)χ(0)i, Eq. (A3) simplifies to Z t0 Z τ Z 0 Z λ2 τ c 0 s 0 dx hχ(x + y)χ(0)i hχ(x)χ(0)i dt dy ds var ξτ (τc ) = 2 τ 0 t0 −τ 0 s0 −τc Z 0 Z Z Z λ2 τc 0 τc 00 τ 0 t + 2 dx hχ(x + s0 )χ(0)i hχ(x − s00 )χ(0)i dt ds ds τ 0 0 t −τ 0 0 Z s0 Z τ Z t0 2 Z τc λ = 2 ds0 dy dt0 dx C(x + y)C(x) τ 0 s0 −τc 0 t0 −τ Z Z 0 Z Z λ2 τc 0 τc 00 τ 0 t dt dx C(x + s0 )C(x − s00 ). ds + 2 ds τ 0 0 t0 −τ 0

(A4)

Taking the derivative of Eq. (A4) with respect to τc yields Z τ Z t0 Z Z Z Z 0 d λ 2 τc 0 τ 0 t λ2 τc 0 dt dx C(x + y)C(x) + 2 var(ξτ (τc )) = 2 dy ds dt dx C(x + s0 − τc )C(x) dτc τ 0 τ 0 0 0 t −τ 0 0 t −τ Z t0 Z t0 Z τ Z τ 2 Z τc 2 Z τc λ λ dx C(x + τc )C(x − y) + 2 dx C(x + y)C(x − τc ). dt0 dt0 dy dy + 2 τ 0 τ 0 t −τ t0 −τ 0 0 0 (A5) Clearly, the slope goes to zero as τc → 0. Now, we examine the limit as τc becomes large relative to the characteristic decay time of the correlation τC . The first two terms on the right-hand side of Eq. (A5) reduce to Z Z τ Z t0 Z Z τ Z t0 λ 2 τc λ2 0 0 0 dt dx C(x + y)C(x) + 2 dt dx C(x + y)C(x) dy dy τ2 0 τ −τc 0 t0 −τ 0 t0 −τ (A6) Z t0 Z t0 Z Z τ Z τ 2 Z τc λ2 τc λ = 2 dy dt0 dx C(x + y)C(x) = 2 dy dt0 dx C(y − x)C(x). τ −τc τ 0 t −τ t0 −τ 0 −τc 0 where the last step was accomplished via switching the limits of y and using the symmetry C(−t) = C(t) of the correlation. Now, applying Fubini’s theorem [58, Section 35.3] and interchanging the two inner integrals on t0 and x, Eq. (A6) can be further reduced to Z Z Z τ Z t0 Z 0 λ2 τc λ2 τ c 0 dy dt dx C(y − x)C(x) = 2 dy dx (τ + x)C(y − x)C(x) τ 2 −τc τ −τc t0 −τ 0 −τ (A7) Z Z τ Z Z τ λ 2 τc λ2 τc + 2 dy dx (τ − x)C(y − x)C(x) = dy dx C(y − x)C(x). τ −τc τ −τc 0 −τ 32

Then, Eq. (A7) can be simplified by assuming that τc and τ are large relative to the characteristic decay time of the correlation τC and applying Fubini’s theorem to the integral of the convolution integral to yield λ2 τ

τc

τ

λ2 dy dx C(y − x)C(x) = τ −τc −τ

Z

Z

Z

τc

Z

τ

dy C(y) −τc

dx C(x) = −τ

4ξ 2 . τ

(A8)

Finally, the third and fourth terms on the right-hand side of Eq. (A5) can be reduced to λ2 2 2 τ

Z

τc

Z

0

τ

dx (τ − x) [C(x − τc )C(x + y) + C(x + τc )C(x − y)]

dy

(A9)

0

by again applying Fubini’s theorem for the inner integrals. By inspection, the second term in Eq. (A9) vanishes since C(x + τc ) vanishes as τc becomes large. Similarly, the first term goes to zero since the product C(x − τc )C(x + y) vanishes as τc becomes large, since the correlations C(x − τc ) and C(x + y) as a function of x are separated by at least τc . From Eq. (A8) and (A9), the derivative of var ξτ (τc ) with respect to τc is d 4ξ 2 , var ξτ (τc ) = dτc τ

(A10)

in the limit τc → ∞, where τ is the total running time. Therefore, for τc  τC , var ξτ (τc ) ∼

4ξ 2 τc . τ

(A11)

Given that the initial slope is 0 and the variance is positive, there is a negative offset from this linear trend which implies var ξτ (τc ) <

4ξ 2 τ. τ c

Appendix B: On-the-fly, adaptive algorithm to calculate the transport coefficient

The basis of the algorithm has been introduced in Section III and the relevant parameters have been defined in Table I and II. In this appendix, we describe the algorithmic details that allow us to estimate the transport coefficient and its error in an on-the-fly and parallel fashion. The primary responsibilities of our algorithm are to: 1. determine that the flux samples generated by MD correspond to a steady state, 2. estimate the maximum significant correlation time, 3. calculate block averages and the variance of the averages, 4. determine when the variances of the averages of various blocks are converged in time, 33

5. estimate the large block limit of the variances, and 6. determine that the integral of correlation is sufficiently converged with respect to averaging and truncation errors. A synopsis of the overall structure of the algorithm is given in Algorithm B.1. Once the trajectories have been initiated on each parallel replica of the material system, the first task is to determine when a steady state has been reached (see Section III A). Once the Hellinger metric kρ(m) (χ) − ρ(m−1) (χ)kH across all replicas falls below the threshold pdf , the algorithm proceeds to estimate the maximum correlation time τc = Nc ∆t, where Nc is defined by Nc =

ξmax λ Cτ (0)∆t

(B1)

and determines the amount of memory needed to calculate the correlation of χ, as described in Section III B. Here, ξmax is a reasonable guess for the transport coefficient ξ times a large multiplier, chosen to control truncation error and limit the sensitivity of the algorithm on the particular a priori guess of ξ. To estimate the averaging error, the algorithm computes the variance in ξτ by dividing the flux data from each replica into blocks. It then calculates the block averages according to Eqs. (22) and (21). To obtain an accurate error estimate, the algorithm needs to efficiently compute the large block limit of the block variance, as described in Section III C. To estimate this limit, each replica r forms sequences of averages for increasing block sizes (1)

(2)

(3)

Na , Na , Na , . . . = 20 Na , 21 Na , 22 Na , . . . using the recursive tree structure (r,2b−1)

ξ2Na

=

 1  (r,b) (r,b−1) , ξNa + ξNa 2

(B2)

which expands to a new block size given enough samples . Here, only two choices need to (1)

be made: the initial block size Na , and the block size multiplier S (in this case S = 2). To (r,b)

(r,b)

update the average ξτ and variance var ξNa given a new block average sample ξNa efficiently, the algorithm employs the well-known Welford-Knuth recurrence relation [59, 60]: Ns = Ns + 1 ∆ξ (r,b) where ∆ξ = ξNa − ξτ Ns   Mξ (r,b) = where Mξ = Mξ + ∆ξ ξNa − ξτ (Ns − 1)

ξτ = ξτ + (r,b)

var ξNa

(B3)

Note that this recurrence relation allows for simple restart from a previous simulation. 34

With a method to calculate the averages and variances for increasing block sizes, the algorithm now evaluates a sequence of convergence tests summarized in Algorithm B.2. The first test is to determine when a given block variance is steady using a tolerance var and a (r,b)

windowed average over the last Nw samples of the running average of var ξNa . This same tolerance is used to determine the large block limit of the variance from Eq. (25) in the second test. If this test is passed, the algorithm determines whether the integral ξτ has a well-defined upper limit. Specifically, the third test compares the coefficient ξτ (Nc ∆t) for the maximum computed upper limit Nc to those obtained from a range smaller upper limits [`Nc , Nc ] using a tolerance ave . Finally, the current coefficient ξτ is compared to its recent preceding values to determine whether the estimate is steady in running time. for each replica r = 1, . . . , Nr initialize momenta {pα } by sampling from the Boltzmann distribution start dynamics to generate {xα (i ∆τ ), pα (i ∆τ )} for i = 0, 1, . . . calculate histogram ρ(m) (χ) of the flux χ & convergence metric kρ(m) (χ) − ρ(m−1) (χ)kH from Eq. (19) if kρ(m) (χ) − ρ(m−1) (χ)kH < pdf then estimate maximum correlation length Nc using Eq. (B1) sample flux χ every ∆t and store χ((i + j) ∆t)), j = 0, . . . , Nc , (r,b)

(r,b)

& calculate correlation CNa and integral ξNa from Eq. (18) (b)

increase maximum block size Na via Eq. (B2) given sufficient block average samples   (r,b)   master replica receive ξNa       (r,b)   update average ξτ & variance var ξNa via Eq. (B3) every Na  sample steps   if (convergence tests (Algorithm B.2) pass ) then terminate.       (r,b)  other replicas send ξN to master a Algorithm B.1: Parallel and adaptive algorithm for estimating the transport coefficient ξ. Communication particular to a replica is highlighted in italics. Note that master replica (r = 1) updates the estimate ξτ using not only its block average samples but also those of all the other replicas. Also, ∆τ is the time-step of the integrator and the sampling interval ∆t is a multiple of ∆τ .

35

(r,b)

determine if the block variance var ξNa is steady in time (r,b)

(r,b)

if (| var ξNa − var ξNa | < var (r,b)

(r,b)

where var ξNa is the average of var ξNa over the last Nw samples then (r,b)

calculate the block limit limNa →∞ Na var ξNa as in Eq. (25) (m)

if (|Na

var ξ

(r,b) (m)

Na

(m+1)

− Na

var ξ

(r,b) (m+1)

Na

| < var for two consecutive block sizes ) then

test whether integral ξτ is constant with respect to its upper limit if (|ξτ (i ∆t) ± ∆ξτ (i ∆t) − ξτ (Nc ∆t)| < ave ) for every i ∈ (`Nc , Nc ) ) then test whether coefficient ξτ is constant in running time τ if (|ξ(k ∆t) − ξτ | < ave for k ∈ (Nτ − Nw , Nτ )) then converged Algorithm B.2: Sequence of convergence tests. Refer to Tables I and II for a summary of notation.

36