# Hierarchy of Exchange-Correlation Functionals in Computing Lattice Thermal
Conductivities of Rocksalt and Zincblende Semiconductors

Jiacheng Wei  Key Laboratory of Advanced Materials and Devices for Post-Moore
Chips, Ministry of Education, University of Science and Technology Beijing,
Beijing 100083, China  School of Mathematics and Physics, University of
Science and Technology Beijing, Beijing 100083, China  Zhonghao Xia  Key
Laboratory of Advanced Materials and Devices for Post-Moore Chips, Ministry of
Education, University of Science and Technology Beijing, Beijing 100083, China
School of Mathematics and Physics, University of Science and Technology
Beijing, Beijing 100083, China  Yi Xia  Department of Mechanical and Materials
Engineering, Portland State University, Portland, OR 97201, USA  Jiangang He
<EMAIL_ADDRESS>Key Laboratory of Advanced Materials and Devices for
Post-Moore Chips, Ministry of Education, University of Science and Technology
Beijing, Beijing 100083, China  School of Mathematics and Physics, University
of Science and Technology Beijing, Beijing 100083, China

###### Abstract

Lattice thermal conductivity ($\kappa_{\rm L}$) is a crucial characteristic of
crystalline solids with significant implications for thermal management,
energy conversion, and thermal barrier coating. The advancement of
computational tools based on density functional theory (DFT) has enabled the
effective utilization of phonon quasi-particle-based approaches to unravel the
underlying physics of various crystalline systems. While the higher order of
anharmonicity is commonly used for explaining extraordinary heat transfer
behaviors in crystals, the impact of exchange-correlation (XC) functionals in
DFT on describing anharmonicity has been largely overlooked. The XC functional
is essential for determining the accuracy of DFT in describing interactions
among electrons/ions in solids and molecules. However, most XC functionals in
solids focus on ground state properties that mainly involve the harmonic
approximation, neglecting temperature effects, and their reliability in
studying anharmonic properties remains insufficiently explored. In this study,
we systematically investigate the room-temperature $\kappa_{\rm L}$ of 16
binary compounds with rocksalt and zincblende structures using various XC
functionals such as local density approximation (LDA), Perdew-Burke-Ernzerhof
(PBE), revised PBE for solid and surface (PBEsol), optimized B86b functional
(optB86b), revised Tao-Perdew-Staroverov-Scuseria (revTPSS), strongly
constrained and appropriately normed functional (SCAN), regularized SCAN
(rSCAN) and regularized-restored SCAN (r2SCAN) in combination with different
perturbation orders, including phonon within harmonic approximation (HA) plus
three-phonon scattering (HA+3ph), phonon calculated using self-consistent
phonon theory (SCPH) plus three-phonon scattering (SCPH+3ph), and SCPH phonon
plus three- and four-phonon scattering (SCPH+3,4ph). Our results show that the
XC functional exhibits strong entanglement with perturbation order and the
mean relative absolute error (MRAE) of the computed $\kappa_{\rm L}$ is
strongly influenced by both the XC functional and perturbation order, leading
to error cancellation or amplification. The minimal (maximal) MRAE is achieved
with revTPSS (rSCAN) at the HA+3ph level, SCAN (r2SCAN) at the SCPH+3ph level,
and PBEsol (rSCAN) at the SCPH+3,4ph level. Among these functionals, PBEsol
exhibits the highest accuracy at the highest perturbation order. The SCAN-
related functionals demonstrate moderate accuracy but are suffer from
numerical instability and high computational costs. Furthermore, the different
impacts of quartic anharmonicity on $\kappa_{\rm L}$ in rocksalt and
zincblende structures are identified by all XC functionals, attributed to the
distinct lattice anharmonicity in these two structures. These findings serve
as a valuable reference for selecting appropriate functionals for describing
anharmonic phonons and offer insights into high-order force constant
calculations that could facilitate the development of more accurate XC
functionals for solid materials.

##  I Introduction

Lattice thermal conductivity is a fundamental characteristic of crystalline
solids and holds significant importance in various modern technologies. It
plays a crucial role in essential processes such as heat dissipation in
integrated circuit chips [1, 2], direct energy conversion between thermal and
electrical energy [3, 4], and safeguarding devices using thermal barrier
coatings [5, 6, 7]. Therefore, the identification of materials with the
desired $\kappa_{\rm L}$ values is of significant interest for various
technical applications. However, the synthesis and experimental measurement of
$\kappa_{\rm L}$ for a specific compound can be a time-consuming and costly
process. With the rapid advancements in computer science and simulation
theory, many simulation methods have become more efficient in predicting
$\kappa_{\rm L}$ of a wide range of crystalline materials. The accuracy of
these predictions is highly dependent on the methodology utilized. Different
approaches may yield conflicting conclusions [8, 9, 10].

Currently, there are primarily two methods used for calculating $\kappa_{\rm
L}$ in literature: molecular dynamics (MD), which includes both equilibrium
and non-equilibrium MD [11, 12], and the Boltzmann’s kinetic approach [13]. MD
accounts for all levels of anharmonic terms and can analyze a wide range of
systems such as amorphous materials, nanostructures, and fluids. However, MD
simulations require extensive computational time and do not consider quantum
effects [14]. On the other hand, Boltzmann’s kinetic approach offers improved
computational efficiency for studying bulk crystals, but its precision depends
on the order of perturbation in phonon-phonon interactions [14]. For instance,
the Peierls-Boltzmann method, which is based on three-phonon interactions,
tends to overestimate $\kappa_{\rm L}$ in materials like BAs [8], with a more
accurate value achieved by incorporating four-phonon scattering [15]. Similar
challenges arise in materials like AgCl and HgTe [10]. However, the
computational expense of including higher-order anharmonic effects increases
substantially [16, 17].

Another critical factor that greatly influences the value of the calculated
$\kappa_{\rm L}$ is the XC functional utilized in DFT. It is widely recognized
that the XC functional is essential in ensuring the precision and reliability
of DFT calculations. It is worth noting that various functionals offer
distinct representations of the gradual behavior with respect to density
gradient, even when operating at the same rung of the Jacob’s ladder [18].
Therefore, it is expected that force constants may be greatly influenced by
the chosen XC functional, ultimately impacting $\kappa_{\rm L}$ substantially
[19]. Nevertheless, the development of XC functionals for all properties poses
a considerable challenge, as most functionals can only yield precise results
for certain properties. With significant dedication to the development of XC
functionals over the course of more than 50 years, more than 200 XC
functionals have been developed to date. These functionals have been tailored
to address a diverse range of properties and scenarios [20]. Within the solid-
state research community, the most commonly utilized XC functionals include
LDA [21], PBE [22], PBEsol [23], and revTPSS [24]. Recently, the SCAN
functional has been developed with meticulous attention to satisfying all 17
known exact constraints [25]. This functional has demonstrated excellent
performance across various properties of solids, including lattice constants,
polymorph stabilities, phonon dispersion, and more [25, 26, 27, 28]. However,
subsequent studies have revealed that the SCAN functional exhibits numerical
instability and has the potential to lead to divergence in self-consistent
calculations in certain instances [29, 30]. To address this issue, the
regularized SCAN (rSCAN) functional was developed, which overcomes the
numerical stability problem but at the cost of breaking certain constraints
[29]. However, the rSCAN functional has been found to exhibit larger errors
compared to SCAN in specific systems. To mitigate these errors, the
regularized-restored SCAN (r2SCAN) functional was introduced, which achieves a
better balance between accuracy and efficiency [30]. While the performance of
these XC functionals has been assessed in our previous work regarding
anharmonicity within the quasi-harmonic approximation [31], their performances
in computing the $\kappa_{\rm L}$ have not been thoroughly tested yet.

The casual observations in literature have already demonstrated the importance
of the XC functional in the calculation of $\kappa_{\rm L}$. For example, Jain
and McGaughey explored the influence of XC functionals on the calculated
$\kappa_{\rm L}$ for pure crystalline silicon and determined that the relative
error in calculated $\kappa_{\rm L}$ can reach up to 17 % [32]. LDA, PBE,
PBEsol, and PW91 [33] underestimate $\kappa_{\rm L}$, while BLYP [34, 35]
overestimates $\kappa_{\rm L}$ by 12 %. Qin et al. investigated the
$\kappa_{\rm L}$ of graphene using ten XC functionals and found that the
$\kappa_{\rm L}$ values ranged from 1396 to 4376 Wm-1K-1, depending on the
specific XC functional employed [32]. Taheri et al. achieved higher values of
$\kappa_{\rm L}$ for graphene, ranging from 5442 to 8677 Wm-1K-1 by utilizing
five different XC functionals and two types of pseudopotentials [36, 37].
Recently, Han and Ruan discovered that four-phonon scattering plays a
significant role in graphene, and the discrepancies found in the literature
regarding the calculated $\kappa_{\rm L}$ are attributed to convergence issues
in the grid sampling of the Brillouin zone and the presence of higher orders
of anharmonicity [38]. In a separate study, Dongre et al. determined that the
calculated $\kappa_{\rm L}$ of GaP using PBE is 1.8 times larger than that
calculated value using LDA [39]. Xia et al. investigated the $\kappa_{\rm L}$
of Tl3VSe4 using the PBEsol functional [40], while Jain conducted similar
calculations for this compound using PBE [41]. The discrepancies in the
calculated $\kappa_{\rm L}$ values obtained using these two XC functionals
have led to different explanations regarding the origin of the ultralow
$\kappa_{\rm L}$ in this compound [40, 41]. It is worth noting that the
interplay between the XC functional and the perturbation order of
anharmonicity complicates the underlying physics of $\kappa_{\rm L}$, as
errors introduced by the XC functional and anharmonicity order may
inadvertently cancel each other out, resulting in better agreement with
experimental values when using low-order anharmonicity approaches, as
demonstrated in previous studies [10].

