Anharmonicity Measure for Materials

Theoretical frameworks used to qualitatively and quantitatively describe nuclear dynamics in solids are often based on the harmonic approximation. However, this approximation is known to become inaccurate or to break down completely in many modern functional materials. Interestingly, there is no reliable measure to quantify anharmonicity so far. Thus, a systematic classification of materials in terms of anharmonicity and a benchmark of methodologies that may be appropriate for different strengths of anharmonicity is currently impossible. In this work, we derive and discuss a statistical measure that reliably classifies compounds across temperature regimes and material classes by their"degree of anharmonicity". This enables us to distinguish"harmonic"materials, for which anharmonic effects constitute a small perturbation on top of the harmonic approximation, from strongly"anharmonic"materials, for which anharmonic effects become significant or even dominant and the treatment of anharmonicity in terms of perturbation theory is more than questionable. We show that the analysis of this measure in real and reciprocal space is able to shed light on the underlying microscopic mechanisms, even at conditions close to, e.g., phase transitions or defect formation. Eventually, we demonstrate that the developed approach is computationally efficient and enables rapid high-throughput searches by scanning over a set of several hundred binary solids. The results show that strong anharmonic effects beyond the perturbative limit are not only active in complex materials or close to phase transitions, but already at moderate temperatures in simple binary compounds.

Theoretical frameworks used to qualitatively and quantitatively describe nuclear dynamics in solids are often based on the harmonic approximation. However, this approximation is known to become inaccurate or to break down completely in many modern functional materials. Interestingly, there is no reliable measure to quantify anharmonicity so far. Thus, a systematic classification of materials in terms of anharmonicity and a benchmark of methodologies that may be appropriate for different strengths of anharmonicity is currently impossible. In this work, we derive and discuss a statistical measure that reliably classifies compounds across temperature regimes and material classes by their "degree of anharmonicity". This enables us to distinguish "harmonic" materials, for which anharmonic effects constitute a small perturbation on top of the harmonic approximation, from strongly "anharmonic" materials, for which anharmonic effects become significant or even dominant and the treatment of anharmonicity in terms of perturbation theory is more than questionable. We show that the analysis of this measure in real and reciprocal space is able to shed light on the underlying microscopic mechanisms, even at conditions close to, e.g., phase transitions or defect formation. Eventually, we demonstrate that the developed approach is computationally efficient and enables rapid high-throughput searches by scanning over a set of several hundred binary solids. The results show that strong anharmonic effects beyond the perturbative limit are not only active in complex materials or close to phase transitions, but already at moderate temperatures in simple binary compounds.

