We present a method, based on the classical Green-Kubo theory of linear response, to compute the heat conductivity of extended systems, leveraging energy-density, rather than energy-current, fluctuations, thus avoiding the need to devise an analytical expression for the macroscopic energy flux. The implementation of this method requires the evaluation of the long-wavelength and low-frequency limits of a suitably defined correlation function, which we perform using a combination of recently-introduced cepstral-analysis and Bayesian extrapolation techniques. Our methodology is demonstrated against standard current-based Green-Kubo results for liquid argon and water, and solid amorphous Silica, and compared with a recently proposed similar technique, which utilizes mass-density, instead of energy-density, fluctuations.

Topics
Heat transfer, Linear response, Signal processing, Bayesian inference, Nonequilibrium statistical mechanics

Heat transport is one of the most fundamental off-equilibrium processes occurring in nature, whose understanding sinks its roots into the classic works of Onsager in the 1930s1 and of Green and Kubo (GK)2,3 in the 1950s. In spite of the strong theoretical foundation that these works provide to the numerical simulation of transport phenomena in condensed matter, until recently widespread misconceptions had hindered the application of standard equilibrium simulation techniques to the numerical modeling of heat transport when atomic forces are derived from quantum-mechanical ab initio methods.4 While the crux of the difficulties was originally ascribed to the impossibility of assigning a well-defined value to the atomic energies that enter the definition of the macroscopic energy current—when atomic forces are derived from quantum mechanical methods—it was later realized that a unique partition of the total energy into local contributions is not possible even when atomic forces are derived from classical potentials. It was then demonstrated that, although the macroscopic current corresponding to different local partitions of the total energy depends on the specific partition being considered, the value of the resulting heat conductivity does not. This independence of the heat conductivity of the local representation of the total energies was dubbed gauge invariance of transport coefficients5–8 and has far-reaching consequences in other branches of transport theory as well, such as, e.g., charge transport in ionic conductors.9,10

Although gauge invariance provides a general framework to compute heat transport coefficients for general atomic forces derived from either electronic-structure theory or arbitrary—possibly machine-trained—inter-atomic potentials, the problem still remains that the expression for the macroscopic energy current for general classical force fields or quantum density functionals may be difficult to compute analytically and to implement numerically. For this reason, it would be desirable to devise a method to compute the heat conductivity not relying on the macroscopic energy current, which is the standard ingredient of the GK formula for the transport coefficient, but rather on the energy density, which—albeit not uniquely defined because of gauge freedom—is easier to compute. A successful attempt in this direction has been recently made by Cheng and Frenkel (CF)11 who have shown how to evaluate the heat conductivity of simple fluids, κ, by analyzing the wave-vector dependence of the widths of the Rayleigh (R) and Brillouin (B) peaks, as estimated from equilibrium molecular dynamics (EMD) simulations:12–14 
(1)
(2)
where DT=κmncp and Dℓ are the heat and longitudinal diffusivities, respectively, RLP=cpcv−1 is the (Landau-Planczek) ratio between the intensities of the Rayleigh and Brillouin peaks, which vanishes for incompressible fluids, cp and cv are the isobaric and isochoric specific heats, and m and n the molecular mass and average number density, respectively. The smallness of the Landau–Planczek ratio in the condensed phases (RLP < 0.01 in liquid water at ambient conditions),15 however, makes the application of this method numerically demanding in liquids and practically prohibitive in solids.

In a pioneering, but not highly cited, 1994 paper, Palmer showed how the heat conductivity can be determined by a direct analysis of the energy-density fluctuations in EMD simulations.16 This method does not rely on the value of the Landau–Planczek ratio, and therefore can be equally applied to liquids and solids. However, its broad applicability is also plagued by the numerical and statistical problems related to the evaluation and long-wavelength extrapolation of the peak of the energy dynamical structure factor (see below for a precise definition) and, on a more conceptual level, by the intrinsic indeterminacy of the energy density at the molecular scale, discussed above.

Building on the work of Palmer,16 in this paper we first reexamine the concept of gauge invariance in the specific context of the energy-density approach to heat transport, and then we show that a combination of recent advances in data analysis and signal processing, based on Bayesian inference17 and cepstral analysis,18,19 can be leveraged to overcome the numerical difficulties that have plagued the application of Palmer’s approach. We argue that the proposed methodology naturally lends itself to a quantitative evaluation of the resulting statistical uncertainties and can be applied using simulation cells of moderate size and EMD runs of moderate length; we finally showcase it with three applications to the paradigmatic cases of liquid water and Argon, and solid amorphous Silica.

