revtex4-2Repair the float

Thermal Conductivity Predictions with Foundation Atomistic Models

Abstract

Recent advances in machine learning have led to foundation models for atomistic materials chemistry, potentially enabling quantum-accurate descriptions of interatomic forces at reduced computational cost. These models are benchmarked by predicting materials’ properties over large databases; however, these computationally intensive tests have been limited to basic quantities related to harmonic phonons, leaving uncertainty about the reliability for complex, technologically and experimentally relevant anharmonic heat-conduction properties. Here we present an automated framework that relies on foundation models to compute microscopic vibrational properties, and employs them within the Wigner formulation of heat transport to predict the macroscopic thermal conductivity in solids with arbitrary composition and structure. We apply this framework with the foundation models M3GNet, CHGNet, MACE-MP-0, and SevenNET to 103 diverse compounds, comparing predictions against first-principles references and introducing a benchmark metric based on conductivity. This framework paves the way for physics-aware, accurate predictions of vibrational and thermal properties, and for uncovering materials that violate semiclassical Boltzmann transport and feature exceptional heat-shielding or thermoelectric performance.

Over the past decades, several research groups have tackled the challenging computational task of fitting the Born-Oppenheimer potential energy surface as a function of atomic coordinates [1, 2, 3, 4, 5, 6, 7]. These efforts resulted in the development of so-called machine-learning potentials (MLPs), which allow us to reproduce first-principles energies and microscopic interatomic forces with nearly the same accuracy and orders-of-magnitude lower computational cost. These developments enable the prediction of macroscopic observables from the integration of atomistic properties at a computational cost much reduced compared to first-principles methods, effectively opening avenues to design materials for target applications from theory. For example, engineering the magnitude of the macroscopic thermal conductivity of a material through changes in its atomistic composition and structure is crucial for neuromorphic computing [8], thermoelectric energy harvesting [9], and aerospace [10] technologies. Major drawbacks of MLPs-based methods are the significant work required to generate a first-principles database used to train and validate MPLs [11], as well as the applicability limited to specific materials’ compositions or structural phases. Past works attempted to bypass these limitations by employing end-to-end machine-learning methods, which predict microscopic [12] or macroscopic [13, 14] materials’ properties very efficiently, but with the compromise of not rigorously resolving the fundamental physics that underlies them. Fittingly, recent work [15] has formally demonstrated the possibility to obtain a complete description (in the mathematical sense) of atomic environments (and thus of forces) by employing Message Passing Neural Networks (MPNNs) models [16] with many-body messages [17] (to be precise, with the MACE architecture [17] this can be achieved using 4-body terms and one message pass). This breakthrough has enabled the development of foundation machine-learning potentials (fMLPs) [18, 19, 20, 21, 22, 23], which are trained across nearly all chemical elements and can be directly combined to describe materials with diverse structures and compositions. Therefore, fMLPs could potentially overcome the problems related to complex training and the very limited transferability of conventional MLPs, ultimately enabling physics-aware predictions of macroscopic observables from the integration of atomistic quantities. Recent research efforts have assessed the accuracy of fMLPs in predicting, e.g., the structural stability of solids [24] or the harmonic vibrational frequencies [25, 26, 27, 28]. However, fMLPs have not yet been tested for predicting vibrational-anharmonic and heat-conduction properties, due to the complexity of the computational framework that relates the interatomic forces to the macroscopic conductivity [29, 30, 31, 32, 33, 34, 35, 36, 37].

Here we present an automated framework that uses fMLPs to compute anharmonic vibrational and thermal properties. We employ this to benchmark against first-principles reference data [30, 31] the predictions obtained from state-of-the-art (SOTA) fMLPs M3GNet [18], CHGNet [19], MACE-MP-0 [20], and SevenNET [21] (all trained on the Materials Project DFT-PBE database [38] and hereafter collectively referred to as ‘mp-fMLPs’). In particular, after detailing the methods to calculate harmonic and anharmonic vibrational properties of solids using fMLPs, we discuss how these quantities determine the thermal conductivity within the recently developed Wigner heat-Transport Equation (WTE) [34, 35], which generalizes the Boltzmann Transport Equation (BTE) [39] accounting not only for heat carried by particle-like propagation of phonons, but also for conduction arising from phonons’ wave-like tunneling. Thus, this framework is employed to predict from DFT-PBE or mp-fMLPs the thermal conductivity of 103 compounds made up of 34 different chemical species and having wurtzite, zincblende, or rocksalt structure. We introduce descriptors to quantify the accuracy of fMLPs in predicting anharmonic and thermal properties, and show that they could potentially be used as a metric to benchmark or fine-tune fMLPs. Finally, we show how the framework introduced paves the way for rigorous, physics-aware predictions of vibrational and thermal properties in materials with arbitrary composition and structure, highlighting its potential to find materials that violate semiclassical Boltzmann conduction relevant for thermal-insulation or thermoelectric applications.