I. INTRODUCTION
In condensed-matter physics and materials science, the dynamics of nuclei plays a decisive role for many materials properties. At the lowest level of theory, these vibrations are commonly described in the harmonic approximation, in which the potential-energy surface V(R) is approximated by V (2) (R), a second-order Taylor expansion in the atomic displacements, {∆R I }, about a minimum energy configuration, {R 0 I }, in terms of force constants, Φ I,J α,β [1]. The dynamical properties are determined by the model Hamiltonian for the nuclear momenta, P I , and positions, R I = R 0 I + ∆R I . Since this Hamiltonian separates into 3N uncoupled harmonic oscillators, called phonons, both the classical equations of motion and the quantum-mechanical Schrödinger equation can be solved analytically. This allows for the computation of thermodynamic properties such as the harmonic free energy and the heat capacity [1,2]. Since computing Φ I,J α,β from first principles is a straightforward and computationally affordable task [3][4][5], the harmonic approximation is a popular tool in materials science.
Material properties and phenomena that are described inaccurately or not at all within such a harmonic model are generally referred to as anharmonic effects. These effects include i) the temperature dependence of equilibrium properties like thermal lattice expansion, ii) thermal shift of vibrational frequencies and linewidth broadening, iii) phase transitions, and iv) heat transport. All of these properties either diverge or vanish in the harmonic approximation [2,[6][7][8].
In recent years, significant advances in modeling anharmonic effects were achieved, especially in the field of thermal transport. To date, anharmonic effects can be addressed either exactly via non-perturbative ab initio Molecular Dynamics (aiMD) [9,10], or, more commonly, with approximate, perturbative models [11][12][13][14][15][16][17][18]. In the latter case, the Taylor expansion of the potential V(R) is extended beyond the second order and the additional terms are treated as a perturbation of the harmonic model. Yet, there is a fundamental question about the applicability of perturbative approaches: The reliability of any perturbation expansion is controlled by the strength of the perturbation, which has to be "small" with respect to the reference [19]. However, no measure exists to date to reliably quantify the strength of anharmonic effects. First, this lack prevents the systematical and quantitative exploration of the applicability limits of perturbative techniques. Second, this hinders a systematic classification of materials across temperature regimes by anharmonicity. Given that anharmonic effects may influence or even determine macroscopic material properties, understanding qualitative trends that govern anharmonicity across material space is a challenge of growing importance.
In this work, we address this open issue by deriving and validating the required anharmonicity measure. As discussed below, it (a) allows for a systematic and quantitative classification of compounds across material space from systems with only mild anharmonic contributions, to strongly anharmonic systems where the phonon picture is invalid, (b) establishes a link between the actuating microscopic mechanisms and macroscopic properties, and (c) requires a fraction of the computational cost of either perturbative or aiMD calculations and thus paves the way for high-throughput anharmonicity classification of materials. After presenting and discussing the underlying theory, definitions, and the numerical aspects for two exemplary materials in Sec. II-V, we show how our anharmonicity measure can be incorporated into an ab initio high-throughput screening of material space in Sec. VI. We discuss our findings and their implications for theoretical solid-state physics and materials science in Sec. VII.

II. DEFINITION OF ANHARMONICITY
In the Born-Oppenheimer approximation [20], the dynamics of a nuclear many-body system is described by the Hamiltonian determined by the potential energy V(R), where R = {R 1 , . . . , R N } denotes the atomic coordinates and P = {P 1 , . . . , P N } the respective momenta.
Treating the nuclei as classical particles, the Hamiltonian generates equations of motion for each atom I, with atomic mass M I , accelerationR I , and the force F I . As mentioned in the introduction, the harmonic approximation replaces the many-body potential V(R) by a second-order Taylor expansion around equilibrium R 0 , in terms of the force constants Accordingly, the forces are approximated as On this approximated potential-energy surface, the equations of motion can be solved analytically [2]. Differences between the actual and the harmonic nuclear dynamics arise from the difference between the exact potential V(R) and the harmonic potential V (2) (R), i. e., from anharmonic contributions. Consequently, we define the anharmonic contribution to the potential as Analogously, the anharmonic contributions to the forces that enter the equations of motion are given by Left: Sketch of a one-dimensional potentialenergy surface V(R) (solid black), its harmonic approximation V (2) (R) (dashed blue), and the anharmonic contribution V A (R) (solid red). Right: The force F (R) given by the derivative of the potential energy V(R) (black), the force F (2) stemming from the harmonic potential V (2) (R) (blue), and the anharmonic contribution with F (2) I obtained from Eq. (6). The construction is qualitatively depicted in Fig. 1. Compared to the potential V(R), working with the forces has the advantage that they naturally decompose into 3N components, hence giving straightforward access to a real-space per-atom and a reciprocal-space per-mode analysis. For the latter purpose, we construct the dynamical matrix Here,Ĩ labels the atoms in the primitive unit cell, L denotes the Bravais lattice points, J(J, L) labels the periodic images ofJ, and q is the phonon wave vector [1]. Diagonalizing D(q) yields the harmonic eigenfrequencies ω n (q) and eigenmodes e n (q). For readability, we use the generalized index s = (q, n) in the following. In this notation, the mode-resolved forces are

III. QUANTIFYING ANHARMONICITY
Two prototypical materials are used as examples in the following to elucidate the concepts developed in this work: Silicon (fcc-diamond) serves as example for a largely harmonic material, whereas the low-temperature structure of the perovskite KCaF 3 (Pnma [22], cf. Fig. 2), is used as example for a complex, strongly anharmonic material. As shown in Fig. 3 and 4, both materials exhibit vibrational spectra of roughly the same frequency range in the harmonic approximation. The perovskite KCaF 3 features an anion octahedral tilting (a − b + a − in the Glazer notation [25,26]) which reduces with temperature, finally leading to a dynamically stabilized cubic crystal at 560 K [26,27]. We will discuss this phase transition and the implications for anharmonicity quantification in more detail in Sec. IV. We investigate both compounds at room temperature via ab initio molecular dynamics (aiMD) simulations [28] at the GGA level of theory using the PBEsol exchangecorrelation functional [29], light default basis sets [30], and a Langevin thermostat [31] to perform canonical sampling in supercells of 216 atoms (Si) and 160 atoms (KCaF 3 ), respectively. The aiMD is performed with a time step of 5 fs for an initial thermalization period of 2 ps and a sampling period of 8 ps. The chosen numerical settings ensure that all quantities of interest, i.e., structural parameters such as lattice constants, dynamical properties such as vibrational frequencies, as well as thermodynamic averages such as the σ A measure introduced below are converged within ±1 %.
All the necessary calculations and tools used for performing the calculations and investigating the results are implemented in our python package FHI-vibes [32]. It builds on top of the Atomistic Simulation Environment (ASE) [33], interfaces with phonopy for building harmonic force constants [4], and integrates tightly with the all-electron, numeric atomic orbitals code FHI-aims for performing ab initio calculations of energy, forces, and stress [30,34].

A. Normalization of Forces
A prerequisite for comparing forces acting in different systems under various thermodynamic conditions is that these forces are normalized. For this purpose, we characterize each force component F I,α (t) observed during the simulation by the probability-density function p I,α (F ), and use the definition of the thermodynamic average to obtain where δ(F ) denotes the delta distribution. To characterize the whole system, we use the mixture probability distribution i. e., the weighted sum of probability distributions for each force component F I,α . Since the average of the individual force components vanishes in the absence of an external force, F I,α = 0, the distribution p(F ) has a mean of zero. However, the width of p(F ) depends on the actual material, as shown in plots of p(F ) for silicon and KCaF 3 in Fig. 5 a) and b). We evaluate the width of the force distribution by computing its standard deviation, with the thermodynamic average Distribution of forces p(F ) observed during aiMD simulations at 300 K before (left) and after normalization (right) for silicon (upper) and KCaF3 (lower row). Dashed lines denote the standard deviation of the distribution.
The analysis reveals that the distribution of forces in silicon exhibts a width of σ[F Si ] = 0.60 eV/Å, while KCaF 3 features a width of σ[F KCaF ] = 0.38 eV/Å. This is consistent with the phonon dispersions in Fig. 3 and 4, since KCaF 3 features more low-energy states, resulting, on average, in smaller restoring forces compared to silicon. We take σ [F ] as a measure for the average magnitude of forces acting in a material in thermodynamic equilibrium, including harmonic and anharmonic contributions. The average force σ[F ] therefore defines a scale in which the forces can be given and compared independent of the system by defining the normalized force As shown in Fig. 5 c) and d), the mixture distributions of these normalized forces, p (F/σ[F ]), exhibit the same unit width for both Si and KCaF 3 . This normalization thus allows to perform a meaningful comparison between the two materials.