In this study, we systematically evaluate the performance of three commonly
used functionals (LDA, PBE, and PBEsol) in solids, one popular van der Waals
density functional (optB86b) [42], and four meta-GGA functionals (revTPSS,
SCAN, rSCAN, and r2SCAN) in predicting the room-temperature $\kappa_{\rm L}$
for 16 binary compounds with rocksalt and zincblende structures at three
perturbation orders, namely HA+3ph, SCPH+3ph, and SCPH+3,4ph. Our results show
that the strong entanglement between the XC functional and the perturbation
order complicates the comparison between the computed and the experimental
$\kappa_{\rm L}$. In other words, the MRAE of the computed $\kappa_{\rm L}$ is
strongly influenced by both the XC functional and the perturbation order,
leading to error cancellation or amplification. The minimal (maximal) MRAE is
achieved with revTPSS (rSCAN) at the HA+3ph level, SCAN (r2SCAN) at the
SCPH+3ph level, and PBEsol (revTPSS) at the SCPH+3,4ph level. Overall, PBEsol
is numerically stable and exhibits the highest accuracy in calculating
$\kappa_{\rm L}$ among these functionals. The SCAN-related functionals
demonstrate moderate accuracy but suffer from numerical instability and higher
computational cost. Furthermore, all XC functionals identify the different
impacts of quartic anharmonicity on $\kappa_{\rm L}$ between rocksalt and
zincblende structures, which can be attributed to their distinct coordination
environment and lattice anharmonicity. These findings serve as a valuable
reference for selecting appropriate functionals for describing anharmonic
phonons and offer insights into high-order force constant calculations that
could facilitate the development of more accurate XC functionals for solid
materials.

##  II COMPUTATION DETAILS

All DFT calculations were performed using the projector-augmented wave
(PAW)[43, 44] method within the Vienna Ab initio Simulation Package (VASP)
[45, 46]. Eight XC functionals were utilized, namely LDA, PBE, PBEsol,
optB86b, revTPSS, SCAN, rSCAN, and r2SCAN. With the exception of LDA, the PBE-
version PAW pseudopotential was employed for all other functionals due to the
unavailability of the PAW pseudopotentials for these functionals in VASP.
While this may introduce some errors in these functionals, previous studies on
solids have predominantly used the PBE pseudopotential, with negligible errors
reported in studies utilizing the other functional [47, 48, 49]. The plane-
wave basis set was used with cutoff energies of 520 eV for LDA, PBE, PBEsol,
and optB86b, 600 eV for revTPSS, and 800 eV for SCAN, rSCAN, and r2SCAN.
Additionally, for SCAN, rSCAN, and r2SCAN, the precision parameter was set to
"accurate" (PREC = ACCURATE), nonspherical contributions within the PAW
spheres were included (LASPH = .TRUE.), and additional grids for the
evaluation of augmentation charges (ADDGRID = .TRUE.) were implemented.

The second-order force constants were calculated using the finite displacement
method as implemented in Phonopy code [50] with 4 $\times$ 4 $\times$ 4
supercell and 2 $\times$ 2 $\times$ 2 $k$-points mesh, and a 0.01
Ådisplacement. During the self-consistent calculations, the energy is
converged to $10^{-8}$ eV. The projection operators were evaluated in
reciprocal space (LREAL = .FALSE.) to more accurately compute forces. Since
all atoms are in high-symmetry Wyckoff positions and the only structure
parameter is the lattice constant, the equilibrium structures of all compounds
studied in this work were generated using the experimental lattice constants
at room temperature to eliminate errors caused by thermal expansion. The
effect of supercell size on harmonic phonons was carefully tested. It was
found that a 5 $\times$ 5 $\times$ 5 supercell can produce artificial
imaginary frequencies for some compounds using certain XC functionals, while
the more numerically stable option for all compounds was the 4 $\times$ 4
$\times$ 4 supercell, as shown in Fig. S1. The nonanalytic correction to
phonon dispersions near the $\Gamma$ point is performed by dipole-dipole
interaction [51, 52], based on the PBEsol functional calculated Born effective
charges and macroscopic static dielectric constants using the density function
perturbation theory implementation in VASP.

The third-order and fourth-order force constants were extracted using the
compressive sensing lattice dynamics (CSLD) [53]. Phonon frequency shifts at
finite temperatures were calculated using the SCPH theory [54, 55, 56]. The
Peierls-Boltzmann transport equation was literately solved by sampling with a
uniform 24 $\times$ 24 $\times$ 24 $q$-point mesh for $\kappa_{\rm L}$
calculations with HA+3ph ($\kappa_{\rm 3ph}^{\rm HA}$) and SCPH+3ph
($\kappa_{\rm 3ph}^{\rm SCPH}$) methods, as well as a sparser 16 $\times$ 16
$\times$ 16 $q$-point mesh for $\kappa_{\rm L}$ with SCPH+3,4ph ($\kappa_{\rm
3,4ph}^{\rm SCPH}$) method. The four-phonon scattering process is accelerated
using the sampling method [57]. All other parameters remain consistent with
the previous study [10].

Elastic constants are determined through the structure deformation method as
implemented in Pymatgen [58]. The mean sound velocity ($v_{\rm m}$) is
obtained as the average of the longitudinal ($v_{\rm L}$) and transverse
($v_{\rm T}$) sound velocities, which are calculated based on the density
($\rho$), bulk ($B$), and shear ($G$) modulus [59].

$v_{m}=[\frac{1}{3}(\frac{1}{v_{\rm L}^{3}}+\frac{2}{v_{\rm
T}^{3}})]^{\frac{-1}{3}}$

$v_{\rm T}=\sqrt{\frac{G}{\rho}}$

$v_{\rm L}=\sqrt{\frac{B+\frac{4}{3}G}{\rho}}$

##  III RESULTS AND DISCUSSION

In the kinetic theory [60], $\kappa_{\rm L}$ can be represented as the product
of the heat capacity ($C_{\rm v}$), the square of the phonon group velocity
($v_{\rm g}$), and the phonon relaxation time ($\tau$), as shown in the
equation $\kappa_{\rm L}=\frac{1}{3}C_{\rm v}v_{\rm g}^{2}\tau$. Therefore,
prior to assessing the accuracy of computing $\kappa_{\rm L}$ using eight XC
functionals at three theoretical perturbation orders, a comparison was
conducted on the precision of these functionals regarding phonon spectra,
$v_{\rm g}$, and the Grüneisen parameter ($\gamma$). The Grüneisen parameter
is commonly utilized to characterize anharmonicity and is closely linked to
the scattering rate ($1/\tau$).

###  III.1 Phonon spectrum

The phonon spectrum contains essential information regarding lattice dynamics,
rendering it indispensable for comprehending $\kappa_{\rm L}$, and is also
significantly influenced by temperature. However, the harmonic phonon,
obtained from the harmonic approximation by diagonalizing the dynamical matrix
formed from harmonic interatomic force constants, fails to consider
temperature effects. The impact of temperature on the phonon spectrum can be
effectively simulated using the SCPH theory [54, 61, 56]. In this study, the
phonon spectrum at room temperature is calculated using the first-order
correction based on quartic anharmonicity [62]. The phonon spectra for all the
compounds are analyzed using both harmonic and quartic anharmonic
approximations, with the results computed using 8 XC functionals presented in
Fig. S2-S7. Additionally, Fig. 1(a) illustrates the discrepancy
($\Delta{\omega}_{\rm HA}$) in the optical phonon frequency at the $\Gamma$
point when calculated using different XC functionals compared to the PBEsol
functional. Note the phonon frequency before nonanalytical correction is
adopted here for the sake of simplicity. It is unexpected that the
$\Delta{\omega}_{\rm HA}$ of the rocksalt compounds is greater than that of
the zincblende, even though the rocksalt compounds have a considerably lower
optical frequency at the $\Gamma$ point, indicating that the phonon
frequencies of the zincblende compounds are less sensitive to the XC
functional than rocksalt compounds. Out of all these functionals, only LDA
exhibits negative $\Delta{\omega}_{\rm HA}$ for all the compounds due to the
fact that LDA calculated chemical bonds are elongated (bonding interaction is
weakened) when the experimental volumes utilized in this study. When comparing
the phonon frequencies calculated using the HA with those obtained through the
SCPH method at 300 K, it is observed that the rocksalt structure exhibits a
slight increase in frequency, particularly in the case of three optical
branches. Conversely, the phonon frequency of the zincblende structure remains
relatively unchanged or even decreases slightly in certain instances. In order
to quantify this, we show the frequency shift of the phonon at the $\Gamma$
point from the SCPH to HA method ($\Delta{\omega}_{\rm SCPH-HA}$ =
$\omega_{\rm SCPH}$ \- $\omega_{\rm HA}$) in Fig. 1(b). All the
$\Delta{\omega}_{\rm SCPH-HA}$ of the rocksalt compounds are large and
positive except that of CaO calculated by SCAN, while the $\Delta{\omega}_{\rm
SCPH-HA}$ of the zincblende compounds are generally small or even negative by
some of these functionals. Exceptions are predominantly observed in four
metaGGA functionals: InAs by rSCAN, SiC by revTPSS, rSCAN, and r2SCAN, ZnSe by
revTPSS, SCAN, and rSCAN, ZnTe by r2SCAN, and GaAs by SCAN. This highlights
the differences of the role played by quartic anharmonicity between rocksalt
and zincblende structures. All functionals predict a large positive
$\Delta{\omega}_{\rm SCPH-HA}$ in BaO, while revTPSS is the only functional
predicting a large $\Delta{\omega}_{\rm SCPH-HA}$ in SiC and SCAN is the only
functional having a large and negative $\Delta{\omega}_{\rm SCPH-HA}$ for CaO.
These observations underscore the large difference between metaGGA and
semilocal functionals in characterizing anharmonicity.

