## Thermal Conductivity Predictions with Foundation Atomistic Models

Balázs Póta, 1 Paramvir Ahlawat, 1 Gábor Csányi, 2 and Michele Simoncelli 1, ∗

1 Theory of Condensed Matter Group of the Cavendish Laboratory, University of Cambridge, Cambridge, UK 2 Engineering Laboratory, University of Cambridge, Cambridge, UK

Advances in machine learning have led to the development of foundation models for atomistic materials chemistry, enabling quantum-accurate descriptions of interatomic forces across diverse compounds at reduced computational cost. Hitherto, these models have been benchmarked relying on descriptors based on atoms' interaction energies or harmonic vibrations; their accuracy and efficiency in predicting observable and technologically relevant heat-conduction properties remains unknown. Here, we introduce a framework that leverages foundation models and the Wigner formulation of heat transport to overcome the major bottlenecks of current methods for designing heat-management materials: high cost, limited transferability, or lack of physics awareness. We present the standards needed to achieve first-principles accuracy in conductivity predictions through model's fine-tuning, discussing benchmark metrics and precision/cost trade-offs. We apply our framework to a database of solids with diverse compositions and structures, demonstrating its potential to discover materials for next-gen technologies ranging from thermal insulation to neuromorphic computing.

Over the past decades, intense research efforts have been dedicated to tackle the challenging task of fitting the Born-Oppenheimer potential energy surface as a function of atomic coordinates 1-10 . These resulted in the development of machine-learning potentials (MLPs), which allow us to reproduce atoms' interaction energies and forces with nearly quantum accuracy and computational cost orders-of-magnitude lower than first-principles methods. These advances have enabled scientists to make progress in understanding how atomistic structure and composition influence macroscopic observables - an example is the thermal conductivity, a property determined by atomic structure and forces, and relevant for technologies ranging from energy harvesting 11 to efficient neuromorphic computing 12 . Major drawbacks of methods based on MLPs are the significant work required to generate a first-principles database needed to train and validate MLPs, 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 13,14 or macroscopic 15-17 materials' properties very efficiently, but without necessarily accounting for the fundamental physics that determines them. Fittingly, recent work 18 has formally demonstrated the possibility to obtain a physics-aware and mathematically complete description of atomic environments (and thus of forces) by employing Message Passing Networks 19 with many-body messages 20 (with the MACE architecture 20 this can be achieved, e.g., using 4-body terms and one message pass). This breakthrough has enabled the development of foundation machine-learning potentials (fMLPs) 21-33 , which are trained across nearly all chemical elements and can be combined to describe almost arbitrary compounds. Recent works have assessed the accuracy of fMLPs in predicting, e.g., atoms' interaction energies (thermodynamic stability 34 ) or harmonic vibrations 35-37 . However, the accuracy of fMLPs in pre- dicting anharmonic and heat-conduction properties remains unexplored, due to the complexity of the relation between these properties and interatomic forces 38,39 .

Here we introduce a framework that leverages fMLPs and the Wigner formulation of heat transport 38,39 to predict with first-principles accuracy anharmonic vibrational properties and thermal conductivity in solids with arbitrary compositions and structures. This overcomes the major bottlenecks of methods to predict conductivity based on standard MLPs, specifically their lack of transferability across different compositions, and the complex and data-intensive MLP training procedure. In particular, after discussing how to calculate harmonic and anharmonic vibrational properties of solids using fMLPs, we rely on these quantities to determine the thermal conductivity with the recently developed Wigner heatTransport Equation (WTE) 38,39 . The WTE unifies the established theories for heat transport in crystals 40 and glasses 41 , comprehensively describing solids with arbitrary composition or structural disorder. We showcase the capabilities of this combined fMLP-WTE framework in all the state-of-the-art (SOTA) non-proprietary fMLPs trained on the Materials Project database 42 and hereafter collectively referred to as 'mp-fMLPs' - M3GNet 21 , CHGNet 22 , MACE-MP-0 23 , SevenNET 24 , and ORBMPtrj 43 . We benchmark their zero-shot predictions for vibrational and thermal properties against a database of first-principles reference data for 103 chemically and structurally diverse compounds 44,45 . We rely on these benchmarks to introduce precision metrics for anharmonic vibrational and thermal properties, and we present a fine-tuning protocol that achieves first-principles accuracy in materials where fMLPs' zero-shot precision is poor. We demonstrate that our framework can efficiently perform physics-aware and quantum-accurate predictions of heat-conduction properties also in materials that violate semiclassical Boltzmann transport, relevant for energy or information-management technologies.

## Atomic vibrations &amp; thermal conductivity

The recently developed WTE 38,39 offers a comprehensive approach to predict from atomistic vibrational properties the thermal conductivity of anharmonic crystals 46 , disordered glasses 47 , as well as the intermediate case of anharmonic-and-disordered 'complex crystals' 48 . To assess how the accuracy in the conductivity prediction is affected by the precision with which fMLPs describe atomic vibrations, it is sufficient to consider the directionally averaged trace of the conductivity tensor obtained from the WTE solution in the relaxation-time approximation 39

̸

$$4$$