B. Anharmonicity Measure
To compute the anharmonic contribution to the forces, we use the aiMD forces F I (t) and obtain their harmonic contribution F  I (t), as defined in Eq. (8). In close analogy the previous section, we use the probability distribution p I,α (F A ) and the mixture probability distribution p(F A ) to characterize the statistical behavior of F A . Likewise, we normalize F A I,α with respect to the force scale σ[F ]. Accordingly, F A I (t) / σ [F ] describes the anharmonic contribution to the force on atom I with respect to the average magnitude of the total forces F I . For both silicon and KCaF 3 at 300 K, the distribu- Fig. 6 as joint probability plots. Given that the force scale σ [F ] introduced above is used, the normalized forces are similarly distributed in x-direction for both materials as discussed earlier. On the y-axis, however, the distribution of the rescaled anharmonic forces F A /σ[F ] is significantly different for the two materials. In silicon, the distribution is sharply peaked around 0 with a width of 0.15 σ [F Si ]. From the distribution, we can quantify that only 15 % of the forces stem from anharmonic contributions on average and the probability of finding anharmonic force contributions of 0.5 σ [F Si ] or larger is < 0.01 %. This confirms the general understanding that silicon is largely harmonic and strong anharmonic contributions are essentially absent. Conversely, the distribution of anharmonic forces for the perovskite KCaF 3 in Fig. 6 (right) is much broader, featuring a width of 0.36 σ [F KCaF ], i. e., 36 % of the forces acting on the nuclei stem from anharmonic effects on average. More importantly, finding strongly anharmonic force contributions of 0.5 σ [F KCaF ] or larger is 16.5 %, thus several orders of magnitude more probable than in silicon. This means that these contributions are indeed significant in this compound, as argued above in the introduction of Sec. III.
In spirit of this discussion, we define the following measure for the quantitative estimation of the degree of anharmonicity in a material: with the thermodynamic ensemble average · T obtained according to Eq. 15. σ A (T ) measures the standard deviation of the distribution of anharmonic force components F A I,α obtained from the ab initio forces F and their harmonic approximation F (2) according to Eq. (8), normalized by the standard deviation of the ab initio force distribution in the absence of external forces. This is mathematically equivalent to the root mean square error (RMSE) of the harmonic model divided by the standard deviation of the force distribution.