Figure 1: (Color online) (a) The HA frequency difference ($\Delta{\omega}_{\rm
HA}$) of the optical phonon frequency at the $\Gamma$ point calculated using
the other functionals with respect to that obtained with PBEsol. (b) The
frequency shift ($\Delta{\omega}_{\rm SCPH-HA}$) between the phonon calculated
by SCPH and HA for all the compounds calculated with 8 XC functionals.

###  III.2 Sound velocity

Phonon group velocity $v_{\rm g}$ is another pivotal metric encapsulating the
speed at which the main heat carrying phonons propagate through a solid
crystalline, determining its thermal transport behavior [63]. Therefore, it is
widely used in both empirical formula and phonon Boltzmann transport theory
for computing $\kappa_{\rm L}$. Herein, $v_{\rm g}$ is approximated by the
mean sound velocity ($v_{\rm m}$), which is calculated from elastic constants,
see the “Computational Details” section for details. The $v_{\rm m}$ is a
critical parameter that measures the speed of the primary heat-carrying
phonons traveling through a solid crystalline material, influencing its
thermal conductivity behavior [63]. Figure 2 illustrates the comparison
between our calculated $v_{\rm m}$ values of all the compounds with the
experimental values measured at room temperature, utilizing all the XC
functionals examined in this study. It should be noted that the majority of
$v_{\rm m}$ values are derived from direct experimental measurements, while a
small subset are calculated using the experimental elastic constants. The
calculated $v_{\rm m}$ values of nearly all compounds, as determined using all
the XC functionals, are consistently lower than the corresponding experimental
values, with the notable exception of BaO. This discrepancy is unexpected,
considering that our calculations are focused solely on the intrinsic
compounds and do not account for potential factors such as defects and grain
boundaries, which are known to exert a substantial slowing down of $v_{\rm m}$
[63]. The MRAE of these compounds exhibit only slight variations across the XC
functionals. The smallest relative errors are observed at CaO and SiC, while
the largest errors are found at NaF and NaCl. Specifically, the smallest and
largest MRAE values are observed in PBE (7.7 %) and SCAN (12.2 %) functionals,
respectively. It is noteworthy that the calculated $v_{\rm m}$ values for NaCl
and NaF show significant deviations from the experimental values. SCAN, among
all the XC functionals, identifies PbS, ZnSe, and InSb as having the largest
relative errors, leading to the highest MRAE, despite its smallest relative
error in NaCl. In general, the variation in $v_{\rm m}$ across various
functionals is relatively small presumably due to the fact that the elastic
constant is the second derivative of energy with respect to displacement,
which is harmonic and can be well described by all these functionals.

Figure 2: (Color online). The relative error of the calculated mean sound
velocity ($v_{\rm m}^{\rm cal.}$) with respect to experiment values ($v_{\rm
m}^{\rm exp.}$). The $v_{\rm m}^{\rm exp.}$ are either taken from Ref. 64 or
calculated using the experimental elastic constants in Ref. 64.

###  III.3 Grüneisen parameter

The Grüneisen parameter $\gamma$, which is defined as the rate of change of
phonon frequency in response to variations in crystal volume [65, 66], is
commonly utilized to quantify the degree of anharmonicity and can be employed
to compute $\kappa_{\rm L}$ and the thermal expansion coefficient of solids
[67, 67, 68]. In our previous study [31], we evaluated the performance of six
XC functionals on computing the average value of $\gamma$ for 17 binary
compounds with rocksalt, zincblende, and fluorite structures. The average
$\gamma$ was determined using the Phonopy code [50] based on the HA phonon
frequencies of 11 volumes around the equilibrium volume at 0 K. In this study,
the total $\gamma$ is computed as a weighted sum of the mode contributions, as
implemented in the ShengBTE package [69], incorporating both the HA and SCPH
force constants. Consequently, two $\gamma$ values are determined for each
compound, specifically at 0 and 300 K. Figure 3(a) illustrates the comparison
between the experimental and calculated $\gamma$ values for all the rocksalt
and zincblende compounds examined in this work. The experimental $\gamma$
values primarily originate from Ref. 70, 71. For the calculated $\gamma$
values based on the HA force constants, the MRAEs of the LDA, PBE, PBEsol,
optB86b, revTPSS, SCAN, rSCAN, and r2SCAN functionals are 20.18 %, 15.49 %,
18.04 %, 16.46 %, 18.32 %, 22.73 %, 15.92 %, and 14.68 %, respectively. It
should be noted that both SCAN and r2SCAN exhibit significant discrepancies
compared to the previous results [31], which can be attributed to differences
in the dataset as well as the impact of lattice constants, i.e., the use of
DFT-optimized versus experimental lattice constants. Notably, the compounds
that resulted in large relative errors with the SCAN and r2SCAN functionals in
our previous work (AlP and InP) were excluded from this study due to the
presence of imaginary frequencies when employing the experimental lattice
constants. Additionally, the substantial difference observed in InSb may be
attributed to the discrepancy in the employed lattice constants. As depicted
in Fig. 3(b), the replacement of the HA force constant with the SCPH one leads
to a reduction in all MRAEs, particularly for the LDA, PBEsol, and SCAN
functionals, primarily due to improved agreement for PbS. Among all the XC
functionals considered, r2SCAN exhibits the lowest MRAE (10.24 %), followed by
PBEsol (12.08 %) and LDA (12.36 %). With the exception of LDA and rSCAN, all
functionals significantly underestimate the $\gamma$ value of InSb. It is
evident that PBE, PBEsol, optB86b, and revTPSS functionals yield smaller MRAEs
for the $\gamma$ values of rocksalt compounds compared to zincblende ones,
likely due to the generally larger $\gamma$ values observed in rocksalt
structure than zincblende [72].

Figure 3: (Color online). (a) and (b) are the relative errors of the total
Grüneisen parameters ($\gamma$) at 0 K and 300 K, respectively, when compared
to experimental values. Gray color represents the calculation with imaginary
frequency at 0 K (a) and at room temperature (b).

###  III.4 Scattering rates

The main process of phonon scattering in an intrinsic semiconductor is the
phonon-phonon scattering, predominantly involving a three-phonon process [73].
However, four-phonon scattering may also be significant in certain systems
[15, 74]. The scattering rate (1/$\tau$) characterizes the strength of
scattering among phonons and is a key factor in the phonon Boltzmann transport
equation, playing a dominant role in determining $\kappa_{\rm L}$ [75]. In
this study, scattering rates associated with the three-phonon and four-phonon
processes are calculated using the single-mode relaxation approximation [75].
As depicted in Fig. S8-S13, all XC functionals indicate that three-phonon
scattering prevails in the scattering process across all compounds.
Specifically, the three-phonon scattering rates in the low-frequency region, a
key factor in heat resistance, are higher in rocksalt compounds ($\sim$ 10
ps-1) compared to those with the zincblende structure ($\sim$ 1 ps-1),
consistent with previous research and reflecting the stronger anharmonicity
(larger $\gamma$) in rocksalt structure compared to zincblende [10, 76].
Furthermore, all XC functionals suggest that the four-phonon scattering rates
in zincblende compounds are relatively higher than those in rocksalt
compounds, aligning with previous observations based on PBE functional [10],
due to the relatively smaller three-phonon scattering rates in zincblende
compounds.