Results
From vibrational energy to thermal conductivity
The temperature-dependent thermal conductivity ($\kappa(T)$) of a solid is a macroscopic, experimentally and technologically relevant quantity that describes the capability to conduct heat. To predict such quantity, we start from the Born-Oppenheimer Hamiltonian for atomic vibrations [35] expanded up to anharmonic third order in displacements from equilibrium, and we account for kinetic-energy perturbations due to isotopes [40],

$$\begin{split}&\hat{H}=\sum_{\bm{R},b,\alpha}{\frac{\hat{p}^{2}_{\bm{R}b\alpha}}{2M_{b}}}+\frac{1}{2}\!\!\!\sum_{\begin{subarray}{c}\scriptscriptstyle{\bm{R},b,\alpha}\\ \scriptscriptstyle{\bm{R^{\prime}}\!,b^{\prime}\!,\alpha^{\prime}}\end{subarray}}\!\!\frac{\partial^{2}{V}}{\partial{u}_{\bm{R}b\alpha}\partial{u}_{\bm{R^{\prime}}b^{\prime}\alpha^{\prime}}\!}\bigg{\lvert}_{\!\rm eq}\!\!{u}_{\bm{R}b\alpha}{u}_{\bm{R^{\prime}}b^{\prime}\alpha^{\prime}}\\ &+\frac{1}{3!}\!\!\!\sum_{\begin{subarray}{c}\scriptscriptstyle{\bm{R},b,\alpha}\\ \scriptscriptstyle{\bm{R^{\prime}}\!,b^{\prime}\!,\alpha^{\prime}}\end{subarray}}\!\!\frac{\partial^{3}{V}}{\partial{u}_{\bm{R}b\alpha}\partial{u}_{\bm{R^{\prime}}b^{\prime}\alpha^{\prime}}\!\partial{u}_{\bm{R}^{\prime\prime}b^{\prime\prime}\alpha^{\prime\prime}}\!}\bigg{\lvert}_{\!\rm eq}\!\!{u}_{\bm{R}b\alpha}{u}_{\bm{R^{\prime}}b^{\prime}\alpha^{\prime}}{u}_{\bm{R}^{\prime\prime}b^{\prime\prime}\alpha^{\prime\prime}}\,\\ &+\sum_{\scriptscriptstyle{\bm{R},b,\alpha}}\left(\frac{m_{b}}{M_{b}}-1\right){\frac{\hat{p}^{2}_{\bm{R}b\alpha}}{2M_{b}}},\end{split}$$ (1)

where $\hat{p}_{\bm{R}b\alpha}$ and $\hat{u}_{\bm{R}b\alpha}$ are the momentum and positions-displacement operators along the Cartesian direction $\alpha$ for the atom $b$ having isotope-averaged mass $M_{b}$ and position $\bm{R}{+}\bm{\tau}_{b}$ (here, $\bm{R}$ is the Bravais-lattice vector and $\bm{\tau}_{b}$ the position in the crystal’s primitive cell); the last term in Eq. (1) describes kinetic-energy perturbations induced by isotopes [40, 29] ($m_{b}$ is the exact mass of the atom at position $\bm{\tau}_{b}$, which can deviate from the isotopically-averaged mass $M_{b}$). The leading (harmonic) term in such an equation determines the vibrational frequencies; in particular, the Fourier transform of the mass-rescaled second-order derivative of the interatomic potential yields the dynamical matrix at wavevector $\bm{q}$, $\mathsfit{D}(\bm{q})_{b\alpha,b^{\prime}\!\alpha^{\prime}}{=}\sum_{\bm{R}}\frac{\partial^{2}{V}}{\partial{u}_{\bm{R}b\alpha}\partial{u}_{\bm{R}^{\prime}b^{\prime}\alpha^{\prime}}}\Big{\lvert}_{\rm eq}\frac{e^{-i\bm{q}\cdot(\bm{R}{+}\bm{\tau}_{b}{-}\bm{\tau}_{b^{\prime}})}}{\sqrt{M_{b}M_{b^{\prime}}}}$. By diagonalizing the dynamical matrix,

$$\textstyle\sum_{b^{\prime}\alpha^{\prime}}\mathsfit{D}(\bm{q})_{b\alpha,b^{\prime}\!\alpha^{\prime}}\mathcal{E}(\bm{q})_{s,b^{\prime}\!\alpha^{\prime}}=\omega^{2}(\bm{q})_{s}\mathcal{E}(\bm{q})_{s,b\alpha},$$ (2)