C. Atom-and Mode-Resolved Anharmonicity
To estimate the degree of anharmonicity in a specific subset of the available degrees of freedom, denoted by X, we evaluate where X can be, e. g., a specific atom I, a group of atoms, or a vibrational mode s. It is important to note that in Eq. (18), we normalize by the width of the force components of interest, σ [F X ], and not by the force scale σ [F ] = X σ[F X ] 2 that averages over all available components, as in Eq. (17). By this means, we assess the relative importance of anharmonicity in the degree(s) of freedom of interest. As an example, Fig. 7 shows the species-resolved anharmonicities for KCaF 3 . This analysis reveals that the forces acting on the K atoms have the largest anharmonic contribution, while the Ca and F atoms are more harmonic. This is a result of the K atoms moving in a relatively shallow potential that allows them to participate significantly in the octahedral tilt, whereas the calcium atoms occupy the relatively stable vertices of the unit cell and its center, cf. Fig. 2.
To estimate the importance of anharmonicity in a specific vibrational mode s, we evaluate Eq. (18) for the mode resolved forces F s obtained from the eigenvectors e s as given by Eq. (10). The mode-resolved degree of anharmonicity σ A s is plotted in Fig. 8 as a function of the mode frequency ω s for silicon and KCaF 3 at 300 K. Across the whole vibrational spectrum, silicon exhibits almost the same mild anharmonicity of σ A s 0.2. Conversely, KCaF 3 , exhibits larger values of σ A s with a significantly increasing magnitude and width for frequencies below 7 THz. Below frequencies of 5 THz, the anharmonic contributions make up for roughly 50 % of the forces, with several modes approaching or even exceeding σ A s = 1, as one of the Γ-point optical modes with ω s = 2.81 THz. This implies that the harmonic model predicts forces that are not even qualitatively correct. They may have a significantly incorrect value and even a wrong sign. In the joint probability density shown in the inset, the anharmonic distribution p(F A s ) is thus broader than the force distribution for this particular mode.
Qualitative insight into the microscopic mechanism underlying anharmonicity can be obtained by inspecting the displacements associated with the modes of interest. For example, the Γ-point mode with σ A s > 1 highlighted in the inset of Fig. 8 participates in the phase transition to the cubic structure above 550 K [26]. This shows that the system begins to "feel" the onset of the phase transition already at 300 K, well below the actual phase transition temperature. This aspect, i.e., the relation between σ A s and phase transition temperatures, is discussed further in more detail and for a broad set of systems in Secs. IV and VI B. Accordingly, the proposed measure is not only a valuable quantification and classifcation tool, but it is also sheds light on the microscopic mechanisms driving anharmonicity.