By comparing the phonon scattering rates of a compound computed using
different XC functionals, it is evident that SCAN-related functionals have a
higher likelihood of predicting a different scattering rate from the other XC
functionals, as illustrated in Fig. S9-S11. For example, rSCAN and r2SCAN
consistently show significantly lower three-phonon scattering rates for NaCl
and NaF compared to other XC functionals. SCAN predicts notably higher three-
phonon scattering rates for PbS and SiC in the low-frequency region. rSCAN
demonstrates significantly higher three-phonon scattering rates for NaCl, NaF,
and BAs in comparison to other functionals, while r2SCAN predicts a notably
high three-phonon scattering rate for ZnTe in the low-frequency region.
Furthermore, rSCAN and r2SCAN predict exceptionally high four-phonon
scattering rates for NaF in the frequency range of 150 $\sim$ 250 cm-1, with
rSCAN also showing a remarkably high four-phonon scattering rate for NaCl.
These findings underscore the complexity and challenges associated with
computing three- and four-phonon scattering rates, especially for metaGGA
functionals.

###  III.5 Lattice thermal conductivity at room temperature

The room-temperature $\kappa_{\rm L}$ of 16 binary compounds with rocksalt and
zincblende structures are computed using 8 commonly used XC functionals in
solids at three different perturbation orders: HA+3ph, SCPH+3ph, and
SCPH+3,4ph. Currently, the mainstream method of computing $\kappa_{\rm L}$ is
HA+3ph, particularly for complex systems. SCPH+3ph accounts for phonon
frequency shifts due to temperature effect, necessitating the calculation of
higher-order force constants, and is particularly relevant for systems with
multiple low-frequencies phonon branches. SCPH+3,4ph is considered as the most
accurate method in this study, as it incorporates both frequency shifts due to
finite temperature and four-phonon scattering.

Figure 4: (Color online). (a)-(c) are relative errors of the calculated
$\kappa_{\rm L}$ at HA+3ph ($\kappa_{\rm 3ph}^{\rm HA}$), SCPH+3ph
($\kappa_{\rm 3ph}^{\rm SCPH}$), and SCPH+3,4ph ($\kappa_{\rm 3,4ph}^{\rm
SCPH}$) three levels with respect to the experimental values at room
temperature, respectively. Gray color indicates the phonon with imaginary
frequencies at 0 K or 300 K that prohibits $\kappa_{\rm L}$ calculations.
Experimental $\kappa_{\rm L}$ values are adopted from 77, 78, 79, 80, 81, 82,
83, 84, 85, 86.

The relative error of the calculated $\kappa_{\rm L}$ for 16 compounds at the
HA+3ph level with 8 XC functionals is depicted in Fig. 4(a). It is important
to note that some compounds exhibit significant imaginary frequencies either
at the experimentally determined room-temperature lattice constants or due to
artificial numerical issues associated with specific functionals, as
previously observed in our lattice constant study [31]. Some of these
imaginary frequencies cannot be stabilized at room temperature using the SCPH
method either. The presence of imaginary frequencies at 0 K prevents the
calculation of $\kappa_{\rm 3ph}^{\rm HA}$, while those at room temperature
impede the determination of $\kappa_{\rm 3ph}^{\rm SCPH}$ and $\kappa_{\rm
3,4ph}^{\rm SCPH}$. The minimum and maximum MRAEs of the calculated
$\kappa_{\rm L}$ at the HA+3ph level are 29.75 % using revTPSS and 49.18 %
using rSCAN, primarily attributed to the overestimation of $\kappa_{\rm L}$ in
BAs, which is a compound known for strong four-phonon scattering [9]. Among
all the functionals considered in this work, only LDA and SCAN show
significant underestimation of $\kappa_{\rm L}$ for compounds with rocksalt
structure. Additionally, PBE underestimates $\kappa_{\rm L}$ only for CaO and
PbS, which is contrary to the previous finding that PBE underestimates
$\kappa_{\rm 3ph}^{\rm HA}$ for most rocksalt compounds [10]. This discrepancy
arises from the use of experimentally determined room-temperature lattice
constants in this study, whereas the previous work utilized the PBE optimized
lattice constants at 0 K [10].

As illustrated in Fig. S14, the experimental volumes measured at room
temperature are smaller than the PBE optimized ones at 0 K for all compounds,
resulting in shorter bond lengths and stronger bond strength, leading to
enlarged $\kappa_{\rm L}$ values. OptB86b, similarly, overestimates the
$\kappa_{\rm L}$ for all compounds except PbS. However, this trend varies
among different functionals. For instance, while revTPSS exhibits larger
optimized volumes compared to optB86b (see Fig. S14), the overestimation of
$\kappa_{\rm L}$ by optB86b surpasses that of revTPSS. The majority of XC
functionals analyzed in this study have a tendency to overestimate the
$\kappa_{\rm L}$ values for compounds with zincblende structure, with BAs and
ZnTe exhibiting the most significant overestimations. This is consistent with
the results reported by Xia et al. utilizing the PBE functional with optimized
lattice constants [10]. Nonetheless, there are exceptions such as ZnTe with
r2SCAN, BAs with rSCAN and revTPSS. This is particularly surprising as it is
found that the calculated $\kappa_{\rm 3ph}^{\rm HA}$ is typically twice as
large as the experimentally measured $\kappa_{\rm L}$ due to its relatively
strong four-phonon scattering [15, 87]. Furthermore, as depicted in Fig. 5 and
2, the revTPSS and rSCAN functionals exhibit similar phonon spectra and
$v_{\rm m}$ to other functionals, such as PBEsol, but display a stronger
three-phonon scattering rate at low-frequency regions due to their larger
weighted phase space [88]. Interestingly, despite being derived from revTPSS
[25], SCAN and r2SCAN do not underestimate the three-phonon scattering rate of
BAs and exhibit very similar MRAE. This suggests the importance of adhering to
the constraints imposed by SCAN and r2SCAN. Among these functionals, revTPSS
functional has the smallest relative error ($\kappa_{\rm 3ph}^{\rm HA}$) for
all compounds except ZnTe, resulting in the lowest MRAE (29.75 %). Although
r2SCAN exhibits a similar $v_{\rm m}$ as PBEsol for ZnTe, as shown in Fig. 2,
it has much stronger three-phonon scattering in the low-frequency region, see
Fig. 5, resulting in good agreement with the experimental $\kappa_{\rm L}$. It
is worth noting that this agreement is likely a result of error cancellation,
as the phonon frequency shifts and scattering rates caused by quadratic
anharmonicity have significant effects on $\kappa_{\rm L}$, leading to
fortuitous agreement between theory and experiment.

Figure 5: (Color online) (a)-(c) are the three-phonon scattering rate,
weighted phase space, and phonon spectrum of BAs, respectively. (d)-(f) are
the three-phonon scattering rate, weighted phase space, and phonon spectrum of
ZnTe, respectively.

The Fig. 4(b) displays the calculated $\kappa_{\rm L}$ values by eight XC
functionals at the SCPH+3ph level. Among these functionals, the SCAN-related
functionals yield both the highest and lowest MRAEs, with SCAN exhibiting the
smallest MRAE of 30.57 % and r2SCAN exhibiting the largest MRAE of 73.11 %.
Following closely behind is rSCAN with a MRAE of 66.22 %. It is evident that
all XC functionals show a significant increase in MRAE compared to the
$\kappa_{\rm L}$ values calculated using HA+3ph. The largest increase occurs
in rSCAN and r2SCAN, with their MRAEs increasing from 49.18 % to 66.22 % and
from 35.23 % to 73.11 %, respectively. The relative errors of each compound
obtained from rSCAN and r2SCAN are notably similar, with the exception of BAs,
which can be attributed to the strong three-phonon scattering predicted by
rSCAN, see Fig. 5. The increase in MRAE from HA+3ph to SCPH+3ph is mainly
attributed to the shift of the phonon frequency towards higher energy (see
Fig. 1) and the prior-existing overestimation of $\kappa_{\rm L}$ by HA+3ph
for most functionals. On the other hand, LDA provides a clearer demonstration
of the alleviation of the $\kappa_{\rm L}$ underestimation from HA+3ph to
SCPH+3ph. Notably, the overestimation of $\kappa_{\rm L}$ for BAs and ZnTe is
slightly reduced due to the slight decrease in phonon frequency caused by
SCPH, as shown in Fig. 1(b). The effects of temperature on $\kappa_{\rm L}$
can be more accurately assessed for each compound by calculating the ratio
between $\kappa_{\rm 3ph}^{\rm SCPH}$ and $\kappa_{\rm 3ph}^{\rm HA}$. As
illustrated in Fig. 6(a), the majority ratios of $\kappa_{\rm SCPH}^{\rm
3ph}/\kappa_{\rm 3ph}^{\rm HA}$ exceed 1, indicating an enlargement of
$\kappa_{\rm L}$ from HA+3ph to SCPH+3ph. Notably, compounds with the rocksalt
structure exhibit larger ratios compared to those with the zincblende
structure, consistent with previous findings using the PBE functional [10],
attributed to the lower coordination number and weaker anharmonicity of the
zincblende structure [72]. Furthermore, the degree of enhancement varies with
the functional used for a specific compound. For instance, the $\kappa_{\rm
3ph}^{\rm SCPH}/\kappa_{\rm 3ph}^{\rm HA}$ ratios for NaCl are 1.67, 1.0,
1.48, 1.28, 1.68, 2.33, 0.79, 1.25 for LDA, PBE, PBEsol, optB86b, revTPSS,
SCAN, rSCAN, and r2SCAN, respectively.