The standard GK expression for the heat conductivity, κ, reads:12–14,
(3)
(4)
In Eq. (3) V and T are the system’s volume and temperature, kB the Boltzmann’s constant, Ĵ is any Cartesian component of the macroscopic average of the heat current density, Ĵ=1V∫Vȷ̂q(r)dr⁠, a hat indicates an implicit dependence on phase-space variables, Γ—as in Ĵ=J(Γ)—and ⟨⋅⟩ indicates an equilibrium average over the initial conditions of a molecular trajectory. The heat current density is defined as:12–14  ȷ̂q=ȷ̂e−(p+e)ȷ̂n/n⁠, where ȷ̂e and ȷ̂n are the energy and molecular-number density currents, respectively, and p, e are the system’s pressure and average energy, respectively. The explicit expressions of the energy and energy-current densities depend on the specific simulation framework (e.g. force-field vs ab initio), as well as on gauge fixing (see below). A detailed derivation of the relevant formulas is presented, e.g., in Ref. 6.
In a one-component system, the macroscopic average of the molecular number current density is proportional to the total momentum, which, being conserved, can be assumed to vanish. As in this work we will only consider this class of systems, we will not make any difference between heat and energy densities, and we will therefore omit any suffixes aimed at distinguishing one from the other. In Eq. (4) χ̃′′ is the imaginary part of the heat-temperature susceptivity, which, according to the fluctuation-dissipation theorem, is proportional to the energy dynamical structure factor (EDSF), C̃(k,ω)⁠:
(5)
(6)
where ẽ̂(k) is the space Fourier transform of the energy density, ê(r)⁠. The equivalence between Eqs. (3) and (4) follows from the fluctuation-dissipation theorem, Eq. (5), and the continuity equation, which, in the frequency-wavenumber domain, reads ωẽ̂(k,ω)=k⋅ȷ̃̂(k,ω)⁠. Thanks to this relation, the density-density and current-current correlation functions can be expressed in terms of each other as: C̃ϵϵ(k,ω)=k2ω2C̃jj(k,ω)⁠.
By combining the Onsager’s regression hypothesis1 with a hydrodynamic description of the long-time evolution of the long-wavelength components of the energy-density fluctuations, the small-frequency/small-wavevector limit of the EDSF can be shown to read:13,14
(7)
Inspection of Eq. (7) reveals that for large wavelengths the static limit of the EDSF reads:
(8)
We conclude that the heat conductivity can be equivalently expressed as:
(9)
(10)
Note the non-commutativity of the limits in Eq. (4), manifested by Eq. (7). This well-known singularity, which has been extensively discussed in the literature following Kadanoff and Martin’s landmark paper,12 should not come as a surprise, for the zero-frequency limit of a generic response function is a static susceptibility. It is an obvious consequence of time-reversal invariance that, strictly speaking, a static perturbation cannot induce a stationary current. The interpretation of this equation when the limits are performed in the correct order is as follows: While the long-time limit of the current induced by a temperature fluctuation vanishes at all finite wavelengths, as the system approaches equilibrium, at any finite time, the current tends to a finite, time-independent, value as the wavelength of the perturbation grows large. This behavior defines a stationary current in response to the temperature fluctuations.

The energy density, ê(r)⁠, entering Eq. (6) is intrinsically ill-defined, as it is affected by a so-called gauge freedom, deriving from the insensitivity of the total energy—and therefore of all the macroscopic thermal properties—on the addition of the divergence of a bounded vector field to the energy density: e(r) → e(r) + ∇ · p(r). The independence of the heat conductivity on the specific representation of the energy density is discussed, e.g., in Refs. 5, 7, 8, and 18, and demonstrated below in the context of the present approach.