IV. APPLICATION TO DYNAMICALLY STABILIZED SYSTEMS
As mentioned before, KCaF 3 undergoes a second-order phase transition to the cubic aristotype structure above 560 K [26]. This structure, which corresponds to an alignment of the octahedra (see Fig. 2), is not a local minimum of the potential-energy surface, but a saddle point. This can be seen from the respective phonon dispersion, which features several imaginary modes, as shown in Fig. 9. We use this aristotype phase of KCaF 3 to exemplify the meaning of σ A for dynamically stabilized systems, since this is a commonly observed stabilization mechanism in strongly anharmonic materials [35][36][37], especially in perovskites [38,39]. Imaginary phonon frequencies as observed in Fig. 9 imply a breakdown of the harmonic approximation, since displacements from equilibrium are not energetically bounded. Accordingly, imaginary modes cannot be used to assess a physically meaningful dynamics and are thus typically neglected in standard perturbative methods. With respect to the potential-energy surface, however, a parabolic expansion around the cubic equilibrium and the evaluation of force constants is still possible. Accordingly, the proposed definition of the σ A -measure remains meaningful, as demonstrated in the following. For this purpose, we compare σ A cubic (T ), computed with the force constants of the cubic structure, to σ A Pnma (T ), computed using the force constants of the low-temperature orthorombic Pnma phase, as shown in Fig. 10. At low temperatures, at which the Pnma phase is stable, σ A Pnma starts from 0.28 and increases roughly linearly up to σ A Pnma (400 K) = 0.47. Between 400 and 600 K, a superlinear increase of σ A Pnma up to a value around 1 above 600 K is observed. Essentially, this means that forces obtained from the harmonic Pnma model are irrelevant for the nuclear dynamics at T > 600 K. This is in line with the finding that the material undergoes a phase transition to the cubic phase at these temperatures, as also observable in the MD simulations.
For the exact same reason, the σ A cubic values obtained using the harmonic model of the cubic structure are very high > 0.9 below the phase transition temperature, since the orthorombic Pnma structure and not the cubic one is thermodynamically stable at these conditions. The fact that σ A cubic and σ A Pnma cross each other at ≈ 500K (σ A = 0.79) indicates that the measure σ A is not only able to quantify the strong anharmonic effects active under these conditions, but that it is also able to qualitatively capture the underlying phase transition.