one obtains from the eigenvalues the phonon energies $\hbar\omega(\bm{q})_{s}$ of the solid ($s$ is a band index ranging from 1 to 3$N_{\rm at}$, where $N_{\rm at}$ is the number of atoms in the primitive cell), and from the eigenvectors $\mathcal{E}(\bm{q})_{s,b\alpha}$ the displacement patterns of atom $b$ in direction $\alpha$ for the normal mode $s$.

The third derivative in Eq. (1), instead, determines the anharmonic linewidth $\hbar\Gamma_{\rm a}(\bm{q})_{s}$ (energy broadening due to three-phonon interactions [29]) of the phonon $(\bm{q})_{s}$:

$$\begin{split}&\hbar\Gamma_{\!\!\rm a}\!(\bm{q})_{\!s}{=}\!\frac{\pi}{N_{c}^{2}}\!{\sum_{\begin{subarray}{c}\bm{q^{\prime}}\!\bm{q^{\prime}\!{}^{\prime}\!}\\ s^{\prime}\!,s^{\prime\prime}\!\end{subarray}}}\!\Big{\{}\!2\!\big{[}\bar{\mathsfit{N}}(\bm{q^{\prime}\!})_{s^{\prime}\!}{-}\bar{\mathsfit{N}}(\bm{q^{\prime}\!{}^{\prime}\!})_{s^{\prime}\!{}^{\prime}\!}\big{]}\delta\big{(}\omega(\bm{q})_{s}\!{+}\omega(\bm{q^{\prime}\!})_{s^{\prime}\!}{-}\omega(\bm{q^{\prime}\!{}^{\prime}\!})_{s^{\prime}\!{}^{\prime}\!}\big{)}\\ &+\big{[}1{+}\bar{\mathsfit{N}}(\bm{q^{\prime}\!})_{s^{\prime}\!}{+}\bar{\mathsfit{N}}(\bm{q^{\prime}\!{}^{\prime}\!})_{s^{\prime}\!{}^{\prime}\!}\big{]}\delta\big{[}\omega(\bm{q})_{s}{-}\omega(\bm{q^{\prime}\!})_{s^{\prime}\!}{-}\omega(\bm{q^{\prime}\!{}^{\prime}\!})_{s^{\prime}\!{}^{\prime}\!}\big{]}\Big{\}}\times\\ &\Bigg{|}\sum_{\begin{subarray}{c}\alpha,\alpha^{\prime}\!,\alpha^{\prime}\!{}^{\prime}\\ b,b^{\prime}\!,b^{\prime}\!{}^{\prime}\!,\bm{R^{\prime}}\!\bm{R^{\prime}\!{}^{\prime}}\end{subarray}}\frac{\partial^{3}V}{\partial u_{\bm{0}b\alpha}\partial u_{\bm{R^{\prime}\!}b^{\prime}\!\alpha^{\prime}\!}\partial u_{\bm{R^{\prime}\!{}^{\prime}}\!b^{\prime}\!{}^{\prime}\!\alpha^{\prime}\!{}^{\prime}}}\mathcal{E}(\bm{q})_{s,b\alpha}\mathcal{E}(\bm{q^{\prime}})_{s^{\prime}\!,b^{\prime}\!\alpha^{\prime}}\mathcal{E}(\bm{q^{\prime}\!{}^{\prime}})_{s^{\prime}\!{}^{\prime}\!,b^{\prime}\!{}^{\prime}\!\alpha^{\prime}\!{}^{\prime}\!}\\ &\sqrt{\frac{\hbar^{3}}{8}}\frac{\Delta(\bm{q}{+}\bm{q^{\prime}\!}{+}\bm{q^{\prime}\!{}^{\prime}\!})e^{-i[\bm{q}\cdot\bm{\tau}_{b}{+}\bm{q^{\prime}\!}\cdot(\bm{R^{\prime}\!}{+}\bm{\tau}_{b^{\prime}\!}){+}\bm{q^{\prime}\!{}^{\prime}\!}\cdot(\bm{R^{\prime}\!{}^{\prime}\!}{+}\bm{\tau}_{b^{\prime}\!{}^{\prime}\!})]}}{\sqrt{M_{b}M_{b^{\prime}\!}M_{b^{\prime}\!{}^{\prime}\!}\omega(\bm{q})_{s}\omega(\bm{q^{\prime}\!})_{s^{\prime}}\omega(\bm{q^{\prime}\!{}^{\prime}\!})_{s^{\prime}\!{}^{\prime}\!}}}\Bigg{|}^{2},\end{split}$$ (3)