In order to prove the gauge invariance of Eq. (9), we examine how a gauge transformation, e′(r) = e(r) + ∇ ⋅p(r), affects the zero-frequency value of the EDSF, C̃(k,ω=0)⁠. Under such a transformation, the integrand in Eq. (6) would go into:
(11)
According to the Wiener-Kintchine theorem,20,21 the Fourier transform of the time auto-correlation function of a zero-mean stationary stochastic process, X̂(t)⁠, in the large-time limit is proportional to the variance of the Fourier transform of the process, X̃̂τ(ω)=∫−τ2τ2X̂(t)eiωtdt⁠:
(12)
The ω = 0 limit of this relation reads:
(13)
When X̂(t) is a conserved flux, the left-hand side of Eq. (13) is the Green-Kubo expression of a transport coefficient. This equation is a generalization of Einstein’s celebrated relation between a particle’s diffusivity and its velocity auto-correlation function.22 Einstein’s formula was extended to generic transport coefficients by Helfand in 1960.23 
According to these considerations, the ω = 0 value of the transformed EDSF can be computed as
(14)
The second term of Eq. (14) is bounded from above by
(15)
and therefore, since p(r) has a finite variance, it vanishes in the k → 0 limit. Since in k → 0 the limit both ẽ̂(k,t) and p̃̂(−k,t′) are real, the third term in Eq. (14) vanishes as well. While at finite k it is to be expected that κ(k) in Eq. (10) may be affected by the specific choices of p(r), its long-wavelength limit does not depend on it, thus providing an independent demonstration of the gauge invariance of the heat conductivity, particularly adapted to the present framework.
The method embodied in Eq. (9) has been benchmarked against standard GK calculations based on Eq. (3) in the cases of liquid water, liquid Argon and solid amorphous Silica. Water was modeled with a sample of 216 molecules interacting through the SPC/E force field24 at T = 300 K and p = 0 bar. Argon was modeled with a sample of 864 atoms interacting through a Lennard-Jones potential25 at T = 100 K and p = 0 bar. Amorphous Silica was modeled with a sample of 648 atoms interacting via the BKS potential26 at 500 K. Long EMD simulations were run in the micro-canonical (NVE) ensemble after careful equilibration performed in a number of different ensembles. In all cases periodic boundary conditions (PBC) were used. All the technical details of the simulations are reported in  Appendix A. In order to estimate the ω = 0 value of the ESDF in Eq. (10), the time series of the Fourier transform of the energy density, ẽ(kℓ,t)=∑nϵn(t)e−ikℓ⋅Rn(t)⁠, was first collected along the EMD trajectory at regular time steps. In the preceding expression ϵn(t) and Rn(t) are the local energy and position of nth atom, respectively, at time t. Leveraging the Wiener-Kintchine theorem, the ESDF is then evaluated as:
(16)
where τ is the length of the EMD run and the integral is actually evaluated as a discrete Fourier transform, as explained, e.g., in Refs. 18 and 19. In the signal-processing literature, a sample of a power spectral density (PSD) drawn from a finite and discrete time series is often referred to as a periodogram. In Eq. (16) and in the rest of this section a hat indicates that the quantity it is put on (the periodogram and the energy density, in the present case) is actually a sample drawn from a stochastic process. It can be demonstrated that the periodogram constitutes a set of stochastic variables that are uncorrelated for different frequencies and that, for a given frequency, ωj=2πjτ⁠, are distributed as:
(17)
where the ξ̂’s are independent χ22 variates. From Eq. (17) it follows that the (energy-density) periodogram is an unbiased estimator of the PSD (the ESDF, in this case), C̃̂(k,ωj)=C̃(k,ωj)⁠. This estimator is not consistent, though, because the variance of any of its values is independent of the length of the EMD run. In order to get such a consistent estimator (at the price of introducing some bias), it is expedient to make use of the so called cepstral analysis method.18,19,27 The cepstrum is defined as the (inverse) discrete Fourier transform of the logarithm of the periodogram. The logarithm turns the multiplicative noise affecting the periodogram into an additive noise affecting the cepstrum; furthermore, the PSD, as well as its logarithm, are much smoother than the periodogram and the cepstrum itself. Therefore, a consistent estimate of the log-PSD can be obtained by applying a low-pass filter to the cepstrum. The optimal value of the low-pass quefrency27—i.e. the number of inverse Fourier coefficients of the cepstrum—is determined in such a way that the estimated bias introduced by the filter is of the order of the resulting statistical noise. This number can be obtained by any model-selection method, and in our applications we use the Akaike information criterion.18,19,28
With the aim of evaluating Eq. (9), the EDSF was first sampled on a homogeneous grid of wave-vectors, {kℓ}, with magnitude smaller than a preassigned cutoff, |kℓ| ≤ kmax, compatible with the chosen PBCs. Because of spherical symmetry, the EDSF only depends on the modulus of its wave-vector argument and can therefore be averaged accordingly. In practice, PBCs reduce the symmetry from spherical to cubic, thus making some PBC-allowed wave-vectors non-equivalent to others with the same magnitude. This artifact is eliminated by averaging the estimates of the EDSF corresponding to different wave-vectors of the same magnitude, thus further reducing the variance of any individual value. The set of wave-vectors with a same magnitude, ki, will be referred to as the ki-star of wave-vectors. In Fig. 1 we display the EDSF of liquid Argon at the k-star of minimum magnitude compatible with the chosen PBCs, kmin = 1.77 nm−1, and the one of largest magnitude utilized in our Bayesian analysis (see below), kmax = 11.3 nm−1. The increasing sharpness of the peaks for increasing wavelength (decreasing wave-vector) is a signature of their hydrodynamic nature, as manifested by Eq. (7), and makes cepstral analysis increasingly difficult to perform, thus giving rise to a larger statistical incertitude. Because of isotropy, the wave-vector dependence of the conductivity, Eq. (10), is an even function of the squared modulus of k, which we assume to be well approximated by a polynomial of order M:
(18)
In order to determine the coefficients of the polynomial fit in Eq. (18), w = {w0⋯wM}, we relied on a Bayesian inference method, by seeking to maximize their posterior probability, given the occurrence of the measured data, D⁠: pw|D=pD|wp(w)/p(D)⁠, where p(D) and p(w) are the data and model-parameter priors, respectively. The data are the conductivities measured at each k-star, which are independent normal variates whose expectations, κi, and variances, σi2 are estimated by the cepstral analysis procedure outlined above. The normality of the distribution of the measured values of κ(k) is showcased in  Appendix B. The probability of the data conditional to the values of the coefficients of the polynomial is thus given by pD|w∝e−L(w)⁠, where L(w)=∑iκi−w⋅Φi22σi2⁠, is the likelihood function and Φi is an array of monomials of even degree up to order 2M, evaluated at the k-stars, ki, for which the k-dependent conductivity has been estimated: Φi=Φ(ki)≡{1,ki2…ki2M}⁠. We now make the further assumption that the prior distribution of the parameters is normal, p(w)=α2πM2⁡exp−αw2⁠, thus allowing the extrapolation to prevent over-fitting, which, in ordinary linear regression, would be equivalent to introducing a regularization term, αw2⁠. By leveraging again Bayes theorem, the optimal number of parameters, M, and regularization parameter, α, are determined as those that maximize the probability of their occurrence, conditionally to actual observation of the data set:
(19)
The actual procedure that we followed to determine M, α, w, and, therefore, κ = w0 is described in Refs. 17 and 30.
FIG. 1.
Upper panel: EDSF of liquid Argon computed at PBC-compatible k-vectors of minimum magnitude, kmin = 1.77 nm−1. Lower panel: same, but for the k-vectors, kmax = 11.3 nm−1. The window-filtered line refers to a moving average29 performed width a window with of Δ = 0.2 THz. The shaded red area indicates the statistical incertitude.
View largeDownload slide

Upper panel: EDSF of liquid Argon computed at PBC-compatible k-vectors of minimum magnitude, kmin = 1.77 nm−1. Lower panel: same, but for the k-vectors, kmax = 11.3 nm−1. The window-filtered line refers to a moving average29 performed width a window with of Δ = 0.2 THz. The shaded red area indicates the statistical incertitude.

FIG. 1.
Upper panel: EDSF of liquid Argon computed at PBC-compatible k-vectors of minimum magnitude, kmin = 1.77 nm−1. Lower panel: same, but for the k-vectors, kmax = 11.3 nm−1. The window-filtered line refers to a moving average29 performed width a window with of Δ = 0.2 THz. The shaded red area indicates the statistical incertitude.
View largeDownload slide