Notably, all of these functionals exhibit a large $\kappa_{\rm 3ph}^{\rm
SCPH}/\kappa_{\rm 3ph}^{\rm HA}$ ratio for BaO, but a smaller ratio for other
alkaline earth metal oxides $M$O ($M$ = Mg, Sr, and Ca). This difference is
likely attributed to the greater phonon frequency shift from SCPH to HA in BaO
compared to the other $M$O compounds, as shown in Fig. 1(b). Except SCAN-
related functionals, all the other functionals demonstrate a similar trend in
the $\kappa_{\rm 3ph}^{\rm SCPH}/\kappa_{\rm 3ph}^{\rm HA}$ ratio for
compounds with the zincblende structure, where most compounds have a ratio of
around 1. Conversely, SCAN-related functionals display more variability in
this ratio. Specifically, SCAN and rSCAN indicate a ratio greater than 1 for
ZnSe, while rSCAN and r2SCAN show ratios greater than 1 for BAs and ZnTe,
respectively. It is worth noting that rSCAN deviates from the trend observed
in other functionals when calculating the $\kappa_{\rm 3ph}^{\rm
SCPH}/\kappa_{\rm 3ph}^{\rm HA}$ ratio of NaCl in rocksalt compounds: while
all other functionals have a ratio of around 1, rSCAN yields a ratio of 0.79.
This suggests that SCAN-related functionals can not provide consistent results
when calculating $\kappa_{\rm 3ph}^{\rm HA}$ and $\kappa_{\rm 3ph}^{\rm
SCPH}$.

Figure 6: (Color online). Comparison of $\kappa_{\rm L}$ calculated using
three methods for the 16 binary with rocksalts and zincblende structure at 300
K. (a)-(c) are the $\kappa_{\rm 3ph}^{\rm SCPH}/\kappa_{\rm 3ph}^{\rm HA}$,
$\kappa_{\rm 3,4ph}^{\rm SCPH}/\kappa_{\rm 3ph}^{\rm SCPH}$, and $\kappa_{\rm
3,4ph}^{\rm SCPH}/\kappa_{\rm 3ph}^{\rm HA}$, respectively. Gray color
represents the calculation that has imaginary frequency at 0 K or 300 K, and
hence the $\kappa_{\rm L}$ calculation is prohibited.

Fig. 4(c) displays the calculated $\kappa_{\rm L}$ at the SCPH+3,4ph level,
which represents the most accurate method utilized in this study. The maximum
and minimum MRAEs are 35.06 % with rSCAN and 20.03 % with PBEsol,
respectively. A comparison with Fig. 4(b) reveals a notable decrease in MRAEs
attributed to additional phonon-phonon scattering from four-phonon
interactions, which is particularly significant for systems with relatively
weak three-phonon scattering like BAs [15]. It is evident that four-phonon
scattering significantly impacts the calculated $\kappa_{\rm L}$ of the
compounds with the zincblende structure, contrary to the effect of SCPH as
discussed earlier. In stark contrast to the results from HA+3ph and SCPH+3ph,
the LDA-calculated $\kappa_{\rm 3,4ph}^{\rm SCPH}$ for all the examined
compounds in this study are lower than the experimental $\kappa_{\rm L}$, with
larger relative errors for the compounds with the rocksalt structure. This
discrepancy is attributed to the overbinding issue of LDA, which results in
smaller volumes (shorter bond lengths) compared to these at room temperature,
see Fig. S14, while the experimental lattice constants at room temperature are
employed in this work, which leads to weaker chemical bonds in our LDA
calculations. A similar trend is observed in PBEsol, which exhibits slightly
larger equilibrium volumes than LDA at 0 K but still smaller than room-
temperature values for most compounds, see Fig. S14. The agreement with
experimental lattice constants will be significantly improved if the thermal
expansion is considered [31]. Interestingly, despite PBE-optimized volumes at
0 K being larger than experimental volumes at room temperature (see Fig. S14),
resulting in enhanced chemical bonds when utilizing room-temperature lattice
constants, PBE only overestimates the $\kappa_{\rm L}$ of three compounds
(BaO, MgO, and PbTe). This implies a more complex relationship between bond
length and $\kappa_{\rm L}$ across various XC functionals. Among these
functionals, only LDA, PBEsol, and SCAN barely address the issue of
overestimation for BaO, showing the best agreement with experimental values
among these functional [77]. While PBE, optB86b, and revTPSS largely
overestimate the $\kappa_{\rm 3,4ph}^{\rm SCPH}$ of PbTe, other functionals
achieve good agreement with the experimental value if no imaginary frequencies
are present at the SCPH level. Additionally, we observed that certain
functionals yield significantly lower $\kappa_{\rm 3,4ph}^{\rm SCPH}$ values
for specific compounds compared to other XC functionals. Examples include SCAN
for CaO, LDA for PbS, revTPSS and rSCAN for BAs, SCAN for ZnSe, and r2SCAN for
ZnTe. Conversely, all functionals largely underestimate the $\kappa_{\rm L}$
of InAs and InSb except SCAN, suggesting potential inaccuracies in the
experimental values or the influence of other factors such as lattice
constants. Overall, the performance of the PBEsol functional is notable in
most of these compounds, with a slightly smaller relative error observed for
compounds with a rocksalt structure compared to those with a zincblende
structure.

The effect of four-phonon scattering on $\kappa_{\rm L}$ can be assessed by
examining the ratio of $\kappa_{\rm 3,4ph}^{\rm SCPH}/\kappa_{\rm 3ph}^{\rm
SCPH}$. Analysis presented in Fig. 6(b) reveals that the compounds with the
zincblende structure are more significantly affected by the four-phonon
scattering compared to those with the rocksalt one, with the exception of SiC,
BP, and GaAs across all functionals. This observation suggests that four-
phonon scattering is relatively weaker in rocksalt compounds because it has
strong three-phonon scattering. Surprisingly, all the functionals find that
the ratios of $\kappa_{\rm 3,4ph}^{\rm SCPH}/\kappa_{\rm 3ph}^{\rm SCPH}$ for
InSb, InAs, and ZnTe are comparable or even smaller than that of BAs, a well-
known compound with strong four-phonon scattering [15], indicating a more
pronounced four-phonon interaction effect. In contrast to other functionals,
rSCAN predicts that the rocksalt compounds NaF and NaCl also exhibit strong
four-phonon scattering, resulting in a significant underestimation of
$\kappa_{\rm L}$. A similar result is predicted by r2SCAN for NaF, but not for
NaCl. The four-phonon scattering rates of these two compounds as predicted by
rSCAN are illustrated in Fig. 7. When compared to the PBEsol functional, rSCAN
demonstrates significantly stronger four-phonon scattering in the low-
frequency region for both NaF and NaCl. In contrast, r2SCAN only predicts high
four-phonon scattering in NaF, aligning with the observed trend in
$\kappa_{\rm L}$.

Figure 7: (Color online) (a) and (b) are the four-phonon scattering rates
calculated using PBEsol, rSCAN, and r2SCAN for NaF and NaCl, respectively.

The $\kappa_{\rm 3,4ph}^{\rm SCPH}/\kappa_{\rm 3ph}^{\rm HA}$ ratio provides
valuable insight into the relative contributions of phonon renormalization at
room temperature and the four-phonon scattering in the calculation of
$\kappa_{\rm L}$. A $\kappa_{\rm 3,4ph}^{\rm SCPH}/\kappa_{\rm 3ph}^{\rm HA}$
ratio of 1 indicates a perfect cancellation between these two effects, while a
ratio of less than 1 suggests that the reduction of $\kappa_{\rm L}$ due to
four-phonon scattering outweighs the enhancement from phonon renormalization,
although the latter can occasionally lead to a decrease in $\kappa_{\rm L}$ as
well. Analysis of the data presented in Fig. 6(c) indicates that all
zincblende compounds demonstrate a $\kappa_{\rm 3,4ph}^{\rm SCPH}/\kappa_{\rm
3ph}^{\rm HA}$ ratio below 1 across all functionals, which is attributed to
the smaller or negative phonon frequency shifts at room temperature and the
heightened four-phonon scattering characteristic of the zincblende structure,
aligning with previous findings [10]. Among these compounds, SiC and BP stand
out as the only zincblende compounds with a $\kappa_{\rm 3,4ph}^{\rm
SCPH}/\kappa_{\rm 3ph}^{\rm HA}$ ratio approaching 1. Conversely, all rocksalt
compounds exhibit a $\kappa_{\rm 3,4ph}^{\rm SCPH}/\kappa_{\rm 3ph}^{\rm HA}$
ratio close to or greater than 1. Notable exceptions to this trend are
primarily observed in NaF and NaCl across nearly all functionals. At the same
time, these two compounds are predicted to have abnormally smaller
$\kappa_{\rm 3,4ph}^{\rm SCPH}/\kappa_{\rm 3ph}^{\rm HA}$ ratio than other
functionals by rSCAN. It is evident that the rocksalt compounds generally
display similar values for $\kappa_{\rm 3,4ph}^{\rm SCPH}/\kappa_{\rm
3ph}^{\rm HA}$ ratio, whereas the zincblende compounds show large frustration
across these compounds. BaO is identified as a notable case where the
inclusion of four-phonon scattering leads to a significant reduction in
$\kappa_{\rm L}$, surpassing the overestimation produced by SCPH+3ph.

