Towards the high-throughput prediction of finite-temperature properties using the quasi-harmonic approximation

Lattice dynamics calculations within the quasi-harmonic approximation (QHA) provide an infrastructure for modelling the finite-temperature properties of periodic solids at a modest computational cost. With the recent widespread interest in materials discovery by data mining, a database of computed finite-temperature properties would be highly desirable. In this work we provide a first step toward this goal with a comparative study of the accuracy of five exchange-correlation functionals, spanning the local density approximation (LDA), generalised-gradient approximation (GGA) and meta-GGA levels of theory, for predicting the properties of ten Group 1, 2 and 12 binary metal oxides. We find that the predictions are bounded by the LDA, which tends to underestimate lattice parameters and cell volumes relative to experiments, but yields the most accurate results for bulk moduli, expansion coefficients and Grüneisen parameters, and the PBE GGA, which shows the opposite behaviour. The PBEsol GGA gives the best overall predictions of the lattice parameters and volumes whilst also giving relatively reliable results for other properties. Our results demonstrate that, given a suitable choice of functional, a variety of finite-temperature properties can be predicted with useful accuracy, and hence that high-throughout QHA calculations are technically feasible.


Introduction
Accelerating the design and optimisation of functional materials is key to the development of important contemporary technologies such as those underpinning sustainable energy.Computational modelling is a powerful tool for materials Original Content from this work may be used under the terms of the Creative Commons Attribution 4.0 licence.Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.discovery, offering a relatively fast and accurate means to predict material stability and properties.This capability is bolstered by the continued investment in contemporary supercomputers which, together with developments in theory and numerical algorithms, are enabling ever more realistic simulations to be performed on increasingly complex materials.Firstprinciples methods such as density-functional theory (DFT) in particular can often predict material properties with experimental accuracy and are thus powerful tools for materials discovery, characterisation and optimisation [1,2].
In recent years the identification of new materials has increasingly been facilitated by high-throughput (HT) modelling studies.Such studies entail the computational characterisation of thousands or even tens of thousands of materials using automated workflows [2,3].DFT is among the most widely used and successful theoretical techniques for HT studies, owing to its excellent balance of accuracy and computational cost, and modern DFT codes have the ability to predict a wide range of structural, mechanical, electronic and optical properties [4].There is also an increasing drive towards providing access to the computed properties via online databases, and such datasets can serve as a basis for further screening studies and to support experimental investigations.This is the spirit of the materials genome initiative [5], which aims to use computational techniques to predict, screen, and optimise material properties [1].Among several notable recent efforts in this area is the materials project (MP), which provides a consistent set of calculated properties for the majority of the published structures in the Inorganic Crystal Structure Database (ICSD) together with a powerful infrastructure for searching, downloading and exploiting this dataset [1].
To date, there have been a number of examples applying HT screening to the discovery of new functional materials, and many of the resulting predictions have subsequently been validated in experiments [6].A good example is the computational screening for new, more efficient and less expensive thermoelectric materials for direct heat-to-electricity power conversion [2,[7][8][9], as DFT calculations can readily access the electronic and heat transport properties that form the basis of the thermoelectric figure of merit ZT. [10][11][12][13][14] HT screening has also been applied to the rational design of organic functional materials with tuneable optoelectronic properties for applications in organic photovoltaics and organic lightemitting diodes.[2,15,16] Other noteworthy examples include several studies targeting the discovery and optimisation of new and improved battery materials [2,[17][18][19][20].
HT screening has the potential to become even more powerful as the datasets they generate provide a foundation for socalled 'materials informatics', [21] which aims to accelerate materials discovery using data mining, statistical analysis and machine learning methods [2].Such is the potential of materials informatics that data-driven discovery has been referred to as the fourth paradigm of science [22], with the possibility of uncovering new materials processing-structure-propertyperformance relationships that could completely revolutionise materials science [2,23].The MP database in particular has formed the basis for a number of data-mining studies with the goal of identifying new trends and predictive empirical rules [24].For instance, Hautier et al [25] performed a HT computational search over the thousands of binary and ternary oxides available in the MP to identify compounds that can be used as p-type transparent conducting oxides (TCOs).TCOs are characterised by high electrical conductivity and transparency to visible light, and their development is essential to technologies including transparent electronics and solar cells [24][25][26].This data mining study led to new design principles for the production of p-type TCOs with lower effective masses and highlighted promising avenues for future investigation in this area [24,25].
Despite the many successes of DFT, a major drawback is that most routine studies work from frozen snapshots of the atomic structure and only take into consideration the average positions of the atoms, neglecting thermal atomic motion.The neglect of temperature effects leads to a gap between theory and experiments, and highlights the importance of developing new HT theoretical methods that can go beyond 'static' calculations.
There are several methods of incorporating temperature effects into calculations, among which lattice-dynamics calculations offer a good balance of computational cost and accuracy [27][28][29].The theory of lattice dynamics builds on the quantum harmonic oscillator model to provide an infrastructure for modelling the atomic vibrations (phonons) in periodic solids.Within the harmonic approximation (HA), the atoms are modelled as classical masses connected by harmonic springs, such that when they move within the crystal lattice they exert forces on neighbouring atoms, possibly in neighbouring unit cells.These vibrations are assumed to be small in amplitude, and the Taylor expansion of the crystal potential energy φ with respect to the atomic positions is truncated to second order [27].The basic HA allows for the prediction of the harmonic phonon dispersion and density of states (DoS) curves together with thermal properties at constant volume including and the volumetric heat capacity C V , the phonon entropy S V , and the Helmholtz free energy F [27,30].However, the HA assumes the average bond lengths do not change from equilibrium and therefore cannot predict changes in volume with temperature.The quasi-harmonic approximation (QHA) addresses this by taking into account the volume dependence of the phonon frequencies [30], and extends the basic HA to enable modelling of the temperature dependence of physical properties including the unit-cell volume V, the thermal expansion coefficient α, the bulk modulus B, the constantpressure heat capacity C p and the thermodynamic Grüneisen parameter γ.
Given the wide range of properties calculable using the QHA, it is of interest to explore the feasibility of HT QHA calculations.Of the large number of lattice-dynamics studies reported in the literature [31][32][33][34], relatively few have employed the QHA.A notable implementation of HT QHA calculations is that of Nath and coworkers [35].Their automated framework was able to predict finite-temperature properties such as the C p , B and γ with good accuracy, and in the latter case with considerably better accuracy than the more approximate quasi-harmonic Debye model [36].This model was, however, less successful in predicting the α, which may be related to the (somewhat arbitrary) choice of the PBE functional.Indeed, a comparative study of the performance of seven exchange-correlation functionals for QHA calculations on the IV-VI semiconductors PbS, PbTe, ZnS, and ZnTe demonstrated that PBE was among the worst-performing for predicting finite-temperature properties using the QHA [28].Similarly, Erba et al compared five functionals for QHA calculations on MgO and CaO and found considerable variation among different methods [37].Further systematic studies on the performance of different exchange-correlation treatments for calculating the finite-temperature properties of a wider range of materials is therefore an important step toward establishing the technical feasibility of HT QHA calculations.
Binary oxides are ideal candidates for systematic quasiharmonic studies for several reasons.Firstly, they are twoelement compounds and mostly adopt crystal structures with high-symmetry cubic space groups.Secondly, many of these materials have been well characterised and have experimental measurements of structural, mechanical and thermal properties available.Finally, oxides are important materials for a wide range of potential applications.Superionic oxides, such as Li 2 O, Na 2 O and K 2 O, which exhibit high ionic conductivity, can be used in photocathodes [38,39], catalysis [40], solid-state batteries [41] and in fuel cells [42].The alkali metal oxides MgO, CaO, SrO and BaO have applications in catalysis [43], microelectronics [43], refractory systems and plasma displays [44].The Group 12 oxides ZnO and CdO are promising TCOs and as such have recently attracted a large amount of scientific attention [45,46].ZnO is non-toxic and relatively abundant, is characterised by a large band gap and low resistivity [45], and has outstanding optical and piezoelectronic properties that make it an ideal candidate for use in visible and ultraviolet light emitters, laser diodes and transparent field-effect transistors [47][48][49][50].Similarly, CdO has found applications in solar cells, gas sensors, phototransistors and many other electronic and optoelectronic devices [46,51].
In this study, we present a comprehensive assessment of the accuracy of five exchange-correlation functionals for predicting the finite-temperature properties of ten cubic binary oxides using the QHA.We examine a selection of popular functionals used in solid-state physics, viz. the localdensity approximation (LDA), the PBE [52] and PBEsol [53] generalised-gradient approximation (GGA) functionals, PBE with the semi-empirical DFT-D3 dispersion correction [54] (i.e.PBE+D3), and the SCAN meta-GGA [55].We test against a selection of binary metal oxides drawn from Groups 1, 2 and 12 of the periodic table in the cubic Fm 3m structure, viz.Li 2 O, Na 2 O, K 2 O, Rb 2 O, MgO, CaO, SrO, BaO, ZnO and CdO.Theoretical predictions are compared quantitatively against available experimental data, based on which we establish trends in the performance of the different functionals and identify the best choice for potential HT QHA calculations.

Ab Initio Lattice dynamics calculations
2.1.The HA Within the HA, the second-order interatomic force-constant matrices Φ(lκ, l ′ κ ′ ) are computed as: where φ is the potential energy of the crystal, u (lκ) = r (lκ) − r 0 (lκ) are the displacements of the κth atoms in the lth unit cells from their equilibrium positions r 0 , F (lκ) are the restoring forces, and the labels α and β denote the Cartesian directions x, y and z. [30] The Φ are used to construct the dynamical matrix D(q) for a given phonon wavevector q according to: where m κ are the atomic masses.Diagonalisation of D(q) yields the 3n a phonon frequencies ω qj and corresponding displacement vectors W qj at the wavevector q, where j is the band index [30].
The harmonic phonon frequencies can then be used to compute the phonon DoS g(ω) and phonon dispersion curves ω(q).g(ω) shows the number of modes as a function of ω integrated over the entire space of phonon wavevectors defined in the Brillouin zone, and is given by: [27,30] where N is the number of unit cells in the crystal, equivalent to the number of wavevectors included in the summation.The frequency dispersion ω(q) on the other hand shows how each of the 3n a frequencies evolve with q over a set path, typically passing through the high-symmetry wavevectors in the irreducible Brillouin zone.Knowledge of the phonon frequencies allows for a variety of thermodynamic quantities to be computed.The zeropoint energies and the thermal population of phonon energy levels at finite temperature contribute to the constant-volume (Helmholtz) free energy F according to: where U L is the lattice internal energy, equivalent to the φ in equation ( 1), F V is the phonon free energy and U V and S V are the corresponding internal energy and entropy.F V can be obtained from the vibrational partition function Z V (T) defined as: where k B is the Boltzmann constant.
The Z V (T) is then used to obtain F V via the bridge relation: The phonon internal energy U V can be calculated separately as: where n(T) is the phonon occupation number (energy level) determined by the Bose-Einstein distribution.The phonon entropy S V is given by: [27, 30] Finally, the heat capacity at constant volume C V (T) can also be computed from the phonon frequencies according to: [27,29] We note that in these expressions the F V , U V , S V and C V are computed for all the unit cells in the crystal and are typically normalised by 1/N to obtain them per (primitive) unit cell.

The QHA
The HA assumes that the bond lengths oscillate about the equilibrium values and thus cannot describe volume-dependent properties.The QHA improves on the HA by considering the volume dependence of the phonon frequencies.This is done by performing harmonic calculations over a series of volumes about the equilibrium [30].The constantpressure (Gibbs) free energy G(T, p) is found by minimising the Helmholtz energy F(V, T) with respect to the volume V at a target temperature T: [30] The minimisation is usually performed by fitting the energy-volume curves to an equation of state (EoS) such as the Birch-Murnaghan equation [56], which yields the minimum energy F 0 , the corresponding equilibrium volume V 0 and also the bulk modulus B 0 and pressure derivative B ′ 0 : Once the volume as a function of temperature has been determined, the volumetric thermal expansion coefficient α V , which describes the rate at which the volume increases with temperature, can be calculated as: The constant-pressure heat capacity C p (T) can be computed from the G(T) as: [27] Finally, the QHA also allows derivation of the thermodynamic Grüneisen parameter γ, which describes how the phonon frequencies change with volume and thus indirectly how temperature affects the lattice dynamics: γ is the average of the (microscopic) mode Grüneisen parameters γ qj , which are given by: [30] The γ qj are also calculable from the third-order force constants Φ(lκ, l ′ κ ′ , κ ′ ′ l ′ ′ ) and are therefore used as a measure of anharmonicity.

Computational methods
In this work, first-principles calculations were performed using pseudopotential plane-wave DFT as implemented in the Vienna ab initio simulation package (VASP) code [57].
Electron exchange and correlation were modelled using the LDA, the PBE [52] and PBEsol [53] GGA functionals, PBE with the DFT-D dispersion correction (PBE+D3), [54] and the SCAN meta-GGA [55].For consistency, we tried to keep the key technical parameters similar between different sets of compounds.Core electrons were treated using the projector augmented-wave (PAW) pseudopotential method [58,59], with the Li 2s, Na and Mg 3s/2p, K and Ca 4s, Rb and Sr 5s, Ba 6s, Zn 4s/3d, Cd 5s/4d and O 2s/2p electrons included in the valence shells.The valence Kohn-Sham orbitals were represented in a plane-wave basis with a cut-off of 700 or 800 eV and the k-point sampling shown in table 1.Initial models of the ten compounds studied were obtained from the MP database [1] and fully optimised to strict energy convergence criteria of 10 −8 eV on the electronic total energy and 10 −2 eV Å −1 on the forces.
Quasi-harmonic lattice dynamics calculations were performed on the optimised structures using the Phonopy software package [30] with VASP as the force calculator.The supercell expansions and k-point sampling used for the supercell finite-displacement calculations are listed in table 1.The force constants of Li 2 O, Na 2 O, K 2 O and Rb 2 O were computed using 4 × 4 × 4 expansions of the three-atom primitive cells with 192 atoms.The force constants of MgO, CaO, SrO, BaO, ZnO and CdO were calculated using 3 × 3 × 3 expansions of Table 1.Parameters used for the geometry optimisations and phonon calculations performed in this work: phonon supercell expansion (number of atoms), k-point sampling used for the optimisation and phonon supercells, and plane-wave cut-off.All k-point meshes are given as the number of subdivisions on Γ-centred Monkhorst-Pack meshes [60].

Group
Compound Phonon SC (no.atoms) Optimisation Phonon Cut-off (eV) the eight-atom conventional cells with 216 atoms and the following transformation matrix applied to map to the primitive cells during post processing:   0.0 0.5 0.5 0.5 0.0 0.5 0.5 0.5 0.0 The phonon DoS curves were produced by interpolating the phonon frequencies onto q-point meshes with 48 × 48 × 48 subdivisions and using the linear tetrahedron method for Brillouin-zone integration.The same q-point grids were used to calculate the thermodynamic properties.In order to generate the phonon band structures, the phonon frequencies were interpolated along lines of q-points passing through the highsymmetry points in the Brillouin zones of the primitive cells.
Since the oxides are polar solids, there is substantial splitting between the longitudinal and transverse optic modes (LO/TO splitting) in the dispersions [27], and a non-analytical correction needs to be applied to the dynamical matrices at q → 0 to account for this.For calculations with the LDA, PBE, PBEsol and PBE+D3 functionals the required Born effective-charge tensors and high-frequency dielectric constants were computed using the density-functional perturbation theory routines in VASP [69].For SCAN, these parameters were computed from the self-consistent response to finite electric fields according to the PEAD approach [70,71] with a field strength of 10 −2 eV Å −1 .
QHA calculations were initially performed by running phonon calculations on optimised structures with up to 5% compressions and expansions of the cell volume about the athermal equilibrium V 0 , in steps of 1%.We discovered post factum that this range was insufficient for many of the compounds, and so for the LDA, PBE, PBEsol and PBE+D3 calculations we extended it to ±15%, with steps of 3% outside the initial ±5% range.We were not able to do this for the SCAN calculations due to a combination of the higher computational cost and also to some convergence issues outlined in the discussion below.Therefore, for the LDA and GGA calculations a total of 19 volumes were included in the QHA calculations covering a range of ±15% of V 0 , while for the SCAN meta-GGA calculations 11 volumes were included covering a range of ±5% of V 0 .The QHA post-processing was performed without an applied pressure (i.e. the pV term in equation ( 10) was set to zero).

Athermal properties: optimised lattice parameters
The optimised equilibrium lattice constants a 0 of the ten binary oxides calculated using the five exchange-correlation treatments are listed together with the % differences to literature values in table 2.
All functionals yield reasonably accurate predictions, with the largest absolute error being less than 4%.Comparing the average % differences and standard deviations, some general trends can be observed.With the exception of ZnO, the LDA consistently underestimates the lattice parameters (−1.95 ± 1.27%), which is expected given the documented tendency of LDA to overbind [53].With the exception of MgO, BaO and ZnO, PBE yields better results for most compounds than the LDA, but consistently overestimates the lattice constants (+1.08 ± 1.04%), which is again anticipated from its known tendency to do so.The PBEsol and SCAN functionals both give relatively good predictions, with average errors of −0.39 ± 0.98 and −0.39 ± 1.03% respectively.Interestingly, although one would expect predominantly electrostatic bonding in the oxides, the dispersion-corrected PBE+D3 functional appears to give the best overall agreement with experiments with an average deviation of +0.24 ± 1.28%.
We should note however that the DFT values are 'athermal' and do not account for temperature effects, whereas the experimental measurements can be expected to include some degree of thermal expansion, the magnitude of which will depend on the temperatures at which the structures were collected.Assuming positive thermal expansion, it is likely that Table 2. Optimised equilibrium lattice constants a 0 of the ten binary oxides examined in this work obtained using the LDA, PBE, PBEsol, PBE+D3 and SCAN functionals.The experimental lattice constants and the % difference between the calculated and measured values are shown for comparison, and the final two rows of the table show the average and standard deviation of the % differences predicted by each functional.the functionals that predict smaller lattice constants than the measured values (here LDA, PBEsol and SCAN) can be considered to provide a better reproduction of experiments.

Harmonic lattice dynamics
We now compare the harmonic phonon spectra and constantvolume thermodynamic functions.It should be noted that the comparisons are performed over a temperature range from 0 K up to 2/3 the melting temperature T m , above which anharmonic effects typically cause the HA to deviate significantly from experiments.For many of the compounds, calculations performed using the SCAN meta-GGA functional showed significant artefacts in the spectra.While in some cases this could be rectified by using finer charge-density grids, in others the required grids appeared to be so large that they resulted in excessive memory requirements.The numerical sensitivity of meta-GGA functionals to the grid density arises because of the dependence on the second derivative of the density, and has been documented previously [28].As a result, the SCAN calculations on Na 2 O, K 2 O, CaO and ZnO are excluded from this discussion (see figures A1-A4 in the supplementary material).

Harmonic phonon dispersion and DoS.
The harmonic phonon DoS and phonon dispersion curves of the Group 1 (Li 2 O, Na 2 O, K 2 O and Rb 2 O) and Group 2 oxides (MgO, CaO, SrO, BaO) computed using the five functionals are compared in figures 1 and 2 respectively.Where available we have compared the calculated dispersions to neutronscattering measurements.Similar comparisons for the Group 12 oxides ZnO and CdO can be found in figure A5 (supplementary material).
As noted above the SCAN phonon spectra of Na 2 O, K 2 O, CaO and ZnO are not included in these comparisons.We also note that the SCAN phonon calculations on Li 2 O required an expanded charge-density grid with ∼30% more points in order to remove spurious imaginary modes from the spectra (see figure A6, supplementary material).
Leaving aside the technical issues with SCAN, as expected there are no imaginary modes present in any of the calculated dispersions, which indicates that all the functionals predict these materials to be dynamically stable in the cubic structure.As shown in equations ( 1) and ( 2), the phonon frequencies are obtained from the restoring forces on the atoms in response to small-amplitude displacements.The phonon frequencies are generally strongly dependent on the cell volume, so the primary effect of the functional is likely to be through differences in the predicted equilibrium volume (cf table 2), as was found to be the leading effect in phonon calculations on several IV-VI binary chalcogenides [28].
As for the lattice constants, we find that figures 1 and 2 reveal some clear trends.Across all the spectra PBE predicts the lowest frequencies, while LDA/SCAN predict the highest, which is consistent with calculations on Fe 2 VAl and Fe 2 TiSn [77].In the Group 1 and 12 oxides the differences between the functionals are most evident in the high-frequency optic branches of the dispersions, whereas the optic modes in the Group 2 oxides appear to be much less sensitive to the exchange-correlation treatment.Despite these differences, all of the calculated phonon dispersions of Li 2 O and the Group 2 oxides show good agreement with experimental measurements, and since the differences between functionals are subtle it is not possible to identify a best-performing one based solely on an assessment of the calculated dispersions.

Vibrational entropy.
The vibrational entropy of many of the oxides has been measured and provides an additional set of data for comparison.
Figure 3 compares the calculated S V (T) curves for all ten oxides obtained with each of the five functionals to experimental data where available.This thermodynamic quantity results from the occupation of the phonon modes at finite temperature, and thus functionals that predict stiffer frequencies are expected to predict  [73], CaO [74], SrO [75], and BaO [76] captured using the engauge digitizer software [72] are overlaid as black markers for comparison.
lower entropy and vice versa.We indeed find that the trends in the calculated phonon spectra are reflected in the vibrational entropies, with LDA and SCAN typically predicting the Experimental data obtained by Chase [78], Robie et al [79] and Chopelas [80] are overlaid as black markers for comparison.Note that there are three sets of experimental for MgO in (e), of which the data from Chase and Chopelas et al [78,80] are measured to 1100 K and that from Robie et al [79] to 1800 K. Despite these trends, the entropy appears to be less dependent on the exchange-correlation functional than the phonon frequencies themselves, especially for the Group 2 oxides.This can be explained by the variation in the predicted phonon frequencies being larger for the higher-frequency modes, which, due to the Bose-Einstein weighting in the phonon occupation number, are less heavily populated and therefore have a smaller impact on the entropy.The relatively small differences in the vibrational entropy calculated with different functionals are consistent with the results of a similar study by Erba et al [37] The computed S V (T) curves in figure 3 also compare favourably to experiments.Some minor deviation occurs at elevated temperature, which can be explained by quasiharmonic effects (i.e.thermal expansion) and/or increasingly prominent higher-order anharmonicity [27].

Constant-volume heat capacity. An additional thermodynamic property that can be obtained through the HA is the constant-volume heat capacity C V (T).
The C V (T) curves for MgO predicted using the five functionals are shown in figure 4 and compared to experimental measurements from Chopelas and Barron et al [80,81] Data for all ten oxides is presented in figure A7 (supplementary material).
Analysis of figures 4 and A7 suggests that, as for the vibrational entropy, the volumetric heat capacity appears to be largely independent of the exchange-correlation treatment.This can be accounted for using the same reasoning as for the entropy, since the modal contributions to C V (T) are weighted by exp(hω qj /k B T) (cf equation ( 9)), and our findings are again in line with the conclusions of Erba et al [37].
Although experimental data was limited to MgO, the calculations with all five functionals appear to be in excellent agreement with experimental measurements.

Quasi-harmonic lattice dynamics
The QHA improves on the HA to take into account the volume dependence of the lattice energy U L (i.e. the crystal potential energy φ) and the phonon frequencies to enable a wider range of properties to be computed.
To perform a QHA calculation, the lattice energies and phonon spectra are computed over a series of expansions and compressions of the unit cell about the calculated equilibrium volume V 0 .Figure 5 shows the phonon dispersion curves of MgO from Γ to L computed for up to 2% compression and expansion around V 0 using each of the five functionals.
Analysis of the band dispersions suggests that the choice of functional has the same relative effect on the phonon spectra at compressed and expanded volumes.The frequencies generally increase (stiffen) at compressed volumes and decrease (soften) under expansion.The LDA consistently predicts the highest frequencies while PBE consistently predicts the lowest, as observed in the computed phonon spectra of the equilibrium structures.This again highlights the strong dependence of the phonon frequencies on the unit-cell volume and lends support to the idea that the primary differences between functionals most likely arises from differences in the predicted equilibrium volumes as found previously [28].We note also that all five functionals predict MgO to remain dynamically stable under the ±2% volume changes, as evidenced by the absence of imaginary modes in the spectra.
In this section we compare QHA property calculations with each of the functionals over a temperature range of 0 K to up to 2/3 T m .To aid the discussion, a table comparing the calculated QHA parameters to experimental data at room temperature (∼300 K), where available, is provided as table A1 (supplementary material).
We note that the difficulties we encountered with the harmonic calculations using SCAN were amplified in the QHA calculations.As a case in point, we have omitted the SCAN dispersions for the 99% and 102% volume scaling in figure 5 due to the presence of spurious imaginary modes (cf figure A8, supplementary material).For this reason we only present SCAN QHA results for Rb 2 O, SrO and BaO in this section.

Unit-cell volume.
The QHA allows the prediction of the cell volume as a function of temperature, which in principle allows for better comparison to experimental measurements than the athermal DFT lattice constants discussed above.
The V(T) curves are obtained by fitting a free-energy EoS, formed from the harmonic Helmholtz free energies as a function of volume at a specified temperature, to a model function, in this case the Birch-Murnaghan EoS in equation (11).[56] Free energy EoS fits for the ten oxides over a range of temperatures are shown in figures A9-A18 (supplementary material).We found that the Birch-Murnaghan EoS fits for BaO and ZnO were questionable, in particular at elevated temperature (figures A16 and A17, supplementary material).We also identified several cases where the predicted thermal expansion at elevated temperature was outside the range of volume Figure 6 compares the calculated V(T) curves obtained for the Group 2 oxides with the five different functionals to experimental measurements, and data for the Group 1 and 12 oxides is provided in figures A19 and A20 respectively (supplementary material).In addition, figure 7 illustrates the percentage errors in the predicted room-temperature (300 K) cell volumes with respect to experimental measurements for all ten compounds and all five functionals, with the exception of some missing SCAN results as outlined above.
In contrast to the harmonic properties we find that the choice of functional has a significant impact on the predicted V(T) curves.Comparison of the five functionals reveals broadly the same trends as were found for the athermal optimised lattice constants.In particular, LDA consistently predicts the lowest volumes and PBE the highest, which matches the findings of Erba et al [37].
The accuracy of the calculations relative to experiments varies considerably for the different exchange-correlation functionals.In particular, the LDA and PBE functionals produce the least accurate volume predictions, with average deviations of −3.20 and +6.60%, respectively, at 300 K.For most of the oxides PBEsol shows the best agreement with measurements across the full temperature range, with an average deviation of +1.82% at room temperature.
The fact that the deviations from experiment are generally larger than the errors in the athermal lattice constants (cf table 2) suggests that the QHA 'amplifies' the tendency of the functionals to under-or over-estimate cell volumes.Given this observation, the issues with the PBE EoS fits for Li 2 O, CaO and SrO can be attributed to an anomalously-large predicted thermal expansion from an exaggeration of the documented underbinding behaviour.(The issues with the three SCAN calculations are simply due to the higher computational cost of these calculations, due to which we only performed calculations over a ±5% range of volume changes rather than the larger ±15% range for the other functionals; see Computational Methods in Section 3).
For most of the compounds the five functionals predict the V(T) to have similar curvature, indicating that the main difference between them is the predicted equilibrium volume.A notable exception is SrO, for which PBE predicts a markedly different expansion to the other functionals at higher temperature.This can be explained by the ±15% range of volumes included in the calculations being insufficient.The five functionals also predict significantly different curvatures for BaO, which we attributed to the questionable EoS fits.Interestingly, however, although the EoS fits for ZnO had similar issues, all five functionals predict similar curvatures of the V(T) (figure A20, supplementary material).
Finally, we note in passing that for cubic systems the lattice constant can be trivially calculated from the volume as a = V 1/3 , whereas for a general unit cell with two or more independent lattice parameters this is not the case.For such systems one would in principle need to minimise a multidimensional free-energy EoS over all of the independent parameters (cf equation ( 10)), although an alternative, more practical, method is to neglect the effect of 'thermal pressure' on the individual lattice parameters and perform an athermal constant-volume cell optimisation at the QHApredicted volumes of interest.

Thermal expansion. Figure 8 compares the volumetric expansion coefficients as a function of temperature α V (T)
for MgO and CaO obtained with LDA, PBE, PBE+D3 and PBEsol.Data for the all ten oxides is provided in figures A21-A23 (supplementary material).
α V (T) is obtained from the volume and the gradient of the V(T) curves (cf equation ( 12)).Since it depends on V it is similarly sensitive to the choice of functional, as was observed by Erba et al [37] In some cases, we found that the predicted α V (T) showed unphysical oscillations at elevated temperature, which we attributed to numerical noise in the derivatives, and in view of this we limited the range of temperatures computed for SrO with PBE and BaO with SCAN to 1600 and 900 K respectively.We were only able to find experimental measurements of α V (T) for Li 2 O, MgO, CaO and BaO, plus a single measurement for SrO at room temperature.Bearing this in mind, across all the QHA properties examined in the present study the α V appear to be predicted least accurately.It can be seen from table A1 that at 300 K the best-performing functional is the LDA, but even this has an average error of 16.18%, while for the GGA functionals, the errors range from 20 to Percentage errors in the 300 K unit-cell volumes of the ten binary oxides predicted using the five exchange-correlation functionals tested in this study.We note that data for some materials using SCAN are not included due to issues with these calculations (see text).The experimental data used for comparison are taken from references [61,62,64,68,[82][83][84][85][86][87][88].30%.As discussed above we were only able to perform SCAN calculations on two of these compounds, and over a limited range of volume changes, so we cannot reasonably compare it to the other functionals.The predicted expansion coefficients also tend to show a marked deviation from experiment at higher temperatures, which could be attributed to the neglect of higher-order anharmonicity in the QHA.This conclusion is supported by recent work from Masuki et al who proposed a non-perturbative method to correct the QHA for intrinsic anharmonic effects that resulted in much better estimates of the experimental α of MgO [92].It is perhaps worth noting also that one might expect the expansion to be more sensitive to subtle numerical errors and/or to anharmonic effects than the volume itself due to the differentiation.Despite these issues, a similar trend to the V(T) predictions is observed among the different functionals whereby for a given material the LDA predicts the lowest expansion coefficients, PBE predicts the highest, and PBEsol and PBE+D3 predict intermediate values.

Bulk Moduli.
Figure 9 compares the predicted temperature dependence of the bulk moduli B(T) of the Group 2 oxides to experimental measurements, and data for the Group 1 and 12 materials can be found in figures A24 and A25 respectively (supplementary material).
As for the unit-cell volume we find that the choice of functional has a significant impact on the predicted B(T).Notably, PBE consistently underestimates the moduli and yields the least accurate predictions relative to experiments, whereas LDA produces the largest estimates and gives the best agreement for all ten oxides.This is clearly shown in a comparison of the percentage differences to experiments at room temperature (figure 10, table A1 in the supplementary material).The average percentage difference for the LDA is 3.52%, while that for PBE is almost 5 × higher at 16.61%.The corresponding errors for PBEsol and SCAN are on the order of 10%, which in context might be considered acceptable, although the errors increase at higher temperature.As for the thermal expansion coefficients, the larger errors in the predictions may be attributable to the neglect of higher-order anharmonic effects in the QHA.
For most of the oxides the different functionals predict a similar-shaped temperature variation, albeit in some cases with different slopes, but the predictions generally deviate significantly from experiments.This suggests the QHA may simply be unable accurately to describe the B(T), which has also been concluded in similar studies [37,98].Two notable exceptions to the trend of similar-shaped predicted B(T) curves are BaO and Rb 2 O, for which the different functionals predict qualitatively different shapes with both negative and positive gradients at some temperatures.The former might be attributable to the relatively poor EoS fit discussed above (cf figure A16, supplementary material), although we saw no obvious issues in the EoS fits for Rb 2 O (figure A12, supplementary material), and the issues with the fits for ZnO again do not seem to have as large an impact on the calculations as they do for BaO.We note that the bulk modulus is derived from the curvature of the free-energy EoS (cf equation ( 11)), so it is likely to be inherently more sensitive to small errors than the V(T).
The tendency of the LDA to underestimate the lattice constants and overestimate the bulk moduli and phonon frequencies is well documented [99,100].Furthermore, the improved prediction of bulk moduli with the LDA relative to experiments has also been observed in other studies on MgO and similar materials [99].On the other hand, while Erba et al found that the predicted bulk moduli were sensitive to the choice of DFT functional, they also found that the LDA tended to overestimate measurements [37].
Moreover, we note that for some compounds we were unable to find experimental data to compare to, or were only able to find data at a single temperature, so our assessment of the performance of the different functionals is less robust than for some of the other QHA properties.With this in mind, it is not clear which of the functionals tested would give the best results beyond the materials examined in the present work.

Constant-pressure heat capacities.
The constantpressure heat capacity C p (T) can be obtained from the QHA and is arguably a more experimentally-relevant quantity than the constant-volume analogue C V (T). Figure 11 compares the C p (T) curves predicted for the Group 1 and 2 oxides with the five functionals to experimental data where available (plots for the Group 12 oxides can be found in figure A26 in the supplementary material).
Following the trends observed for other properties, we find that LDA tends to yield the lowest heat capacities and PBE the highest.The behaviour of the different functionals at low-tomoderate and high temperatures can be discussed separately.Over the former range the choice of functional does not appear to heavily influence the predicted C p (T), as was observed for the C V (T).Comparing the percentage differences to experiments at 300 K (table A1, supplementary material), the differences between the functionals are negligible and the average absolute errors are all less than 5%, suggesting that all five produce relatively accurate results at room temperature.In contrast, larger deviations between the functionals are seen at higher temperatures.Across all ten oxides the C p (T) curves predicted using the LDA and PBEsol remain closest to experimental measurements, but with larger relative errors.Despite this, the heat capacities appear to be predicted more accurately than the bulk moduli.
As for the volumetric expansion coefficient, in some calculations the C p (T) curves displayed unphysical oscillations at high temperature.We found this issue to be more prevalent in the calculated C p (T) than in the calculated α V (T), and since the heat capacities are computed from the second derivative of the Gibbs energies (cf equation ( 13)) we tentatively attribute this to amplification of numerical noise.For this reason, we only examined the C p (T) curve for SrO with SCAN up to T = 1200 K and the curves for BaO with PBE+D3 and SCAN up to 1100 and 600 K respectively.

Grüneisen parameter.
Finally, figure 12 compares the predicted temperature dependence of the thermodynamic Grüneisen parameters γ of MgO and CaO to experimental measurements.Data for all ten oxides is presented in figures A27-A29 respectively (supplementary material).
In practice, γ is a combination of several of the other QHA properties, viz. the volume, thermal expansion coefficient, bulk modulus and constant-volume heat capacity (cf equation ( 14)).The result is a comparable trend to that seen for many of the other QHA properties, with LDA predicting the lowest values and PBE the highest.
Although we were only able to find experimental data for MgO, CaO and SrO, the percentage differences listed in table A1 (supplementary material) suggest that LDA and PBEsol give the best predictions at 300 K, with absolute errors of 3.63 and 3.79% respectively.The other functionals also give relatively good predictions with average errors below 8%.However, the predictions become less accurate at elevated temperature, and when taking into account the entirety of the temperature ranges covered by the calculations LDA and PBE appear to give respectively the most and least accurate predictions.

Conclusions
This study has provided a quantitative assessment of the performance of quasi-harmonic lattice-dynamics calculations with five different exchange-correlation functionals for predicting the finite-temperature properties of ten binary metal oxides from Groups 1, 2 and 12 of the periodic table.Our findings offer an insight into the capabilities of the QHA and reveal some systematic trends in the performance of a selection of widely-used exchange-correlation functionals, making an important first step toward establishing the technical feasibility of HT QHA calculations.
Overall, our results indicate that both the HA and QHA are effective methods for predicting the finite-temperature properties of binary oxides.The harmonic phonon spectra, entropies and constant-volume heat capacities calculated at the equilibrium volume are in excellent agreement with experiments even at elevated temperature.Comparison to experimental measurements suggests the QHA can accurately predict the temperature dependence of the lattice volume at temperatures up to at least 2/3 the melting temperature.The constant-pressure heat capacities and Grüneisen parameters can also be computed with reasonable accuracy, albeit with some discrepancies at higher temperature.The bulk moduli can be estimated within a reasonable margin of error but the temperature dependence is not well described, as observed in previous studies [37,98], and the predicted thermal-expansion coefficients also show large deviations from experiment at higher temperature.
In the QHA, it is natural to attribute errors in predicted properties to higher-order anharmonic effects becoming increasingly prominent at elevated temperature.The central assumption of the QHA is that the HA remains valid at each of the volume expansions and compressions included in the calculations, which effectively correspond to different pressures and 'lattice temperatures'.This assumption may fail obviously in cases such as large expansions introducing imaginary modes, but may also fail more subtly in 'soft' systems where the frequencies remain real but the HA cannot accurately describe the potential-energy surface.We speculate that this may be the root cause of the generally poor predictions of the bulk moduli, which are inherently more sensitive to the free-energy EoS fits than the volume and Gibbs free energy, and of the poor predictions of properties such as the thermal expansion coefficients that are derivatives of other quantities.
Indeed, the percentage differences in the predicted expansion coefficients to experiments were found to be high even at 300 K, which has also been observed in other studies [92].The idea that the neglect of higher-order anharmonicity is the leading cause of error in the QHA calculations is also evident in the properties of some of the heavier oxides being less well predicted than their lighter counterparts.For example, the properties of CaO were generally less well predicted than those of MgO, which could be explained by the more polarisable cation in the former leading to 'softer' chemical bonding and stronger phonon anharmonicity.
Turning now to the relative performance of different exchange-correlation functionals, our results reveal some interesting systematic trends.Of the five functionals studied the predictions were invariably bounded by the LDA and PBE.The LDA tends to predict the lowest athermal equilibrium lattice constants/cell volumes, vibrational entropies, heat capacities, expansion coefficients and Grüneisen parameters and the highest phonon frequencies and bulk moduli.PBE generally shows the opposite behaviour, and PBEsol and PBE+D3 predict values intermediate between these extremes.While technical issues meant that we were unable to obtain SCAN results for all ten oxides, the results we were able to obtain suggest that this functional predicts harmonic phonon frequencies and thermal properties similar to the LDA and QHA properties intermediate between LDA and PBE.
Our results also highlight that some properties depend more strongly on the choice of exchange-correlation functional than others.All five of the functionals tested predicted very similar harmonic phonon spectra, with the largest variation seen in the high-frequency optic modes of the Group 1 and Group 12 oxides.The harmonic thermal properties were found to be largely independent of the functional, especially at lower temperature.Properties predicted using the QHA were generally more sensitive to the choice of functional, and particularly so at elevated temperature, with the volume, bulk modulus and thermal expansion coefficients most dependent on the exchange-correlation treatment.
When choosing an appropriate functional for HT QHA calculations, it is necessary to balance accuracy against computational cost, possibly taking into consideration the chemical nature of the materials being investigated.
Based on the findings from this study we would consider PBE and SCAN to be poor choices.PBE consistently yields the least accurate predictions compared to experiments, while SCAN is more computationally demanding than the LDA and GGA functionals and in our experience can be difficult to converge.The poor performance of PBE can be traced back to its well-known tendency to under-bind.This results in overpredicted equilibrium volumes and softer phonon frequencies, and the latter may lead to more anharmonic lattice dynamics that become particularly pronounced when calculations are performed on expanded volumes; assuming positive thermal expansion, this would be expected to produce anomalous behaviour at high temperature, which is evident in several of the calculations in this study.The requirement for dense grids to avoid spurious imaginary modes with SCAN can be explained by numerical noise being amplified in the second derivatives of the density, and is an issue we have previously observed with other meta-GGAs [28].
PBE+D3 gives better results than PBE, but in our opinion this is most likely because it corrects for its tendency to underbind.While the present study suggests that the dispersion correction improves upon PBE and gives good results even in systems with predominantly ionic bonding, as was also found by Reckien et al [102] there is no guarantee that it might not lead to unpredictable errors in other systems.
Given that the LDA and PBEsol also show better agreement with experiments, these likely represent better choices than PBE+D3.The simpler LDA functional tends to significantly underestimate lattice constants and equilibrium volumes, but gives the most accurate results for bulk moduli, volumetric expansion coefficients and Grüneisen parameters, and also has a marginally lower computational cost than the GGAs.On the other hand, PBEsol is optimised for solids and gives the most accurate predictions of the lattice constants and volumes together with relatively good predictions of other properties.Based on these considerations, either of the LDA or PBEsol could potentially be candidates for HT calculations.
However, the LDA is known to benefit from a significant cancellation errors [28] and can describe the properties of 'harder' materials like metal oxides well, but would not necessarily perform as well across systems with a wider range of chemistries.For this reason, we propose PBEsol as the best option for HT QHA calculations.

Figure 4 .
Figure 4. Constant-volume heat capacity as a function of temperature C V (T) for the equilibrium structure of MgO calculated using the LDA (purple), PBE (pink), PBE+D3 (yellow), PBEsol (blue) and SCAN (cyan) functionals, with experimental data overlaid as black markers for comparison [79-81].

Figure 5 .
Figure 5. Phonon dispersion curves of MgO from Γ to L calculated at a series of ± 2% expansions and compressions of the equilibrium volume V 0 using LDA (a), PBE (b), PBE+D3 (c), PBEsol (d) and SCAN (e).As noted in the text, the SCAN dispersions obtained for the 99% and 102% volume scaling have been omitted due to the presence of spurious imaginary modes, and a comparison including these volumes can be found in figure A8 (supplementary material).

Figure 7 .
Figure 7.Percentage errors in the 300 K unit-cell volumes of the ten binary oxides predicted using the five exchange-correlation functionals tested in this study.We note that data for some materials using SCAN are not included due to issues with these calculations (see text).The experimental data used for comparison are taken from references[61,62,64,68,[82][83][84][85][86][87][88].

Figure 10 .
Figure10.Percentage errors in the 300 K bulk moduli of the four Group 2 oxides predicted using the five exchange-correlation functionals tested in this study.Note that SCAN data for MgO and CaO are not included due to issues with these calculations (see text).The experimental data for comparison are taken from references[89,[93][94][95][96][97].

Figure 12 .
Figure 12.Thermodynamic Grüneisen parameter as a function of temperature γ(T) for MgO and CaO calculated using the quasi-harmonic approximation with the LDA (purple), PBE (pink), PBE+D3 (yellow) and PBEsol (blue) functionals.Experimental data are overlaid as black markers where available for comparison [91, 101].