$$\text{WTE solution in the relaxation-time approximation} ^ { 3 9 } & \text{ and } \Gamma \\ \kappa ( T ) & = \kappa _ { P } ( T ) + \kappa _ { C } ( T ) = \frac { 1 } { \mathcal { V } N _ { c } } \sum _ { q, s } \left [ \begin{array} { c } & \text{and } \Gamma \\ & \text{tional} \\ & & \text{relying} \\ & & C ( q ) _ { s } \frac { \| \mathbf v ( q ) _ { s, s } \| ^ { 2 } } { 3 } [ \Gamma ( q ) _ { s } ] ^ { - 1 } & ( 1 ) & \text{Errze} \\ & + \sum _ { s ^ { \prime } \neq s } \frac { C ( q ) _ { s } & \frac { \omega ( q ) _ { s } + \omega ( q ) _ { s ^ { \prime } } } { 2 } \left ( \frac { C ( q _ { s } } { \omega ( q ) _ { s } } + \frac { C ( q _ { s ^ { \prime } } } { \omega ( q ) _ { s ^ { \prime } } } \right ) } { 2 } \left ( \frac { C ( q _ { s } } { \omega ( q ) _ { s } } + \frac { C ( q _ { s ^ { \prime } } } { \omega ( q ) _ { s ^ { \prime } } } \right ) } { \quad \text{these} \colon } & \text{to call} \\ & \times \frac { \| \mathbf v ( q ) _ { s, s ^ { \prime } } \| ^ { 2 } } { 3 } \frac { \frac { 1 } { 2 } \left [ \Gamma ( q ) _ { s } + \Gamma ( q ) _ { s ^ { \prime } } \right ] } { \left [ \omega ( q ) _ { s ^ { \prime } - \omega ( q ) _ { s } \right ] ^ { 2 } + \frac { 1 } { 4 } \left [ \Gamma ( q ) _ { s } + \Gamma ( q ) _ { s ^ { \prime } } \right ] ^ { 2 } } \right ]. & \text{pound} \\ \text{Hero} & \text{$C(\omega\colon = \frac{h^{2} \omega^{2} ( q ) _ { s } \bar { \mathcal { U } } ( \omega \wedge \, \left \lceil \bar { \mathcal { U } } ( \omega \, + 1 ) \text{ is the non-offic has} } \end{array}$$

Here, C ( q ) s = ℏ 2 ω 2 ( q ) s k B T 2 ¯ ( N q ) s [ ¯ ( N q ) s +1 ] is the specific heat at temperature T of the vibration having: wavevector q ; mode s ; energy ω ( q ) s ; population given by BoseEinstein distribution ¯ ( N q ) s =[exp( ℏ ω ( q ) s /k B T ) -1] -1 ; total linewidth Γ( q ) s =Γ ( a q ) s +Γ ( i q ) s , where Γ ( a q ) s depends on T and is determined by anharmonic third-order derivatives of the interatomic potential, while Γ ( i q ) s is T -independent and regulated by concentration of isotopes (see Methods for details). v β ( q ) s,s ′ is the velocity operator obtained from wavevector differentiation of the solid's dynamical matrix 39 , its diagonal elements s = s ′ are the phonon group velocities and its off-diagonal elements describe couplings between modes s and s ′ at the same q ; N c is the number of q -points sampling the Brillouin zone and V is the crystal's unit-cell volume.

Eq. (1) shows in the first line that the total macroscopic conductivity is determined by two contributions: the particle-like conductivity ( κ P ), and the wave-like 'coherence' conductivity ( κ C ). In such an equation, the term inside the square brackets is the single-phonon contribution to the total conductivity, K ( q ) s , which, as we will show later, is useful for benchmarking fMLPs. The second line shows that κ P originates from phonons that carry heat C ( q ) s by propagating particle-like with velocity v ( q ) s,s over lifetime [Γ( q ) s ] -1 ; it can be shown 38 that κ P is equivalent to the conductivity described by the phonon Boltzmann Transport Equation (BTE). The term on the third and fourth lines, instead, shows that κ C describes conduction through wave-like tunneling between pairs of phonons s, s ′ at wavevector q . It has been shown 38,39 that in simple crystals particlelike propagation dominates over wave-like tunneling ( κ P ≫ κ C ), in complex crystals both these mechanism co-exist ( κ P ∼ κ C ), and in strongly disordered glasses tunneling dominates ( κ P ≪ κ C ) 47 .

## Vibrational and thermal properties from fMLPs

In general, κ T ( ) is highly sensitive to both harmonic and anharmonic vibrational properties 39,49-51 . Therefore, as a first benchmark test for fMLPs, it is natural to compare their predictions for κ T ( ) with those obtained from reference first-principles data. To this end, we begin by calculating the harmonic ( ℏ ω ( q ) s , v β ( q ) s,s ′ , and Γ ( i q ) s ) and third-order anharmonic ( Γ ( a q ) s ) vibrational properties that determine κ T ( ) . We do this by: (i) relying on the Density-Functional Theory (DFT) data of Refs. 44,45 , calculated using the Perdew, Burke, and Ernzerhof (PBE) functional; or (ii) using SOTA nonproprietary mp-fMLPs trained on DFT-PBE 42 . Then, these atomistic vibrational properties are used in Eq. (1) to calculate the conductivity of the materials in the phononDB-PBE database 44,45 -103 diverse binary compounds, involving 34 different chemical species and having rocksalt, zincblende, or wurtzite structure.

In Fig. 1 we compare the conductivity predicted

FIG. 1. Thermal conductivity computed at 300 K from DFT-PBE and MACE-MP-0 large, for 103 compounds taken from the phononDB-PBE database, which have rocksalt (green), zincblende (orange), or wurtzite (blue) structure. The solid line indicates perfect agreement and the dashed lines discrepancies within a factor of 2. The color gradient shows the transition from the WTE regime where heat transport occurs through both propagation and tunneling mechanisms (red) to the semiclassical propagation-dominated BTE regime (white). The arrows highlight materials discussed in detail later. Inset, distribution of relative deviations between the conductivities predicted by DFT-PBE and MACE-MP-0.

![image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00755v4-figure-597874a9-2ee4-4118-9893-43ba5355a2c8.jpg)

a) wurtzite BeO

FIG. 2. Phonon bands, specific heat at constant volume, energy-linewidth distributions, and conductivity of wurtzite BeO (panel a ), zincblende BeTe ( b ), and rocksalt LiBr ( c ). DFT-PBE values are computed using the data from Refs. 44,45 , while 'MACE' values are computed using the MACE-MP-0 fMLP 23 . For the DFT-PBE phonons, we show the effect of considering (red) or not (green) the long-range non-analytical contribution (NAC); MACE-MP-0 does not fully account for NAC (see text). In the frequency-linewidths distributions at 300 K, we resolve the anharmonic ( Γ a ) and isotopic ( Γ i ) parts of the linewidths; the markers' areas are proportional to the contribution of the given phonon mode to the conductivity.

![image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00755v4-figure-a266b7d5-5ca3-4547-a373-e369d39269bd.jpg)

from first principles (DFT-PBE) or from the mp-fMLP MACE-MP-0 large (trained on PBE) 23 , underscoring that predictions are generally compatible within a factor of 2. We also highlight that in about 20% of the compounds studied, the semiclassical particle-like BTE fails to fully describe heat transport, and in these cases it is crucial to employ the more general WTE to perform accurate predictions. These cases have very low conductivity ( κ TOT ≲ 2 W/mK), which receives significant contributions from both wave-like tunneling (described by the WTE but missing from the BTE) and particlelike propagation (described by both the WTE and BTE) heat-transport mechanism, see Methods Fig. 6.

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 LiBr - and show in Fig. 2 their phonon band structures, specific heat at constant volume, and macroscopic conductivity resolved in terms of contributions from microscopic phonon modes. Examining the phonon band structures in the first column, we observe that MACEMP-0 tends to underestimate the high-frequency opti- cal vibrational modes compared to DFT-PBE. To further investigate these differences, we plot the DFT-PBE phonons considering or not the non-analytical correction term 52 (NAC). The NAC long-range term is responsible for the energy splitting between the longitudinal-optical and transversal-optical modes in polar dielectrics 52 , and is not fully considered in fMLPs trained using a radial cutoff (e.g., for MACE-MP-0 23 such cutoff is 6 Å, while it is 5 Å for SevenNet 24 , CHGNet 22 , and M3GNet 21 ). This explains why the DFT-PBE phonons without NAC are in closer agreement with the MACE-MP-0 phonons. Importantly, it is noticeable that even without NAC, DFTPBE energies tend to be higher than MACE-MP-0 energies, confirming the general tendency of MACE-MP-0 to underestimate phonon energies 36 . Nevertheless, Fig. 2 shows that considering or not the NAC term has a negligible impact on the specific heat at constant volume and on κ T ( ) . This can be understood from the third column of Fig. 2, where we show that in these materials a significant amount of heat is carried by low-energy phonons that are negligibly affected by NAC. Specifically, we report the energy-linewidth distributions, resolving the anharmonic and isotopic parts of the linewidths ( Γ ( a q ) s

and Γ ( i q ) s , respectively), and quantifying how much each phonon contributes to the conductivity with the singlemode conductivity contribution K ( q ) s defined in Eq. (1).

The energy-linewidth distributions in Fig. 2 a show overall agreement between DFT-PBE and MACE-MP-0 for the phonon modes that are most relevant for conduction, and we see that this implies macroscopic κ T ( ) compatible within 7% . Importantly, we show in the Methods (Fig. 9) that these fMLP conductivity predictions have a computational cost more than 3 orders of magnitude lower than the corresponding DFT predictions.

Fig. 2 b shows that compatibility between κ T ( ) predicted from DFT or fMLP can also result from cancellation of errors. E.g., in BeTe there are visible differences in the energy-linewidth distributions, and these largely cancel out when integrated to determine the conductivity. Finally, Fig. 2 c illustrates that discrepancies in the energy-linewidth distributions do not always compensate, and can translate into significant differences (a factor of 2) on κ T ( ) . We also note that depending on the chemical composition, the anharmonic linewidths at room temperature can dominate over the isotopic linewidth (e.g., in BeO) or not (e.g., in BeTe).

The results above motivate us to investigate whether good agreement between DFT and fMLP conductivities is obtained as a result of accurately described microscopic

FIG. 3. Errors on conductivity: SRE &amp; SRME in 103 chemically and structurally diverse compounds. Large discrepancies between DFT-PBE and mp-fMLP (MACE-MP-0) on microscopic, single-phonon conductivity contributions are described by large SRME (3), which implies a large SRE (2). We highlight how in BeO a low SRME implies a low SRE; LiBr displays high SRME and high SRE; BeTe exhibits high SRME but low SRE because of compensation of microscopic errors. The inset shows materials with very low SRE.

![image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00755v4-figure-ac24b8c0-9077-4e12-92a4-b3d7c30c7a81.jpg)

harmonic and anharmonic vibrational properties, or because of compensating errors. To study this, we quantify the discrepancies in the total macroscopic thermal conductivities using the Symmetric Relative Error (SRE),

$$S R E [ \kappa ] = 2 \frac { | \kappa _ { f M L P } - \kappa _ { D F T } | } { \kappa _ { f M L P } + \kappa _ { D F T } } \,. \quad \quad ( 2 )$$

Next, to quantify the error on the microscopic phonon properties, we introduce the Symmetric Relative Mean Error (SRME) on the single-phonon conductivity:

$$\underset { \text{meta} } { \text{pre} } ^ { \text{pre} } _ { \text{ell} a } \underset { \text{meta} } { \text{sRME} } \left [ \left \{ \mathcal { K } ( \boldsymbol q ) _ { s } \right \} \right ] = & \frac { 2 } { N _ { c } \mathcal { V } } \frac { \sum _ { \boldsymbol q s } | \mathcal { K } _ { \text{fMLP} } ( \boldsymbol q ) _ { s } - \mathcal { K } _ { \text{DFT} } ( \boldsymbol q ) _ { s } | } { \kappa _ { \text{fMLP} } + \kappa _ { \text{DFT} } } \,, ( 3 )$$

where K DFT ( q ) s refers to the single-phonon conductivity (term inside square brackets in Eq. (1)) evaluated using DFT, and K fMLP ( q ) s refers to the same expression evaluated using fMLP. 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. We note that the SRME[ K ( q ) s ] error can stem from discrepancies in the harmonic (second-order) or anharmonic (third-order) force constants. Therefore, in the Methods (Fig. 7) we discuss how to decompose the SRME and SRE into errors on harmonic and anharmonic vibrational properties.

## Accuracy of various SOTA foundation models

Here we utilize the SRE and SRME descriptors to study the accuracy of non-proprietary mp-fMLPs on the structures contained in the phononDB-PBE database 44,45 , comparing in Fig. 4 M3GNet 21 , CHGNet 22 , MACE-MP0 23 , SevenNet 24 , and ORB-MPtrj 43 . To determine if a particular mp-fMLP tends to systematically overestimate or underestimate the conductivity, we rely on the Symmetric Relative Difference (SRD) in the total Wigner conductivity from DFT-PBE or mp-fMLP,

$$S R D [ \kappa ] = 2 \frac { \kappa ^ { \text{fMLP} } - \kappa ^ { \text{DFT} } } { \kappa ^ { \text{fMLP} } + \kappa ^ { \text{DFT} } }, \quad \quad ( 4 )$$

which ranges from -2 to +2, resolving both overestimation and underestimation of the conductivity. Figure 4 a summarizes in violin plots the SRD for the compounds in the phononDB-PBE database. Unstable structures with imaginary phonon frequencies were included considering κ fMLP = 0 , i.e. , SRD= -2 . This analysis reveals that the most accurate mp-fMLPs tend to underestimate the thermal conductivity, as a consequence of underestimating vibrational frequencies 36 , and in some cases overestimating anharmonic linewidths.

To examine whether the SRD in Fig. 4 a is influenced by error compensation, we show in Fig. 4 b violin plots for SRME [ {K ( q ) s } ] . The most accurate mp-fMLPs, MACE-MP-0 23 , yields zero-shot conductivities compati-

FIG. 4. Performance of mp-fMLPs in predicting thermal conductivity. a , violin plots for the Symmetric Relative Difference (SRD) in the total Wigner conductivity (4) for the materials in phononDB-PBE simulated with ORB-MPtraj 43 , SevenNet 24 , MACE-MP-0 23 , CHGNet 22 , M3GNet 21 . b , SRME in the same materials database. 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. The pie charts in c display the number of compounds with non-negative frequencies that retained or not the correct crystal symmetry after unconstrained relaxation (green and blue, respectively); we also show whether structural instabilities (imaginary phonon frequencies) were observed in relaxations performed with (orange) or without (red) enforcing symmetries.

![image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00755v4-figure-97fd3433-eb1c-4d2b-9479-633f540cdad6.jpg)

ble within a factor of two of the DFT-PBE values for 69% of the materials in the phononDB-PBE database 44,45 .

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 (2) - 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 sign-symmetric distribution. Importantly, the mean for SRE[ κ ] and mean for SRME[ {K ( q ) s } ] (Table I) are comparable in the absence of cancellation of microscopic errors, while the mean SRE[ κ ] is significantly lower than mean SRME[ {K ( q ) s } ] when cancellation of microscopic errors occurs.

Importantly, we note that our benchmarks on

TABLE I. Benckmark metrics for fMLPs. The mean of SRE [ κ ] , and SRME [ {K ( q ) s } ] introduced here resolve the accuracy of fMLPs in describing interatomic forces and derived structural and heat-conduction properties. The Matbench F1 34 descriptor is based on atomic energies and resolves thermodynamic stability. Different (similar) values for SRE [ κ ] and SRME [ {K ( q ) s } ] indicate presence (absence) of cancellation of errors on microscopic vibrational properties. The lower SRE [ κ ] , the higher is the accuracy in predicting macroscopic κ (300K) . The lower SRME [ {K ( q ) s } ] , the higher is the accuracy in predicting both microscopic vibrational properties, and macroscopic κ (300K) . The higher F1 score, the higher is the accuracy in describing thermodynamic stability.

|                            |   ORB | SevenNet   | MACE    | CHGN.   | M3GN.   |
|----------------------------|-------|------------|---------|---------|---------|
| mean SRE[ κ ]              | 1.594 | 0 . 597    | 0 . 512 | 1 . 695 | 1 . 397 |
| mean SRME[ { K ( q ) s } ] | 1.808 | 0 . 767    | 0 . 664 | 1 . 717 | 1 . 469 |
| F1 (11/9/2024) 34          | 0.763 | 0.724      | 0.669   | 0.613   | 0.569   |

SRME[ {K ( q ) s } ] and SRE[ κ ] evaluate the accuracy of fMLPs in predicting the experimentally observable conductivity, which is a property derived from interatomic forces. Therefore, they provide information different and complementary to those given by the established tests based on atomic interaction energies reported on the Matbench Discovery leaderboard 34 . Table I highlights how fMLPs that are highly accurate according to tests based on interaction energies can be inaccurate in predicting crystal structures, interatomic forces, and thermal conductivity. Conversely, it worth noting that in the materials studied so far, fMLPs that are highly accurate according to the SRME[ {K ( q ) s } ] metric are also highly ranked by energy-based tests, and Fig. 4 c also shows that low SRME[ {K ( q ) s } ] is correlated with correctly describing crystal symmetries.

## Achieving first-principles accuracy via fine-tuning

In the previous section we discussed the ionic conductor LiBr as paradigmatic example where the zero-shot conductivity predictions from MACE-MP-0 are in strong disagreement-about a factor of two-with DFT-PBE reference data. These conductivity discrepancies might be due to lack of relevant reference data in the MP dataset 42 used to train MACE-MP-0, and in this section we introduce a fine-tuning protocol to correct them and achieve first-principles accuracy.

Employing the same DFT settings used to generate reference DFT data (see Methods for details), we prepared three distinct datasets for fine-tuning. Each dataset contains a set of different DFT-PBE frames - with 'frame' we mean a set comprising forces, total energy, and stress computed from DFT for a set of atomic positions perturbed from equilibrium. Perturbations consists of rattling (random atomic displacements drawn from Gaus-

FIG. 5. Achieving first-principles accuracy on LiBr's vibrational &amp; thermal properties through fine tuning. Panel a top, conductivity as a function of number of DFT-PBE supercell frames (see text) in the fine-tuning dataset, and number of fine-tuning epochs. The MACE zero-shot prediction is purple; performing fine-tuning on one single DFT frame and 100 finetuning epochs yields compatibility within 7% with the reference DFT-PBE value (green). Fine-tuning with datasets containing 3 frames with rattling and volume variations yields compatibility within 2% at 200 epochs. Panel a bottom, SRME [ {K ( q ) s } ] converges to a minimum upon increasing the number of epochs, and larger datasets yield faster convergence. After fine-tuning on 3 frames for 200 epochs, remarkable agreement between fine-tuned MACE and DFT-PBE is observed for: ( b ) phonons; ( c ) microscopic phonon energy-linewidth distributions; ( d ) macroscopic κ T ( ) . We also note that after fine-tuning, our WTE predictions for κ T ( ) agree with experiments 53 , and at 300 K with ab-initio (PBE) molecular-dynamics simulations (aiMD) 50 .

![image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00755v4-figure-d0261fbd-30d8-411a-a78a-a60239d9087b.jpg)

![image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00755v4-figure-c643eccd-4abd-4774-94ba-39b4d3774233.jpg)

![image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00755v4-figure-fafc3e61-a092-4de5-8085-36518a18e448.jpg)

![image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00755v4-figure-fe8ff3dc-d727-4b42-bda5-2d065889bf85.jpg)

sian distribution with standard deviation of 0.03Å) 54 , and in some cases also a few-percent isotropic rescaling of volume. DFT frames have the same size as those used to compute the reference phononDB-PBE phonons, 4x4x4 supercell of LiBr (512 atoms).

We generated three fine-tuning datasets, containing: (i) one single frame at equilibrium volume and with atomic rattling; (ii) 3 frames, out of which 2 feature ±1% isotropic volume rescaling and rattling, plus the frame in (i) at equilibrium volume and with rattling; (iii) 20 frames, out of which 10 are at equilibrium volume with different rattling, and 10 frames with rattling and volume variations of ±0.25 %, ±0.50%, ±1%, ±2%, and ±5%. Using these datasets, we performed fine-tuning training with number of epochs ranging from 1 to 1000. Fine-tuning MACE yields predictions in closer agreement with the reference DFT-PBE value. Specifically, from Fig. 5 a we see that fine-tuning relying on the 3-frame dataset (ii) and for 200 epochs yields κ (300K) compatible within 2% with reference DFT-PBE conductivity. Using the larger fine-tuning dataset (iii) does not show significant improvements, and the single-frame dataset (i) allows to reduce conductivity discrepancies from 47% to 7% in 100 fine-tuning epochs. The excellent agreement between the value of κ (300K) computed from reference DFT-PBE and fine-tuned MACE (3 frames, 200 epochs) originates from low SRME[ {K ( ) q s } ] ; this implies also remarkable compatibility for phonons in Fig. 5 b , microscopic energy-linewidth distribution in c , and macroscopic κ T ( ) in d . It is worth mentioning that conductivity compatibility within 2% corresponds to rootmean-square error (RMSE) on energy &lt;1.0 meV/atom, on forces &lt;5 meV/Å, and on stresses &lt;0.3 meV/Å ; see 3 Methods Fig. 10 for details oh how these errors depend on the fine-tuning dataset's size, and on the number of fine-tuning training epochs.

Finally, we highlight how the predictions for LiBr's κ T ( ) in the temperature range 100-800 K reproduce experimental measurements 53 , in both trend and magnitude (error bars show the compatibility threshold of 20% discussed in Ref. 53 ). We also emphasize that employing the WTE is essential for accurately describing heat transport in LiBr at high temperatures. In fact, upon increasing temperature the tunneling conductivity κ C becomes more important and non-negligible compared to the propagation conductivity κ P , and the total WTE conductivity κ T ( T )= κ P ( T )+ κ C ( T ) is in better agreement with the experimental trend compared to the BTE (propagation-only) conductivity κ P ( T ) . Finally, at 300 K our κ T is compatible also with predictions from ab-initio (PBE) molecular dynamics (aiMD) 50 .

## Conclusions

We introduced a framework that leverages foundation models for atomistic materials chemistry 21-33 and the Wigner formulation of heat transport 38,39 to address the major challenges of current techniques for designing heatmanagement materials: high computational cost, limited transferability, or lack of physics awareness. We demonstrated that our framework can systematically achieve quantum accuracy in conductivity predictions at computational costs up to 3 orders of magnitude lower than current methods, which rely on density-functional theory for direct evaluations or for full training of conventional (not fully composition-transferable) machine-learning potentials. Importantly, our physics-aware framework also overcomes the reliability issues (unknown unknowns) of current end-to-end machine-learning methods for predicting vibrational or thermal properties.

We developed benchmark tests based on technologically relevant observables, such as thermal conductivity, and showed applications to a database 44,45 of chemically

d

DFT

MA

MA

and structurally diverse materials using state-of-the-art foundation models (ORB-MPtrj 43 , SevenNet 24 , MACEMP-0 23 , CHGNet 22 , and M3GNet 21 ). The metrics we introduced test the accuracy of fMLPs in describing interatomic forces and derived vibrational or heat-conduction properties; therefore, they provide information that are complementary to SOTA energy-based descriptors 34 , and allow to test the accuracy of fMLP in describing crystal structures and atomistic dynamics.

Finally, we have shown that our framework enables the screening or prediction of materials that violate semiclassical Boltzmann transport, which are crucial for applications ranging from thermal insulation 11,39 to neuromorphic computing 12 . Ultimately, our framework paves the way for the theory-driven optimization, design, and discovery of materials for next-generation energy and information-management technologies.

- ∗ ms2855@cam.ac.uk
- [1] T. B. Blank, S. D. Brown, A. W. Calhoun, and D. J. Doren, The Journal of Chemical Physics 103 , 4129 (1995).
- [2] J. Behler and M. Parrinello, Physical Review Letters 98 , 146401 (2007).
- [3] A. P. Bartók, M. C. Payne, R. Kondor, and G. Csányi, Physical Review Letters 104 , 136403 (2010).
- [4] M. Rupp, A. Tkatchenko, K.-R. Müller, and O. A. von Lilienfeld, Physical Review Letters 108 , 058301 (2012).
- [5] L. Zhang, J. Han, H. Wang, R. Car, and W. E, Physical Review Letters 120 , 143001 (2018).
- [6] R. Drautz, Physical Review B 99 , 014104 (2019).
- [7] A. Seko, A. Togo, and I. Tanaka, Physical Review B 99 , 214108 (2019).
- [8] V. L. Deringer, A. P. Bartók, N. Bernstein, D. M. Wilkins, M. Ceriotti, and G. Csányi, Chemical Reviews 121 , 10073 (2021).
- [9] S. Batzner, A. Musaelian, L. Sun, M. Geiger, J. P. Mailoa, M. Kornbluth, N. Molinari, T. E. Smidt, and B. Kozinsky, Nature Communications 13 , 2453 (2022).
- [10] E. Kocer, T. W. Ko, and J. Behler, Annual Review of Physical Chemistry 73 , 163 (2022).
- [11] X. Qian, J. Zhou, and G. Chen, Nature Materials 20 , 1188 (2021).
- [12] G. F. Nataf, S. Volz, J. Ordonez-Miranda, J. Íñiguez González, R. Rurali, and B. Dkhil, Nature Reviews Materials , 1 (2024).
- [13] R. Okabe, A. Chotrattanapituk, A. Boonkird, N. Andrejevic, X. Fu, T. S. Jaakkola, Q. Song, T. Nguyen, N. Drucker, S. Mu, Y. Wang, B. Liao, Y. Cheng, and M. Li, Nature Computational Science 4 (2024).
- [14] Z. Guo, P. Roy Chowdhury, Z. Han, Y. Sun, D. Feng, G. Lin, and X. Ruan, npj Computational Materials 9 , 1 (2023).
- [15] J. Ojih, M. Al-Fahdi, Y. Yao, J. Hu, and M. Hu, Journal of Materials Chemistry A 12 , 8502 (2024).
- [16] A. Rodriguez, C. Lin, H. Yang, M. Al-Fahdi, C. Shen, K. Choudhary, Y. Zhao, J. Hu, B. Cao, H. Zhang, et al. , npj Computational Materials 9 , 20 (2023).
- [17] Y. Srivastava and A. Jain, Journal of Applied Physics 134 (2023).
- [18] V. D. Rose, A. Kozachinskiy, C. Rojas, M. Petrache, and P. Barceló, Three iterations of $(d-1)$-WL test distinguish non isometric clouds of $d$-dimensional points (2023).
- [19] J. Gilmer, S. S. Schoenholz, P. F. Riley, O. Vinyals, and G. E. Dahl, in Proceedings of the 34th International Conference on Machine Learning (PMLR, 2017) pp. 1263-1272.
- [20] I. Batatia, D. P. Kovács, G. N. C. Simm, C. Ortner, and G. Csányi, MACE: Higher Order Equivariant Message Passing Neural Networks for Fast and Accurate Force Fields (2023).
- [21] C. Chen and S. P. Ong, Nature Computational Science 2 (2022).
- [22] B. Deng, P. Zhong, K. Jun, J. Riebesell, K. Han, C. J. Bartel, and G. Ceder, Nature Machine Intelligence 5 , 1031 (2023).
- [23] I. Batatia, P. Benner, Y. Chiang, A. M. Elena, D. P. Kovács, J. Riebesell, X. R. Advincula, M. Asta, M. Avaylon, W. J. Baldwin, F. Berger, N. Bernstein, A. Bhowmik, S. M. Blau, V. Cărare, J. P. Darby, S. De, F. Della Pia, V. L. Deringer, R. Elijošius, Z. El-Machachi, F. Falcioni, E. Fako, A. C. Ferrari, A. Genreith-Schriever, J. George, R. E. A. Goodall, C. P. Grey, P. Grigorev, S. Han, W. Handley, H. H. Heenen, K. Hermansson, C. Holm, J. Jaafar, S. Hofmann, K. S. Jakob, H. Jung, V. Kapil, A. D. Kaplan, N. Karimitari, J. R. Kermode, N. Kroupa, J. Kullgren, M. C. Kuner, D. Kuryla, G. Liepuoniute, J. T. Margraf, I.-B. Magdău, A. Michaelides, J. H. Moore, A. A. Naik, S. P. Niblett, S. W. Norwood, N. O'Neill, C. Ortner, K. A. Persson, K. Reuter, A. S. Rosen, L. L. Schaaf, C. Schran, B. X. Shi, E. Sivonxay, T. K. Stenczel, V. Svahn, C. Sutton, T. D. Swinburne, J. Tilly, C. van der Oord, E. Varga-Umbrich, T. Vegge, M. Vondrák, Y. Wang, W. C. Witt, F. Zills, and G. Csányi, A foundation model for atomistic materials chemistry (2024), arXiv:2401.00096.
- [24] Y. Park, J. Kim, S. Hwang, and S. Han, Journal of Chemical Theory and Computation 20 , 4857 (2024).
- [25] H. Yang, C. Hu, Y. Zhou, X. Liu, Y. Shi, J. Li, G. Li, Z. Chen, S. Chen, C. Zeni, M. Horton, R. Pinsler, A. Fowler, D. Zügner, T. Xie, J. Smith, L. Sun, Q. Wang, L. Kong, C. Liu, H. Hao, and Z. Lu, Mattersim: A deep learning atomistic model across elements, temperatures and pressures (2024), arXiv:2405.04967 [cond-mat.mtrl-sci].
- [26] A. Merchant, S. Batzner, S. S. Schoenholz, M. Aykol, G. Cheon, and E. D. Cubuk, Nature 624 , 80 (2023).
- [27] K. Choudhary and B. DeCost, npj Computational Materials 7 , 10.1038/s41524-021-00650-1.
- [28] C. Chen, W. Ye, Y. Zuo, C. Zheng, and S. P. Ong, Chemistry of Materials 31 , 3564 (2019).
- [29] L. Ward, R. Liu, A. Krishna, V. I. Hegde, A. Agrawal, A. Choudhary, and C. Wolverton, Physical Review B 96 , 024104 (2017).
- [30] Y. Zuo, M. Qin, C. Chen, W. Ye, X. Li, J. Luo, and S. P. Ong, Materials Today 51 , 126 (2021).
- [31] T. Xie and J. C. Grossman, Physical Review Letters 120 , 145301 (2018).
- [32] J. Gibson, A. Hire, and R. G. Hennig, npj Computational Materials 8 , 211 (2022).

- [33] R. E. A. Goodall, A. S. Parackal, F. A. Faber, R. Armiento, and A. A. Lee, Science Advances 8 , eabn4117 (2022).
- [34] J. Riebesell, R. E. A. Goodall, P. Benner, Y. Chiang, B. Deng, A. A. Lee, A. Jain, and K. A. Persson, Matbench Discovery - A framework to evaluate machine learning crystal stability predictions (2024).
- [35] H. Yu, M. Giantomassi, G. Materzanini, J. Wang, and G.-M. Rignanese, Systematic assessment of various universal machine-learning interatomic potentials (2024).
- [36] B. Deng, Y. Choi, P. Zhong, J. Riebesell, Z. Li, K. Jun, K. A. Persson, and G. Ceder, Overcoming systematic softening in universal machine learning interatomic potentials by fine-tuning (2024).
- [37] H. Lee, V. I. Hegde, C. Wolverton, and Y. Xia, Accelerating High-Throughput Phonon Calculations via Machine Learning Universal Potentials (2024).
- [38] M. Simoncelli, N. Marzari, and F. Mauri, Nature Physics 15 , 809 (2019).
- [39] M. Simoncelli, N. Marzari, and F. Mauri, Physical Review X 12 , 041011 (2022).
- [40] R. Peierls, Annalen der Physik 395 , 1055 (1929).
- [41] P. B. Allen and J. L. Feldman, Physical Review Letters 62 , 645 (1989).
- [42] A. Jain, S. P. Ong, G. Hautier, W. Chen, W. D. Richards, S. Dacek, S. Cholia, D. Gunter, D. Skinner, G. Ceder, and K. A. Persson, APL Materials 1 , 011002 (2013).
- [43] O. M. Ltd., orb-models, https://github.com/ orbital-materials/orb-models (2024).
- [44] A. Togo, L. Chaput, and I. Tanaka, Phys. Rev. B 91 , 094306 (2015).
- [45] A. Seko, A. Togo, H. Hayashi, K. Tsuda, L. Chaput, and I. Tanaka, Physical Review Letters 115 , 10.1103/PhysRevLett.115.205901.
- [46] M. Simoncelli, N. Marzari, and A. Cepellotti, Phys. Rev. X 10 , 011019 (2020).
- [47] M. Simoncelli, F. Mauri, and N. Marzari, npj Computational Materials 9 , 1 (2023).
- [48] A. Pazhedath, L. Bastonero, N. Marzari, and M. Simoncelli, Phys. Rev. Appl. 22 , 024064 (2024).
- [49] J. J. Plata, P. Nath, D. Usanmaz, J. Carrete, C. Toher, M. de Jong, M. Asta, M. Fornari, M. B. Nardelli, and S. Curtarolo, npj Computational Materials 3 , 1 (2017).
- [50] F. Knoop, T. A. Purcell, M. Scheffler, and C. Carbogno, Physical Review Letters 130 , 236301 (2023).
- [51] S. Bandi, C. Jiang, and C. A. Marianetti, Machine Learning: Science and Technology (2024).
- [52] X. Gonze and C. Lee, Physical Review B 55 , 10355 (1997).
- [53] B. Hakansson and R. G. Ross, Journal of Physics: Condensed Matter 1 , 3977 (1989).
- [54] F. Eriksson, E. Fransson, and P. Erhart, Advanced Theory and Simulations 2 , 1800184 (2019).

## Methods

## Vibrational properties relevant to compute κ T ( )

In this section we discuss the relation between interatomic vibrational energies and the quantities that appear in the WTE thermal conductivity expression (1). The starting point is the Born-Oppenheimer Hamiltonian for atomic vibrations 39 expanded up to anharmonic third order 49,55-61 in the atomic displacements from equilibrium, and accounting for energy perturbations due to presence of isotopes 62-66

$$\underset { \text{and} } { \overset { \dots } { \text{lat} } } } \overset { \dots } { \text{line} } = \sum _ { R, b, \alpha } \frac { \hat { p } _ { R b \alpha } ^ { 2 } } { 2 M _ { b } } + \frac { 1 } { 2 } \sum _ { \underset { \alpha } { R ^ { \prime }, b ^ { \prime }, \alpha ^ { \prime } } } \frac { \partial ^ { 2 } V } { \partial u _ { R b \alpha } \partial u _ { R ^ { \prime } b ^ { \prime } \alpha ^ { \prime } } } \Big | \hat { u } _ { R b \alpha } \hat { u } _ { R ^ { \prime } b ^ { \prime } \alpha ^ { \prime } }$$

$$\text{and} \\ \text{uni} \\ \text{tice} \\ \text{atice}$$

$$\text{iac} \, \text{po} ^ { \quad } + \sum _ { R, b, \alpha } \left ( \frac { M _ { b } } { m _ { b } } - 1 \right ) \frac { \hat { p } _ { R b \alpha } ^ { 2 } } { 2 M _ { b } }.$$

Here, ˆ p R bα and ˆ u R bα are the momentum and positiondisplacement operators for the atom b having isotopeaveraged mass M b , position R τ + b ( R is the Bravaislattice vector and τ b the position in the crystal's unit cell), and along the Cartesian direction α ; the last term describes kinetic-energy perturbations induced by isotopes ( m b is the exact mass of atom at τ b , we adopt the common approximation 59,62 that the mass m b depends only on the index b and not on the lattice vector R ). The leading (harmonic) term in Eq. (5) determines the vibrational frequencies. In particular, the Fourier transform of the mass-rescaled Hessian of the interatomic potential yields the dynamical matrix at wavevector q ,

$$\text{and} _ { \text{hy-} } \quad D ( q ) _ { b \alpha, b ^ { \prime } \alpha ^ { \prime } } = \sum _ { R } \frac { \partial ^ { 2 } V } { \partial u _ { R b \alpha } \partial u _ { R \Phi ^ { \prime } \alpha ^ { \prime } } } \Big | _ { \text{eq} } \frac { e ^ { - i q \cdot ( R + \tau _ { b } - \tau _ { b ^ { \prime } } ) } } { \sqrt { M _ { b } M _ { b ^ { \prime } } } }, \ \ ( 6 )$$

and by diagonalizing the dynamical matrix,

$$\sum _ { b ^ { \prime } \alpha ^ { \prime } } D ( q ) _ { b \alpha, b ^ { \prime } \alpha ^ { \prime } } \mathcal { E } ( q ) _ { s, b ^ { \prime } \alpha ^ { \prime } } = \omega ^ { 2 } ( q ) _ { s } \mathcal { E } ( q ) _ { s, b \alpha }, \ \ ( 7 )$$

one obtains from the eigenvalues the phonon energies of the solid ℏ ω ( q ) s ( s is a band index ranging from 1 to 3 N at , where N at is the number of atoms in the unit cell), and from the eigenvectors E ( q ) s,bα the displacement patterns of atom b in direction α for the phonon having wavevector q and mode s .

As discusses in detail in Ref. 39 , the velocity operator appearing in the WTE conductivity expression (1) is determined by the wavevector derivative of the square root of the dynamical matrix:

$$\cdot \cdot \cdot \cdot \cdot \cdot \cdot v ^ { \beta } ( q ) _ { s, s ^ { \prime } } = \sum _ { b, \alpha, b ^ { \prime }, \alpha ^ { \prime } } \mathcal { E } ^ { * } ( q ) _ { s, b \alpha } [ \nabla _ { q } ^ { \beta } \sqrt { D ( q ) } _ { b \alpha, b ^ { \prime } \alpha ^ { \prime } } ] \mathcal { E } ( q ) _ { s ^ { \prime }, b ^ { \prime } \alpha ^ { \prime } }. ( 8 )$$

The anharmonic linewidths ℏ Γ ( a q ) s , which contribute to the total linewidths ℏ Γ( q ) s = ℏ Γ ( a q ) s + ℏ Γ ( i q ) s appearing in Eq. (1), are the energy broadenings due to anharmonic three-phonon interactions 55,59,67-69 . It can be shown 58 that these are determined by the third derivative

![image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00755v4-figure-48572b9f-725c-4eca-833d-e3daef77e3dc.jpg)

P

C

FIG. 6. Failures of semiclassical Boltzmann transport in binary crystals with low conductivity and having rocksalt (green), zincblende (orange), and wurtzite (blue) structure. These occur when the conductivity described by the BTE in terms of particle-like propagation of phonons ( κ P ) does not dominate over the coherences conductivity due to phonons' wave-like tunneling ( κ C ). The WTE 39 accounts for both κ P and κ C , and failures from semiclassical Boltzmann transport are highlighted by a relative strength of particle-wave transport ( κ P -κ C ) / κ ( P + κ C ) appreciably smaller than one (more intense red background, as in Fig. 1). These conductivities were obtained from the force constants in the phononDB-PBE dataset 44,45 .

of the interatomic potential (see Eq. (5)),

$$\hbar { \Gamma } _ { a } ( q ) _ { s } = & \frac { \pi } { N _ { c } ^ { 2 } } \underset { \substack { q ^ { \prime } q ^ { \prime \prime } \\ s ^ { \prime }, s ^ { \prime \prime } } } { \sum } \left [ 2 [ \bar { N } ( q ^ { \prime } ) _ { s ^ { \prime } } - \bar { N } ( q ^ { \prime \prime } ) _ { s ^ { \prime \prime } } ] \delta ( \omega ( q ) _ { s } + \omega ( q ^ { \prime } ) _ { s ^ { \prime } } - \omega ( q ^ { \prime \prime } ) _ { s ^ { \prime \prime } } ) \begin{array} { c } \text{particle-1} \\ \text{both the} \\ \text{in E. } ( 1 \\ \text{mechanis} \end{array} }.$$

$$+ \left [ 1 + \bar { N } ( q ^ { \prime } ) _ { s ^ { \prime } } + \bar { N } ( q ^ { \prime \prime } ) _ { s ^ { \prime \prime } } \right ] \delta \left [ \omega ( q ) _ { s } - \omega ( q ^ { \prime } ) _ { s ^ { \prime } } - \omega ( q ^ { \prime \prime } ) _ { s ^ { \prime \prime } } \right ] \right \} \times \quad \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \$$

$$\left | \sum _ { \substack { \alpha, \alpha ^ { \prime \prime }, \alpha ^ { \prime \prime } \\ b, b ^ { \prime }, b ^ { \prime \prime } R ^ { \prime \prime } } } } \frac { \alpha ^ { 3 } V } { \partial u _ { 0 b \alpha } \partial u _ { R ^ { \prime } b ^ { \prime } \alpha ^ { \prime } } \partial u _ { R ^ { \prime \prime } b ^ { \prime \prime } \alpha ^ { \prime \prime } } } \mathcal { E } ( q ) _ { s, b \alpha } \mathcal { E } ( q ^ { \prime } ) _ { s ^ { \prime }, b ^ { \prime } \alpha ^ { \prime } } \mathcal { E } ( q ^ { \prime \prime } ) _ { s ^ { \prime \prime }, b ^ { \prime \prime } \alpha ^ { \prime \prime } } \quad \text{anham} \\ b \text{where} \\ \text{$\text{order}$} \\ \text{$\text{frame}$}$$

$$\sqrt { \frac { \hbar { ^ { 3 } } } { 8 } } \frac { \Delta ( q + q ^ { \prime } + q ^ { \prime \prime } ) e ^ { - i [ q \cdot \tau _ { b } + q ^ { \prime } \cdot ( R ^ { \prime } + \tau _ { b ^ { \prime } } ) + q ^ { \prime \prime } \cdot ( R ^ { \prime \prime } + \tau _ { b ^ { \prime } } ) ] } } { \sqrt { M _ { b } M _ { b ^ { \prime } } M _ { b ^ { \prime \prime } } \omega ( q ) _ { s } \omega ( q ^ { \prime } ) _ { s ^ { \prime \prime } } } } \Big | ^ { 2 }, \quad ( 9 ) \quad \text{mater}$$

where ¯ ( N q ) s =[exp( ℏ ω ( q ) s /k B T ) -1] -1 is the BoseEinstein distribution, ∆( + q q ′ + q ′ ′ ) is the Kronecker delta (equal to 1 if q + + q ′ q ′ ′ is a reciprocal lattice vector, 0 otherwise), δ is the Dirac delta. The linewidth due to isotopic-mass disorder is determined exclusively by harmonic properties and concentration of isotopes 62 :

$$\hbar { \Gamma } _ { i } ( q ) _ { s } { = } & \frac { \hbar { \pi } } { 2 N _ { c } } [ \omega ( q ) _ { s } ] ^ { 2 } \sum _ { q ^ { \prime }, s ^ { \prime } } \delta \left [ \omega ( q ) _ { s } { - } \omega ( q ^ { \prime } ) _ { s ^ { \prime } } \right ] \quad \text{(\sum_{nof\colon})} } ^ { \text{ as un} } \\ & \times \sum _ { b } g _ { 2 } ^ { b } \left | \sum _ { \alpha } \mathcal { E } ( q ) _ { s, b \alpha } ^ { * } \mathcal { E } ( q ^ { \prime } ) _ { s ^ { \prime }, b \alpha } \right | ^ { 2 }, \quad \text{(\sum_{nof\colon})} } ^ { \text{ as un} }$$

semiclassical particle-like BTE fails. This happens when particle-like propagation mechanisms - described by both the BTE and WTE 38,46,70 , and determining κ P in Eq. (1) - do not dominate over wave-like tunneling mechanisms - missing from the BTE but described by the WTE and determining κ C in Eq. (1). Failures of the BTE have been discussed to appear in, e.g., strongly anharmonic complex crystals for energy harvesting 71-83 , thermal barrier coatings 48,84-86 , as well as in several disordered functional materials 87-92 . Our automated framework computes the total Wigner thermal conductivity κ P + κ C 39 , and therefore can be used to find materials that violate semiclassical Boltzmann thermal transport. In Fig. 6 we show that among the 103 compounds analyzed at 300 K, the BTE fails in materials with low thermal conductivity ( κ TOT ≲ 2 W/mK), such as AgBr and AgCl, since the wave-like tunneling conductivity becomes comparable to the particle-like propagation conductivity. We also note that in simple crystals with large conductivity, the BTE is accurate as the particle-like propagation conductivity dominates over the wave-like tunneling conductivity.

where g b 2 = ∑ i f i,b ( M b -m i,b M b ) 2 describes the variance of 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 = ∑ i f i,b m i,b is the weighted average mass).

## Failures of semiclassical Boltzmann transport

As anticipated in the main text, our automated framework can be used to identify materials in which the

Symmetric Relative Mean Error on harmonic and anharmonic vibrational properties. The analyses reported in the main text discuss SRME [ {K ( q ) s } ] as a descriptor that is informative of the accuracy of fMLPs in describing microscopic harmonic and anharmonic properties. In this section we provide the tools to resolve whether a large SRME [ {K ( q ) s } ] originates from large errors in microscopic harmonic properties, or large errors in anharmonic properties, or from both. To this aim, we employ the expression for the microscopic, singlephonon contributions to the thermal conductivity (see

terms inside the square brackets in Eq. (1)) to define descriptors that quantity the accuracy with which harmonic and anharmonic properties are predicted. In particular, we note that the conductivity contribution of a phonon ( q ) s is a function of: (i) all the phonon frequencies { ω ( q s ′ } at a fixed wavevector q and variable mode s ′ = 1 , . . . , N s ( N s is the number of phonon bands, equal to 3 times the number of atoms in the crystal's unit cell 93 ); (ii) all velocity-operator elements { v ( q ) s,s ′ } at fixed q , s and variable s ′ ; (iii) all isotopic linewidths at fixed q and variable s ′ , { Γ( i q ) s ′ } , which depend solely on harmonic properties, see Eq. 10; (iv) the anharmonic linewidths { Γ ( a q ) s ′ } at fixed q and variable s ′ . This can be summarized using the following notation: K ( q ) s = K ( q ) s [ { ω ( q ) s ′ } { , v ( q ) s,s ′ } { , Γ ( i q ) s ′ } { , Γ ( a q ) s ′ } ] where the curly brackets denote the set of values of a certain quantity over all the bands at fixed q (e.g., { ω ( q ) s ′ } = { ω ( q ) 1 , ω ( q ) 2 , . . . , ω ( q ) N s } ).

To quantify the impact of errors on the harmonic (har) and anharmonic (anh) properties on the single-phonon conductivity contributions, we define the Symmetric Relative Mean Error (SRME) on these properties as follows:

$$\text{SRME} [ \text{har} ] = & 2 \frac { \sum _ { q s } | \mathcal { K } ( q ) _ { s } \left [ \{ \omega _ { \text{MMLP} } ( q ) _ { s ^ { \prime } } \}, \{ \nu _ { \text{MLP} } ( q ) _ { s, s ^ { \prime } } \}, \{ \Gamma _ { i, \text{MLP} } ( q ) _ { s ^ { \prime } } \}, \{ \Gamma _ { a, \text{DFT} } ( q ) _ { s ^ { \prime } } \} \right ] - \mathcal { K } _ { \text{DFT} } ( q ) _ { s } \right ] } { \sum _ { q s } \left ( \mathcal { K } ( q ) _ { s } \left [ \{ \omega _ { \text{MLP} } ( q ) _ { s ^ { \prime } } \}, \{ \nu _ { \text{MLP} } ( q ) _ { s, s ^ { \prime } } \}, \{ \Gamma _ { i, \text{MLP} } ( q ) _ { s ^ { \prime } } \}, \{ \Gamma _ { a, \text{DFT} } ( q ) _ { s ^ { \prime } } \right ] + \mathcal { K } _ { \text{DFT} } ( q ) _ { s } \right ) }, \quad ( 1 1 )$$

$$\text{SRME} [ \tanh]=2 \frac { \sum _ { q s } \left | \mathcal { K } ( q ) _ { s } \left [ \{ \omega _ { D F T } ( q ) _ { s ^ { \prime } } \}, \{ \nu _ { D F T } ( q ) _ { s, s ^ { \prime } } \}, \{ \Gamma _ { i, D F T } ( q ) _ { s ^ { \prime } } \}, \{ \Gamma _ { a, \text{MLP} } ( q ) _ { s ^ { \prime } } \} \right ] - \mathcal { K } _ { D F T } ( q ) _ { s } \right | } { \sum _ { q s } \left ( \mathcal { K } ( q ) _ { s } \left [ \{ \omega _ { D F T } ( q ) _ { s ^ { \prime } } \}, \{ \nu _ { D F T } ( q ) _ { s, s ^ { \prime } } \}, \{ \Gamma _ { i, D F T } ( q ) _ { s ^ { \prime } } \}, \{ \Gamma _ { a, \text{MLP} } ( q ) _ { s ^ { \prime } } \} \right ] + \mathcal { K } _ { D F T } ( q ) _ { s } \right ) }, \quad ( 1 2 )$$

where we have used the shorthand notation K DFT ( q ) s = K DFT ( q ) s [ { ω DFT ( q ) s ′ } { , v DFT ( q ) s,s ′ } { , Γ i , DFT ( q ) s ′ } { , Γ a DFT , ( q ) s ′ } ] . Intuitively, SRME[har] (11) is large when harmonic vibrational properties ( ω ( q ) s , or v ( q ) s,s ′ , or Γ( i q ) s ) differ

FIG. 7. Relation between harmonic and anharmonic errors, and thermal conductivity, calculated for 103 compounds using Eq. (11) and Eq. (12) considering DFT and MACE-MP-0 23 . The area of the marker is proportional to SRE[ κ ], Eq. (2). We highlight how SRME[har] (Eq. (11)), SRME[anh] (Eq. (12)), and SRE[ κ ] are all low (large) in BeO (LiBr); in contrast, BeTe displays large SRME[har], large SRME[anh], but low SRE[ κ ] due to compensation of errors.

![image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00755v4-figure-ac28a842-0af4-4668-81a8-69bd4977e3f5.jpg)

between DFT and fMLP, while SRME[anh] (12) is large when the anharmonic linewidths ( Γ ( a q ) s ) differ between DFT and fMLP. We show in Fig. 7 that small SRME[har] and small SRME[anh] (e.g., as in BeO), imply a small SRE[ κ ] Fig. 7 also shows that often large SRME[har] and large SRME[anh] (e.g. LiBr) translate into large SRE [ κ ] ; however there also cases (e.g. BeTe) in which large SRME[har] and large SRME[anh] can also compensate each other and result in a small SRE [ κ ] . This confirms that SRE [ κ ] is not a reliable descriptor for the capability of fMLPs to capture the harmonic and anharmonic physics underlying heat conduction. In contrast, having both small SRME[har] and small SRME[anh] is a sufficient condition to accurately capture harmonic and anharmonic vibrational properties, as well as conductivity.

Overall, these tests illustrate the importance of benchmarking the accuracy of fMLPs in predicting both microscopic harmonic and anharmonic vibrational properties, further motivating the introduction of the SRME[ {K ( q ) s } ].

Influence of isotopic scattering on conductivity. To analyze the influence of isotope scattering on the conductivity, we compare thermal conductivity values with and without considering the linewidths from isotope mass-disorder (Eq. (10)). The results are summarized in Table II. In wurtzite BeO and rocksalt LiBr, the impact of isotope scattering is minimal or small, respectively. However, in zincblende BeTe isotopic and anharmonic linewidths have comparable values (Fig. 2 b ), and therefore both affect the conductivity. In this material, MACE-MP-0 significantly overestimates some linewidths

compared to DFT-PBE, and this compensates for the overestimation in other regions. Consequently, MACEMP-0 shows a significantly larger error in thermal conductivity when isotope scattering is not considered.

TABLE II. Influence of isotopic scattering on the thermal conductivity at 300 K for BeO, BeTe, and LiBr.

|                     | wurtzite BeO   | wurtzite BeO   | zincblende BeTe   | zincblende BeTe   | rocksalt LiBr   | rocksalt LiBr   |
|---------------------|----------------|----------------|-------------------|-------------------|-----------------|-----------------|
|                     | DFT            | MACE           | DFT               | MACE              | DFT             | MACE            |
| with Γ i ( q ) s    | 286.391        | 259.443        | 78.246            | 69.138            | 1.789           | 0.948           |
| without Γ i ( q ) s | 291.963        | 263.956        | 289.374           | 402.030           | 1.803           | 0.951           |

Automated Wigner Conductivity Workflow. The workflow, outlined in Fig. 8, includes structure relaxation, force constant computation, thermal conductivity evaluation and analysis. To automate the calculation of interatomic forces, we developed an interface between the packages ase 94 (version 3.23.0) and phono3py 44,95 (version 3.2.0). Atomic positions and unit-cell lattice parameters were simultaneously relaxed to minimize energy, forces, and stresses. We note that the determination of irreducible q points needed for thermal conductivity calculations depends on the crystal symmetry, and to evaluate SRME [ {K ( q ) s } ] (Eq. (3)) it is essential to compare structures that have the same crystal symmetry when studied using DFT or fMLPs.

In the second stage, instead, the constraint on the crystal symmetry is removed, to test whether the fMLP correctly preserves the physical crystal symmetry or not. After having removed the crystal-symmetry constraints, the relaxation proceeds similarly to step one: first, atomic positions and cell parameters are relaxed together, and then an additional fixed-cell relaxation is performed. If the symmetry is preserved, the final relaxed structure is used. If symmetry is broken, the structure from the first stage, with enforced symmetry, is retained. We detect symmetries utilizing the spglib package 96 with the default precision parameters. As discussed in Fig. 4 c , in some cases retaining the symmetric structure yielded structural instabilities (imaginary phonon frequencies),

We employed a two-stage relaxation protocol. During the first stage, symmetries are explicitly enforced as constraints while simultaneously relaxing atomic positions and cell parameters. After each step of atomic position relaxation, cell parameters are adjusted if the forces exceed a certain threshold (values are fMLP-specific and discussed later). Finally, an additional relaxation of atomic positions with fixed cell parameters is performed using a stricter threshold, to ensure that the structure reaches a stable energy minimum at a fixed volume, necessary to perform physically sound phonon calculations. The force thresholds employed for each fMLP are discussed below.

FIG. 8. Automated Wigner Conductivity Workflow. The first part of the workflow, involving force evaluation, is performed on a GPU. The input structures undergo relaxation in two stages. In the first stage, initial symmetries are enforced while simultaneously relaxing atomic positions and cell parameters. In the second stage, the structure is further relaxed by removing the symmetry constraint. If the symmetry group of the final structure changes, the structure with enforced symmetry relaxation is used instead. After initializing the displacement structures, the force sets are calculated for each displacement. CPUs handle the remainder of the workflow. Force constants are constructed from the displacement force sets, and thermal conductivity is calculated using phono3py . Structures with imaginary phonon frequencies are discarded from the final analysis (conducted on a PC).

![image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00755v4-figure-ea9cc52e-a3d7-47c9-baa6-d1023d906566.jpg)

which disappeared when the structure was allowed to break the symmetry and relax. While these structurally stable and broken-symmetry structures are not useful for evaluating SRME [ {K ( q ) s } ] , they provide information about shortcomings of fMLPs that could addressed in future studies or used as benchmark metric for accuracy in the description of crystal structures.

After noting that more precise fMLPs allowed to use stricter relaxation thresholds, we tailored these thresh-

olds based on the performance of each fMLP. In particular, M3GNet and ORB-MPtraj reached convergence using force thresholds of 2 0 . × 10 -3 eV/Å and 8 0 . × 10 -4 eV/Åfor the joint cell-position and the additional fixed-cell relaxation respectively. For these two models, using lower thresholds resulted in relaxation oscillating around different configurations without reaching convergence. For CHGNet, MACE-MP-0, and SevenNet, it was possible to reach convergence using force thresholds of 10 -4 eV/Å and 5 0 . × 10 -5 eV/Å for joint cell-position and additional fixed-cell relaxation, respectively. It is worth noting that ORB-MPtraj yielded non-converging relaxation (oscillating between different structures, using up to 500 iterations in stage 1 and 200 iterations in stage 2) in 40 structures, while all other models successfully converged for all materials. The BFGS algorithm was employed for structural relaxation, whereas for CHGNet and M3GNet, an initial relaxation was performed using their built-in protocols, which was then followed by the 2stage relaxation scheme described above. In the supercell force-constant calculations, we used the same parameters used in the DFT reference data 44,45 .

Thermal conductivities were computed using phono3py following Refs. 44,95 . A 19 × 19 × 19 q -mesh was used for rocksalt and zincblende structures, and a 19 × 19 × 15 q -mesh for wurtzite structures. The collision operator was computed using the tetrahedron method. For Fig. 2, linewidths were calculated on a 9 × 9 × 9 q -mesh, for graphical clarity. We utilized version 3 of phono3py for anharmonic linewidth calculations and employed the Wang method 97,98 for the NAC term, as implemented in the phonopy 95,99 package (version 2.26.6). When structures with imaginary phonon frequencies were found with mp-fMLPs, these were discarded from subsequent thermal-conductivity analysis.

DFT and zero-shot mp-fMLP calculations. Computational details of the DFT calculations underlying the phononDB-PBE database can be found in Refs. 44,45 . These employed the Perdew, Burke, and Ernzerhof (PBE) exchange-correlation functional 100 , consistent with the approach used for generating the MP dataset 42 , used to train mp-fMLPs.

Fig. 9 shows that in wurtzite BeO, the zero-shot MACE-MP-0-L2 calculation has a computational cost that is more than 3 orders of magnitude lower than the conventional DFT calculation, and yields a conductivity compatible within 7% with DFT reference values (Fig. 2 a ). Importantly, we note that speedups much stronger than those highlighted in Fig. 9 are expected in disordered materials without symmetries 47,87 . In fact, BeO has a simple unit cell containing just 4 atoms and with 12 crystal symmetries 102 , which allow to reduce the number of displacements 95 . In contrast, in disordered materials symmetries are broken and consequently a larger number of finite-displacement calculations are needed. Moreover, standard DFT approaches scale with the cube of the number of atoms, while fMLPs have a linear scaling. Based on these considerations, we expect even stronger speedups in disordered materials 103-108 .

We note that the DFT parameters used to generate the phononDB-PBE database are similar but not always exactly equal to those used for DFT calculations in the MP database. In the regime of computational convergence, variations of the DFT parameters at fixed functional are expected to yield unimportant effects on the conductivity. In particular, Ref. 109 shows that for the prototypical case of silicon, computing force constants using DFT parameters within the range of 'computationally converged values' found in the literature produces conductivity variations that are within 7 % , and also shows that in ideal cases varying DFT parameters within the domain of computational convergence should yield conductivity variations smaller than 2 % . As such, the significant discrepancies shown in Fig. 3 are expected to be not caused by possible small differences in the DFT parameters.

The computational details on the mp-fMLPs can be found in the following references: M3GNet 21 , CHGNet 22 , MACE-MP-0 23 , SevenNet 24 , and ORBMPtraj 43 . M3GNet was used through the matgl 110 (version 1.1.2) package with the M3GNet-MP-2021.2.8PES model. The CHGNet version 0.3.8 was used. The MACE-MP-0 2024-01-07-mace-128-L2 model was run through LAMMPS 111 . SevenNet was used through the sevenn package (version 0.9.2) with the SevenNet-0\_11July2024 model. ORB-MPtraj was used through

FIG. 9. Computational cost: DFT vs fMLP MACEMP-0-L2 for wurtzite BeO. Orange and red show the cost for computing harmonic and anharmonic force constants, respectively; green is the cost to obtain the conductivity from the force constants through the solution of the Wigner transport equation. The force constants were computed on CPU for DFT, and on GPU for MACE (see text for details); these costs were compared using the conversion 1 GPUh = 55 CPUh 101 .

![image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00755v4-figure-a07ec2de-d9e1-4072-a798-3bb443d5e7a6.jpg)

the orb-models package (version 0.3.0, commit d142854) with the orb-mptraj-only-v1 model.The details on the violin plots presented in Fig. 4 can be found in Ref. 112 .

The ORB-MPtraj, SevenNet, MACE, and CHGNet calculations were performed on an Nvidia Tesla A100 SXM4 80GB GPU, with the relaxation and force calculations for all 103 structures costing approximately ∼ 3 GPU hours per model. M3GNet was executed on AMD EPYC 7702 CPUs, with relaxation and force calculations costing around ∼ 30 CPU hours. Thermal conductivity calculations for all 103 structures required about ∼ 50 CPU hours per model, and similar resources were needed for evaluating thermal conductivity from the DFT force constants. Overall, the analysis of these five models required approximately ∼ 12 GPU hours and ∼ 330 CPU hours.

Details on fine-tuning. In this section we discuss how the size of the fine-tuning dataset and the number of finetuning epochs influence the accuracy with which energies, forces, and stresses are described. Results for LiBr are shown in Fig. 10; we underscore that to achieve 2% accuracy on the conductivity it is necessary to have RMSE on energy smaller than 1.0 meV/atom, RMSE on forces smaller than 5 meV/Å (about 2% relative error), and RMSE on stresses smaller than 0.3 meV/Å . 3

The three DFT datasets used for fine-tuning are made up of DFT calculations in 4x4x4 supercells of LiBr (same size used for the phonon calculation in the phononDBPBE 44,45 database) perturbed from equilibrium with atomic rattling and in some cases also volume variations from 0.25% to 5%. All DFT calculations are performed employing the Perdew, Burke, and Ernzerhof (PBE) exchange-correlation functional 100 , and the projector augmented wave (PAW) formalism for electronic minimization, implemented within the VASP (Vienna ab initio simulation package) code 113,114 . Following Ref. 45 , we use a 600 eV plane-wave-basis cutoff (20% higher than the default one), and the Brillouin zone is sampled at the Gamma point. Using these DFT settings, the LiBr structure in the phononDB-PBE database shows a residual negative stress of 960 Bar. Therefore, to make the most accurate possible comparison between fine-tuned MACE-MP-0 and reference DFT data, we relaxed the fine-tuned MACE-MP-0 structure exactly to the same residual stress present in the reference DFT-PBE data.

The MACE-MP-0 2024-01-07-mace-128-L2 model https://github.com/ACEsuit/mace-mp/ commits/mace\_mp\_0 was fine-tuned using the mace\_run\_train script in the main branch of https://github.com/ACEsuit/mace . The training started from the last checkpoint of MACE-MP-0 L2. The three datasets discussed above were used to perform

FIG. 10. Improving fMLP accuracy in LiBr via finetuning. Upon increasing the number of fine-tuning epochs, the root mean square errors on energies (panel a ), forces ( b and c ), and stresses ( d ) decrease. First-principles accuracy is reached faster using larger fine-tuning datasets. The green region corresponds to compatibility within 2% with reference DFT values for κ .

![image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00755v4-figure-4b07e86e-67a5-4be1-80cd-b88f7583f37e.jpg)

three independent trainings. For trainings on the 1-frame and 3-frame datasets we employed a batch size of 1, and for trainings on the 20-frame dataset we used a batch of size 4. Since our fine-tuning datasets were computed with the DFT-PBE functional, which was also employed to generate the MP database 42 used to train MACE-MP-0, in the fine-tuning we kept the E0s fixed at their original MACE-MP-0 value. Representative examples of training scripts are available at: https: //github.com/MSimoncelli/fMLP\_conductivity.git . The fine-tuning was performed on Nvidia Tesla A100

Data availability. The Matbench Discovery results are available at https://matbench-discovery. materialsproject.org/ . The phononDB-PBE dataset including the displacements generated by phono3py and the corresponding force sets are available https://github.com/atztogo/phonondb 44,45 . The remaining dataset will be made available in the future.

Code availability. The phono3py and phonopy packages are available at https://github.com/phonopy/ ; the ase package is available at https://gitlab. com/ase/ase ; the spglib package is available at https://github.com/spglib/spglib . matgl containing the M3GNet model is available at https://github. com/materialsvirtuallab/matgl ; the CHGNet model and package is available at https://github.com/ CederGroupHub/chgnet ; the SevenNet model and package is available at https://github.com/MDIL-SNU/ SevenNet ; the ORB-MPtraj model and package is available at https://github.com/orbital-materials/ orb-models . The MACE used through LAMMPS 111 is available at https://github.com/ACEsuit/lammps , and the MACE-MP-0 model is available at https://github. com/ACEsuit/mace-mp . The Hiphive package is available at https://gitlab.com/materials-modeling/ hiphive .

- ∗ ms2855@cam.ac.uk
- [55] J. Carrete, B. Vermeersch, A. Katre, A. van Roekeghem, T. Wang, G. K. Madsen, and N. Mingo, Computer Physics Communications 220 , 351 (2017).
- [56] T. Tadano, Y. Gohda, and S. Tsuneyuki, J. Phys. Condens. Matter 26 , 225402 (2014).
- [57] G. Barbalinardo, Z. Chen, N. W. Lundgren, and D. Donadio, Journal of Applied Physics 128 , 135104 (2020).
- [58] L. Paulatto, F. Mauri, and M. Lazzeri, Physical Review B 87 , 214303 (2013).
- [59] G. Fugallo, M. Lazzeri, L. Paulatto, and F. Mauri, Phys. Rev. B 88 , 045430 (2013).
- [60] A. Cepellotti, J. Coulter, A. Johansson, N. S. Fedorova, and B. Kozinsky, Journal of Physics: Materials 5 , 035003 (2022).
- [61] I. Savić, D. Donadio, F. Gygi, and G. Galli, Applied Physics Letters 102 (2013).
- [62] S.-i. Tamura, Phys. Rev. B 27 , 858 (1983).
- [63] C. J. Foss and Z. Aksamija, Physical Review Materials 4 , 124006 (2020).
- [64] K. R. Hahn, C. Melis, F. Bernardini, and L. Colombo, Frontiers in Mechanical Engineering 7 , 712989 (2021).
- [65] O. Abou El Kheir, L. Bonati, M. Parrinello, and M. Bernasconi, npj Computational Materials 10 , 33 (2024).
- [66] D. A. Stewart, I. Savić, and N. Mingo, Nano Letters 9 ,
- 81 (2009).
- [67] A. Togo and A. Seko, The Journal of Chemical Physics 160 , 211001 (2024).
- [68] T. Tadano, Y. Gohda, and S. Tsuneyuki, Journal of Physics: Condensed Matter 26 , 225402 (2014).
- [69] W. Li, J. Carrete, N. A. Katcho, and N. Mingo, Computer Physics Communications 185 , 1747 (2014).
- [70] J. Dragašević and M. Simoncelli, arXiv:2303.12777 (2023).
- [71] E. Di Lucente, M. Simoncelli, and N. Marzari, Physical Review Research 5 , 033125 (2023).
- [72] A. Jain, Phys. Rev. B 102 , 201201 (2020).
- [73] Y. Xia, V. I. Hegde, K. Pal, X. Hua, D. Gaines, S. Patel, J. He, M. Aykol, and C. Wolverton, Phys. Rev. X 10 , 041029 (2020).
- [74] X. Shen, N. Ouyang, Y. Huang, Y.-H. Tung, C.-C. Yang, M. Faizan, N. Perez, R. He, A. Sotnikov, K. Willa, C. Wang, Y. Chen, and E. Guilmeau, Advanced Science 11 , 10.1002/advs.202400258 (2024).
- [75] T. Pandey, M.-H. Du, D. S. Parker, and L. Lindsay, Materials Today Physics 28 , 10.1016/j.mtphys.2022.100881.
- [76] A. Fiorentino and S. Baroni, Physical Review B 107 , 054311 (2023).
- [77] S. Zhou, E. Xiao, H. Ma, K. Gofryk, C. Jiang, M. E. Manley, D. H. Hurley, and C. A. Marianetti, Physical Review Letters 132 , 106502 (2024).
- [78] J. Zheng, C. Lin, C. Lin, G. Hautier, R. Guo, and B. Huang, npj Computational Materials 10 (2024).
- [79] F. Jia, S. Zhao, J. Wu, L. Chen, T.-H. Liu, and L.-M. Wu, Angewandte Chemie 135 , e202315642 (2023).
- [80] T. Tadano and W. A. Saidi, Physical Review Letters 129 , 185901 (2022).
- [81] T. Bernges, R. Hanus, B. Wankmiller, K. Imasato, S. Lin, M. Ghidiu, M. Gerlitz, M. Peterlechner, S. Graham, G. Hautier, Y. Pei, M. R. Hansen, G. Wilde, G. J. Snyder, J. George, M. T. Agne, and W. G. Zeier, Advanced Energy Materials 12 , 2200717 (2022).
- [82] X. Yang, J. Tiwari, and T. Feng, Materials Today Physics 24 , 100689 (2022).
- [83] Z. Zeng, X. Shen, R. Cheng, O. Perez, N. Ouyang, Z. Fan, P. Lemoine, B. Raveau, E. Guilmeau, and Y. Chen, Nature Communications 15 (2024).
- [84] G. Caldarelli, M. Simoncelli, N. Marzari, F. Mauri, and L. Benfatto, Physical Review B 106 , 024312 (2022).
- [85] Y. Xia, D. Gaines, J. He, K. Pal, Z. Li, M. G. Kanatzidis, V. Ozolinš, and C. Wolverton, Proceedings , of the National Academy of Sciences 120 , e2302541120 (2023).
- [86] Y. Luo, X. Yang, T. Feng, J. Wang, and X. Ruan, Nature Communications 11 , 2554 (2020).
- [87] A. F. Harper, K. Iwanowski, W. C. Witt, M. C. Payne, and M. Simoncelli, Physical Review Materials 8 , 043601 (2024).
- [88] Y. Liu, H. Liang, L. Yang, G. Yang, H. Yang, S. Song, Z. Mei, G. Csányi, and B. Cao, Advanced Materials 35 , 2210873 (2023).
- [89] M. Ndour, P. Jund, and L. Chaput, Journal of NonCrystalline Solids 621 , 122618 (2023).
- [90] N. W. Lundgren, G. Barbalinardo, and D. Donadio, Physical Review B 103 , 024204 (2021).
- [91] S. Kielar, C. Li, H. Huang, R. Hu, C. Slebodnick, A. Alatas, and Z. Tian, Nature Communications 15 , 6981 (2024).

- [92] M. Simoncelli, D. Fournier, M. Marangolo, E. Balan, K. Béneut, B. Baptiste, B. Doisneau, N. Marzari, and F. Mauri, Temperature-invariant heat conductivity from compensating crystalline and glassy transport: from the Steinbach meteorite to furnace bricks (2024).
- [93] J. M. Ziman, Electrons and phonons: the theory of transport phenomena in solids (Oxford university press, 1960).
- [94] A. H. Larsen, J. J. Mortensen, J. Blomqvist, I. E. Castelli, R. Christensen, M. Dułak, J. Friis, M. N. Groves, B. Hammer, C. Hargus, E. D. Hermes, P. C. Jennings, P. B. Jensen, J. Kermode, J. R. Kitchin, E. L. Kolsbjerg, J. Kubal, K. Kaasbjerg, S. Lysgaard, J. B. Maronsson, T. Maxson, T. Olsen, L. Pastewka, A. Peterson, C. Rostgaard, J. Schiøtz, O. Schütt, M. Strange, K. S. Thygesen, T. Vegge, L. Vilhelmsen, M. Walter, Z. Zeng, and K. W. Jacobsen, Journal of Physics: Condensed Matter 29 , 273002 (2017).
- [95] A. Togo, L. Chaput, T. Tadano, and I. Tanaka, Journal of Physics: Condensed Matter 35 , 353001 (2023).
- [96] A. Togo and I. Tanaka, $ { Spglib}$: a software library for crystal symmetry search (2018).
- [97] Y. Wang, J. J. Wang, W. Y. Wang, Z. G. Mei, S. L. Shang, L. Q. Chen, and Z. K. Liu, Journal of Physics: Condensed Matter 22 , 202201 (2010).
- [98] X. Gonze and C. Lee, Physical Review B 55 , 10355 (1997).
- [99] A. Togo, Journal of the Physical Society of Japan 92 , 012001 (2023).
- [100] J. P. Perdew, K. Burke, and M. Ernzerhof, Physical Review Letters 77 , 3865 (1996).
- [101] The Cambridge Service for Data Driven Discovery (CSD3), cost of CPU and GPU hours (2024-09-02).
- [102] A. Togo, K. Shinohara, and I. Tanaka, Science and Technology of Advanced Materials: Methods 0 , 2384822.
- [103] L. Reicht, L. Legenstein, S. Wieser, and E. Zojer, Molecules 29 , 3724 (2024).
- [104] I. Kriuchevskyi, V. V. Palyulin, R. Milkus, R. M. Elder, T. W. Sirk, and A. Zaccone, Phys. Rev. B 102 , 024108 (2020).
- [105] H. de Araujo Oliveira, Z. Fan, A. Harju, and L. F. C. Pereira, ACS Applied Nano Materials (2024).
- [106] E. Bosoni, D. Campi, D. Donadio, G. C. Sosso, J. Behler, and M. Bernasconi, Journal of Physics D: Applied Physics 53 , 054001 (2020).
- [107] F. DeAngelis, M. G. Muraleedharan, J. Moon, H. R. Seyf, A. J. Minnich, A. J. H. McGaughey, and A. Henry, Nanoscale and Microscale Thermophysical Engineering 23 , 81 (2019).
- [108] M. Guerboub, S. D. Wansi Wendji, C. Massobrio, A. Bouzid, M. Boero, G. Ori, and E. Martin, The Journal of Chemical Physics 158 (2023).
- [109] A. Jain and A. J. McGaughey, Computational Materials Science 110 , 115 (2015).
- [110] T. W. Ko, M. Nassar, S. Miret, E. Liu, and S. P. Ong, Zenodo (2023), https://doi.org/10.5281/zenodo. 8025189 .
- [111] S. Plimpton, Journal of Computational Physics 117 , 1 (1995).
- [112] J. L. Hintze and R. D. Nelson, The American Statistician 52 , 181 (1998).
- [113] G. Kresse and J. Hafner, Physical Review B 47 , 558 (1993).
- [114] G. Kresse and J. Furthmüller, Computational Materials Science 6 , 15 (1996).

Acknowledgements. The computational resources were provided by: (i) the Kelvin2 HPC platform at the NI-HPC Centre (funded by EPSRC and jointly managed by Queen's University Belfast and Ulster University); (ii) the UK National Supercomputing Service ARCHER2, for which access was obtained via the UKCP consortium and funded by EPSRC [EP/X035891/1]. M. S. acknowledges support from Gonville and Caius College. P.A. acknowledges support from SNSF through Post.Doc mobility fellowship P500PN\_206693. We thank William C. Witt and Ilyes Batatia for the useful discussions, and we gratefully acknowledge Atsushi Togo for having provided harmonic and anharmonic force constants computed using density functional theory (PBE functional) from the phononDB-PBE database 44,45 .

Author contributions. M.S. conceived and supervised the project. B.P. developed the automated computational framework and performed the numerical calculations with inputs from P.A., G.C., and M.S. The results were analyzed and organized by B.P. and M.S., with inputs from P.A. All authors contributed to discussing the results, writing and editing the manuscript.

Competing interests. C.G. has equity interest of Symmetric Group LLP that licenses force fields commercially and also in Ångstrom AI. The other authors declare that they have no competing interest.
Balázs Póta Paramvir Ahlawat Gábor Csányi Michele Simoncelli
Thermal Conductivity Predictions with Foundation Atomistic Models
Advances in machine learning have led to the development of foundation models
for atomistic materials chemistry, enabling quantum-accurate descriptions of
interatomic forces across diverse compounds at reduced computational cost.
Hitherto, these models have been benchmarked relying on descriptors based on
atoms' interaction energies or harmonic vibrations; their accuracy and
efficiency in predicting observable and technologically relevant
heat-conduction properties remains unknown. Here, we introduce a framework that
leverages foundation models and the Wigner formulation of heat transport to
overcome the major bottlenecks of current methods for designing heat-management
materials: high cost, limited transferability, or lack of physics awareness. We
present the standards needed to achieve first-principles accuracy in
conductivity predictions through model's fine-tuning, discussing benchmark
metrics and precision/cost trade-offs. We apply our framework to a database of
solids with diverse compositions and structures, demonstrating its potential to
discover materials for next-gen technologies ranging from thermal insulation to
neuromorphic computing.