where $\bar{\mathsfit{N}}(\bm{q})_{s}{=}[\exp(\hbar\omega(\bm{q})_{s}/k_{\rm B}T){-}1]^{-1}$ is the Bose-Einstein distribution at temperature $T$, $\Delta(\bm{q}{+}\bm{q^{\prime}\!}{+}\bm{q^{\prime}\!{}^{\prime}\!})$ is the Kronecker delta (equal to 1 if $\bm{q}{+}\bm{q^{\prime}\!}{+}\bm{q^{\prime}\!{}^{\prime}\!}$ is a reciprocal lattice vector, zero otherwise), $\delta$ is the Dirac delta. The last line in Eq. (1) accounts for the presence of isotopic-mass disorder and yields the following linewidth [40]

$$\begin{split}\hbar\Gamma_{\rm{i}}(\bm{q})_{s}{=}&\frac{\hbar\pi}{2N_{c}}[\omega(\bm{q})_{s}]^{2}{\textstyle\sum_{\bm{q^{\prime}},s^{\prime}}}\delta\big{[}\omega(\bm{q})_{s}{-}\omega(\bm{q^{\prime}})_{s^{\prime}}\big{]}\\ &\times\textstyle\sum_{b}g_{2}^{b}\Big{|}\textstyle\sum_{\alpha}\mathcal{E}(\bm{q})^{\star}_{s,b\alpha}\mathcal{E}(\bm{q^{\prime}})_{s^{\prime},b\alpha}\Big{|}^{2},\end{split}$$ (4)

where $g_{2}^{b}=\sum_{i}f_{i,b}\big{(}\frac{M_{b}-m_{i,b}}{M_{b}}\big{)}^{2}$ describes the variance the isotopic masses of atom $b$ ($f_{i,b}$ and $m_{i,b}$ are the mole fraction and mass, respectively, of the $i$th isotope of atom $b$; $M_{b}=\sum_{i}f_{i,b}m_{i,b}$ is the weighted average mass).

The recently developed WTE [34, 35] allows to predict the thermal conductivity of solids accounting for the interplay between structural disorder, anharmonicity, and Bose-Einstein statistics of vibrations. This offers a comprehensive approach to describe ordered ‘simple crystals’ having phonon interband spacings much larger than the linewidths [41, 42], completely disordered glasses [43, 36, 44], as well as the intermediate regime of ‘complex crystals’ with interband spacings smaller than the linewidths [35, 45, 46]. To assess how the accuracy in the prediction of the conductivity is affected by the precision with which fMLPs describe harmonic and anharmonic (third-order) force constants in Eq. (1), it is sufficient to consider the conductivity obtained from the WTE solved in the single-mode relaxation-time approximation (SMA)

$$\begin{split}&\kappa(T)=\frac{1}{\mathcal{V}N_{c}}\sum_{\bm{q},s}C(\bm{q})_{s}\frac{|\!|\bm{\mathsfit{v}}(\bm{q})_{s,s}|\!|^{2}}{3}[\Gamma(\bm{q})_{s}]^{-1}\\ &{+}\frac{1}{\mathcal{V}{N_{\rm c}}}{\sum_{\bm{q},s\neq s^{\prime}}}\frac{\omega(\bm{q})_{s}{+}\omega(\bm{q})_{s^{\prime}}}{4}\!\left[\frac{C(\bm{q})_{s}}{\omega(\bm{q})_{s}}{+}\frac{C(\bm{q})_{s^{\prime}\!}}{\omega(\bm{q})_{s^{\prime}\!}}\right]\frac{|\!|\bm{\mathsfit{v}}(\bm{q})_{s,s^{\prime}}|\!|^{2}}{3}\\ &\hskip 42.67912pt\times\frac{\frac{1}{2}\big{[}\Gamma(\bm{q})_{s}{+}\Gamma(\bm{q})_{s^{\prime}}\big{]}}{\big{[}\omega(\bm{q})_{s}{-}\omega(\bm{q})_{s^{\prime}}\big{]}^{2}+\frac{1}{4}\big{[}\Gamma(\bm{q})_{s}{+}\Gamma(\bm{q})_{s^{\prime}}\big{]}^{2}}\;;\end{split}$$ (5)