V. ACCELERATING ANHARMONICITY QUANTIFICATION
In the previous sections, accurate but computationally involved aiMD simulations were used to explore the potential-energy surface for performing anharmonicity quantification. For scanning through material space in a high-throughput fashion, it is desirable to obtain reliable estimates for σ A at a more moderated cost, i.e., by a fast computation in order to decide whether a full aiMD calculation is necessary to model the nuclear dynamics, or if the harmonic approximation (with or without further perturbative corrections) might suffice.
For this purpose, we note that the thermodynamic averages entering the anharmonicity metric given in Eq. (17) and (18) can be evaluated approximately by sampling with the harmonic Hamiltonian, i. e., ≈ O with β = 1/k B T . For the practical evaluation of Eq. (20), we generate atomic configurations via [40] ∆R in which e s are the harmonic eigenvectors, A s = √ 2k B T /ω s is the mean mode amplitude in the classical limit [2], and ζ s is a normally distributed random number [40]. To estimate the convergence of σ A with respect to the number of samples generated by Eq. (21), we evaluate σ A as defined in Eq. (17) for each individual sample R n and denote this value by σ A [R n ]. The values are plotted in Fig. 11 for a total of 30 samples. For silicon, each of the individual samples is sufficient to obtain an estimation of σ A within more than 99 % accuracy, given that the harmonic approximation in Eq. (21) holds in this case. In the case of KCaF 3 , where the harmonic approximation is not expected to yield a reliable dynamics, the described approach yields a value of σ A = 0.38 that differs from the one obtained by aiMD (σ A = 0.36) by 5 %, as shown in Fig. 11. Furthermore, we find that a good estimate of σ A = 0.39 ± 0.09 can be obtained with less than 10 samples, thus with a computational cost that is reduced by up to two orders of magnitude with respect to a full aiMD. Note that the estimated value of σ A obtained via sampling also allows to judge how reliable the sampling itself is, as reflected by the fact aiMD and harmonic sampling coincide for Si, but differ in the case of KCaF 3 .
Exploiting this rational allows to speed up this statistical approach even further by assigning ζ s = (−1) s−1 , which generates a single, deterministic sample from the most probable part of the random distribution [41]. In the classical limit explored in this work, this essentially implies evaluating σ A (T ) at the turning point of the oscillation estimated by the harmonic model. This allows one to reliably single out very harmonic materials (σ A ≤ 0.2) within a single force evaluation, thus saving a further order of magnitude in computational cost. Most importantly, one can exclude that highly-anharmonic materials are misclassified when σ A OS ≤ 0.2, given that the anharmonicity is that small even at the turning point. We compare the temperature dependence of σ A (T ) for this one-shot sampling technique and aiMD in Fig. 12. For silicon, the one-shot sampling is found to be sufficient to obtain σ A , as expected. Interestingly, also for KCaF 3 , the one-shot sampling is able to estimate σ A at low temperatures, but no longer yields quantitatively reliable estimates at elevated temperatures, at which the actual dynamical behavior of KCaF 3 deviates significantly from the harmonic reference as discussed in the previous section. Over the whole temperature range, the one-shot approach does however detect that strongly anharmonic effects are active and that aiMD simulations are necessary to reliably treat its dynamics. This is further substantiated in Fig. 13 for a wider set of materials discussed in more detail in the following section. Here, we compare the values of σ A at 300 K obtained by using the one-shot and molecular dynamics approaches with each other for 63 materials. As expected, σ A OS and σ A MD are in good agreement with each other for Si like materials with a σ A < 0.2, while deviations are observed for larger values of σ A OS , whereby the errors are more pronounced the larger σ A OS becomes. In particular, we observe that σ A OS can yield qualitatively wrong results for σ A OS > 0.4. For this reason, σ A OS is used in the following to pre-screen the materials and to single out materials with σ A < 0.2, whereas aiMD is used to obtain reliable values of σ A whenever σ A OS > 0.2.

VI. APPLICATION TO MATERIAL SPACE
To substantiate and generalize the insights obtained for the two example materials in the previous section, we compute σ A for two distinct groups of materials at multiple temperatures: simple binary compounds (rock salts, zincblende and wurtzites) and perovskites. For both cases, we perform symmetry-preserving geometry optimization for the structures using parametric constraints in FHI-aims until all forces are converged to a numerical precision better than 10 −3 eV/Å [42]. From there, we calculate a converged harmonic model of each material's vibrational properties and then generate thermally displaced supercells using either molecular dynamics or the one-shot approach according to Eq. (21). All calculations use the PBEsol functional to calculate the exchange-correlation energy and an SCF convergence criteria of 10 −6 eV/Å and 5 × 10 −4 eV/Å for the density and forces, respectively. Relativistic effects are included in terms of the scalar atomic ZORA approach and all other settings are taken to be the default in FHI-aims. For all calculations we use the light basis sets and numerical settings in FHI-aims. These settings ensure a convergence in lattice constants of ±0.1Å and a relative accuracy in phonon frequencies of 3%.

