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.
I. INTRODUCTION
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
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.
II. THEORY
A. Gauge invariance
The energy density, , 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.
III. COMPUTER SIMULATIONS
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.
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.
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 ; corresponding to 216, 1000 and 1728 molecules respectively.
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.
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.
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.
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.
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.
IV. CONCLUSIONS
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.
ACKNOWLEDGMENTS
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).
AUTHOR DECLARATIONS
Conflict of Interest
The authors have no conflicts to disclose.
Author Contributions
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).
DATA AVAILABILITY
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.
APPENDIX A: COMPUTATIONAL DETAILS
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 . 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.
APPENDIX B: STATISTICAL PROPERTIES OF THE DATASET
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.
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.
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.
APPENDIX C: FINITE-SIZE EFFECTS
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 ; 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).
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 ; 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.
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 ; 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.