Upper panel: EDSF of liquid Argon computed at PBC-compatible k-vectors of minimum magnitude, kmin = 1.77 nm−1. Lower panel: same, but for the k-vectors, kmax = 11.3 nm−1. The window-filtered line refers to a moving average29 performed width a window with of Δ = 0.2 THz. The shaded red area indicates the statistical incertitude.

Close modal

In Fig. 2 we display the wave-vector dependence of Eq. (10), κ(k), computed for SPC/E water, liquid Argon and solid amorphous Silica via the cepstral analysis, and the k → 0 extrapolation of Eq. (9), κ, obtained from the Bayesian regression method. The results are in good accordance with the GK integral estimated via the cepstral analysis on the same trajectory, also reported in the figures. In  Appendix C we study the finite size effects on the computation of Eq. (10). We studied the size convergence of Eq. (10) for SPC/E water modeled with samples of different linear box size, L: 18.5, 31.1 and 37.4Å⁠; corresponding to 216, 1000 and 1728 molecules respectively.

FIG. 2.
Upper panel: Wave-vector dependence of the heat conductivity, κ(k), Eq. (10), of SPC/E water at T = 300 K and 0 bar, as evaluated from the EMD simulation, kmin = 3.40 nm−1. Central panel: same for liquid Argon at T = 100 K and 0 bar, kmin = 1.77 nm−1. Lower panel: same for solid amorphous Silica at T = 500 K, kmin = 3.04 nm−1. The data are averages of k-vectors with a same magnitude, compatible with the PBCs in use. Cubic average refers to an average performed over k-vectors that are strictly equivalent by cubic symmetry, i.e. excluding vectors that are equivalent by spherical, but not cubic symmetry. The red line indicates the Bayesian fit with the associated estimated uncertainty. At k = 0 we report the Bayesian extrapolation and the GK estimate, on the same EMD trajectory, of the heat conductivity.
View largeDownload slide

Upper panel: Wave-vector dependence of the heat conductivity, κ(k), Eq. (10), of SPC/E water at T = 300 K and 0 bar, as evaluated from the EMD simulation, kmin = 3.40 nm−1. Central panel: same for liquid Argon at T = 100 K and 0 bar, kmin = 1.77 nm−1. Lower panel: same for solid amorphous Silica at T = 500 K, kmin = 3.04 nm−1. The data are averages of k-vectors with a same magnitude, compatible with the PBCs in use. Cubic average refers to an average performed over k-vectors that are strictly equivalent by cubic symmetry, i.e. excluding vectors that are equivalent by spherical, but not cubic symmetry. The red line indicates the Bayesian fit with the associated estimated uncertainty. At k = 0 we report the Bayesian extrapolation and the GK estimate, on the same EMD trajectory, of the heat conductivity.

FIG. 2.
Upper panel: Wave-vector dependence of the heat conductivity, κ(k), Eq. (10), of SPC/E water at T = 300 K and 0 bar, as evaluated from the EMD simulation, kmin = 3.40 nm−1. Central panel: same for liquid Argon at T = 100 K and 0 bar, kmin = 1.77 nm−1. Lower panel: same for solid amorphous Silica at T = 500 K, kmin = 3.04 nm−1. The data are averages of k-vectors with a same magnitude, compatible with the PBCs in use. Cubic average refers to an average performed over k-vectors that are strictly equivalent by cubic symmetry, i.e. excluding vectors that are equivalent by spherical, but not cubic symmetry. The red line indicates the Bayesian fit with the associated estimated uncertainty. At k = 0 we report the Bayesian extrapolation and the GK estimate, on the same EMD trajectory, of the heat conductivity.
View largeDownload slide

Upper panel: Wave-vector dependence of the heat conductivity, κ(k), Eq. (10), of SPC/E water at T = 300 K and 0 bar, as evaluated from the EMD simulation, kmin = 3.40 nm−1. Central panel: same for liquid Argon at T = 100 K and 0 bar, kmin = 1.77 nm−1. Lower panel: same for solid amorphous Silica at T = 500 K, kmin = 3.04 nm−1. The data are averages of k-vectors with a same magnitude, compatible with the PBCs in use. Cubic average refers to an average performed over k-vectors that are strictly equivalent by cubic symmetry, i.e. excluding vectors that are equivalent by spherical, but not cubic symmetry. The red line indicates the Bayesian fit with the associated estimated uncertainty. At k = 0 we report the Bayesian extrapolation and the GK estimate, on the same EMD trajectory, of the heat conductivity.

Close modal

In Fig. 3 we display the Bayesian extrapolation (k = 0) of the heat conductivity of liquid Argon as a function of the maximum magnitude of the wave-vectors considered in the fit, kc. The Bayesian extrapolation is robust and converges to the GK result computed on the same trajectory, also reported in the figures, when we enlarge the data set, thus demonstrating the numerical stability and statistical consistency of our method.

FIG. 3.
Upper panel: Bayesian prediction of liquid Argon at 100 K and 0 bar, at k = 0 as a function of the maximum k-vector included in the fit, kc, confronted with the GK integral computed on the same trajectory with the cepstral analysis. Lower panel: Number of parameters of the Bayesian regression fit as a function of the maximum k-vector included in the fit.
View largeDownload slide

Upper panel: Bayesian prediction of liquid Argon at 100 K and 0 bar, at k = 0 as a function of the maximum k-vector included in the fit, kc, confronted with the GK integral computed on the same trajectory with the cepstral analysis. Lower panel: Number of parameters of the Bayesian regression fit as a function of the maximum k-vector included in the fit.