A. Rock salts, Zincblende, and Wurtzites
To understand how prevalent harmonic, Si-like materials are across a broader chemical space, we use σ A OS to screen over an initial test set that includes 97 rock salt (RS), 67 zincblende (ZB), and 45 wurtzite (WZ) binary and elemental solids, as summarized in Figure 14.
At 300 K only 35% of all 209 materials tested can be classified as highly harmonic with a σ A OS < 0.2. At elevated temperatures anharmonic effects get significantly stronger. At 500 K only 10% of the materials have a σ A OS < 0.2, while 34% feature a σ A OS value > 0.4. For both temperatures, zincblende and wurtzite materials are more harmonic on average, while the majority of rock salts has a σ A OS value indicative of more complicated dynamical processes. These results are in line with experimental and theoretical studies that show materials with a higher coordination number have longer bond lengths and softer lattices leading to stronger anharmonic interactions [43,44]. This screening demonstrates that anharmonic effects are more prevalent in material space than previously thought, particularly at technologically relevant temperatures.
One of the most anharmonic subclasses in the initial set of materials is noble metal halides, seven of which are among the eleven most anharmonic materials at 300 K, with only one of the remaining three materials having a σ A OS < 0.4. In order to analyze the suspected high anharmonicity of this class, we use aiMD 1 to accurately calculate σ A for all stable noble metal halide zincblende and rock salt materials at 300 K, as summarized in Table I. In all cases, the aiMD calculations confirm (σ A MD ≥ σ A OS ) the strong anharmonicity indicated qualitatively by the one-shot approach. Quantitatively, the aiMD reveals the the actual values of σ A can be even substantially higher. A closer inspection of the dynamics reveal that this is related to spontaneous defect formation, which we observe for every material except zincblende CuI. In these cases, the noble metal ions move into the the metastable interstitial sites of the lattice forming metallic clusters within the material, as illustrated in the inset of Figure 15b for rock salt AgBr. Similar effects have been observed in earlier aiMD studies of cuprous halides [45,46] and have been debated extensively in literature [47][48][49]. Al- though a discussion of these defect-related aspects goes beyond the scope of this work, we note that large jumps in σ A (i.e. σ A approaches or exceeds 1.0) are correlated with the formation of defects. The jumps are a result of the harmonic model obviously failing to account for the occupation of interstitial sites, as can be seen in Figure 15a. The actual magnitude of these σ A -jumps is determined by the number of defects formed in the structure and by the magnitude of the distortion from the pristine structure, whereby the occurrence of such jumps leads to strong fluctuations in σ A , as also summarized in Tab. I. This data reveals that these defects are particularly pronounced for zincblende CuCl, CuBr, and AgBr, leading to σ A values much greater than 1.0. In this context, we would like to stress that neither the employed trajectory length nor the used supercell size is sufficient for an accurate description of defect formation in thermodynamic equilibrium. Nevertheless, the metric σ A is a useful indicator for where interesting, highly-anharmonic lattice dynamical phenomena occur.

B. Perovskites
As discussed in Section IV and shown in Figure 10 as KCaF 3 approached the cubic phase transition σ A quickly rose to ∼ 1.0. To see how general this behavior is we calculate σ A at 300 K and 600 K for a set of ten orthorhombic perovskites with phase transition temperatures span- ning three orders of magnitude and summarize the results in Table II. As the data in Table II illustrates once the perovskites are around or above their respective transition temperatures, σ A tends to go above 0.45, with its overall magnitude determined by the extent of deformation away from the harmonic reference structure. These results combined with the data from the previous sections indicate that σ A could be useful as an indicator for phase transitions or defect formation in material. To get a better understanding of how individual modes behave at different points near transition temperatures we illustrate the mode resolved σ A s values for CaZrO 3 , NaTaO 3 , KCdF 3 , and CsSnI 3 at 300 and 600 K using violin plots in Figure 16. When the material is far from the phase-transition temperature the mode projection is similar to what was seen for Si in Figure 8 with the relative anharmonicity grouped together in a band centered around the average σ A value for the material, as seen for CaZrO 3 . As a material approaches a transition (e.g. NaTaO 3 at 600 K or KCdF 3 and CsSnI 3 at 300 K) the cloud broadens with a few highly anharmonic modes present. Finally as the temperature further increases, many more modes become anharmonic as the harmonic model fails to even qualitatively describe the system. These results demonstrate the potential of this measure to not only predict when a potential phase transition will happen, but also which vibrational modes are responsible for that transition.

