Evaluating the harmonic approximation for the prediction of thermodynamic formation properties of solids

Thermodynamic stability is regularly quantified using the enthalpy of formation ( 𝛥 𝑓 𝐻 ) and the Gibbs free energy of formation ( 𝛥 𝑓 𝐺 ). Although knowledge of these properties is crucial for many applications, literature data are often missing for many solids. We evaluate the harmonic approximation (HA) for the prediction of 𝛥 𝑓 𝐻 and 𝛥 𝑓 𝐺 at constant volume and pressure, based on phonon calculations. Using density functional theory to carry out phonon calculations, we show that the HA excellently describes the temperature-dependence of 𝛥 𝑓 𝐻 and 𝛥 𝑓 𝐺 for 14 compounds; mean absolute error (MAE) of 2.1 kJ ⋅ mol − 1 and 1.1 kJ ⋅ mol − 1 , respectively, in the temperature interval 0–2500 K. Moreover, the performance of the HA was evaluated using computational data from the Materials Project database for 69 compounds yielding an MAE of 1.1 kJ ⋅ mol − 1 and 1.1 kJ ⋅ mol − 1 , respectively, in the temperature interval 100–800 K. Very good performance is also observed for a number of additional compounds, including several hydrated salts, at 298.15 K. The model is subsequently applied to a number of phase-transition phenomena that demonstrate the strengths and weaknesses of the model. In addition, it is demonstrated that the HA model can provide quantitative performance that rivals that of the conventional quasi-harmonic approximation for the prediction of formation properties, at a significantly reduced computational effort ( ∼ 5–10 times faster).


Introduction
The enthalpy of formation (  ) and the Gibbs free energy of formation (  ) are crucial properties for quantifying thermodynamic stability and equilibria in the context of chemical reactions.However, literature data for thermodynamic formation properties are incomplete or unavailable for many solids [1][2][3].As a result, being able to predict the value of these properties is highly relevant for many applications involving solids and electrolytes, including carbon capture and storage [4][5][6], pharmaceuticals [7,8], solar cells [9], fuel cells [10] and inorganic waste recycling [11].
Experiments for determining thermodynamic properties are challenging to carry out, especially across broad temperature and pressure ranges [12].Consequently, it is to be expected that experiments cannot match the rate at which new materials are discovered or proposed.Correspondingly, alternative approaches like empirical and semi-empirical models to predict   ,    and related properties of solids have been developed and used for several decades [13][14][15][16][17].However, many of these models have only been developed and tested for small groups of compounds or narrow temperature ranges, which limits their applicability.Primitive phonon models like the Einstein [18] and Debye [19] models can provide a good description of thermodynamic properties for elements [20], but they require accurate input parameters and they do not suffice for chemical compounds.Modifications have been introduced to extend the applicability and accuracy of these models [21,22], but they do not consistently provide chemical accuracy for the thermodynamic properties.Recently, approaches based on machine learning (ML) and data science have been shown capable of providing quite accurate predictions of   ( ) and   ( ) for large sets of compounds [23][24][25][26][27]. Nonetheless, these models are typically trained on experimental data and are, thus, inherently limited by the availability and accuracy of this data.In addition, these methods may only be applicable in certain temperature and pressure intervals.
In recent years, the growth of computer-power has enabled extensive development and application of electronic structure models and methods that are able to predict the temperature-and pressuredependence of thermodynamic properties of solids [28].The quasiharmonic approximation (QHA) is the most prominent of these methods for predicting isobaric thermodynamic properties.In the QHA, phonon https://doi.org/10.1016/j.commatsci.2023 calculations are carried out at several cell volumes around the equilibrium volume in the harmonic approximation (HA), followed by the application of an equation of state [29].The phonon calculations are typically performed using the finite displacement method (FDM) [30] or density functional perturbation theory (DFPT) [31].The QHA recovers a significant part of the anharmonicity of the real system at a moderate level of computational complexity compared to more exact treatments such as ab initio molecular dynamics [32][33][34][35], Monte Carlo simulation [36][37][38] or self-consistent phonon theory [39][40][41].Some developments on the QHA have successfully reduced its cost and complexity while maintaining decent accuracy [42][43][44][45].Even so, QHA calculations are often associated with considerable computational cost and/or complexity.As a result, only rarely are attempts made to predict the even more laborious formation properties using QHA [29,46].Moreover, extensive work has been carried out in recent years on first-principles, high-throughput lattice dynamics calculations in combination with materials databases, namely AFLOW [47], PhononDB [3] and Materials Project (MP) [48].Accordingly, temperature-dependent data for the constant volume thermodynamic properties calculated in the harmonic approximation (internal energy  , entropy , Helmholtz free energy  and isochoric heat capacity   ) are now available for thousands of solids.These properties are significantly more straightforward to calculate (e.g. the isochoric properties have been seen to be less sensitive to the exchange-correlation functional used [49]), and considerably less computationally demanding, than the commonly applied QHA calculations used to calculate the isobaric properties (enthalpy , Gibbs free energy , isobaric heat capacity   ).The AFLOW and PhononDB databases contain QHA calculations for numerous solids.Despite this, the calculations and methodologies around these databases were not (directly) intended for application to the calculation of the thermodynamic formation properties.That is, there are no first-principles-based frameworks for computationally inexpensive and reliable predictions of thermodynamic formation properties of solids; specifically, the temperature-dependence of said properties.Such a framework/methodology should permit high-throughput calculations of    and    and provide, even to non-experts, the ability to obtain fast and reliable predictions of the formation properties when little or no experimental data are available.
It is long-established that the specific heat capacity of solids at elevated temperatures is inadequately described by the harmonic model at constant volume [50].And, since thermal expansion is not considered in the harmonic model, properties like   ,  and  are not considered, conventionally.In this work, we evaluate the use of the harmonic approximation at constant volume for calculating  and , as a means to predict    and   .Since phonon calculations need only be carried out for a single structure with this approach, this significantly reduces computational time and complexity compared to the QHA.By comparison to tabulated values from the NIST-JANAF and NBS databases, we show that, despite its limitations, this model accurately predicts the temperature-dependence of    and    at 1 bar for several solids.