FIG. 3.
Upper panel: Bayesian prediction of liquid Argon at 100 K and 0 bar, at k = 0 as a function of the maximum k-vector included in the fit, kc, confronted with the GK integral computed on the same trajectory with the cepstral analysis. Lower panel: Number of parameters of the Bayesian regression fit as a function of the maximum k-vector included in the fit.
View largeDownload slide

Upper panel: Bayesian prediction of liquid Argon at 100 K and 0 bar, at k = 0 as a function of the maximum k-vector included in the fit, kc, confronted with the GK integral computed on the same trajectory with the cepstral analysis. Lower panel: Number of parameters of the Bayesian regression fit as a function of the maximum k-vector included in the fit.

Close modal

In this work we have integrated a method originally proposed by Palmer in 1994,16 to compute the heat conductivity in extended systems from energy-density fluctuations, with modern data and signal analysis techniques that considerably enhance its applicability. As a downside, as already anticipated by Palmer, the numerical efficiency of the methodology presented here appears to be poorer than that of the standard GK procedure based on the fluctuations of the energy flux. In the specific case of Argon, we estimated that, in order to achieve the same statistical accuracy of the standard GK procedure, one would need EMD trajectories one-order of magnitude as long. The enhanced statistical performance of current-based methods for estimating transport coefficients, as opposed to density-based methods, can be ascribed to the necessity of performing long-wavelength extrapolations when utilizing the latter. Avoiding any such extrapolation would significantly enhance the performance of any statistical analysis, no matter how smart this may be. Of course, this numerical superiority may, or may not, be tarnished by the need to devise and implement an explicit expression for the macroscopic energy flux.

In particular, it may be challenging, if not impossible, to devise an analytical form of the energy flux when advanced energy functionals, such as hybrid or meta-GGA, are used for ab initio MD. We believe that a density-based approach, such as the one examined in the present paper or in CF’s,11 would be mandatory in these cases. The relevance of these cases is however increasingly limited by the emergence of artificial-intelligence methods to train classical force fields with quantum mechanical accuracy,31 which allow for a fast and accurate evaluation of various transport coefficients using the standard GK approach, including the heat conductivity.32,33 Furthermore, all the density-based methods proposed so far are at present limited to one-component systems, whereas a systematic extension of the standard GK one has been recently established.34 In summary, we conclude that standard current-based methods are preferable in all those cases where an explicit expression for the energy flux is available. In the rare cases where it is not (e.g. for ab initio MD with advanced functionals), density-based methods may be the way to go. The method discussed here and the one proposed by CF11 are conceptually similar, but the latter is in practice only applicable when the Landau-Planczek ratio is sizeable, i.e. for compressible fluids, whereas the former can be applied to solids and incompressible fluids as well. We also hint at the superiority of the data analysis technique introduced in Ref. 30 and utilized here, which allows for an easy control and optimal estimate of the statistical accuracy of the transport coefficients. This technique would readily integrate into CF’s approach,11 as well as in any methodology relying on the long-wavelength extrapolation of correlation functions estimated for finite, periodic models.

The authors are grateful to Riccardo Bertossa, Alfredo Fiorentino, Federico Grasselli, Paolo Pegolo, and Davide Tisi for many insightful suggestions and valuable discussions. This work was partially supported by the European Commission through the MaX Centre of Excellence for supercomputing applications (Grant No. 101093374) and by the Italian MUR, through the PRIN project FERMAT (Grant No. 2017KFY7XF) and the Italian National Centre from HPC, Big Data, and Quantum Computing (Grant No. CN00000013).

The authors have no conflicts to disclose.

ED, MGI, and SB conceived the study and provided the overall conceptual framework. ED and MGI conducted the analytical work and developed the necessary computer codes. ED run and analyzed the computer simulations. SB provided supervision and guidance throughout the project. The manuscript was collaboratively written by ED and SB.

Enrico Drigo: Conceptualization (equal); Data curation (equal); Formal analysis (equal); Investigation (equal); Methodology (equal); Project administration (equal); Software (equal); Writing – original draft (equal); Writing – review & editing (equal). Maria Grazia Izzo: Conceptualization (equal); Data curation (equal); Formal analysis (equal); Investigation (equal); Supervision (equal); Writing – original draft (equal); Writing – review & editing (equal). Stefano Baroni: Conceptualization (equal); Formal analysis (equal); Funding acquisition (equal); Investigation (equal); Methodology (equal); Project administration (equal); Supervision (equal); Writing – original draft (equal); Writing – review & editing (equal).

The data and scripts that support the plots and relevant results within this paper are available on the Materials Cloud platform.35 See https://doi.org/10.24435/materialscloud:ft-b3, Ref. 40.

We modeled water with the SPC/E force field that assumes the intra-molecular bond lengths and angles to be constrained to their equilibrium values, while inter-molecular atomic forces are described by Coulomb-plus-Lennard-Jones potentials whose parameters are fitted to experiment.24 The resulting agreement between computed and measured static and dynamic properties, such as the radial distribution function ad diffusivity, is very good. The long-range Coulomb interaction was treated via the Ewald summation method with a 9 Å  cutoff. We have performed NVE EMD simulations of 216 SPC/E water molecules at T = 300 K and p = 0 bar, using PBCs. After a careful equilibration in the NPT and NVT ensembles for 125 ps, we sampled the energy density from a 1-ns NVE trajectory with a time step of 0.25 fs. The energy density was computed and dumped at intervals of 2.5 fs. Liquid Argon was simulated at T = 100 K and p = 0 bar using a sample of 864 atoms interacting through an iter-atomic Lennard-Jones potential whose parameters are fitted to experiments.25 After a preliminary equilibration in the NPT ensemble, we sampled the system from a 2-ns NVE trajectory using a time step of 1 fs. In this case, the energy density was computed and dumped at intervals of 5 fs. Solid amorphous Silica was simulated at 500 K using a sample of 648 atoms modeled via the BKS potential.26 The initial configuration was taken from Ref. 36 and optimized so as to make atomic forces smaller than a preassigned threshold of 10−8eV/Å⁠. After careful equilibration for 125 ps in the NVT ensemble, we simulated the system for 1 ns in the NVE ensemble with a timestep of 0.25 fs. The energy density was computed and dumped at intervals of 5 fs. All simulations were performed with the parallel code LAMMPS37,38 using the velocity-Verlet algorithm.