VII. CONCLUSION AND OUTLOOK
In this work, we present a measure for the degree of anharmonicity, σ A (T ) that quantifies the importance of anharmonic effects in a crystalline material. In practice, this is done by statistically analyzing the ab initio interactions in thermodynamic equilibrium. This measure allows for a rapid scan through material space for the purpose of ranking materials by anharmonicity. Our results indicate that materials whose properties are significantly affected by anharmonic effects are not uncommon at all. Rather, largely harmonic materials like silicon and diamond with σ A < 0.2 are the exception to the rule. In fact, only 35% of the binary compounds and none of the perovskites we screened over were that harmonic at room temperature. At more elevated temperatures, anharmonic effects are even more prevalent.
The proposed metric gives access to key aspects of anharmonicity itself, since σ A (T ) is rigorously based on the actual interactions driving the dynamics in a material. As demonstrated for the phase transitions occurring in perovskites, analyzing the per-mode contributions to σ A sheds light on the microscopic origin of anharmonic effects that determine the macroscopic properties of materials in thermodynamic equilibrium.
As an outlook, we present an additional finding in Fig. 17, in which the experimental lattice thermal conductivity κ L at 300 K is plotted against σ A (300 K) for those RS/ZB/WZ compounds with reliable measurements of κ L [60].
Fitting the data on log-log scale to a linear model, we get a slope of -4.79, implying an inverse power law between κ L and σ A , with an average factor difference (AFD) of 1.48. AFD was introduced by Miller and coworkers to measure the accuracy of a model via |log κ L,exp − log κ L,model | , (22) where N is the number of samples in the test set [43]. The strong inverse correlation between these two properties shown in Fig. 17 illustrates that σ A is a good descriptor for anharmonicity by itself because as a material's vibrational properties become more anharmonic, its phonon lifetimes, and therefore κ L , decrease. It is remarkable that even without explicitly including any of the other material properties that influence κ L , such as group velocities or heat capacities, we get a similar AFD as other semi-empirical models [43,60]. Clearly, this calls for future, more extensive research on more exhaustive datasets. Along these lines, we like to stress again that the scope of the presented method goes beyond thermal transport and phase transition mechanisms, and potentially applies to any phenomenon governed by anharmonic effects such as free energies [61], thermal stability [62], thermal expansion [63], defect formation [64,65], ferroelectricity [66], and electron-phonon coupling [67]. Besides these practical aspects, these findings call for a systematic analysis and scrutiny of the validity of perturbative techniques, which are commonly used, e. g., for computing lattice thermal conductivities. In these approaches, anharmonic effects are treated by perturbation theory, starting from the harmonic approximation and assuming the perturbation to be small V A V. This assumption seems to be justified for materials like silicon with σ A < 0.2, in which anharmonic effects are responsible for 20% of the interatomic interactions at most. However, this assumption becomes highly questionable for the majority of materials, which, as shown in this work, exhibit σ A > 0.2. Especially thermal insulators generally feature large values of σ A , cf. Fig. 17. The formalism developed in this work is ideally suited to single out and classify well-defined test systems with different anharmonic strength and character, ranging from simple harmonic materials (σ A ≤ 0.20) up to complex materials featuring phase transitions (σ A ≥ 1). Across this anharmonicity range, non-perturbative methodologies such as non-equilibrium techniques [68][69][70][71] and equilibrium Green-Kubo approaches [9,10] can provide reliable benchmarks, against which the various perturba-tive techniques at different degrees of sophistication [11][12][13][14][15][16][17][18][72][73][74][75][76][77][78] need to be validated. This comparison will allow to identify up to which strength of anharmonicity σ A these different techniques work reliably and above which threshold of σ A perturbation theory breaks down completely.