where $C(\bm{q})_{s}{=}\frac{\hbar^{2}\omega^{2}(\bm{q})_{s}}{k_{\rm B}T^{2}}\bar{\mathsfit{N}}(\bm{q})_{s}\big{[}\bar{\mathsfit{N}}(\bm{q})_{s}{+}1\big{]}$ is the specific heat of the vibration with energy $\omega(\bm{q})_{s}$ and total linewidth $\Gamma(\bm{q})_{s}{=}\Gamma_{\!\!\rm a}(\bm{q})_{s}{+}\Gamma_{\!\!\rm i}(\bm{q})_{s}$, $\mathsfit{v}^{\beta}(\bm{q})_{s,s^{\prime}}{=}\!\sum_{\begin{subarray}{c}b,\alpha,b^{\prime}\!,\alpha^{\prime}\end{subarray}}\!\mathcal{E}^{\star}(\bm{q})_{{s},b\alpha}[{{\nabla^{\beta}_{\bm{q}}\sqrt{\mathsfit{D}(\bm{q})}}}_{b\alpha,b^{\prime}\alpha^{\prime}}]\mathcal{E}(\bm{q})_{s^{\prime},b^{\prime}\alpha^{\prime}}$ is the velocity operator coupling eigenstates $s$ and $s^{\prime}$ at the same wavevector $\bm{q}$ (its diagonal elements $s=s^{\prime}$ are the usual phonon group velocities) [35], $N_{\rm c}$ is the number of $\bm{q}$-points used to sample the Brillouin zone and $\mathcal{V}$ is the crystal’s primitive-cell volume. The first line on the right-hand side of Eq (5) describes a conduction mechanism in which vibrations carry the heat $C(\bm{q})_{s}$ by propagating particle-like with velocity ${\bm{\mathsfit{v}}(\bm{q})_{s,s}}$ over the lifetime $[\Gamma(\bm{q})_{s}]^{-1}$. It can be rigorously shown [34, 35] that it coincides with the conductivity emerging from the Peierls-Boltzmann equation [39], and will be henceforth referred to as $\kappa_{\rm P}$. The term on the second and third lines of Eq. (5) accounts for conduction through a wave-like tunneling mechanism between pairs of vibrational eigenstates. This term arises from the coherence between two different phonon modes $s,s^{\prime}$ at the same wavevector $\bm{q}$ (i.e., it becomes more signicant as their frequency difference $\omega(\bm{q})_{s}-\omega(\bm{q})_{s^{\prime}}$ becomes smaller) and is therefore referred to as ‘coherences conductivity’, $\kappa_{\rm C}$. It has been shown in Refs. [34, 35, 47] that in simple crystals $\kappa_{{\rm P}}{\gg}\kappa_{{\rm C}}$ because particle-like propagation dominates over the wave-like tunneling. In contrast, in complex crystals both these mechanisms co-exist and can have comparable magnitude, implying $\kappa_{{}_{\rm P}}{\sim}\kappa_{{}_{\rm C}}$. Finally, Refs. [36, 44] have shown that in strongly disordered oxide glasses $\kappa_{{\rm P}}{\ll}\kappa_{{\rm C}}$.

Vibrational and thermal properties from fMLPs
In general, $\kappa(T)$ is highly sensitive to both harmonic and anharmonic vibrational properties [35, 45, 46, 10, 47, 48, 49, 37, 50]; therefore, the accuracy of fMLPs in describing these properties can be quantified by comparing their predictions for $\kappa(T)$ with those obtained from reference first principles data. To accomplish this goal, we begin by calculating the harmonic and third-order anharmonic force constants that determine the solid’s vibrational energy (Eq. (1)). We do this by using either: (i) Density-Functional Theory (DFT) with the Perdew, Burke, and Ernzerhof (PBE) functional (see Refs. [30, 31] for details); or (ii) SOTA non-proprietary mp-fMLPs trained on DFT-PBE. Then, we employ these force constants to determine harmonic vibrational properties (frequencies, velocity operators, and isotopic linewidths, see Eq. (2) and Eq. (4)) and anharmonic linewidths (Eq. (3)). Subsequently, these atomistic vibrational properties are used in Eq. (5) to calculate the thermal conductivity of the 103 diverse compounds in the phononDB-PBE database. The database contains rocksalt, zincblende, and wurtzite binary compunds and involves 34 chemical species, including alkali metals (Cs, K, Li, Na, Rb), alkaline earth metals (Be, Mg, Ca, Sr, Ba), transition metals (Ag, Cu, Zn, Cd), post-transition metals (Al, Ga, In, Pb), metalloids (As, Sb, Si, Te, B), nonmetals (H, C, N, O, P, S, Se), and halogens (F, Cl, Br, I). Computational details are provided in the Methods.

Figure 1: Thermal conductivity computed at 300 K from DFT-PBE or MACE-MP-0, for 103 binary compounds taken from the phononDB-PBE database, which have rocksalt (green), zincblende (orange), and wurtzite (blue) structure. Solid lines indicate perfect agreement, while dashed lines represent discrepancies of up to a factor of 2. The arrows highlight three selected materials, which will be analyzed in detail in terms of their microscopic vibrational properties later in the manuscript. The inset displays the distribution of relative deviations between the conductivities predicted by DFT-PBE and MACE-MP-0.