Our Bayesian regression procedure is based on the assumption that the data [the values of the wavevector-dependent heat conductivities, κ(k), Eq. (10)] are normal variates at any value of their argument, k. In order to evaluate the soundness of this hypothesis, we have examined the statistics of the values of the κ’s collected from 200-ps segments extracted from our 20-ns long EMD trajectory for liquid water. In Fig. 4 we display the distribution of the data thus collected, corresponding to the wavevectors of minimum magnitute, kmin = 3.40 nm−1 and of maximum magnitude, kmax = 9.63 nm−1 utilized in our Bayesian procedure. In both cases, the data pass the Shapiro-Wilk (SW) normality test33,39 with respect to a standard significance level of α = 0.05.

FIG. 4.
Distribution of the wave-vector dependence of the heat conductivity, κ(k), Eq. (10), of SPC/E water at T = 300 K and 0 bar as evaluated from multiple 200 ps-long EMD segments extracted from a 20 ns long trajectory and the p-value of the Shapiro-Wilk normality test. The solid lines represent a Gaussian curve fitted on the distribution of the heat conductivities at the wavevector presented.
View largeDownload slide

Distribution of the wave-vector dependence of the heat conductivity, κ(k), Eq. (10), of SPC/E water at T = 300 K and 0 bar as evaluated from multiple 200 ps-long EMD segments extracted from a 20 ns long trajectory and the p-value of the Shapiro-Wilk normality test. The solid lines represent a Gaussian curve fitted on the distribution of the heat conductivities at the wavevector presented.

FIG. 4.
Distribution of the wave-vector dependence of the heat conductivity, κ(k), Eq. (10), of SPC/E water at T = 300 K and 0 bar as evaluated from multiple 200 ps-long EMD segments extracted from a 20 ns long trajectory and the p-value of the Shapiro-Wilk normality test. The solid lines represent a Gaussian curve fitted on the distribution of the heat conductivities at the wavevector presented.
View largeDownload slide

Distribution of the wave-vector dependence of the heat conductivity, κ(k), Eq. (10), of SPC/E water at T = 300 K and 0 bar as evaluated from multiple 200 ps-long EMD segments extracted from a 20 ns long trajectory and the p-value of the Shapiro-Wilk normality test. The solid lines represent a Gaussian curve fitted on the distribution of the heat conductivities at the wavevector presented.

Close modal

In order to study the size-convergence of the method, we computed the wave-vector dependence of the heat conductivity, κ(k), Eq. (10), for SPC/E water at 300 K and 0 bar, for three different linear box sizes, L: 18.5, 31.1 and 37.4Å⁠; corresponding to 216, 1000 and 1728 water molecules respectively.  Appendix A describes the simulation protocol that we employed for the simulations of SPC/E water at each size.

In Fig. 5 we display the size-convergence of Eq. (10), for SPC/E water. Increasing the number of molecules, e.g. augmenting the box size at fixed mass density, allows to investigate further the k → 0 limit of κ(k), computing smaller module wave-vectors. Focusing on the wave-vectors that are compatible with all the sizes we considered, e.g. approximatively 3.40 nm−1, compared to the simulation of 1728 molecules, already with 216 molecules we were able to obtain a well converged estimate of Eq. (10).

FIG. 5.
Wave-vector dependence of the heat conductivity, κ(k), Eq. (10), of SPC/E water at 300 K and 0 bar, as evaluated from the EMD simulations for three different linear box sizes, L: 18.5, 31.1 and 37.4Å; corresponding to 216, 1000 and 1728 water molecules respectively. The data are averages of k-vectors with a same magnitude, compatible with the PBCs in use. Cubic averages, represented as crosses, refer to the averages performed over k-vectors that are strictly equivalent by cubic symmetry, i.e. excluding vectors that are equivalent by spherical, but not cubic symmetry. At k = 0 we report, as dots, the Bayesian extrapolations of the heat conductivity at each size.
View largeDownload slide

Wave-vector dependence of the heat conductivity, κ(k), Eq. (10), of SPC/E water at 300 K and 0 bar, as evaluated from the EMD simulations for three different linear box sizes, L: 18.5, 31.1 and 37.4Å⁠; corresponding to 216, 1000 and 1728 water molecules respectively. The data are averages of k-vectors with a same magnitude, compatible with the PBCs in use. Cubic averages, represented as crosses, refer to the averages performed over k-vectors that are strictly equivalent by cubic symmetry, i.e. excluding vectors that are equivalent by spherical, but not cubic symmetry. At k = 0 we report, as dots, the Bayesian extrapolations of the heat conductivity at each size.