###  III.6 Dependence of lattice thermal conductivity on volume

Our previous study demonstrates that there is a significant disparity between
the computed room-temperature lattice constants using QHA and the
corresponding experimental values for many of these XC functionals [31].
However, due to the considerable cost associated with computing thermal
expansion through anharmonic approximation for such an extensive dataset, this
work employs the room-temperature experimental lattice constants to mitigate
the influence of lattice thermal expansion on $\kappa_{\rm L}$. As illustrated
in Fig. S14, most compounds exhibit significant discrepancies between the
optimized volumes at 0 K and the experimental ones at room temperature.
Consequently, the chemical bonding of these compounds may be either
strengthened or weakened by the smaller or larger volumes in our calculations.
A stronger (weaker) chemical bond could result in a higher (lower) $v_{\rm g}$
and, consequently, an increased (decreased) $\kappa_{\rm L}$ [89]. To further
elucidate the impact of volume on the $\kappa_{\rm L}$ of these compounds, we
present the correlation between $\kappa_{\rm L}$ and volume for two
represented compounds ZnTe and PbTe in Fig. 8. As demonstrated above, PBEsol
exhibits exceptional performance in computing $\kappa_{\rm L}$ across all
perturbation orders. Consequently, PBEsol is utilized to explore the
dependence of $\kappa_{\rm L}$ on volume. ZnTe is a zincblende compound, and
our calculated $\kappa_{\rm L}$ is notably higher than that reported in a
previous study [10]. On the other hand, PbTe serves as a representative of the
rocksalt compound with significant anharmonicity. It is observed that the
$\kappa_{\rm L}$ increases monotonically as the volume decreases
(corresponding to an increase in pressure) for both ZnTe and PbTe.
Interestingly, while the $\kappa_{\rm 3,4ph}^{\rm SCPH}$ of ZnTe is nearly
linearly dependent on volume, the changes of $\kappa_{\rm 3ph}^{\rm HA}$ and
$\kappa_{\rm 3ph}^{\rm SCPH}$ with volume is nonlinear. Additionally, both
$\kappa_{\rm 3ph}^{\rm HA}$ and $\kappa_{\rm 3ph}^{\rm SCPH}$ are much higher
than $\kappa_{\rm 3,4ph}^{\rm SCPH}$, and the difference between them widens
as the volume decreases, suggesting a significant enhancement in four-phonon
scattering. Notably, a crossover is observed between $\kappa_{\rm 3ph}^{\rm
HA}$ and $\kappa_{\rm 3ph}^{\rm SCPH}$ at a volume close to the experimental
one. When the volume is smaller than this critical volume, $\kappa_{\rm
3ph}^{\rm SCPH}$ is lower than $\kappa_{\rm 3ph}^{\rm HA}$. In contrast to
ZnTe, the behavior of PbTe is distinct, with $\kappa_{\rm 3ph}^{\rm HA}$,
$\kappa_{\rm 3ph}^{\rm SCPH}$, and $\kappa_{\rm 3,4ph}^{\rm SCPH}$ exhibiting
a linear increase as the volume decreases, albeit with slightly different
slopes. Notably, for any given volume, $\kappa_{\rm 3ph}^{\rm SCPH}$
consistently exceeds $\kappa_{\rm 3ph}^{\rm HA}$. These differences stem from
the disparity in the impact of the SCPH effect on rocksalt and zincblende
structures.

Figure 8: (Color online). (a) and (b) are the PBEsol calculated dependence of
$\kappa_{\rm L}$ on the volume for ZnTe and PbTe calculated at three
perturbation orders (HA+3ph, SCPH+3ph, and SCPH+3,4ph), respectively. The
vertical dash line indicates the PBEsol optimized volume.

##  IV Conclusions

In this study, we systematically investigate the performance of eight commonly
used and newly proposed XC functionals in the solid community on lattice
thermal conductivity based on the phonon gas model with three different
perturbation orders. Within the approximation of harmonic phonons and three-
phonon scattering, all the XC functionals tend to overestimate $\kappa_{\rm
L}$ for almost all the compounds. Among these functionals, revTPSS exhibits
the smallest MRAE (29.75 %), followed by SCAN (30.11 %), LDA (34.00 %), and
r2SCAN (35.23 %), whereas PBE (45.51 %) and rSCAN (49.18 %) possess the
largest MREA among these XC functionals. If the phonon spectrum is
renormalized at room temperature using the SCPH method and only include three-
phonon scattering (SCPH+3ph), the $\kappa_{\rm L}$ values are further enlarged
for nearly all the functionals, due to the phonon frequency shift towards
higher energy, resulting in minimal MRAE of 30.57 % by SCAN and maximal MRAE
of 73.11 % by r2SCAN. However, if the four-phonon scattering is included in
the SCPH+3ph as well (SCPH+3,4ph), which is the most accurate method employed
in this work, the MRAEs of all the functionals are significantly reduced,
achieving the minimal MRAE of 20.03 % by PBEsol and maximal MRAE of 35.06 % by
rSCAN. Both rSCAN and r2SCAN have improved numerical stability, but only
r2SCAN has comparable MARE with PBEsol, implying the importance of adhering to
the constraints imposed by SCAN. However, although metaGGA functionals have a
higher hierarchy and are computationally more expensive than LDA and GGAs,
their accuracy in capturing anharmonicity is slightly inferior to the lower
hierarchy functionals, particularly PBEsol. Due to the small difference in
phonon spectra and the smaller MRAE in sound velocity compared to the
Grüneisen parameter, the primary source of error across these functionals
predominantly arises from the difference in dealing with anharmonicity. Our
findings not only offer a comprehensive assessment of XC functionals for
computing $\kappa_{\rm L}$, but also provide a useful guide for selecting
appropriate XC functionals when dealing with anharmonicity. Furthermore, our
results highlight the challenges of accurately describing anharmonicity using
XC functionals.

##  V ACKNOWLEDGMENTS

J.W. Z.X. and J.H. acknowledge the support of the National Science Fundation
of China (12374024). Y. X. acknowledges support from the US National Science
Foundation through award 2317008. The computing resource was supported by USTB
MatCom of Beijing Advanced Innovation Center for Materials Genome Engineering.