Figure 2: Phonon bands, specific heat at constant volume, phonon frequency-linewidth distributions, and macroscopic thermal conductivity, for wurtzite BeO (panel a), zincblende BeTe (panel b), and rocksalt RbH (panel c). DFT-PBE values are computed using the data from Refs. [30, 31], while ‘MACE’ values are computed using the MACE-MP-0 fMLP[20]. For the DFT-PBE calculations, phonons are plotted considering the long-range non-analytical correction term (NAC, red lines or scatter points) or not (dotted green); MACE-MP-0 does not account for long-range contributions to the vibrational energy (see text). In the frequency-linewidths distributions at 300 K, we resolve the anharmonic ($\Gamma_{a}$) and isotopic ($\Gamma_{i}$) part of the linewidths; the areas of the markers are proportional to the contribution of the given phonon mode to the conductivity (see Eq. (6)). Specific heat and thermal conductivity are negligibly affected by considering or not the long-range NAC term.

In Fig. 1 we compare the thermal conductivity predicted from first principles (DFT-PBE) or from the mp-fMLP MACE-MP-0 (trained on DFT-PBE data) [20]. We highlight how predictions from DFT-PBE or MACE-MP-0 are generally consistent within a factor of 2 across most materials, see inset. We note that in about 20% of the compounds studied the semiclassical particle-like BTE fails to fully describe heat transport, and it is crucial to employ the more general WTE, see Fig. Thermal Conductivity Predictions with Foundation Atomistic Models in the Methods. These cases generally have ultra-low conductivity ($\kappa_{\rm TOT}\lesssim 2$ W/mK), and are characterized by having wave-like tunneling conductivity ($\kappa_{C}$, accounted for by the WTE but missing from the BTE) with magnitude comparable to the particle-like propagation conductivity ($\kappa_{P}$, accounted for by both the WTE and BTE).

To understand the microscopic origin of the discrepancies in the macroscopic DFT-PBE and MACE-MP-0 conductivities, we select three representative materials — wurtzite BeO, zincblende BeTe, and rocksalt RbH — and show in Fig. 2 their phonon band structures, specific heat at constant volume, and the macroscopic thermal conductivity resolved in terms of contributions from microscopic phonon modes. Examining the phonon band structures in the first column, we observe that MACE-MP-0 tends to underestimate the high-frequency optical vibrational modes compared to DFT-PBE. To further investigate these differences, we plot the DFT-PBE phonons considering the non-analytical correction term [51] (NAC, red lines or scatter points), or not (dotted green line). This long-range interaction is responsible for the energy splitting between the longitudinal-optical and transversal-optical modes in polar dielectric materials [51], and is not fully considered in fMLPs trained using a radial force cutoff (e.g., for MACE-MP-0 such cutoff is 6 Å, while it is 5 Å  for SevenNet[21], CHGNet[19], and M3GNet[18]). This explains why the DFT-PBE phonons without NAC are in closer agreement with the MACE-MP-0 phonons (blue lines or scatter points). Importantly, it is noticeable that even without NAC, DFT-PBE frequencies tend to be higher than MACE-MP-0 frequencies, confirming the general tendency of MACE-MP-0 to underestimate vibrational frequencies discussed in Ref. [27]. Nevertheless, Fig. 2 shows that considering or not the NAC term has a negligible impact on the specific heat at constant volume and on the temperature-dependent conductivity ($\kappa(T)$) of BeO, BeTe and RbH. This can be understood from the third column of Fig. 2, where we show that thermal transport in these materials is dominated by low-energy (acoustic) phonons, and these are negligibly affected by NAC. Specifically, these plots report the frequency-linewidth distributions, resolving both the anharmonic and isotopic contributions to the linewidths (Eq. 3 and Eq. 4, respectively), and also quantifying how much a single phonon mode contributes to the total conductivity with the following expression:

$$\begin{split}&{\mathcal{K}}(\bm{q})_{s}{=}C(\bm{q})_{s}\frac{|\!|\bm{\mathsfit{v}}(\bm{q})_{s,s}|\!|^{2}}{3}[\Gamma(\bm{q})_{s}]^{-1}\\ &+\sum_{s^{\prime}\neq s}\!\frac{C(\bm{q})_{s}}{C(\bm{q})_{s}{+}C(\bm{q})_{s^{\prime}\!}}\frac{\omega(\bm{q})_{s}{+}\omega(\bm{q})_{s^{\prime}\!}}{2}\!\!\left[\frac{C(\bm{q}_{s}}{\omega(\bm{q})_{s}}{+}\frac{C(\bm{q}_{s^{\prime}\!}}{\omega(\bm{q})_{s^{\prime}\!}}\right]\\ &\times\frac{|\!|\bm{\mathsfit{v}}(\bm{q})_{s,s^{\prime}}|\!|^{2}}{3}\frac{\frac{1}{2}\left[\Gamma(\bm{q})_{s}{+}\Gamma(\bm{q})_{s^{\prime}\!}\right]}{[\omega(\bm{q})_{s^{\prime}\!}{-}\omega(\bm{q})_{s}]^{2}{+}\frac{1}{4}[\Gamma(\bm{q})_{s}{+}\Gamma(\bm{q})_{s^{\prime}\!}]^{2}}.\end{split}$$ (6)

This equation describes how much a single phonon mode $(\bm{q})_{s}$ contributes to heat conduction. Specifically, the term on the first line accounts for the aforementioned particle-like conduction mechanism; the term on the second and third lines, instead, describes the wave-like conduction that originates from tunneling between two non-degenerate phonons $(\bm{q})_{s}$ and $(\bm{q})_{s^{\prime}}$ — here the single-phonon contributions from $(\bm{q})_{s}$ is resolved using the ratio between specific heats $\tfrac{C(\bm{q})_{s}}{C(\bm{q})_{s}+C(\bm{q})_{s^{\prime}}}$ as weight, as discussed by Eq. (E3) in [35].

The frequency-linewidth distributions in Fig. 2a show overall agreement between DFT-PBE and MACE-MP-0 for the phonon modes that are mainly contributing to conduction — the acoustic modes and the optical modes with $\hbar\omega(\bm{q})_{s}<k_{B}T$ — and we see that this implies very similar values for the corresponding macroscopic $\kappa(T)$ (with a relative difference $<10\%$). Importantly, Fig. 2b shows compatibility between $\kappa(T)$ predicted from DFT-PBE or MACE-MP-0 can also result from cancellation of errors; specifically, in BeTe there are visible differences in the frequency-linewidth distributions, and these largely cancel out when integrated to determine the conductivity. Finally, Fig. 2c illustrates that the presence of systematic, non-compensating differences in the microscopic frequency-linewidth distributions can also directly translate into significant discrepancies on the macroscopic conductivities (difference of a factor of 2.7). Additionally, we note that depending on the the chemical composition, the anharmonic linewidths at room temperature can either dominate over the isotopic linewidth (e.g., in BeO) or not (e.g., in BeTe). In the former case, the thermal conductivity is practically unaffected by considering or not the isotopic contributions to the linewidth, while in the latter case considering isotopic scattering yields a reduction of $\kappa(300\,\rm{K})$ of about ${\sim}75\%$ (see Table. 2).

The results above motivate us to investigate when good agreement between DFT and fMLP conductivities is obtained from accurately described microscopic harmonic and anharmonic vibrational properties, or because of compensation of errors. To achieve this goal, we resolve the discrepancies in the total macroscopic thermal conductivities using the Symmetric Relative Error (SRE),

$$\text{SRE}\big{[}\kappa\big{]}=2\frac{\left|{\kappa}_{\rm{fMLP}}-{\kappa}_{\rm{DFT}}\right|}{{\kappa}_{\rm{fMLP}}+{\kappa}_{\rm{DFT}}}\,.$$ (7)

Then, to determine whether a low value for $\text{SRE}\big{[}\kappa\big{]}$ originates from accurately described microscopic vibrational properties, or because of compensation of errors, we introduce the Symmetric Relative Mean Error (SRME) on the single-phonon conductivity contribution $\mathcal{K}(\bm{q})_{s}$:

$$\begin{split}\text{SRME}\big{[}\{\mathcal{K}(\bm{q})_{s}\}\big{]}{=}\frac{2}{N_{c}\mathcal{V}}\frac{\sum_{\bm{q}s}\!\left|{\mathcal{K}}_{\rm{fMLP}}\!(\bm{q})_{s}{-}{\mathcal{K}}_{\rm{DFT}}\!(\bm{q})_{s}\right|}{{\kappa}_{\rm{fMLP}}+{\kappa}_{\rm{DFT}}}\,,\;\;\end{split}$$ (8)

where ${\mathcal{K}}_{\rm{DFT}}\!(\bm{q})_{s}$ refers to Eq. (8) evaluated using DFT, and ${\mathcal{K}}_{\rm{fMLP}}\!(\bm{q})_{s}$ refers to the same equation evaluated using fMLP.

Figure 3: Errors on conductivity: SRE & SRME, computed comparing DFT-PBE and fMLP (MACE-MP-0) predictions in 103 chemically and structurally diverse compounds. Large discrepancies between DFT-PBE and fMLP on microscopic, single-phonon conductivity contributions are described by large SRME (Eq. (8)), which implies a significant discrepancy between the macroscopic DFT-PBE and fMLP conductivities, and thus large SRE (Eq. (7)). Three materials are highlighted: wurtzite BeO shows low SRME and thus low SRE; rocksalt RbH displays high error in both SRME and SRE; zincblende BeTe exhibits high SRME but low SRE because of compensation of microscopic errors.

Fig. 3 illustrates that a large SRME generally implies large SRE. Importantly, knowing both SRE and SRME enables us to identify when microscopic error compensation occurs — this is indicated by a large SRME but small SRE. Three representative cases are highlighted in Fig. 3: wurtzite BeO shows low SRME and thus correspondingly low SRE; rocksalt RbH exhibits high error in both SRME and SRE; zincblende BeTe displays high SRME but low SRE due to compensation of microscopic errors. We note that the SRME[$\mathcal{K}(\bm{q})_{s}$] error can stem from discrepancies in the harmonic (second-order) or anharmonic (third-order) force constants. Therefore, in Fig. 6 in the Methods we discuss how to disentangle the SRME and SRE into errors on harmonic vibrational properties and anharmonic vibrational properties.

Overall, this analysis highlights the importance of benchmarking the accuracy of fMLPs not only at the macroscopic but also at the microscopic level. In particular the macroscopic SRE[$\kappa$] alone is not a reliable descriptor for the fMLPs’ ability to capture the harmonic and anharmonic physics underlying heat conduction. However, achieving both a small macroscopic SRE[$\kappa$] and a small microscopic SRME[$\kappa$] is a sufficient condition for accurately describing both the microscopic physics and the macroscopic thermal conductivity.

Accuracy of various SOTA foundation models
The previous section demonstrated that the SRE and SRME descriptors can effectively measure the accuracy of fMLPs in describing vibrational and thermal properties. Therefore, here we utilize these descriptors to compare the accuracy of different non-proprietary mp-fMLPs on the 103 diverse structures contained in the phononDB-PBE database[30, 31]. Fig. 4 compares the accuracy of the four foundation models M3GNet[18], CHGNet[19], MACE-MP-0[20], and SevenNet[21] mp-fMLPs. To determine if a certain mp-fMLP tend to systematically overestimate or underestimate the conductivity, it is useful to rely on the Symmetric Relative Difference (SRD) in the total Wigner conductivity from DFT-PBE or mp-fMLP,

$$\rm{SRD}\big{[}\kappa\big{]}=2\frac{\kappa^{\rm fMLP}-\kappa^{\rm DFT}}{\kappa^{\rm fMLP}+\kappa^{\rm DFT}},$$ (9)

which ranges from -2 to +2, resolving both overestimation and underestimation of the conductivity. Figure  4a summarizes the SRD for all 103 binary compounds in the phononDB-PBE database, as depicted in the violin plot [52]. In this plot, the width of each violin shape is related to the percentage of materials with SRD values indicated on the y-axis. The median SRD value is marked by a white scatter point, while the black boxes illustrate the interquartile range, and the whiskers show $1.5$ times the interquartile range below and above the first and third quartile points. Unstable structures with negative phonon frequencies were included considering $\kappa^{\rm fMLP}=0$, i.e., SRD=$-2$. This analysis reveals that all the mp-fMLPs assessed in this study tend to underestimate the thermal conductivity. This finding is consistent with the observed systematic underestimation of vibrational frequencies noted in Ref. [27], and might also derived from overestimation of anharmonic linewidths.

Figure 4: Performance of mp-fMLPs in calculating thermal conductivity. Panel a presents violin plots depicting the Symmetric Relative Difference (SRD) in the total Wigner thermal conductivity (Eq. (9)) for mp-fMLPs SevenNet[21], MACE-MP-0[20], CHGNet[19], M3GNet[18]. Panel b shows the SRME for mode-resolved thermal conductivity. In both panels, the medians are marked with white scatter points, the widths of the boxes represent the interquartile range, and the whiskers show the range of data points without outliers.

To examine whether the SRD shown in Fig. 4a is influenced by error compensation, we analyze in Fig. 4b the violin plots for $\text{SRME}\big{[}\{\mathcal{K}(\bm{q})_{s}\}\big{]}$. We see that the SRME distribution for MACE-MP-0 is the closest to zero, followed by SevenNet, M3GNet, and CHGNet. In addition, to quantify the overall accuracy of a fMLPs in predicting the macroscopic conductivity over a materials’ database, without resolving the possible compensation of microscopic errors, it is informative to consider the mean of the modulus of the deviations, i.e., the mean of the distribution of SRE (7) — we prefer this over the mean of the SRD distribution, as the latter can be close to zero in the presence of very broad but symmetric distribution. Importantly, we note that the mean for SRE[$\kappa$] and mean for SRME[$\mathcal{K}(\bm{q})_{s}$] are expected to be comparable in the absence of compensation of microscopic errors. In contrast, in the presence of compensation of microscopic errors, the mean SRE[$\kappa$] is significantly lower than mean SRME[$\mathcal{K}(\bm{q})_{s}$]. The mean values for SRME[$\mathcal{K}(\bm{q})_{s}$] and SRE[$\kappa$] are reported in Table Thermal Conductivity Predictions with Foundation Atomistic Models.