FIG. 5.
Wave-vector dependence of the heat conductivity, κ(k), Eq. (10), of SPC/E water at 300 K and 0 bar, as evaluated from the EMD simulations for three different linear box sizes, L: 18.5, 31.1 and 37.4Å; corresponding to 216, 1000 and 1728 water molecules respectively. The data are averages of k-vectors with a same magnitude, compatible with the PBCs in use. Cubic averages, represented as crosses, refer to the averages performed over k-vectors that are strictly equivalent by cubic symmetry, i.e. excluding vectors that are equivalent by spherical, but not cubic symmetry. At k = 0 we report, as dots, the Bayesian extrapolations of the heat conductivity at each size.
View largeDownload slide

Wave-vector dependence of the heat conductivity, κ(k), Eq. (10), of SPC/E water at 300 K and 0 bar, as evaluated from the EMD simulations for three different linear box sizes, L: 18.5, 31.1 and 37.4Å⁠; corresponding to 216, 1000 and 1728 water molecules respectively. The data are averages of k-vectors with a same magnitude, compatible with the PBCs in use. Cubic averages, represented as crosses, refer to the averages performed over k-vectors that are strictly equivalent by cubic symmetry, i.e. excluding vectors that are equivalent by spherical, but not cubic symmetry. At k = 0 we report, as dots, the Bayesian extrapolations of the heat conductivity at each size.

Close modal
2.
M. S.
Green
,
J. Chem. Phys.
20
,
1281
(
1952
);
https://doi.org/10.1063/1.1700722
Crossref
Search ADS
Google Scholar
 
M. S.
Green
J. Chem. Phys.
22
,
398
(
1954
).
https://doi.org/10.1063/1.1740082
Crossref
Search ADS
Google Scholar
 
3.
R.
Kubo
,
J. Phys. Soc. Jpn.
12
,
570
(
1957
);
https://doi.org/10.1143/JPSJ.12.570
Crossref
Search ADS
Google Scholar
 
R.
Kubo
,
M.
Yokota
, and
S.
Nakajima
,
J. Phys. Soc. Jpn.
12
,
1203
(
1957
).
https://doi.org/10.1143/JPSJ.12.1203
Crossref
Search ADS
Google Scholar
 
4.
S.
Stackhouse
,
L.
Stixrude
, and
B. B.
Karki
,
Phys. Rev. Lett.
104
,
208501
(
2010
).
https://doi.org/10.1103/physrevlett.104.208501
Google Scholar
Crossref
Search ADS
PubMed
 
5.
A.
Marcolongo
,
P.
Umari
, and
S.
Baroni
,
Nat. Phys.
12
,
80
–
84
(
2016
).
https://doi.org/10.1038/nphys3509
Google Scholar
Crossref
Search ADS
 
6.
L.
Ercole
,
A.
Marcolongo
,
P.
Umari
, and
S.
Baroni
,
J. Low Temp. Phys.
185
,
79
(
2016
).
https://doi.org/10.1007/s10909-016-1617-6
Google Scholar
Crossref
Search ADS
 
7.
F.
Grasselli
 and
S.
Baroni
,
Eur. Phys. J. B
94
,
160
(
2021
).
https://doi.org/10.1140/epjb/s10051-021-00152-5
Google Scholar
Crossref
Search ADS
PubMed
 
8.
A.
Marcolongo
,
L.
Ercole
, and
S.
Baroni
,
J. Chem. Theory Comput.
16
,
3352
(
2020
).
https://doi.org/10.1021/acs.jctc.9b01174
Google Scholar
Crossref
Search ADS
PubMed
 
9.
F.
Grasselli
 and
S.
Baroni
,
Nat. Phys.
15
,
967
(
2019
).
https://doi.org/10.1038/s41567-019-0562-0
Google Scholar
Crossref
Search ADS
 
10.
P.
Pegolo
,
F.
Grasselli
, and
S.
Baroni
,
Phys. Rev. X
10
,
041031
(
2020
).
https://doi.org/10.1103/physrevx.10.041031
Google Scholar
Crossref
Search ADS
 
11.
B.
Cheng
 and
D.
Frenkel
,
Phys. Rev. Lett.
125
,
130602
(
2020
).
https://doi.org/10.1103/physrevlett.125.130602
Google Scholar
Crossref
Search ADS
PubMed
 
12.
L. P.
Kadanoff
 and
P. C.
Martin
,
Ann. Phys.
24
,
419
(
1963
).
https://doi.org/10.1016/0003-4916(63)90078-2
Google Scholar
Crossref
Search ADS
 
13.
D.
Forster
,
Hydrodynamic Fluctuations, Broken Symmetry, and Correlation Functions
(
CRC Press
,
2018
).
Google Scholar
Crossref
Search ADS
 
14.
J.-P.
Hansen
 and
I. R.
McDonald
,
Theory of Simple Liquids
, 4th ed. (
Academic Press
,
Oxford
,
2013
), p.
i
.
Google Scholar
 
15.
Y. A.
Cengel
 and
M. A.
Boles
,
Thermodynamics: An Engineering Approach
(
McGraw-Hill Education
,
2011
).
Google Scholar
 
16.
B. J.
Palmer
,
Phys. Rev. E
49
,
2049
(
1994
).
https://doi.org/10.1103/physreve.49.2049
Google Scholar
Crossref
Search ADS
 
17.
C. M.
Bishop
,
Pattern Recognition and Machine Learning
, 1st ed. (
Springer
,
2006
).
Google Scholar
 
18.
L.
Ercole
,
A.
Marcolongo
, and
S.
Baroni
,
Sci. Rep.
7
,
15835
(
2017
).
https://doi.org/10.1038/s41598-017-15843-2
Google Scholar
Crossref
Search ADS
PubMed
 