## References

  * Moore and Shi [2014] A. L. Moore and L. Shi, Emerging challenges and materials for thermal management of electronics, Mater. Today 17, 163 (2014).
  * Zhang  _et al._ [2021] Z. Zhang, X. Wang, and Y. Yan, A review of the state-of-the-art in electronic cooling, e-Prime 1, 100009 (2021).
  * Rowe [1995] D. E. Rowe,  _CRC Handbook of Thermoelectrics (1st ed.)_ (CRC Press, 1995).
  * Bell [2008] L. E. Bell, Cooling, heating, generating power, and recovering waste heat with thermoelectric systems, Science 321, 1457 (2008).
  * Mondal  _et al._ [2021] K. Mondal, L. I. Nuñez, C. M. Downey, and I. J. van Rooyen, Thermal barrier coatings overview: Design, manufacturing, and applications in high-temperature industries, Ind. Eng. Chem. Res. 60, 6061 (2021).
  * Padture  _et al._ [2002] N. P. Padture, M. Gell, and E. H. Jordan, Thermal barrier coatings for gas-turbine engine applications, Science 296, 280 (2002).
  * Clarke and Phillpot [2005] D. R. Clarke and S. R. Phillpot, Thermal barrier coating materials, Mater. Today 8, 22 (2005).
  * Lindsay  _et al._ [2013] L. Lindsay, D. A. Broido, and T. L. Reinecke, First-principles determination of ultrahigh thermal conductivity of boron arsenide: A competitor for diamond?, Phys. Rev. Lett. 111, 025901 (2013).
  * Yang  _et al._ [2019] X. Yang, T. Feng, J. Li, and X. Ruan, Stronger role of four-phonon scattering than three-phonon scattering in thermal conductivity of iii-v semiconductors at room temperature, Phys. Rev. B 100, 245203 (2019).
  * Xia  _et al._ [2020a] Y. Xia, V. I. Hegde, K. Pal, X. Hua, D. Gaines, S. Patel, J. He, M. Aykol, and C. Wolverton, High-throughput study of lattice thermal conductivity in binary rocksalt and zinc blende compounds including higher-order anharmonicity, Phys. Rev. X 10, 041029 (2020a).
  * Green [1954] M. S. Green, Markoff Random Processes and the Statistical Mechanics of Time-Dependent Phenomena. II. Irreversible Processes in Fluids, J. Chem. Phys. 22, 398 (1954).
  * Evans and Morriss [2008] D. J. Evans and G. Morriss,  _Statistical Mechanics of Nonequilibrium Liquids_, 2nd ed. (Cambridge University Press, 2008).
  * Klemens [1958] P. G. Klemens, Thermal conductivity and lattice vibrational modes, J. Phys. C: Solid State Phys. 7, 1 (1958).
  * Puligheddu  _et al._ [2019] M. Puligheddu, Y. Xia, M. Chan, and G. Galli, Computational prediction of lattice thermal conductivity: A comparison of molecular dynamics and boltzmann transport approaches, Phys. Rev. Mater. 3, 085401 (2019).
  * Feng  _et al._ [2017] T. Feng, L. Lindsay, and X. Ruan, Four-phonon scattering significantly reduces intrinsic thermal conductivity of solids, Phys. Rev. B 96, 161201 (2017).
  * Han  _et al._ [2022] Z. Han, X. Yang, W. Li, T. Feng, and X. Ruan, Fourphonon: An extension module to shengbte for computing four-phonon scattering rates and thermal conductivity, Comput. Phys. Commun. 270, 108179 (2022).
  * Xia [2018] Y. Xia, Revisiting lattice thermal transport in PbTe: The crucial role of quartic anharmonicity, Appl. Phys. Lett. 113, 073901 (2018).
  * Perdew and Schmidt [2001] J. P. Perdew and K. Schmidt, Jacob’s ladder of density functional approximations for the exchange-correlation energy, AIP Conf. Proc. 577, 1 (2001).
  * Zhou  _et al._ [2023] H. Zhou, S. Zhou, Z. Hua, K. Bawane, and T. Feng, Extreme sensitivity of higher-order interatomic force constants and thermal conductivity to the energy surface roughness of exchange-correlation functionals, Appl. Phys. Lett. 123, 192201 (2023).
  * Mardirossian and Head-Gordon [2017] N. Mardirossian and M. Head-Gordon, Thirty years of density functional theory in computational chemistry: an overview and extensive assessment of 200 density functionals, Mol. Phys. 115, 2315 (2017).
  * Kohn and Sham [1965] W. Kohn and L. J. Sham, Self-consistent equations including exchange and correlation effects, Phys. Rev. 140, A1133 (1965).
  * Perdew  _et al._ [1996] J. P. Perdew, K. Burke, and M. Ernzerhof, Generalized gradient approximation made simple, Phys. Rev. Lett. 77, 3865 (1996).
  * Perdew  _et al._ [2008] J. P. Perdew, A. Ruzsinszky, G. I. Csonka, O. A. Vydrov, G. E. Scuseria, L. A. Constantin, X. Zhou, and K. Burke, Restoring the density-gradient expansion for exchange in solids and surfaces, Phys. Rev. Lett. 100, 136406 (2008).
  * Perdew  _et al._ [2009] J. P. Perdew, A. Ruzsinszky, G. I. Csonka, L. A. Constantin, and J. Sun, Workhorse semilocal density functional for condensed matter physics and quantum chemistry, Phys. Rev. Lett. 103, 026403 (2009).
  * Sun  _et al._ [2015] J. Sun, A. Ruzsinszky, and J. P. Perdew, Strongly constrained and appropriately normed semilocal density functional, Phys. Rev. Lett. 115, 036402 (2015).
  * Sun  _et al._ [2016a] J. Sun, R. C. Remsing, Y. Zhang, Z. Sun, A. Ruzsinszky, H. Peng, Z. Yang, A. Paul, U. Waghmare, X. Wu, _et al._ , Accurate first-principles structures and energies of diversely bonded systems from an efficient density functional, Nat. Chem. 8, 831 (2016a).
  * Kitchaev  _et al._ [2016] D. A. Kitchaev, H. Peng, Y. Liu, J. Sun, J. P. Perdew, and G. Ceder, Energetics of ${\mathrm{mno}}_{2}$ polymorphs in density functional theory, Phys. Rev. B 93, 045132 (2016).
  * Zhang  _et al._ [2017] Y. Zhang, J. Sun, J. P. Perdew, and X. Wu, Comparative first-principles studies of prototypical ferroelectric materials by lda, gga, and scan meta-gga, Phys. Rev. B 96, 035143 (2017).
  * Bartók and Yates [2019] A. P. Bartók and J. R. Yates, Regularized scan functional, J. Chem. Phys. 150, 161101 (2019).
  * Furness  _et al._ [2020] J. W. Furness, A. D. Kaplan, J. Ning, J. P. Perdew, and J. Sun, Accurate and numerically efficient r2scan meta-generalized gradient approximation, J. Phys. Chem. Lett. 11, 8208 (2020).
  * Shao  _et al._ [2023] X. Shao, P. Liu, C. Franchini, Y. Xia, and J. He, Assessing the performance of exchange-correlation functionals on lattice constants of binary solids at room temperature within the quasiharmonic approximation, Phys. Rev. B 108, 024306 (2023).
  * Jain and McGaughey [2015] A. Jain and A. J. H. McGaughey, Effect of exchange–correlation on first-principles-driven lattice thermal conductivity predictions of crystalline silicon, Comput. Mater. Sci. 110, 115 (2015).
  * Perdew  _et al._ [1992] J. P. Perdew, J. A. Chevary, S. H. Vosko, K. A. Jackson, M. R. Pederson, D. J. Singh, and C. Fiolhais, Atoms, molecules, solids, and surfaces: Applications of the generalized gradient approximation for exchange and correlation, Phys. Rev. B 46, 6671 (1992).
  * Becke [1988] A. D. Becke, Density-functional exchange-energy approximation with correct asymptotic behavior, Phys. Rev. A 38, 3098 (1988).
  * Lee  _et al._ [1988] C. Lee, W. Yang, and R. G. Parr, Development of the colle-salvetti correlation-energy formula into a functional of the electron density, Phys. Rev. B 37, 785 (1988).
  * Qin  _et al._ [2018] G. Qin, Z. Qin, H. Wang, and M. Hu, On the diversity in the thermal transport properties of graphene: A first-principles-benchmark study testing different exchange-correlation functionals, Comput. Mater. Sci. 151, 153 (2018).
  * Taheri  _et al._ [2018] A. Taheri, C. Da Silva, and C. H. Amon, First-principles phonon thermal transport in graphene: Effects of exchange-correlation and type of pseudopotential, J. Appl. Phys. 123 (2018).
  * Han and Ruan [2023] Z. Han and X. Ruan, Thermal conductivity of monolayer graphene: Convergent and lower than diamond, Phys. Rev. B 108, L121412 (2023).
  * Dongre  _et al._ [2022] B. Dongre, J. Carrete, N. Mingo, and G. K. H. Madsen, Thermal conductivity of group-iii phosphides: The special case of gap Phys. Rev. B 106, 205202 (2022).
  * Xia  _et al._ [2020b] Y. Xia, K. Pal, J. He, V. Ozoliņš, and C. Wolverton, Particlelike phonon propagation dominates ultralow lattice thermal conductivity in crystalline tl3vse4, Phys. Rev. Lett. 124, 065901 (2020b).
  * Jain [2020] A. Jain, Multichannel thermal transport in crystalline tl 3 vse 4, Phys. Rev. B 102, 201201 (2020).
  * Klimeš  _et al._ [2011] J. c. v. Klimeš, D. R. Bowler, and A. Michaelides, Van der waals density functionals applied to solids, Phys. Rev. B 83, 195131 (2011).
  * Blöchl [1994] P. E. Blöchl, Projector augmented-wave method, Phys. Rev. B 50, 17953 (1994).
  * Kresse and Joubert [1999] G. Kresse and D. Joubert, From ultrasoft pseudopotentials to the projector augmented-wave method, Phys. Rev. B 59, 1758 (1999).
  * Kresse and Furthmüller [1996] G. Kresse and J. Furthmüller, Efficient iterative schemes for ab-initio total-energy calculations using a plane-wave basis set, Phys. Rev. B 54, 11169 (1996).
  * Kresse and Furthmüller [1996] G. Kresse and J. Furthmüller, Efficiency of ab-initio total energy calculations for metals and semiconductors using a plane-wave basis set, Comput. Mater. Sci. 6, 15 (1996).
  * Sun  _et al._ [2016b] J. Sun, R. C. Remsing, Y. Zhang, Z. Sun, A. Ruzsinszky, H. Peng, Z. Yang, A. Paul, U. Waghmare, X. Wu, M. L. Klein, and J. P. Perdew, Accurate first-principles structures and energies of diversely bonded systems from an efficient density functional, Nat Chem 8, 831 (2016b).
  * Fu and Singh [2018] Y. Fu and D. J. Singh, Applicability of the strongly constrained and appropriately normed density functional to transition-metal magnetism, Phys Rev Lett 121, 207201 (2018).
  * Mejia-Rodriguez and Trickey [2019] D. Mejia-Rodriguez and S. B. Trickey, Comment on "regularized scan functional" [j. chem. phys. 150, 161101 (2019)], J Chem Phys 151, 207101 (2019).
  * Togo and Tanaka [2015] A. Togo and I. Tanaka, First principles phonon calculations in materials science, Scr. Mater. 108, 1 (2015).
  * Gonze  _et al._ [1994] X. Gonze, J.-C. Charlier, D. Allan, and M. Teter, Interatomic force constants from first principles: The case of $\alpha$-quartz, Phys. Rev. B 50, 13035 (1994).
  * Gonze and Lee [1997] X. Gonze and C. Lee, Dynamical matrices, born effective charges, dielectric permittivity tensors, and interatomic force constants from density-functional perturbation theory, Phys. Rev. B 55, 10355 (1997).
  * Zhou  _et al._ [2014] F. Zhou, W. Nielson, Y. Xia, and V. Ozolins, Lattice anharmonicity and thermal conductivity from compressive sensing of first-principles calculations, Phys Rev Lett 113, 185501 (2014).
  * Hooton [1955] D. Hooton, Li. a new treatment of anharmonicity in lattice thermodynamics: I, Philos. Mag. J. Sci. 46, 422 (1955).
  * Pawley  _et al._ [1966] G. S. Pawley, W. Cochran, R. A. Cowley, and G. Dolling, Diatomic ferroelectrics, Phys. Rev. Lett. 17, 753 (1966).
  * Werthamer [1970] N. R. Werthamer, Self-consistent phonon formulation of anharmonic lattice dynamics, Phys. Rev. B 1, 572 (1970).
  * Guo  _et al._ [2024] Z. Guo, Z. Han, D. Feng, G. Lin, and X. Ruan, Sampling-accelerated prediction of phonon scattering rates for converged thermal conductivity and radiative properties, npj Comput. Mater. 10, 31 (2024).
  * Ong  _et al._ [2013] S. P. Ong, W. D. Richards, A. Jain, G. Hautier, M. Kocher, S. Cholia, D. Gunter, V. L. Chevrier, K. A. Persson, and G. Ceder, Python materials genomics (pymatgen): A robust, open-source python library for materials analysis, Comput. Mater. Sci. 68, 314 (2013).
  * Hill [1952] R. Hill, The elastic behaviour of a crystalline aggregate, Proc. Phys. Soc. A 65, 349 (1952).
  * Tritt [2010] T. M. Tritt,  _Thermal Conductivity: Theory, Properties, and Applications_ (2010).
  * Koehler [1966] T. R. Koehler, Theory of the self-consistent harmonic approximation with application to solid neon, Phys. Rev. Lett. 17, 89 (1966).
  * Xia  _et al._ [2020c] Y. Xia, V. Ozoliņš, and C. Wolverton, Microscopic mechanisms of glasslike lattice thermal transport in cubic ${\mathrm{cu}}_{12}{\mathrm{sb}}_{4}{\mathrm{s}}_{13}$ tetrahedrites, Phys. Rev. Lett. 125, 085901 (2020c).
  * Hanus  _et al._ [2019] R. Hanus, M. T. Agne, A. J. E. Rettie, Z. Chen, G. Tan, D. Y. Chung, M. G. Kanatzidis, Y. Pei, P. W. Voorhees, and G. J. Snyder, Lattice softening significantly reduces thermal conductivity and leads to high thermoelectric efficiency, Adv. Mater. 31, 1900108 (2019).
  * Chen  _et al._ [2018] Z. Chen, X. Zhang, S. Lin, L. Chen, and Y. Pei, Rationalizing phonon dispersion for lattice thermal conductivity of solids, Natl. Sci. Rev. 5, 888 (2018).
  * Grüneisen [1908] E. Grüneisen, Über die thermische ausdehnung und die spezifische wärme der metalle, Ann. Phys. 331, 211 (1908).
  * Grüneisen [1912] E. Grüneisen, Theorie des festen zustandes einatomiger elemente, Ann. Phys. 344, 257 (1912).
  * Morelli and Slack [2006] D. T. Morelli and G. A. Slack, High lattice thermal conductivity solids, in  _High thermal conductivity materials_ (Springer, 2006) pp. 37–68.
  * Ritz  _et al._ [2019] E. T. Ritz, S. J. Li, and N. A. Benedek, Thermal expansion in insulating solids from first principles, J. Appl. Phys. 126 (2019).
  * Li  _et al._ [2014] W. Li, J. Carrete, N. A. Katcho, and N. Mingo, ShengBTE: a solver of the Boltzmann transport equation for phonons, Comp. Phys. Commun. 185, 1747–1758 (2014).
  * Shindé and Goela [2006] S. L. Shindé and J. Goela,  _High thermal conductivity materials_, Vol. 91 (Springer, 2006).
  * Slack [1979] G. A. Slack, The thermal conductivity of nonmetallic crystals, SOLID STATE PHYS 34, 1 (1979).
  * Miller  _et al._ [2017] S. A. Miller, P. Gorai, B. R. Ortiz, A. Goyal, D. Gao, S. A. Barnett, T. O. Mason, G. J. Snyder, Q. Lv, V. Stevanović, and E. S. Toberer, Capturing anharmonicity in a lattice thermal conductivity model for high-throughput predictions, Chem. Mater. 29, 2494 (2017).
  * Peierls [1996] R. E. Peierls,  _Quantum theory of solids_ (Clarendon Press, 1996).
  * Xie  _et al._ [2020] L. Xie, J. H. Feng, R. Li, and J. Q. He, First-principles study of anharmonic lattice dynamics in low thermal conductivity ${\mathrm{agcrse}}_{2}$: Evidence for a large resonant four-phonon scattering, Phys. Rev. Lett. 125, 245901 (2020).
  * Feng and Ruan [2016] T. Feng and X. Ruan, Quantum mechanical prediction of four-phonon scattering rates and reduced thermal conductivity of solids, Phys. Rev. B 93, 045202 (2016).
  * Zeier  _et al._ [2016] W. G. Zeier, A. Zevalkink, Z. M. Gibbs, G. Hautier, M. G. Kanatzidis, and G. J. Snyder, Thinking like a chemist: Intuition in thermoelectric materials, Angew. Chem. Int. Ed. 55, 6826 (2016).
  * Simoncelli  _et al._ [2019] M. Simoncelli, N. Marzari, and F. Mauri, Unified theory of thermal transport in crystals and glasses, Nat. Phys. 15, 809 (2019).
  * El-Sharkawy  _et al._ [1983] A. El-Sharkawy, A. Abou El-Azm, M. Kenawy, A. Hillal, and H. Abu-Basha, Thermophysical properties of polycrystalline pbs, pbse, and pbte in the temperature range 300–700 k, Int. J. Thermophys. 4, 261 (1983).
  * Akhmedova and Abdinov [2009] G. Akhmedova and D. S. Abdinov, Effect of thallium doping on the thermal conductivity of pbte single crystals, Inorg. Mater. 45, 854 (2009).
  * Håkansson and Ross [1985] B. Håkansson and R. Ross, Thermal conductivity of solid naf under high pressure, Int. J. Thermophys. 6, 353 (1985).
  * Håkansson and Andersson [1986] B. Håkansson and P. Andersson, Thermal conductivity and heat capacity of solid nacl and nai under pressure, J Phys Chem Solids 47, 355 (1986).
  * Steigmeier and Kudman [1963] E. Steigmeier and I. Kudman, Thermal conductivity of iii-v compounds at high temperatures, Phys. Rev. 132, 508 (1963).
  * Madelung [2004] O. Madelung,  _Semiconductors: data handbook_ (Springer Science & Business Media, 2004).
  * Kang  _et al._ [2018] J. S. Kang, M. Li, H. Wu, H. Nguyen, and Y. Hu, Experimental observation of high thermal conductivity in boron arsenide, Science 361, 575 (2018).
  * Kumashiro  _et al._ [1989] Y. Kumashiro, T. Mitsuhashi, S. Okaya, F. Muta, T. Koshiro, Y. Takahashi, and M. Mirabayashi, Thermal conductivity of a boron phosphide single-crystal wafer up to high temperature, J. Appl. Phys. 65, 2147 (1989).
  * Slack [1972] G. A. Slack, Thermal conductivity of ii-vi compounds and phonon scattering by fe 2+ impurities, Phys. Rev. B 6, 3791 (1972).
  * Chen  _et al._ [2020] K. Chen, B. Song, N. K. Ravichandran, Q. Zheng, X. Chen, H. Lee, H. Sun, S. Li, G. A. G. U. Gamage, F. Tian, Z. Ding, Q. Song, A. Rai, H. Wu, P. Koirala, A. J. Schmidt, K. Watanabe, B. Lv, Z. Ren, L. Shi, D. G. Cahill, T. Taniguchi, D. Broido, and G. Chen, Ultrahigh thermal conductivity in isotope-enriched cubic boron nitride, Science 367, 555 (2020).
  * Li and Mingo [2014] W. Li and N. Mingo, Thermal conductivity of fully filled skutterudites: Role of the filler, Phys. Rev. B 89, 184304 (2014).
  * He  _et al._ [2022] J. He, Y. Xia, W. Lin, K. Pal, Y. Zhu, M. G. Kanatzidis, and C. Wolverton, Accelerated discovery and design of ultralow lattice thermal conductivity materials using chemical bonding principles, Adv. Funct. Mater. 32, 2108532 (2022).