Theoretical methodology
In this section, we present the harmonic approximation that is used to treat the thermodynamic properties of solids.From the Legendre transforms of Gibbs we know that the Gibbs free energy may be written in terms of the Helmholtz free energy  as: where  and  are pressure and volume, respectively.We proceed, by assuming that the volume is constant with temperature, i.e. that there is no thermal expansion.In addition, contributions to the thermodynamic properties from temperature-dependent electron excitation and configurational effects are considered negligible.At temperatures below the melting point and for most systems/applications, this assumption is valid [12,51].There are examples of systems where configurational [52], electronic [53] and magnetic [54,55] contributions to the thermodynamic properties significantly contribute to thermodynamic stability, even at room-temperature, necessitating explicit calculation of these contributions.However, for high-throughput applications these effects are not generally significant enough to warrant the additional computational complexity required to model them.In summary, we consider a strictly harmonic solid.This renders the method unsuitable for studying parameters such as bulk moduli, thermal expansion coefficients or Grüneisen parameters which are proportional to the volume derivative Appendix B. Conversely, the thermodynamic properties are well-approximated by the harmonic term [56].Furthermore, when considering formation properties, we can expect some degree of error-cancellation, primarily of errors pertaining to the neglect of anharmonicity.For solids, the Helmholtz free energy in the HA is given by: where   is the zero-motion, electronic internal energy and  ℎ is the phonon Helmholtz energy.In the HA,  ℎ may be expressed in terms of the phonon frequencies via the canonical ensemble [56].We can then write Eq. ( 1) as: Similarly, the enthalpy becomes: ℎ ( ) and  ℎ ( ) are defined in Eqs.(A.3) and (A.1), respectively.At low pressures, the contribution from the  -term is small in a chemical context, e.g. for MgCO 3 ,  = 0.0028 kJ mol −1 at 1 bar, and this term is further reduced by cancellation with the solid constituent atoms when formation properties are considered.In this isochoric model,   must be equal to   .In addition, if we disregard  , we recover  and  from  and , respectively.Hence, Eqs.(3) and ( 4) are almost the same as the well-known expressions for the Helmholtz free energy and the internal energy, respectively, in solids [56].Despite its small magnitude, we include the  -term when calculating  and  as it does not require additional computational effort.A typical QHA calculation requires phonon calculations for ∼ 7 volumes beyond the equilibrium volume [29,49], consequently the HA will be ∼7 times faster, disregarding the additional time required for post-processing of the thermodynamic properties in the QHA, such as choosing and applying an equation of state [57].This represents a significant cost-reduction, especially considering how resource-intensive phonon calculations can be.At high pressures it will most likely be necessary to account for the effect of pressure on the equilibrium cell geometry in order to get accurate formation properties.A pressureconstrained geometry relaxation [58] could be employed in such a case.More sophisticated structure-searching methods beyond the HA may be required for complex cases where several competing phases exist [59,60], in order to determine the stable structure(s) at a given  and  .These approaches are not considered presently, as the focus of this work is kept on the prediction of the thermodynamic formation properties of stable structures at 1 bar.
The thermodynamic property of formation of a substance A n 1 B n 2 … refers to the change in a thermodynamic property associated with the reaction: where A, B, . . .are the constituent elements of A n 1 B n 2 … in their standard state.For a thermodynamic property , the thermodynamic property of formation (  ) is determined as: The enthalpy of formation and Gibbs free energy of formation are denoted    • and    • , respectively.The symbol • is used to denote that the respective property of formation is evaluated at 1 bar.
It is well-established that generalized gradient approximation (GGA) functionals struggle in predicting    for solids, since the electronic energy of formation     is especially inaccurate when the phases of the solid and its constituent species differ [24].Our intention is to demonstrate the performance of the HA in modeling the temperaturedependence of the thermodynamic formation properties.Thus, to eliminate the errors associated with the PBEsol-calculated     , we introduce a correction whereby the enthalpy is corrected by reference to an experimental value of the formation enthalpy at 298.15 K,    •  (298.15K): And the corrected Gibbs free energy is given by: The corresponding formation properties (i.e.   •  ( ) and    •  ( )) can then be derived from Eq. ( 6), using either elemental data derived from theory or from experiment/literature.In the former case, for a compound containing  elements,  + 1 phonon calculations are required to calculate the formation properties.The reference value    •  (298.15K) was chosen because it is reported for many substances.However, if no experimental data for this quantity exists, alternative sources for the reference value may be necessary.In recent years, a number of methods have been developed to predict    that significantly improve upon the local density approximation (LDA) and GGA functionals.These include: correction scheme methods [24,[61][62][63], new exchange correlation functionals [25,[64][65][66][67] and machine learning assisted methods [68,69].Inclusion of these methods is outside the scope of this paper, but enthalpies predicted using these approaches can straightforwardly be applied instead of the experimental reference value which is used here.Moreover, if    •  (298.15K) is available, it could be used to obtain a correction to  • (298.15K) instead of using    •  (298.15K), and/or the correction may be ''centered'' at a different temperature if experimental/theoretical values of the thermodynamic formation properties at other temperatures are available.
Prior to the phonon calculations, the unit cell geometry is relaxed.Input structures are given by experimentally determined structures from literature of the species in their stable state and structure in the relevant temperature range and at 1 bar, see Table S1.Magnetization is assumed to be negligible for the ground-state structures considered in this paper, we expect this to be a good assumption.Even if a magnetic transition is expected to occur, its contribution to the thermodynamic formation properties is expected to be very small and this is thus not considered.A convergence threshold of 10 −9 a.u.(atomic units) was used for the energy, in all calculations, and, in addition, a threshold of 10 −5 a.u. and 0.5 ⋅ 10 −3 kbar was used for the forces and pressure, respectively, in the relaxation calculations.No pressure is employed in the relaxation procedure and its inclusion in the relaxation is assumed to have a negligible effect on the thermodynamic properties at 1 bar.The Monkhorst-Pack scheme [74] was applied for Brillouin zone sampling.The grid density was set to approximately 1500-2000 points per reciprocal atom.This density was chosen based on its reliable performance in high-throughput phonon calculations [48,49].For the plane-wave kinetic energy cutoff and corresponding charge density cutoff, the values suggested in the SSSP library are used.
For metals and semi-metals, Methfessel-Paxton smearing [75] is used.The default value suggested by the QE input generator [73] is used for the smearing width ().  is found by the formula [76]: For calculations of molecular gases, each calculation is carried out at the  -point, the molecule is centered in a cubic cell with side-lengths = 30 Bohr and the Martyna-Tuckerman correction [77] is applied to the total energy.
Phonon calculations of the solids were carried out on equilibrium structures with the finite displacement method, in the HA, using QE and Phonopy [56].A 2 × 2 × 2 supercell is used for all the phonon calculations.This was chosen as it has been shown in literature to be a robust choice for a number of systems [49,78].
For the gaseous molecules, the phonon module of QE is used for calculation of vibrational frequencies.A rigid-rotor harmonic oscillator approximation was assumed and the thermodynamic properties were determined in the ideal-gas approximation using the Atomic Simulation Environment [79].
All thermodynamic quantities were calculated in steps of 2 K.At very low temperatures (≤10 K) the step-size was increased for some systems to avoid numerical instabilities.All thermodynamic formation properties considered in this work are evaluated at a pressure of 1 bar.The temperature interval considered for each compound is determined by experimental data for phase change temperatures of all formation reaction species at 1 bar (Table S1).Only the stable phase at 298.15 K and 1 bar is considered for all species.The effect of pressure on the unit cell's equilibrium structure was not considered.The reference value    See definitions of the mean absolute error (MAE) and mean absolute relative error (MARE) in the supporting material.

Results and discussion
Thermodynamic properties for the elements have been determined from data tabulated in the NIST-JANAF database [2].This approach was chosen based on property comparisons of DFT-based calculations of the elemental thermodynamic properties to corresponding calculations based on experimental elemental data (see Fig. S1).The two approaches yield quite similar results; their performance is briefly compared and covered in the discussion.Because the use of DFT-computed elements is limited by the phase of the elemental constituents and the additional computational effort required to obtain this data, the use of experimental data for the elements (which are widely available) was adopted in the present paper in order to extend the data-sets.As a result, for compounds, differences of the type The former notation is used for the sake of consistency and ease of interpretation.

𝛥 𝑓 𝐻 • (𝑇 ): Comparison to the conventional approach
Data pertaining to the formation enthalpies calculated via Eq.( 7) are visualized in Fig. 1a for 14 compounds.Since the reference value used is the experimental enthalpy of formation at 298.15 K, the presented data reflect how well the temperature-dependence of    • is modeled by the constant-volume HA.All the data have an error associated with the approximations of the DFT electronic structure calculation used.The PBEsol functional has shown good performance for the entropy [80] and geometry [81] of solids, which has motivated its use in this paper.The magnitude of the absolute errors roughly correlate with system size which is expected behavior as the magnitude of the property increases accordingly.Another expected behavior which can be observed is the systematic underestimation of    • with temperature as the effects of thermal expansion (and other temperature-dependent effects) become more pronounced; a result of the overestimation of phonon-frequencies associated with assuming a harmonic potential.Since  is only weakly temperature-dependent, the errors associated with modeling its temperature-dependence are generally quite low.Despite this, the errors obtained with the HA are still significantly lower than if    • ( ) is assumed to be constant with temperature (Fig. 1b) as is often the case for estimations of    • in literature.The data confirm the consensus that this assumption is good in the intermediate temperature-interval around the    •  (298.15K) reference used here, but also that it becomes quite erroneous at more extreme temperatures.Moreover, it should be noted that if/when the elemental constituents undergo phase changes (which were avoided in the present data) the error associated with this assumption can be significant if left unaccounted for.From Fig. 1 it can be concluded that, for this data-set, both approaches perform within chemical accuracy (i.e.MAE ≤ 4.18 kJ mol −1 ) for    • ( ) but carrying out HAcalculations provides a notable improvement over assuming a constant value for    • ( ).

𝛥 𝑓 𝐻 • (𝑇 ): Extending the temperature range
Unlike Figs.1a, 2a contains a comparison of    •  ( ) to    •  ( ) in the entire temperature-interval of the compounds in the structure and state at which they are stable at 298.15 K and 1 bar.Here, the MAE (=2.12 kJ mol −1 ) is still well within chemical accuracy.For MgO and ZrC, which have melting points ≥ 2500 K, the errors become noticeably large at the highest temperatures.Clearly, some care should be taken if such high temperatures are to be considered in the HA framework, emphasizing that the model is most useful at lower temperatures (e.g.≤ 1500 K).

𝛥 𝑓 𝐻 • (𝑇 ): Extending the compound data-set
In order to investigate how the HA performs for a greater variety of crystal chemistries, we have leveraged the thermodynamic property data of the Materials Project (MP) database [48] (Fig. 3).MP contains constant-volume thermodynamic properties computed using PBEsol DFT and DFPT [48] for more than 1000 compounds in the temperatureinterval 5 K to 800 K. Since their methodology is similar to that used here, we expect similar performance and this is supported by comparison of the MP entropies to our entropies (see Fig. S2).Fig. 3a shows a parity plot of the corresponding simulated and experimental enthalpy for 69 compounds that are present in both the MP thermodynamics database and the NIST-JANAF database.Although the temperatures are not as extreme as those seen in the previous data, the agreement with experiment is nonetheless excellent (MAE = 1.05 kJ mol −1 , MARE = 0.29%, r 2 = 1.00).From Fig. 3b we can see the familiar increase of error with temperature.However, the center of the error remains low even at 800 K.This demonstrates well the usefulness and applicability of the HA thermodynamic properties in high-throughput applications.

𝐺(𝑇 ): Comparing theoretical calculations
In order to evaluate how well the constant volume HA performs in comparison to other theoretical approaches, the temperature dependence of ( ) relative to experiment has been determined from the QHA data of PhononDB [3] and from the quasi-harmonic Debye data of AFLOWs Automatic Gibbs Library (AGL) [82] for a number of compounds.The AGL model is a computationally cheaper alternative to the QHA which enables high-throughput calculations of isobaric thermodynamic properties.The results are presented in Fig. 4. The HA data is shown in Fig. 4a.The temperature-dependence of the HA error of ( ) has a distinct shape: underestimation at low temperatures and overestimation at high temperatures.This can be attributed to the contributions of the  and −  terms to the Gibbs function .As seen previously,  tends to be systematically underestimated in the HA (Fig. 2a); a negative contribution to  relative to the true system.With increasing temperature,  also tends to be underestimated [83].Consequently, with increasing temperature, the −  term imparts a positive contribution to  relative to the true system; becoming increasingly dominant as the temperature increases.This behavior is fortunate for the quantitative determination of ( ), since it results in increased tendency for error cancellation; especially at lower temperatures.The result, evident in Fig. 4a, is very low errors in the HA-determined temperature-dependence of ( ) for the investigated data-set.In comparison to the HA, the AGL (Fig. 4b) and QHA (Fig. 4c) exhibit notably different behavior for the temperature-dependence of ( ).Here, the errors behave more linearly, an indication that the qualitative description of the temperature-dependence is more apt.However, the quantitative performance is significantly worse than in the HA.Presumably, this can, in part, be explained by the lessening of error-cancellation (as previously discussed) leading to increased exposition of the errors associated with the shortcomings both in these models and the underlying electronic structure calculations.An investigation of where and how exactly errors in the thermodynamic properties arise across models and approximation is very relevant but outside the scope of the present study.It should be noted that limitations of the electronic structure calculations used in the AGL and PhononDB datasets also contribute to the error.Overall, Fig. 4 helps highlight key strengths and weaknesses of the HA model. ( ) to    •  ( ) in the entire temperature-interval of 14 compounds in the structure and state at which they are stable at 298.15 K and 1 bar.The trend and performance of the HA for    • ( ) here resemble that seen in Fig. 4 for  • ( ).This is understandable, since the thermodynamic properties of the constituent elements come from experimental data.This can be contrasted to calculations in which these data come from DFT in which case the qualitative description of the formation properties (interestingly) Fig. 4. Errors in temperature dependence of the Gibbs free energy as a function of temperature ( ), evaluated as ( ) = (( ) − (0 K)) − (  ( ) −   (0 K)) for 9 solid compounds.a ( ) evaluated in the harmonic approximation (this paper), b ( ) evaluated using data from AFLOWs Automatic Gibbs Library (AGL) [82], c ( ) evaluated using QHA data from the PhononDB [3].Experimental data from NIST-JANAF [2].The data-points for each compound have been offset based on temperature as a visual aid.Temperature step = 100 K. Data at  = 0 K are excluded from the figure data.
appears to be better, presumably owing to error-cancellation between the compound and its constituent elements (see Fig. S1).The nature and importance of this behavior for further application of the formation thermodynamic properties from the HA is not entirely clear and this should be considered in future works on this topic.Nonetheless, the quantitative performance of the present implementation is very good and the MAE of    •  ( ) is even lower (1.07 kJ mol −1 ) than that of    •  ( ) for this data-set.For comparison, the ML-based   ( )descriptor proposed by Bartel et al. [23], which requires only the compound structure and an estimate of the formation energy/enthalpy, yields a MAE = 8.27 kJ mol −1 for the same compound-set in the temperature-interval 300-2000 K when    •  (298.15K) is used as the reference formation energy (see Fig. S3).For the same temperatureinterval, the HA-model yields MAE = 1.03 kJ mol −1 .Furthermore, analogously to Fig. 1b, we have evaluated how well    •  (298.15)predicts    •  ( ) (see Fig. S4).Clearly, at low temperatures this is a decent assumption but the error increases rapidly with temperature for most of the systems (MAE = 53.4kJ mol −1 ), showing the importance of taking the entropy into account.
Table 1 shows   -data at 298.15 K, computed using reference data from the NBS tables.This enables    •   to be calculated for some additional compounds that are not present in NIST-JANAF.Among these are four MgCl 2 hydrates: the mono-, di-, tetra-and hexahydrate.The agreement with experiment is generally good.The hydrates exhibit somewhat larger absolute errors which can, in part, be explained by the presence of more atoms and increased complexity of these compounds.However, the errors could also be related to the presence of water which may be inadequately described by the PBEsol functional [84] or due to the increased amount of gas-phase constituent species which may affect cancellation of errors.For now, more work is needed to ascertain the exact reason for this behavior and whether this applies generally.
Another hydrate of interest is the dihydrate of sodium chloride, NaCl⋅2H 2 O.In a 1992 paper [85], through a literature review of data related to the NaCl/water system, D.G.  1.The simulated data slightly overestimate Archer's estimates.Nonetheless, the agreement is good, which serves to mutually verify the results, providing further reliability for research pertaining to NaCl⋅2H 2 Ocontaining systems.Overall, Table 1 echoes the good performance of the HA method seen in the prior figures.

𝛥 𝑓 𝐺 • (𝑇 ): Extending the compound data-set
Fig. 3c shows a parity plot of the corresponding simulated and experimental Gibbs free energy for 69 compounds that are present in both the MP thermodynamics database and the NIST-JANAF database.Like seen with the formation enthalpy, the agreement with experiment is again excellent (MAE = 1.05 kJ mol −1 , MARE = 14.0%, r 2 = 1.00).The relatively large MARE is due to the low absolute    • ( ) values of a subset of the compounds.As observed for the corresponding formation enthalpy data, Fig. 3d shows a slight but distinct increase in the error with temperature which nonetheless remains low even at 800 K.This shows that the Gibbs free energy is similarly well predicted for this larger set of compound chemistries.

Case studies
In recent decades, bischofite and other MgCl 2 hydrates have received attention for their heat storage capabilities which make them attractive as depositories for excess energy produced by dynamic and irregular power sources such as solar energy [87][88][89].Knowledge of the possible reactions which the hydrates may undergo, and the conditions at which they are favorable, is important for designing thermal energy storage systems.Fig. 5 shows the equilibrium temperature for the dehydration of MgCl 2 ⋅6H 2 O at a given partial pressure of water, simulated via the HA (Eq.( 8)), see also Fig. S5 and S6 for the dehydration reactions of MgCl 2 ⋅4H 2 O and MgCl 2 ⋅2H 2 O, respectively.This is determined by the temperature at which the Gibbs free energy of the dehydration reaction is equal to zero.This is straightforward to find once the   values of the reaction species are known.The simulated data agrees quite well with the range of temperatures/vapor pressures and slope found by Smeets et al. [90] in a DFT study where the authors modeled the hydrate as a gas-phase molecule.There is also decent agreement   with the estimates based on experiments proposed by Kipouros et al. [91].Transition temperatures proposed by others [92][93][94] also coincide with the simulated data.This illustrates the strength of the HA model beyond its ability to predict thermodynamic formation properties.
Operating within the harmonic approximation does impose some limits on the phenomena which the model can describe.This is demonstrated in Fig. 6, for the  ↔  transition of SiO 2 (quartz), a geologically [96] and technologically [97] important material.Note that no correction is applied in this case because of a lack of room-temperature experimental data for -quartz; this should not affect the conclusion since we are considering a relative quantity between two very similar solids.It is well-documented that the  ↔  transition occurs at ≈ 846 K [98,99], but evidently the HA model fails to predict this.It has been shown that the transition is governed by the effects of thermal expansion [100,101].And this has been shown to be true of other polymorphic transitions as well [102].Naturally, this makes the HA unsuited for modeling this transition since it does not account for thermal expansion.Accordingly, for some phenomena, caution should be exercised when applying the harmonic model [103].Especially at temperatures near the melting point, methods beyond the HA and QHA that are able to recover the anharmonic potential [104][105][106] are likely required for adequate prediction of the stable phase.

Conclusion
We have evaluated the performance of the isochoric, harmonic approximation model for the prediction of thermodynamic properties of solids.In its current form the model is computationally significantly more economical than current, conventional methods and by optimizing the input parameters and/or employing other ab initio or semi-empirical methods, the calculations could be sped up beyond what has been demonstrated here, enabling fast and straightforward calculations even for relatively large and/or complex compounds such as molecular crystals.When applied to the prediction of thermodynamic formation properties, including enthalpy of formation and Gibbs free energy of formation, the model accurately predicts their temperature dependence for a variety of solids.This was substantiated by evaluating the formation properties of 69 compounds of varied chemistry using data sourced from the Materials Project database which showed the HAs excellent performance compared to experiment in the 100-800 K temperature range.Furthermore, we have shown that the harmonic approximation is able to compete with the conventional quasi-harmonic approximation for the prediction of thermodynamic formation properties as well as state-of-the-art algorithms based on machine learning.More thorough investigations of larger datasets (higher temperatures/pressures and more, varied compounds) are required to establish the applicability of the model.Nonetheless, while the harmonic model exhibits limitations in some aspects, its speed and relative accuracy makes it a promising candidate for use in high-throughput predictions of thermodynamic formation properties of solids.

Fig. 1 .
Fig. 1. a Error of the simulated standard enthalpy of formation    •  ( ) as a function of temperature compared to experimental data    •  ( ) for 14 solid compounds, b Error of the experimental enthalpy of formation at 298.15 K,    •  (298.15K), compared to experimental data as a function of temperature    •  ( ) for 14 compounds.Only the state in which the compound and constituent elements are stable at 298.15 K and 1 bar are considered.Experimental data from NIST-JANAF [2].The data-points for each compound have been offset based on temperature as a visual aid.Temperature step = 100 K. Data at  = 300 K are excluded from the figure data.

Fig. 2 .
Fig. 2. a Error of the simulated standard enthalpy of formation as a function of temperature    •  ( ) compared to experimental data    •  ( ) for 14 solid compounds (analogous to 1a with an extended temperature-range), b Error of the simulated Gibbs free energy of formation    •  ( ) (different to Fig. 1b) as a function of temperature compared to experiment    •  ( ) for 14 compounds.Experimental data from NIST-JANAF [2].The data-points for each compound have been offset based on temperature as a visual aid.Temperature step = 100 K.For    • , data at  = 300 K are excluded from the figure data.

Fig. 3 .
Fig.3.Simulated thermodynamic properties using data from Materials Project[48] for 69 species (see TableS4) compared to experimental data, as a function of temperature.a Parity plot of the simulated standard enthalpy of formation    •  ( ) and the experimental enthalpy of formation    •  ( ) as a function of temperature, b kernel density plot of the absolute error in the MP-simulated enthalpy, as a function of temperature, c Parity plot of the simulated standard Gibbs free energy of formation    •  ( ) and the experimental Gibbs free energy of formation    •  ( ) as a function of temperature, d kernel density plot of the absolute error in the MP-simulated Gibbs free energy, as a function of temperature.Experimental data from NIST-JANAF [2].Temperature step = 100 K.For    • , data at  = 300 K are excluded from the figure data.
Fig.3.Simulated thermodynamic properties using data from Materials Project[48] for 69 species (see TableS4) compared to experimental data, as a function of temperature.a Parity plot of the simulated standard enthalpy of formation    •  ( ) and the experimental enthalpy of formation    •  ( ) as a function of temperature, b kernel density plot of the absolute error in the MP-simulated enthalpy, as a function of temperature, c Parity plot of the simulated standard Gibbs free energy of formation    •  ( ) and the experimental Gibbs free energy of formation    •  ( ) as a function of temperature, d kernel density plot of the absolute error in the MP-simulated Gibbs free energy, as a function of temperature.Experimental data from NIST-JANAF [2].Temperature step = 100 K.For    • , data at  = 300 K are excluded from the figure data.

3. 5 .
Fig. 2b contains a comparison of    • ( ) to    •  ( ) in the entire temperature-interval of 14 compounds in the structure and state at which they are stable at 298.15 K and 1 bar.The trend and performance of the HA for    • ( ) here resemble that seen in Fig.4for  • ( ).This is understandable, since the thermodynamic properties of the constituent elements come from experimental data.This can be contrasted to calculations in which these data come from DFT in which case the qualitative description of the formation properties (interestingly)

Fig. 5 .
Fig. 5. Simulated equilibrium temperature at different water partial pressures p H 2 O compared to literature [90,91] for the dehydration reaction of MgCl 2 ⋅6H 2 O. Reference  •  (298.15K) values from the NBS tables [86].To construct this figure, data points from figures in literature have been processed using the Webplotdigitizer tool [95].

Fig. 6 .
Fig. 6.Calculated reaction Gibbs free energy    •  for the quartz  ↔  polymorphic transition as a function of temperature.

Table 1
Archer introduced new equations for thermodynamic properties of NaCl⋅2H 2 O, reporting estimates of    •  (298.15K) and    •  (298.15K) at 298.15 K.    •  (298.15Gibbs free energy of formation    • in kJ mol −1 from simulation using Eq.(8) (sim) and from experiment (exp), including the difference between the two, at 298.15 K. Experimental and reference data taken from the NBS tables [86] unless otherwise stated.