19.
L.
Ercole
,
R.
Bertossa
,
S.
Bisacchi
, and
S.
Baroni
,
Comput. Phys. Commun.
280
,
108470
(
2022
).
https://doi.org/10.1016/j.cpc.2022.108470
Google Scholar
Crossref
Search ADS
 
21.
A.
Khintchine
,
Math. Ann.
109
,
604
(
1934
).
https://doi.org/10.1007/bf01449156
Google Scholar
Crossref
Search ADS
 
24.
H. J. C.
Berendsen
,
J. R.
Grigera
, and
T. P.
Straatsma
,
J. Phys. Chem.
91
,
6269
(
1987
).
https://doi.org/10.1021/j100308a038
Google Scholar
Crossref
Search ADS
 
25.
L.
Rowley
,
D.
Nicholson
, and
N.
Parsonage
,
J. Comput. Phys.
17
,
401
(
1975
).
https://doi.org/10.1016/0021-9991(75)90042-x
Google Scholar
Crossref
Search ADS
 
26.
B. W. H.
van Beest
,
G. J.
Kramer
, and
R. A.
van Santen
,
Phys. Rev. Lett.
64
,
1955
(
1990
).
https://doi.org/10.1103/physrevlett.64.1955
Google Scholar
Crossref
Search ADS
PubMed
 
27.
B. P.
Bogert
,
J. R.
Healy
, and
J. W.
Tukey
,
Proceedings of the Symposium on Time Series Analysis, Brown University, June 11–14th, 1962
,
SIAM Series in Applied Mathematics
(
Wiley
,
Hoboken, NJ
,
1963
), pp.
209
–
243
.
Google Scholar
 
28.
H.
Akaike
,
IEEE Trans. Autom. Control
19
,
716
(
1974
).
https://doi.org/10.1109/tac.1974.1100705
Google Scholar
Crossref
Search ADS
 
29.
E. W.
Weisstein
, “
Moving average. From MathWorld–A Wolfram web resource
,” https://mathworld.wolfram.com/MovingAverage.html
30.
E.
Drigo
 and
S.
Baroni
, “
Seebeck coefficient of liquid water from equilibrium molecular dynamics
,”
J. Chem. Theory Comput.
(published online,
2023
).
https://doi.org/10.1021/acs.jctc.3c00760
Google Scholar
 
32.
D.
Tisi
,
L.
Zhang
,
R.
Bertossa
,
H.
Wang
,
R.
Car
, and
S.
Baroni
,
Phys. Rev. B
104
,
224202
(
2021
).
https://doi.org/10.1103/physrevb.104.224202
Google Scholar
Crossref
Search ADS
 
33.
C.
Malosso
,
L.
Zhang
,
R.
Car
,
S.
Baroni
, and
D.
Tisi
,
npj Comput. Mater.
8
,
139
(
2022
).
https://doi.org/10.1038/s41524-022-00830-7
Google Scholar
Crossref
Search ADS
 
34.
R.
Bertossa
,
F.
Grasselli
,
L.
Ercole
, and
S.
Baroni
,
Phys. Rev. Lett.
122
,
255901
(
2019
).
https://doi.org/10.1103/physrevlett.122.255901
Google Scholar
Crossref
Search ADS
PubMed
 
35.
L.
Talirz
,
S.
Kumbhar
,
E.
Passaro
,
A. V.
Yakutovich
,
V.
Granata
,
F.
Gargiulo
,
M.
Borelli
,
M.
Uhrin
,
S. P.
Huber
,
S.
Zoupanos
,
C. S.
Adorf
,
C. W.
Andersen
,
O.
Schütt
,
C. A.
Pignedoli
,
D.
Passerone
,
J.
VandeVondele
,
T. C.
Schulthess
,
B.
Smit
,
G.
Pizzi
, and
N.
Marzari
,
Sci. Data
7
,
299
(
2020
).
https://doi.org/10.1038/s41597-020-00637-5
Google Scholar
Crossref
Search ADS
PubMed
 
36.
A.
Fiorentino
,
P.
Pegolo
, and
S.
Baroni
,
npj Comput. Mater.
9
,
157
(
2023
).
https://doi.org/10.1038/s41524-023-01116-2
Google Scholar
Crossref
Search ADS
 
37.
S.
Plimpton
,
J. Comput. Phys.
117
,
1
(
1995
).
https://doi.org/10.1006/jcph.1995.1039
Google Scholar
Crossref
Search ADS
 
38.
A. P.
Thompson
,
H. M.
Aktulga
,
R.
Berger
,
D. S.
Bolintineanu
,
W. M.
Brown
,
P. S.
Crozier
,
P. J.
in 't Veld
,
A.
Kohlmeyer
,
S. G.
Moore
,
T. D.
Nguyen
,
R.
Shan
,
M. J.
Stevens
,
J.
Tranchida
,
C.
Trott
,
S. J.
Plimpton
, and
S. J.
Plimpton
,
Comput. Phys. Commun.
271
,
108171
(
2022
).
https://doi.org/10.1016/j.cpc.2021.108171
Google Scholar
Crossref
Search ADS
 
39.
S. S.
Shapiro
 and
M. B.
Wilk
,
Biometrika
52
,
591
(
1965
).
https://doi.org/10.2307/2333709
Google Scholar
Crossref
Search ADS
 
40.
Enrico
Drigo
,
Maria
Grazia Izzo
, and
Stefano
Baroni
 (
2023
). “
Heat conductivity from energy-density fluctuations
,” Materials Cloud. https://doi.org/10.24435/materialscloud:ft-b3.
Google Scholar
 
© 2023 Author(s). All article content, except where otherwise noted, is licensed under a Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).
Published open access through an agreement with Scuola Internazionale Superiore di Studi Avanzati