Abstract
The quasiharmonic approximation (QHA) in combination with density-functional theory is the main computational method used to calculate thermodynamic properties under arbitrary temperature and pressure conditions. QHA can predict thermodynamic phase diagrams, elastic properties and temperature- and pressure-dependent equilibrium geometries, all of which are important in various fields of knowledge. The main drawbacks of QHA are that it makes spurious predictions for the volume and other properties in the high temperature limit due to its approximate treatment of anharmonicity, and that it is unable to model dynamically stabilized structures. In this work, we propose an extension to QHA that fixes these problems. Our approach is based on four ingredients: (i) the calculation of the n-th order force constants using randomly displaced configurations and regularized regression, (ii) the calculation of temperature-dependent effective harmonic frequencies ω(V, T) within the self-consistent harmonic approximation (SCHA), (iii) Allen’s quasiparticle (QP) theory, which allows the calculation of the anharmonic entropy from the effective frequencies, and (iv) a simple Debye-like numerical model that enables the calculation of all other thermodynamic properties from the QP entropies. The proposed method is conceptually simple, with a computational complexity similar to QHA but requiring more supercell calculations. It allows incorporating anharmonic effects to any order. The predictions of the new method coincide with QHA in the low-temperature limit and eliminate the QHA blowout at high temperature, recovering the experimentally observed behavior of all thermodynamic properties tested. The performance of the new method is demonstrated by calculating the thermodynamic properties of geologically relevant minerals MgO and CaO. In addition, using cubic SrTiO3 as an example, we show that, unlike QHA, our method can also predict thermodynamic properties of dynamically stabilized phases. We expect this new method to be an important tool in geochemistry and materials discovery.
Similar content being viewed by others
Introduction
The accurate calculation of thermodynamic properties of solids is important in many fields of science and technology. Computational methods are used whenever determining these properties experimentally is expensive, impractical, or unfeasible. The prediction of thermodynamic properties is relevant in geochemistry and geophysics, where they are used to interpret seismic data, inform models of the Earth’s interior, and make predictions about its composition and behavior1,2,3,4,5. The inherent complexity involved in experimentally reproducing the conditions that prevail in a planetary interior, as well as the relative scarcity of information that can be gleaned from the analysis of meteorites and seismic waves, make the accurate prediction of thermodynamic properties a valuable tool in this field. Likewise, the accurate computational prediction of thermodynamic properties is important in the design of new materials, in the development of databases, which enable the high-throughput screening for materials with desirable properties6,7,8,9, and in materials modeling under extreme conditions10. Calculable thermodynamic properties of solids include the Gibbs energy, which gives access to the thermodynamic phase diagram, and the temperature- and pressure- dependence of the density (the equation of state), elastic moduli, and other properties of technological and fundamental relevance such as thermal conductivities.
The workhorse for the computational prediction of thermodynamic properties in solids is density-functional theory (DFT), which provides a convenient balance between accuracy and computational speed5,11. Unlike experimental methods, the application of an arbitrary pressure or stress on a material within DFT is straightforward, whereas incorporating temperature effects is much more challenging. The effect of temperature on the thermodynamic properties of a solid manifests primarily through atomic vibrations. The most straightforward method for incorporating the vibrational contribution to the thermodynamic properties is the harmonic approximation12,13,14,15,16,17 (HA), which assumes atoms vibrate under a harmonic potential. In the HA, atomic forces depend linearly on the displacements via the second-order force constants (FC2)13,14,15:
where E is the system energy and \({u}_{i}^{\alpha }\) is the displacement of atom i in the α Cartesian direction. Second-order force constants, as well as phonon frequencies and normal modes, are routinely calculated in DFT, either via the direct (also known as frozen-phonon) method17,18,19 or using perturbation theory20.
In the HA, crystal vibrations behave as an ideal Bose gas, and consequently the harmonic vibrational contribution to the thermodynamic properties of a material can be calculated exactly. However, within the HA it is not possible to predict any thermodynamic property that involves the volume-dependence of the vibrational frequencies (e.g., the thermal expansion coefficient) or transport properties such as thermal conductivities17,21. The most popular approach for going beyond the HA is the quasi-harmonic approximation (QHA), in which anharmonicity is partially incorporated by making the harmonic vibrational frequencies depend on the volume explicitly (ω(V))5,16,22,23,24,25,26. In a typical QHA calculation, a grid of volumes is sampled and constant-volume relaxations are carried out, followed by harmonic frequency calculations. Standard software packages are available for obtaining the QHA thermodynamic properties from the energy and phonon data on the volume grid24,25,27,28,29,30.
The quasi-harmonic approximation has had great success in the prediction of material properties under various pressure and temperature conditions31,32,33,34. The main drawbacks of QHA are that it produces spurious results at sufficiently high temperature, particularly at low pressure22,31,32,35, and that it is unable to model dynamically stabilized phases. The crux of the problem is that the temperature dependence of the vibrational frequencies in QHA is included indirectly through the volume expansion (ω(V(T))), which is correct only to first order36. The QHA validity temperature limit depends on the particular system and on pressure, but it is expected that systems with soft (low-frequency) phonons or light atoms are particularly affected37. Correcting the blowout in the QHA thermodynamic predictions at high temperature requires a proper description of anharmonicity, which necessitates the inclusion of the higher-order terms in the vibrational Hamiltonian17,36,38. It also requires a method for deriving the thermodynamic properties, since the statistical treatment beyond the HA no longer has an exact solution. One popular way of overcoming the limitations of QHA is using molecular dynamics (MD) simulations at classical or ab initio (AIMD) level39,40,41,42,43,44,45,46, but carrying out accurate MD simulations is complex and requires a substantial investment of time and effort (although machine learning methods can be used to alleviate part of the computational cost).
In this article, we propose a computationally efficient alternative method for correcting the spurious behavior of QHA at high temperature based on Allen’s quasiparticle (QP) theory36,38. In this theory, the temperature-dependent frequencies (ω(V, T)) are known, and the entropy is approximated by assuming a system of non-interacting bosonic quasiparticles with energy levels given by those frequencies. QP theory is combined with n-th order force constant matrices (FCn) determined from DFT forces using machine learning regression techniques47, from where effective temperature-dependent harmonic frequencies are obtained via a self-consistent harmonic approach47,48. All properties other than the entropy are calculated from QP in a thermodynamically consistent way using a numerical method based on the Debye model that simplifies the practical implementation of the QP theory. The proposed method has a computational complexity similar to QHA but requires more supercell calculations. It incorporates the effect of the n-th order force constants, and coincides with the QHA results at low to moderate temperatures by construction. In addition, the FCn determined during the procedure enable the calculation of other properties inaccessible or difficult to model accurately in plain QHA such as the thermal conductivity, temperature-dependent frequencies, and temperature-induced second-order phase transitions. With two illustrative examples (MgO and CaO), we show that our implementation of QP theory prevents the QHA blowout and substantially improves the accuracy of the predicted thermodynamic properties at high temperature. Another example (cubic SrTiO3) shows how our new method can be used to model dynamically stabilized structures, which present imaginary frequencies at zero temperature. We believe this method is an important step towards the easy calculation of free energies, thermal conductivities, elastic constants, and other properties in large systems under arbitrary pressure and temperature.
Results
Magnesium oxide
To illustrate our new method, QHA and QP calculations were run on MgO using the calculation parameters indicated in “Computational details”. We used a volume grid containing 18 points between 12.26 Å3 and 24.29 Å3 per formula unit. At each volume the SCF calculation was carried out followed by a DFPT calculation for the Born charges and the dielectric tensor. Then, a single structure with randomly displaced atoms was calculated to generate the harmonic FC2, which was then used to generate 10 more structures by phonon rattle “Calculation of the force constants”. The calculated forces at these rattled structures were used to construct the FCn, followed by the generation of the SQP at 13 temperatures between 50 K and 3000 K. The resulting entropies were used to calculate the Debye model parameters “Calculation of thermodynamic properties”, and the thermodynamic properties were computed with gibbs2.
Figure 1 shows a selection of four calculated thermodynamic properties in MgO as a function of temperature at zero pressure: volume per formula unit (V), adiabatic bulk modulus (BS), constant-pressure heat capacity (Cp), and thermal expansion coefficient (α). The figure compares the QP and QHA results against reported experimental and calculated data. The plot shows clearly the spurious behavior of QHA at high temperature. In agreement with the literature31,36, the volumes predicted by QHA start to diverge at around 1000 K and the divergence is so large at 1500 K that the equilibrium volume exceeds the volume grid used for the calculation (and therefore cannot be predicted reliably and is not represented). The QHA divergence at high temperature extends to all other calculated thermodynamic properties as well. As predicted by Allen36,38, quasiparticle theory corrects the spurious QHA behavior, and all thermodynamic properties show trends that are consistent with experimental behavior at high temperature, up to the melting point. In the case of α and Cp, the calculated data are very close to the experimental values; for V and BS the disagreement between experimental and calculated results may be ascribed to the shortcomings of the exchange-correlation functional (PBE), as shown below. As expected from the construction of our method, QP and QHA coincide at low temperature. Figure 1 also compares the calculated α with the results obtained by Kwon et al.49, who used CSLD to obtain up to fourth-order force constants followed by an anharmonic free energy calculation using the improved self-consistent phonon (ISCPH) theory50,51, that incorporates explicit terms for the third-order and fourth-order anharmonic contributions. The agreement with Kwon et al.’s thermal expansion data is remarkable. The mismatch in the calculated volumes in the QHA and anharmonic results can be attributed to the different functionals used (PBEsol in Kwon et al., PBE in this work).
The plots show the volume (a), adiabatic bulk modulus (b), constant-pressure heat capacity (c), and volumetric thermal expansion coefficient (d) as a function of temperature using QHA (dashed lines) and QP (full lines) theories. Experimental data95,96,113,114,115 are shown as discrete points (circles and triangles). Our QHA/QP results are shown in red and the results of Kwon et al.49 using compressive sensing lattice dynamics (CSLD)44,92 and an anharmonic free energy calculation with the improved self-consistent phonon (ISCPH) theory50,51 are shown in blue.
Figure 2a shows the 300 K and 3000 K isotherms calculated with QHA and QP, using the PBE functional. As before, our implementation of QP theory coincides with QHA at low temperature, which makes the 300 K isotherms essentially indistinguishable. However, the 3000 K isotherm shows significant deviations between both theories below 50 GPa. At low pressure, QHA predicts equilibrium volumes so large that they exceed our calculation grid (and therefore are not represented). QP corrects this behavior and restores the correct low-pressure behavior of the high-temperature isotherm. The predicted volumes are slightly higher than experiment but this is again caused by shortcomings of the exchange-correlation functional. The QHA and anharmonic results from Kwon et al.49 (using ISCPH) are similar, although shifted to lower volumes due to their use of PBEsol instead of PBE.
a Volume as a function of pressure calculated with quasiharmonic (QHA, dashed/dotted lines) and quasiparticle (QP, full lines) theories at two temperatures, compared against experimental data for the 300 K isotherm (Duffy et al.93 and Jacobsen et al.116), and the 3000 K isotherm (Tange et al.97) shown as discrete points, and compared also with the results reported by Kwon et al.49 (purple lines). b Volume as a function of temperature using various functionals: PBE (red, blue)84, PBEsol (black)117, B86bPBE (green, orange)84,118, with and without the exchange-hole dipole moment (XDM) model dispersion correction119,120,121. Experimental data are shown for comparison as discrete points (circles and triangles)95,113.
Figure 2b shows the calculated volume as a function of temperature at zero pressure using various functionals, compared against experimental data. There are significant deviations between the predictions of different functionals both in the equilibrium zero-temperature volume as well as at higher temperatures. As noted in the literature, the inclusion of dispersion interactions unexpectedly seems to improve agreement with experiment significantly in ionic solids such as MgO52. The figure evidences that using QP instead of QHA makes the choice of functional the leading contribution to error at high temperature.
Another advantage of using QP for the calculation of thermodynamic properties is that the FCn can be used for obtaining quantities that are not accessible within the harmonic theory such as, for instance, the vibrational (lattice) contribution to the thermal conductivity (κl)26,53,54. The thermal conductivity is a key physical property used to understand the behavior of the Earth’s mantle55,56,57,58. It can be obtained by solving the Boltzmann transport equation (BTE) using various models59,60,61,62,63,64,65. In this work, we solved the linearized phonon Boltzmann equation including two- and three-phonon scattering processes63,64,65 for MgO, which assumes a small temperature gradient and requires the FC2 and FC3. Both are extracted from the FCn model used in the thermodynamic QP calculation at no additional cost. The linearized BTE is then solved using the ShengBTE code65 to calculate κl, with a converged 15 × 15 × 15 reciprocal-space-point grid, and a “scalebroad” parameter of 0.5. The BTE is solved iteratively until a convergence of 1 × 10−6 W m−1 K−1 is obtained in κl. The Born effective charges and the dielectric tensor, also obtained in the QP thermodynamic calculation, were used to introduce long range corrections65.
Figure 3 shows the calculated κl using QHA and QP (with the PBE functional), and without considering the effect of thermal expansion, as a function of temperature at zero pressure. In both cases, the κl is determined from the FCn obtained with hiphive, but the volume dependence with temperature is determined with either of the two approximations, which indirectly affects the predicted temperature dependence of κ. Even at low temperatures, failing to take into account the effect of thermal expansion on the thermal conductivity results in a severe overestimation of the experimental values. At low temperature, the QHA and QP results agree with each other and with the experimental results (Hofmeister et al.66 reported κl for an MgO single-crystal and Slifka et al.67 for a polycrystalline sample, hence the disagreement). At high temperature, however, the QHA conductivity drops significantly and, because of the spurious volume expansion, it cannot be calculated above 1700 K. In contrast, the QP result is very close to the experimental data even at temperatures close to the melting point, and similar to the results obtained by Kwon et al.49. As in the case of the other thermodynamic properties, the disagreement with the experimental κl may be ascribed to errors coming from our choice of exchange-correlation functional. Figure S4 in the SI shows that different functionals display different behaviors in the calculated κl as a function of temperature. Although none of them diverge, the choice of functional is clearly the leading contribution to error at high temperature.
Lastly, we consider the question of anharmonic effects on zero-point energies. Allen’s QP theory can be used with effective frequencies ω(V, T) that incorporate anharmonicity in the T → 0 limit. However, our method to obtain the effective frequencies causes ω(V, T) to converge to the harmonic values at zero temperature, which makes the zero-point contribution to the free energy identical to the harmonic value. Incorporating anharmonic effects as well as correcting for systematic errors from the exchange-correlation functional is important in the calculation of thermodynamic properties in molecules, and this is most easily done by applying a scaling factor to the resulting harmonic frequencies, a method that is widely used in the field68. To examine the effect of introducing scale-factor-style anharmonic corrections to the zero-point contribution, we calculated the PBE/QP thermodynamic properties with the zero-point energy (F0 in Eq. (17)) scaled by factors between 0.8 and 1.2. The V(T) curves obtained from this change are shown in Fig. S5 of the SI. Clearly, at least in MgO, the effect of incorporating anharmonic effects into the zero-point energy is minimal, and reduces to a small shift in the equilibrium volume with essentially no change in slope. Hence, we are well justified to use the harmonic F0 for the calculation of thermodynamic properties in this system. However, in more complex systems featuring light atoms or strong anharmonicity, it may not be the case that anharmonic effects on zero-point energies can be neglected.
Calcium oxide
As a further example, we now consider calcium oxide (CaO), which, like MgO, also presents a high-temperature melting point (2886 K69). In CaO, the B1 (rock-salt) phase is stable in a wide temperature and pressure range. It transitions to the B2 structure (CsCl-like) at around 60 GPa70,71. The calculation details are given in “Computational details”, and the other QP parameters (number of structures, number of volumes, cluster configuration) are equivalent to those used for MgO.
Figure 4a–d shows the calculated thermodynamic properties as a function of pressure using QHA and QP (with the PBE functional) compared against the experimental data where available. The QHA results are accurate at low temperature but, same as in the case of MgO, there is a divergence in the QHA properties at high temperature. The onset temperature for this divergence is lower than in MgO (around 500 K) indicating a greater importance of anharmonicity in the description of CaO, probably caused by the increase in mass of the cation and the resulting decrease in vibrational frequencies. The QP results coincide with QHA at low temperature, but they correct the spurious behavior of QHA past 500 K, resulting in an excellent agreement with the experimental data in the entire temperature range.
The plots show the volume (a), adiabatic bulk modulus (b), constant-pressure heat capacity (c), and volumetric thermal expansion coefficient (d) as a function of temperature at zero pressure, as well as volume as a function of pressure (e). QP (full lines) and QHA (dashed) are compared against experimental data113,122 (circles). In plot (e), two isotherms are shown: 300 K (green, the QHA and QP isotherms overlap) and 2400 K (red).
The room-temperature and 2400 K isotherms for CaO are shown in Fig. 4e, and the former is compared with experimental data. As in the case of MgO, the QHA and QP results are essentially equivalent at room temperature, and their agreement with the experimental isotherm is excellent. However, the QHA isotherm diverges significantly at low pressure close to the melting point. As for the other thermodynamic properties, QP corrects this behavior and restores the proper trend of the isotherm in the entire pressure range.
CaO has a phase transition between the ambient-pressure B1 (NaCl-like) structure and the high-pressure B2 (CsCl-like) phase, experimentally determined to be in the 60–65 GPa range72,73,74. We calculated the thermodynamic properties of the B2 structure of CaO with the same procedure and calculation level as for the B1 phase “Computational details” and computed the transition pressure as a function of temperature and the phase diagram of this system. Both are shown in Fig. 5. The phase diagram shows that, at low temperature, QHA and QP coincide in their prediction of the transition pressure and are both in very good agreement with experiment. We could find no experimental data for the transition at higher temperatures, but QP shows a weaker dependence of the transition pressure with temperature, related to the better behavior of the equilibrium volumes that enter the transition ΔV in the Clapeyron equation.
Strontium titanate
As a further illustration of our methodology we consider now a phase that is dynamically stabilized: cubic strontium titanate (SrTiO3). The cubic phase of SrTiO3 has a perovskite-type structure and is stable only at relatively high temperature75. Below the critical temperature (105 K at ambient pressure), cubic SrTiO3 undergoes a second-order displacive phase transition resulting in a tetragonal structure, in which TiO6 octahedra rotate relative to each other75. The tetragonal to cubic temperature-induced second-order phase transition in SrTiO3 has been extensively studied in the literature75,76,77, and the rotation of the TiO6 octahedra is responsible for the anomalous dielectric properties in this material at low temperature77. Zero-point vibrations ("quantum effects”) in SrTiO3 are enough to make it a paraelectric down to absolute zero, and it is often described as an incipient ferroelectric77.
Cubic SrTiO3 is not modelable using QHA (or perturbative approaches based on QHA78) because the cubic structure is dynamically unstable according to DFT, and the existence of imaginary phonon frequencies precludes the calculation of the free energy and other thermodynamic properties within the harmonic approximation. We calculated cubic SrTiO3 using the PAW method and the PBE functional at 10 volumes spanning the 55 Å3 to 64 Å3 range (roughly up to 20 GPa). We used a 80 Ry wavefunction cutoff, 800 Ry density cutoff, and a 5 × 5 × 5 k-point grid. For the phonon calculations, we employed supercells 36 times larger than the cubic primitive cell (with lattice vectors a = (− 3, − 1, − 1), b = (− 1, − 1, 3), and c = (− 1, 3, 1) relative to the primitive vectors). Cluster configurations comparable to those for MgO and CaO were used: up to sixth order with cutoffs 6.24 Å (2nd order), 5.00 Å (3rd), 4.00 Å (4th), and 2.50 Å (5th and 6th). The calculation of the harmonic force constants required 2 single-point calculations in the supercell structures (with 180 atoms each) per volume, and the phonon rattle was carried out at 500 K, generating 20 structures at each volume. The training of the n-th order force constants, self-consistent phonon calculation, and calculation of the thermodynamic properties proceeded as described in “Computational details”, except a damping value of α = 0.1 in SCHA was used for increased smoothness in the transition temperatures.
At all volumes examined, imaginary frequencies are observed in the harmonic calculation, correctly indicating the instability of the cubic SrTiO3 phase. Application of the SCHA procedure to calculate the effective frequencies ω(V, T) shows that, at a certain volume-dependent temperature, the imaginary frequencies disappear, indicating a temperature-induced stabilization of the cubic phase. This is consistent with previous computational studies using similar methods for computing temperature-dependent effective frequencies78,79. Figure 6a shows the S(T) curves at three grid volumes obtained from our workflow. The temperatures at which imaginary frequencies are present are marked with a lighter color, and can be considered the transition temperature for the tetragonal-cubic second-order phase transition predicted by SCHA. Below the transition temperature, the imaginary frequencies are not used in the entropy calculation, and therefore the S shown in the figure in that region corresponds to a constrained system forced to adopt cubic symmetry. Interestingly, the S(T) curves are perfectly smooth through the transition, indicating that QP theory allows modeling the thermodynamics of the (constrained) cubic phase in the temperature region where it is dynamically unstable.
a S(T) curves predicted for cubic SrTiO3 with our method at various volumes (V0 is the equilibrium volume at ambient conditions). Light color indicates the presence of imaginary frequencies. b Transition temperature for the second-order tetragonal-to-cubic phase transition in SrTiO3 calculated with quasiparticle (QP) theory (blue) against experiment (green)80.
Figure 6b displays our calculated transition temperatures as a function of pressure compared with experimental data80. The transition temperature at each volume is determined as the temperature at which imaginary frequencies disappear in the converged effective FC2 resulting from the SCHA procedure. Because of the stochastic nature of SCHA, there is some variability in the effective FC2, and this results in an uncertainty in the calculated thermodynamic properties and in the transition temperature. While for the former the noise introduced by SCHA is negligible (in the range of 0.01 J K−1 mol−1 to 0.04 J K−1 mol−1), as evidenced by the smoothness of the S(T) curve in Fig. 6a, the predicted transition temperatures carry a more noticeable uncertainty, in the range of a few tens of K (Fig. 6b). The transition temperature at zero pressure is slightly overestimated (slightly under 175 K compared with the experimental 105 K), but otherwise the agreement with experiment is excellent, and our predictions correctly capture the increase in transition temperature with pressure. The deviations from experimental behavior could again be ascribed to the small free energy differences between phases involved in a temperature-induced phase transition, and to idiosyncrasies of the chosen exchange-correlation functional. A strong dependence of the transition temperatures on the choice of exchange-correlation functional has been reported previously78.
Discussion
This article presents a new method for overcoming the limitations of the quasiharmonic harmonic theory (QHA) in the calculation of thermodynamic properties of solids, particularly at high temperature. Our main objectives in the design of the new approach have been computational and conceptual simplicity.
Same as QHA, the proposed method operates on a volume grid spanning the volumes corresponding to the temperature and pressure ranges of interest. At each volume, first the second-order force constant matrix (FC2) is calculated using the DFT forces calculated at randomly displaced atomic configurations, plus regularized regression methods (the hiphive method). Once the FC2 is determined, additional random configurations are generated according to a Boltzmann distribution, and further DFT calculations are carried out to determine the corresponding atomic forces. These forces are then used to fit the n-th order force constant matrices (FCn) up to the order desired by the user.
The FCn obtained from the previous step are then used in a self-consistent harmonic approximation (SCHA) method, in which effective frequencies ω(V, T) are derived for a given temperature T. The SCHA approach used in this work is a modified version of the one used in the hiphive program, in which the long-range non-analytic contribution is calculated beforehand and care is taken to fit the effective FC2 only to the short-range component of the calculated forces. Using QP theory, the quasiparticle entropies SQP are calculated on a temperature grid at every volume.
The third ingredient in the proposed method is a Debye-like numerical model used to calculate all thermodynamic properties from the QP S(V, T) data obtained in the previous step. This model, inspired by similar models used to fit experimental equations of state, is defined only in terms of four parameters per volume: the Debye temperature (ΘD) and three anharmonicity coefficients (ak). These quantities are fitted to S(V, T) and, through them, the free energy, the constant-volume heat capacities, and the rest of the thermodynamic properties are obtained.
The use of the new approach is demonstrated by calculating the thermodynamic properties of MgO and CaO in a wide pressure range up to the zero-pressure melting point. It is shown that the new method corrects the spurious behavior of QHA at high temperature for all thermodynamic properties and recovers the correct experimental trends, to the point where the leading source of error in the calculation of thermodynamic properties at very high temperature is the choice of exchange-correlation functional. Another nice feature of our method is that the FCn is available at every point in the volume grid, which enables the calculation of properties beyond the harmonic approximation. We demonstrate this by predicting the lattice thermal conductivity as a function of temperature in MgO, and showing how the new method achieves excellent agreement with experiment up to the melting point of the material. It was also shown that the new method can be used to model the thermodynamics of dynamically stabilized phases using cubic SrTiO3 as an illustrative example. Lastly, we propose a simple way for incorporating anharmonicity in zero-point energies via scaling factors.
There are several limitations of the present method, some of which will be addressed in future work. First, there is no treatment made of contributions to the free energy from configurational disorder or from low-frequency vibrational modes akin to hindered rotations (e.g. methyl groups in molecular crystals). Second, significant deviations from the static constant-volume equilibrium geometries at high temperature can be problematic since QP, same as QHA, works in the statically constrained approach. Third, there is no way of directly incorporating anharmonicity into the zero-point energies.
In summary, the new method is a simple, computationally efficient, and accurate way of incorporating anharmonic effects up to arbitrary order into thermodynamic calculations in periodic solids that coincides with QHA at low temperature but corrects its spurious blowout in the high temperature range. We expect the new method to be important for the calculation of thermodynamic properties at very high temperature, particularly in the field of geochemistry and high-throughput materials discovery.
Methods
Computational details
All DFT calculations were carried out using the Quantum ESPRESSO package, version 6.581. The projector-augmented wave (PAW) method82 was used with datasets from the pslibrary83. The Perdew-Burke-Ernzerhof (PBE) exchange-correlation functional84 was used for the parameter exploration; other functionals were used for comparison where necessary. The cutoffs for the plane-wave and density expansions are 100 Ry (1361 eV) and 1000 Ry (13610 eV) respectively. The k-point meshes were 5 × 5 × 5 for the primitive cell in MgO and 9 × 9 × 9 for both phases of CaO (B1 and B2). Smaller k-point samplings were used for the corresponding supercells; in particular, 2 × 2 × 2 k-point meshes were used for the 3 × 3 × 3 supercells used in the frequency calculations. Tight criteria (1 × 10−10 Ry) were used for the convergence of the self-consistent field (SCF) cycle. The Born charges required for the long-range non-analytic correction to the force constants were calculated using density-functional perturbation theory (DFPT)20.
The second and n-th order force constant matrices were calculated using random atomic displacements as implemented in the hiphive package47. For the force constant model regressions, we used a least-squares fit and a shuffle-split cross validation with 5 folds (training and testing set sizes of 80% and 20% of the available data, respectively). Given the relatively low number of parameters in our models (because of the high symmetry of MgO and CaO) relative to the number of DFT forces available, we chose a simple least-squares approach instead of more sophisticated regularized and machine-learning-style regression methods43,44,47,85,86,87,88,89.
When necessary for comparison, second-order force constant matrices were obtained using the direct method with phonopy16,17. Phonopy was also used for the calculation of harmonic free energies, entropies, and constant-volume heat capacities using FC2 obtained with either QHA and QP, employing a Fourier-reinterpolation mesh of 20 × 20 × 20. The rest of the thermodynamic properties were obtained either in the QHA or QP approaches using a modified version of the gibbs2 program25,27.
It is important to note that MgO and CaO have no atomic or cell degrees of freedom, i.e. the volume completely determines the system geometry. In more complex systems, the calculation of the force constants, and of the subsequent thermodynamic properties, would require a constant-volume geometry optimization at each point in the grid. This “statically constrained” approach90 is routinely used in QHA studies, and can also be applied in our method in an analogous way.
Calculation of the force constants
Our proposed method starts in the same way as a typical QHA calculation. First, an evenly spaced grid of volumes is set up spanning the volume range of interest. As in QHA, the lowest volume in the grid determines the maximum pressure and the highest volume determines the maximum temperature at which thermodynamic properties may be calculated. Then, constant-volume geometry relaxations are carried out at each of those volumes with tight thresholds, to make atomic forces negligible.
At each volume, the FCn are calculated. In general, the FCn cannot be determined efficiently using either the direct method (because of poor scaling) or the linear response method (because higher-order derivatives need to be computed)17,20,43,44,81,91. Instead, we use a relatively new approach based on computing the DFT forces at randomly perturbed atomic configurations, followed by fitting a set of force constants against those forces. The energy of the crystal in the vicinity of the equilibrium geometry is given by:
where E0 is the equilibrium energy, \({u}_{i}^{\alpha }\) is the displacement of atom i in the α Cartesian direction relative to its equilibrium position, and ϕ are the force constants. The first-derivative term is missing because the atomic forces (\(-{\phi }_{i}^{\alpha }\)) are zero at the equilibrium geometry. Differentiation of this expression gives the ith component of the force on atom α:
If one chooses a set of force constants determined by the maximum order and an order-dependent cutoff beyond which their value is assumed to be zero (the “force constant model”), Eq. (3) allows a determination of the force constants from the calculated DFT forces and displacements using a simple linear least-squares fit. This method, which may use more sophisticated techniques such as regularized regression to simplify the force constant model, is implemented in the hiphive package47. It is similar to other approaches proposed in the literature such as compressive sensing lattice dynamics (CSLD)44,92 but avoids the use of MD trajectories. Compared to the direct approach, it has the advantage that each randomly perturbed configuration provides 3N displacements and forces (N is the number of atoms in the unit cell), and also that the force constant model may contain force constants up to any arbitrary order, provided enough forces are available for the fit. The relative abundance of forces per calculated structure makes this approach advantageous not only in the calculation of FCn, but also for the calculation of the harmonic force constants in large or low-symmetry systems47. Our FCn calculation at the constant-volume equilibrium geometry starts with the generation of an initial set of FC2. This is achieved by first generating a number of structures with randomly displaced atoms and calculating the atomic forces.
In this work, we use periclase (MgO) as an illustrative example for our new approach. MgO has a number of advantages as test case: it has a simple structure, its B1 phase (NaCl-like) has a very wide stability range both in temperature and pressure93,94, its thermodynamic properties, including the equation of state, have been thoroughly studied experimentally94,95,96,97, and it has a very high melting point69 (around 3000 K). The latter is important because it allows observing the breakdown of QHA at high temperature (at around 1000 K31,36 with the PBE functional), something that may not be the case for other systems with lower melting point, even close to the melting temperature5,31,36,70,98. In addition, MgO is a component of the Earth’s mantle, which in combination with the above makes it a commonly studied model system in computational geochemistry23,31,33,34,49,99,100.
For the generation of the initial FC2 in MgO, the atomic displacements have amplitudes of 0.0106 Å (0.02 bohr, a value similar to the 0.01 Å used in phonopy) and a random direction16,17. The number of structures required to converge the FC2 calculation, as well as the effect of other parameters on the calculation, were examined using the harmonic entropy at room temperature (\({S}_{{\rm{HA}}}^{300\,{\rm{K}}}\)) and at 3000 K (\({S}_{{\rm{HA}}}^{3000\,{\rm{K}}}\)), calculated at the equilibrium geometry. The advantage of using the entropy for the convergence tests is that it is calculated with the same formula (and therefore with the same, unmodified software tools) both in the harmonic and quasiparticle approaches.
Regarding the number of structures, because of its very high symmetry (space group \({\rm{Fm}}\bar{3}{\rm{m}}\)), the direct method (phonopy) requires only 2 structures to obtain the full set of FC2 in MgO. Using a 3 × 3 × 3 supercell, the random displacement (hiphive) approach yields results essentially coincident with phonopy with just one structure, and a small cross-validation root-mean-square error (RMSE) of 9.84 meV Å−1 (c.f. the average force is 204.6 meV Å−1). The differences per formula unit between hiphive and phonopy are: \({S}_{{\rm{HA}}}^{300\,{\rm{K}}}\) and \({S}_{{\rm{HA}}}^{3000\,{\rm{K}}}\), in the order of 0.1 J K−1 mol−1; \({F}_{{\rm{HA}}}^{300\,{\rm{K}}}\), 0.03 kJ mol−1; and \({F}_{{\rm{HA}}}^{3000\,{\rm{K}}}\), 0.3 kJ mol−1. The full S(T) curves predicted by phonopy and hiphive are essentially identical—they are shown in Fig. S1 of the Supporting Information (SI). The difference in the number of structures required to converge the thermodynamic properties is more substantial for larger systems, resulting in a reduced computational cost relative to the direct method47.
Another aspect of the initial FC2 calculation is the size of the supercell used, which determines the maximum range of the second-order force constants. For the harmonic calculation, we found that entropies are converged within 0.1 J K−1 mol−1 for a 2 × 2 × 2 supercell with either phonopy or hiphive, corresponding to a maximum distance between interacting atoms of 4.25 Å at the MgO equilibrium geometry. As we shall see below, converging the anharmonic entropies requires a larger supercell, so we decided to use a 3 × 3 × 3 cell for consistency with later calculations. This supercell contains 216 atoms and has an interatomic interaction cutoff distance of 6.38 Å at the equilibrium geometry.
The FC2 calculated in this first step are then used to generate the randomly displaced structures with which we calculate the FCn using hiphive. For the parameter exploration regarding the FCn calculation, we use the quasiparticle entropy at room temperature (\({S}_{{\rm{QP}}}^{300\,{\rm{K}}}\)) and at 3000 K (\({S}_{{\rm{QP}}}^{3000\,{\rm{K}}}\)) as proxies for the convergence of the calculated thermodynamic properties. The entropies are computed as described below. The random configurations are generated from the FC2 according to a Boltzmann distribution with a fixed temperature (the “phonon rattle” method)47,101. Exploration of the possible phonon-rattle generation temperatures reveals that a too low temperature slows convergence of SQP with the number of structures, whereas a temperature that is too high results in large DFT forces and, consequently, larger RMSE from the model fitting with the same FCn model (although not necessarily a poor SQP). In the particular case of MgO, there is a wide temperature plateau in which the SQP are constant within a few tenths of a J K−1 mol−1 regardless of the generation temperature, as seen in Fig. S2 of the SI. We chose 1000 K for the generation temperature.
It is important to note a this point that the generation temperature for the phonon rattle and the maximum temperature accessible by the resulting force constant model are not the same. Table S1 in the SI shows that a force constant model with up to 4th-order constants is enough to converge \({S}_{{\rm{QP}}}^{3000\,{\rm{K}}}\), the entropy at the highest temperature examined, but a much higher order is required to converge the RMSE of the calculated forces generated at a lower temperature (1000 K). For this reason, our recommended procedure is to converge thermodynamic properties like the entropy with respect to the force constant model (as in Table S1 of the SI), and use a suitable phonon-rattle temperature to reliably calculate all energy derivatives in the model.
Two other crucial aspects of the FCn calculation using the hiphive method are the supercell size and the choice of FCn model (the “cluster configuration”43,47,89). The latter refers to the maximum force constant order (n) in the FCn model as well as the cutoff distance for each of the orders. The cluster configuration determines the number of parameters in the FCn model, and therefore indirectly the number of randomly displaced structures required and the overall cost of the calculation. Naturally, the supercell size determines the maximum allowable cutoff distances in the FCn model, which are given by the radius of the largest sphere that can be inscribed in the chosen supercell. Contrary to what happens with the harmonic calculations, a supercell > 2 × 2 × 2 is required for converging the quasiparticle entropies as a function of the order n and the corresponding cutoff distances.
An important point is that, in highly ionic solids like MgO, it is essential to include the long-range non-analytic contribution to the forces in order to converge the FCn with supercell size. This non-analytic contribution, which arises from the macroscopic polarization of the crystal caused by atomic displacements, is responsible for the longitudinal-optical transverse-optical (LO-TO) splitting of the vibrational frequencies at Γ, and can be calculated from the Born charges (Z*) and the dielectric constant (ϵ) at an arbitrary point in the first Brillouin zone102,103. The atomic forces calculated with DFT can be split into a short-range (SR) contribution from the atoms in the supercell, and a long-range (LR) contribution from the macroscopic dipole created by the atomic displacements17,102,103,104:
In the same way, the force constants can be split into a second-order long-range contribution (ϕLR), which can be obtained from Gonze et al.’s formula directly102,103, and a short-range contribution (ϕSR) that contains the rest of the second-order force constants plus all other force constants up to n-th order:
Because only the short-range contribution can be calculated from the displacements, it is necessary to subtract the long-range contribution from the atomic forces prior to the FCn model fitting, and calculate ϕLR separately. This approach is shown in Fig. 7. We show in Fig. S3 of the SI that, by separating the long-range contribution from the forces prior to the model fitting, the calculated entropies converge to within 0.1 J K−1 mol−1 per formula unit with a supercell size of only 3 × 3 × 3, and the RMSE of the forces is consistently smaller (at 3000 K, 62.1 meV Å−1 with the long-range forces separated prior to the fit, compared to 67.8 eV Å−1 without). This approach is used on all FCn models in this article.
The LR contribution is subtracted from the calculated forces and the resulting short-range (SR) forces are fitted with the n-th order force constant (FCn) model. This results in a short-range FCn model that, together with the second-order long-range contribution, generates the complete FCn model.
After a thorough examination of the convergence of \({S}_{{\rm{QP}}}^{3000\,{\rm{K}}}\) with the cluster configuration, we opted for using an FCn model with up to 6-th order interactions and cutoffs corresponding to 8 atomic shells around the considered atom for the n = 2 (second-order) force constants (cutoff distance = 6.38 Å at the equilibrium geometry), 4 shells for n = 3 (4.21 Å), 3 shells for n = 4 (3.64 Å), and 2 shells for n = 5 and n = 6 (2.96 Å). Interestingly, as shown in Table S1 of the SI, \({S}_{{\rm{QP}}}^{3000\,{\rm{K}}}\) can be converged to within a few tenths of a J K−1 mol−1 by any cluster configuration containing all second-order FCs up to the maximum cutoff allowed by the supercell (8 atomic shells) plus a few shells of third- and fourth-order FCs. However, in order to keep the cross-validation RMSE low (indicating a better fit to the DFT forces) and also a reasonable number of parameters, we decided to use the FCn model indicated above. Nonetheless, these observations may be relevant when applying the proposed method to a new, unknown material—the fact that few parameters are required in the FCn model may reduce the computational effort for larger systems.
Lastly, we consider how many structures are required for convergence of the entropy for the chosen supercell size and cluster configuration. Figure 8 shows the convergence in the QP entropy and free energy at 3000 K as a function of the number of structures. The QP entropy is essentially converged to about 0.2 J K−1 mol−1, and the free energy to less than 1.0 kJ mol−1, with only 4 structures, although converging the RMSE requires around 10 structures, in contrast with the single structure required for the harmonic entropies. This is not surprising since the FCn has many more parameters than the FC2 used in the harmonic approximation. Although fewer structures could have been used, we opted for calculating 10 structures at each volume for the FCn fit. It is important to note that in larger and more complicated systems, the difference in the number of structures required to run a QP calculation would be reduced relative to those required by the QHA approach in combination with the direct method.
\({S}_{{\rm{QP}}}^{3000\,{\rm{K}}}\) (a) and \({F}_{{\rm{QP}}}^{3000\,{\rm{K}}}\) (b) are calculated with the effective temperature-dependent harmonic frequencies. c RMSE for the least-squares fit in the training (red circles) and testing (blue stars) sets.
As a final note, alternative methods for calculating DFT forces and n-th order force constants, such as trained machine learning methods, could greatly facilitate the calculation of the FCn required for our method.
Calculation of the temperature-dependent harmonic frequencies
Now we obtain the harmonic frequencies renormalized at a given temperature T from the FCn calculated in the previous step (ω(V, T)). For this purpose, we employ the self-consistent harmonic approximation (SCHA) method implemented in hiphive47,105, modified as described below. In SCHA, the FC2 that minimizes the free energy at temperature T is found by iteration to self-consistence from a trial FC2, employing a variational principle stemming from the Gibbs-Bogoliubov inequality48. The self-consistent FC2 is, in general, different for every temperature and different from the static FC2105,106. The SCHA method, in particular the stochastic self-consistent harmonic approximation (SSCHA) proposed by Errea et al.48,107, has become very popular recently as an excellent way to study strongly anharmonic systems where perturbation theory is inapplicable. One difference between SSCHA and the SCHA proposed here is that our SCHA obtains the forces at the randomly displaced configurations from the FCn fitted in “Calculation of the force constants”, whereas SSCHA calculates them from DFT, which is computationally much more expensive (although SSCHA could in principle also be used with the fitted FCn). In addition, our SCHA is limited in that it relaxes only the degrees of freedom related to the force-constant matrix, not the crystal geometry. Similar SCHA approaches can be developed that rely on configuration sampling using molecular dynamics43,46,48,108,109,110, which our version of SCHA avoids.
In our modified hiphive SCHA approach at temperature T, we start from a trial FC2, which is usually the harmonic FC2. As in “Calculation of the force constants”, any FC2 in the SCHA process is always composed of a long-range part, calculated from the Born charges and independent of temperature, and a short-range part corresponding to the interactions with nearby atoms that is adjusted in the SCHA procedure. The SCHA steps are:
-
1.
A number Nc of randomly displaced atomic configurations are generated according to a Boltzmann distribution given by the current FC2 (a “phonon rattle”), comprising both the short-range and long-range parts.
-
2.
The atomic forces at each configuration are calculated using the FCn model obtained in “Calculation of the force constants”, and the long-range FC2 non-analytic contribution is subtracted.
-
3.
A new short-range FC2 is fitted to the calculated short-range forces. The new FC2 is then obtained by combining the resulting short-range FC2 with the temperature-independent long-range FC2.
These three steps are repeated until convergence of the FC2. Two important aspects are essential to stabilize the convergence of the SCHA process. First, subtracting the long-range contribution in step 2, but not doing the same for the structure generation in step 1, is crucial to prevent interference from the long-range contribution. Second, a damping strategy is used to stabilize the convergence of the SCHA:
where FC2fit is the FC2 obtained in step 3 above. A value of α = 0.2 gives a suitable balance between speed and stability.
In quasiparticle theory, entropies are calculated from the effective harmonic frequencies at temperature T using the harmonic expression. These entropies are the fundamental result that allows the calculation of quasiparticle thermodynamic properties. For this reason, we use the quasiparticle entropies at temperature T, calculated by diagonalizing the corresponding effective FC2, as a way to track the convergence of the SCHA process. In MgO at the equilibrium geometry and 3000 K, the SCHA converges the calculated entropy to within a few hundredths of a J K−1 mol−1, which is accurate enough for our purposes. The number of configurations Nc used in SCHA is also an important quantity. In MgO, we used Nc = 50 because the cost of a high Nc is reasonably low, given that the atomic forces are generated from an already-known FCn (rather than calculated with DFT as in SSCHA48). Exploratory tests reveal that much fewer configurations (Nc ≈ 10) are enough to converge the quasiparticle entropy.
Quasiparticle theory
Allen’s quasiparticle theory36,38 is a simple way of including next-order anharmonic effects into the quasiharmonic approximation that avoids molecular dynamics simulations. In QP theory, the vibrations in the system behave as a non-interacting boson gas at the considered temperature, with vibrational energy levels given by the frequencies ωqn(V, T), which are assumed to be known. The QP thermodynamic properties of the system at V and T are then calculated from the vibrational frequencies, which may originate from experiment or calculations. QP theory therefore assumes that the occupation of the different vibrational modes is uncorrelated.
The entropy in QP is calculated using the same expression as the harmonic entropy13,14,15,16,22,38, which ensures consistency with the anharmonic corrections from perturbation theory:
where fqn are the equilibrium occupations in a Bose-Einstein distribution:
The rationale for QP theory is the observation that the temperature-dependence of the vibrational frequencies predicted by quasiharmonic theory do not match experimental data38, and it was shown that QP theory corrects the quasiharmonic approximation to next order in the calculated thermodynamic properties38. Here, we employ the effective vibrational frequencies obtained from the SCHA procedure “Calculation of the temperature-dependent harmonic frequencies” as quasiparticle energies for QP theory.
The calculation of the QP entropy is straightforward and it does not require implementing new code. Allen provided expressions for other thermodynamic properties36,38, but these are more complicated. As an example, the QP free energy cannot be obtained by merely using the temperature-dependent effective frequencies in the harmonic free energy formula because, in order for the theory to be thermodynamically consistent, SQP must be \(-{\left(\partial {F}_{{\rm{QP}}}/\partial T\right)}_{V}\), and the frequencies themselves have an explicit dependence on temperature. Therefore, the QP free energy is38:
where Esta(V) is the static (DFT) energy, \({F}_{{\rm{vib}}}^{{\rm{HA}}}(V,T)\) is the free energy calculated with the harmonic expression, \(\Delta {F}_{{\rm{vib}}}^{{\rm{anh}}}\) is the anharmonic correction term, ωqn(V) are the static harmonic frequencies, and ωqn(V, T) are the effective frequencies at temperature T. More complex expressions have been given for other thermodynamic properties, and more details can be found in Allen’s original work36,38.
Calculation of thermodynamic properties
The practical implementation of Eqns. (10) to (12) is problematic for two reasons. First, a correct use of Eq. (12) requires matching the q-vector and band indices from the static harmonic frequencies with their effective counterparts coming out of the SCHA procedure, which may not be straightforward. This is particularly true at high temperature where frequencies may shift in position significantly, or if some of the harmonic frequencies are imaginary. Second, expressions would have to be derived and implemented for all other thermodynamic properties we wish to calculate.
Instead, we opt for a simpler alternative inspired by the thermal equations of state used to fit experimental V(p, T) data22,111. We define a model based on Debye theory containing a few parameters and fit the SQP(T) data calculated using QP theory at each volume. For this purpose, we first define the auxiliary function YQP and the corresponding Debye function YDebye as:
where \({F}_{{\rm{vib}}}^{{\rm{QP}}}\) is the QP free energy, \({F}_{{\rm{vib}}}^{{\rm{Debye}}}\) is the free energy in the Debye model, F0 is the vibrational zero-point energy, and \({F}_{0}^{{\rm{Debye}}}\) is the zero-point energy in the Debye model, N is the number of atoms per cell, D3 is the Debye function:
and x = T/ΘD, with ΘD the Debye temperature. The YDebye function (Eq. (14)) contains a single parameter (the Debye temperature) and we expect its shape closely matches that of \({{Y}^{\rm{QP}}}\) (Eq. (13)) at low temperature, where anharmonic effects are usually less important. The deviation between the two is expected to increase with temperature, as the anharmonic effects increase in magnitude. Therefore, it is reasonable to define our model for YQP as the Debye value plus a polynomial in x:
There is only a handful of parameters in this model: the Debye temperature (ΘD) and the coefficients ak. No more than three of the latter (a0, a1, and a2) are required to obtain an excellent fit to the quasiparticle entropies (Fig. 9), so n = 2 is used this article. Substituting in Equation (16) and reordering:
It is straightforward to show that:
and therefore:
where22:
The figure shows the calculated quasiparticle entropy at the highest volume (-20 GPa, red stars), at zero pressure (yellow circles), and at the lowest volume (169 GPa, blue squares) in the grid. The corresponding entropies from the four-parameter Debye model are also shown (full line).
This model is used as a way of calculating all thermodynamic properties as follows. At every volume, we choose a uniformly spaced temperature grid up to the temperature of interest, and calculate the effective harmonic frequencies (ω(V, T), “Calculation of the temperature-dependent harmonic frequencies”) and, with them, the quasiparticle entropies (Eq. (8)). The resulting S(T; V) data is fitted using Eq. (19) to determine the parameters in our model: ΘD and the ak. Figure 9 shows the result of the model fit to the calculated SQP in MgO at the equilibrium geometry. The agreement with the calculated QP entropies is excellent up to the highest temperature considered at all pressures. Clearly, the model is able to represent the quasiparticle thermodynamic behavior accurately.
The model parameters determined in the fit provide a way of calculating F(T; V) (Eq. (17)) and CV(T; V) (Eq. (20)) at any temperature. The only unknown is the vibrational zero-point energy (F0 in Eq. (17)). The harmonic and anharmonic F0 are different in general, and there is no way to calculate the anharmonic F0 from the SCHA approach. One simple approach to recover the anharmonicity in this term would be to scale the harmonic F0 by a functional-dependent constant factor. This approach is used in molecular calculations68. However, as we shall see, applying scaling factors to the F0 has barely any noticeable effect on the calculated thermodynamic properties of the systems examined. Therefore, the harmonic F0 are used in this work. The calculation of the thermodynamic properties other than F, S, and CV follows the usual procedure, employing a strain-polynomial equation of state for the required volume derivatives, as described in our previous work24,25,27,112. This methodology has been implemented in the open-source gibbs2 program25,27, which is used in this work.
The new workflow is summarized in Fig. 10. At a given volume, first the harmonic force constant matrix (FC2) is determined using hiphive. Then, this FC2 is used to generate randomly displaced structures with which the n-th order force constant matrix (FCn) is obtained “Calculation of the force constants”. Using a self-consistent harmonic approximation (SCHA), we determine the effective harmonic frequencies at all temperatures on a uniformly spaced grid that spans the temperature range of interest (ω(V, T)), “Calculation of the temperature-dependent harmonic frequencies”. These frequencies are then used to calculate the quasiparticle entropies (SQP, “Quasiparticle theory”), which are then fitted with a Debye-based model to determine the value of the four parameters (ΘD, a0, a1, a2) corresponding to that volume “Calculation of thermodynamic properties”. The set of four parameters as a function of volume is then used to calculate the free energy, constant-volume heat capacity and, with them, the rest of the thermodynamic properties25,27.
Workflow summary of our new methodology.
Data availability
The datasets used and/or analysed during the current study available from the corresponding author on reasonable request.
Code availability
The underlying code for this study is available in github and can be accessed via these links: https://github.com/aoterodelaroza/gibbs2, https://github.com/aoterodelaroza/hiphive-scripts.
References
Anderson, D. L. Theory of the Earth (Blackwell Scientific Publications, 2023).
Poirier, J. P. Introduction to the Physics of the Earth’s Interior 2nd edn, Vol. 328 (Cambridge University Press, 2000).
Stixrude, L. & Lithgow-Bertelloni, C. Thermodynamics of mantle minerals—I. Physical properties. Geophys. J. Int. 162, 610–632 (2005).
Stixrude, L. & Lithgow-Bertelloni, C. Thermodynamics of mantle minerals—II. Phase equilibria. Geophys. J. Int. 184, 1180–1213 (2011).
Wentzcovitch, R., Yu, Y. & Wu, Z. Thermodynamic properties and phase relations in mantle minerals investigated by first principles quasiharmonic theory. Rev. Mineral. Geochem. 71, 59–98 (2010).
Jain, A. et al. Commentary: The materials project: a materials genome approach to accelerating materials innovation. APL Mater. 1, 011002 (2013).
Talirz, L. et al. Materials cloud, a platform for open computational science. Sci. Data 7, 299 (2020).
Sarikurt, S., Kocabaş, T. & Sevik, C. High-throughput computational screening of 2D materials for thermoelectrics. J. Mater. Chem. A 8, 19674–19683 (2020).
Kirklin, S., Meredig, B. & Wolverton, C. High-throughput computational screening of new Li-ion battery anode materials. Adv. Energy Mater. 3, 252–262 (2013).
Eswarappa Prameela, S. et al. Materials for extreme environments. Nature Rev. Mater. 8, 81–88 (2023).
Bartel, C. J. Review of computational approaches to predict the thermodynamic stability of inorganic solids. J. Mater. Sci. 57, 10475–10498 (2022).
Born, M. & Huang, K. Dynamical Theory of Crystal Lattices Revised ed. edition, Vol. 432 (Oxford University Press, 1996).
Dove, M. T. Introduction to Lattice Dynamics. Cambridge Topics in Mineral Physics and Chemistry 1st edn, Vol. 280 (Cambridge University Press, 1993).
Dove, M. T. Structure and dynamics—an atomic view of materials. Mater. Today 6, 59 (2003).
Wallace, D. Thermodynamics of Crystals, Vol. 484 (Dover Publications, 1998).
Togo, A. & Tanaka, I. First principles phonon calculations in materials science. Scr. Mater. 108, 1–5 (2015).
Togo, A. First-principles phonon calculations with phonopy and phono3py. J. Phys. Soc. Jpn. 92, 012001 (2023).
Parlinski, K., Li, Z. & Kawazoe, Y. First-principles determination of the soft mode in cubic ZrO2. Phys. Rev. Lett. 78, 4063 (1997).
Togo, A., Oba, F. & Tanaka, I. First-principles calculations of the ferroelastic transition between rutile-type and CaCl2-type SiO2 at high pressures. Phys. Rev. B 78, 134106 (2008).
Baroni, S., De Gironcoli, S., Dal Corso, A. & Giannozzi, P. Phonons and related crystal properties from density-functional perturbation theory. Rev. Mod. Phys. 73, 515 (2001).
Ashcroft, N. W. & Mermin, N. D. Solid State Physics, Vol. 311 (Holt-Saunders, 1976).
Anderson, O. L. Equations of State for Solids in Geophysics and Ceramic Science. (Oxford Monographs on Geology and Geophysics) 1st edn, Vol. 432 (Oxford University Press, USA, 1995).
Karki, B., Wentzcovitch, R. D., De Gironcoli, S. & Baroni, S. First-principles determination of elastic anisotropy and wave velocities of mgo at lower mantle conditions. Science 286, 1705–1707 (1999).
Blanco, M., Francisco, E. & Luaña, V. Gibbs: isothermal-isobaric thermodynamics of solids from energy curves using a quasi-harmonic debye model. Comput. Phys. Commun. 158, 57–72 (2004).
Otero-de-la Roza, A., Abbasi-Pérez, D. & Luaña, V. Gibbs2: A new version of the quasiharmonic model code. II. models for solid-state thermodynamics, features and implementation. Comput. Phys. Commun. 182, 2232–2248 (2011).
Santana-Andreo, J. et al. High-throughput prediction of the thermal and electronic transport properties of large physical and chemical spaces accelerated by machine learning: charting the ZT of binary skutterudites. ACS Appl. Mater. Interfaces 16, 4606–4617 (2024).
Otero-de-la Roza, A. & Luaña, V. Gibbs2: A new version of the quasi-harmonic model code. I. robust treatment of the static data. Comput. Phys. Commun. 182, 1708–1720 (2011).
Nath, P. et al. High-throughput prediction of finite-temperature properties using the quasi-harmonic approximation. Comput. Mater. Sci. 125, 82–91 (2016).
Qin, T., Zhang, Q., Wentzcovitch, R. M. & Umemoto, K. Qha: A python package for quasiharmonic free energy calculation for multi-configuration systems. Comput. Phys. Commun. 237, 199–207 (2019).
Jofré, J., Gheribi, A. E. & Harvey, J.-P. Development of a flexible quasi-harmonic-based approach for fast generation of self-consistent thermodynamic properties used in computational thermochemistry. Calphad 83, 102624 (2023).
Karki, B. B., Wentzcovitch, R. M., de Gironcoli, S. & Baroni, S. High-pressure lattice dynamics and thermoelasticity of MgO. Phys. Rev. B 61, 8793–8800 (2000).
Carrier, P., Wentzcovitch, R. & Tsuchiya, J. First-principles prediction of crystal structures at high temperatures using the quasiharmonic approximation. Phys. Rev. B 76, 064116 (2007).
Erba, A., Shahrokhi, M., Moradian, R. & Dovesi, R. On how differently the quasi-harmonic approximation works for two isostructural crystals: thermal properties of periclase and lime. J. Chem. Phys. 142, 044114 (2015).
Belmonte, D. First principles thermodynamics of minerals at HP–HT conditions: MgO as a prototypical material. Minerals 7, 183 (2017).
Wentzcovitch, R., Karki, B., Cococcioni, M. & De Gironcoli, S. Thermoelastic properties of MgSiO3-perovskite: insights on the nature of the Earth’s lower mantle. Phys. Rev. Lett. 92, 018501 (2004).
Allen, P. B. Theory of thermal expansion: quasi-harmonic approximation and corrections from quasi-particle renormalization. Mod. Phys. Lett. B 34, 2050025 (2020).
Monserrat, B., Drummond, N. & Needs, R. Anharmonic vibrational properties in periodic systems: energy, electron-phonon coupling, and stress. Phys. Rev. B 87, 144302 (2013).
Allen, P. B. Anharmonic phonon quasiparticle theory of zero-point and thermal shifts in insulators: Heat capacity, bulk modulus, and thermal expansion. Phys. Rev. B 92, 064106 (2015).
Sun, T., Zhang, D.-B. & Wentzcovitch, R. M. Dynamic stabilization of cubic CaSiO3 perovskite at high temperatures and pressures from ab initio molecular dynamics. Phys. Rev. B 89, 094109 (2014).
Kawai, K. & Tsuchiya, T. P-V-T equation of state of cubic CaSiO3 perovskite from first-principles computation. J. Geophys. Res. Solid Earth 119, 2801–2809 (2014).
Zhang, D.-B., Sun, T. & Wentzcovitch, R. M. Phonon quasiparticles and anharmonic free energy in complex systems. Phys. Rev. Lett. 112, 058501 (2014).
Wang, H., Zhuang, J., Zhang, Z., Zhang, Q. & Wentzcovitch, R. M. Pgm: A python package for free energy calculations within the phonon gas model. Comput. Phys. Commun. 291, 108845 (2023).
Zhou, F., Nielson, W., Xia, Y. & Ozoliņš, V. Lattice anharmonicity and thermal conductivity from compressive sensing of first-principles calculations. Phys. Rev. Lett. 113, 185501 (2014).
Zhou, F., Nielson, W., Xia, Y. & Ozoliņš, V. Compressive sensing lattice dynamics. I. General formalism. Phys. Rev. B 100, 184308 (2019).
Zhang, Z. & Wentzcovitch, R. M. Ab initio anharmonic thermodynamic properties of cubic CaSiO3 perovskite. Phys. Rev. B 103, 104108 (2021).
Hellman, O., Abrikosov, I. A. & Simak, S. I. Lattice dynamics of anharmonic solids from first principles. Phys. Rev. B 84, 180301 (2011).
Eriksson, F., Fransson, E. & Erhart, P. The hiphive package for the extraction of high-order force constants by machine learning. Adv. Theory Simul. 2, 1800184 (2019).
Errea, I., Calandra, M. & Mauri, F. Anharmonic free energies and phonon dispersions from the stochastic self-consistent harmonic approximation: application to platinum and palladium hydrides. Phys. Rev. B 89, 064302 (2014).
Kwon, C., Xia, Y., Zhou, F. & Han, B. Dominant effect of anharmonicity on the equation of state and thermal conductivity of MgO under extreme conditions. Phys. Rev. B 102, 184309 (2020).
Xia, Y. Revisiting lattice thermal transport in PbTe: the crucial role of quartic anharmonicity. Appl. Phys. Lett. 113, 073901 (2018).
Goldman, V., Horton, G. & Klein, M. An improved self-consistent phonon approximation. Phys. Rev. Lett. 21, 1527 (1968).
Otero-de-la-Roza, A. & Johnson, E. R. Application of xdm to ionic solids: the importance of dispersion for bulk moduli and crystal geometries. J. Chem. Phys. 153, 054121 (2020).
Plata, J. J., Posligua, V., Márquez, A. M., Fernandez Sanz, J. & Grau-Crespo, R. Charting the lattice thermal conductivities of I–III–VI2 chalcopyrite semiconductors. Chem. Mater. 34, 2833–2841 (2022).
Blancas, E. J. et al. Unraveling the role of chemical composition in the lattice thermal conductivity of oxychalcogenides as thermoelectric materials. J. Mater. Chem. A 10, 19941–19952 (2022).
de Koker, N. Thermal conductivity of MgO periclase from equilibrium first principles molecular dynamics. Phys. Rev. Lett. 103, 125902 (2009).
de Koker, N. Thermal conductivity of MgO periclase at high pressure: Implications for the D” region. Earth Planet. Sci. Lett. 292, 392–398 (2010).
Stackhouse, S., Stixrude, L. & Karki, B. B. Thermal conductivity of periclase (MgO) from first principles. Phys. Rev. Lett. 104, 208501 (2010).
Manthilake, G. M., de Koker, N., Frost, D. J. & McCammon, C. A. Lattice thermal conductivity of lower mantle minerals and heat flux from Earth’s core. Proc. Natl. Acad. Sci. 108, 17901–17904 (2011).
Callaway, J. Model for lattice thermal conductivity at low temperatures. Phys. Rev. 113, 1046–1051 (1959).
Omini, M. & Sparavigna, A. An iterative approach to the phonon Boltzmann equation in the theory of thermal conductivity. Physica B: Condens. Matter 212, 101–112 (1995).
Tang, X. & Dong, J. Lattice thermal conductivity of mgo at conditions of earth’s interior. Proc. Natl. Acad. Sci. 107, 4539–4543 (2010).
Tian, Z. et al. Phonon conduction in PbSe, PbTe, and PbTe1−xSex from first-principles calculations. Phys. Rev. B 85, 184303 (2012).
Allen, P. B. Improved Callaway model for lattice thermal conductivity. Phys. Rev. B 88, 144302 (2013).
Fugallo, G., Lazzeri, M., Paulatto, L. & Mauri, F. Ab initio variational approach for evaluating lattice thermal conductivity. Phys. Rev. B 88, 045430 (2013).
Li, W., Carrete, J., Katcho, N. A. & Mingo, N. ShengBTE: a solver of the Boltzmann transport equation for phonons. Comp. Phys. Commun. 185, 1747–1758 (2014).
Hofmeister, A. M. Thermal diffusivity and thermal conductivity of single-crystal MgO and Al2O3 and related compounds as a function of temperature. Phys. Chem. Miner. 41, 361–371 (2014).
Slifka, A. J., Filla, B. J. & Phelps, J. Thermal conductivity of magnesium oxide from absolute, steady-state measurements. J. Res. Natl. Inst. Stan. Technol. 103, 357 (1998).
Alecu, I., Zheng, J., Zhao, Y. & Truhlar, D. G. Computational thermochemistry: scale factor databases and scale factors for vibrational frequencies obtained from electronic model chemistries. J. Chem. Theory Comput. 6, 2872–2887 (2010).
Haynes, W. M. CRC Handbook of Chemistry and Physics 96th edition, Vol. 2677 (CRC Press, 2015).
Karki, B. B. & Wentzcovitch, R. M. Vibrational and quasiharmonic thermal properties of CaO under pressure. Phys. Rev. B 68, 224304 (2003).
Liu, Z.-J. et al. Thermoelasticity of CaO from first principles. Chin. Phys. 16, 499 (2007).
Richet, P., Mao, H.-K. & Bell, P. M. Static compression and equation of state of CaO to 1.35 Mbar. J. Geophys. Res. Solid Earth 93, 15279–15288 (1988).
Mammone, J. F., Mao, H. K. & Bell, P. M. Equations of state of CaO under static pressure conditions. Geophys. Res. Lett. 8, 140–142 (1981).
Jeanloz, R., Ahrens, T. J., Mao, H. K. & Bell, P. M. B1-B2 transition in calcium oxide from shock-wave and diamond-cell experiments. Science 206, 829–830 (1979).
Cowley, R. A. The phase transition of strontium titanate. Philos. Trans. R. Soc. A 354, 2799–2814 (1996).
Slonczewski, J. & Thomas, H. Interaction of elastic strain with the structural transition of strontium titanate. Phys. Rev. B 1, 3599 (1970).
Kvyatkovskii, O. Quantum effects in incipient and low-temperature ferroelectrics (a review). Phys. Solid State 43, 1401–1419 (2001).
Tadano, T. & Tsuneyuki, S. Self-consistent phonon calculations of lattice dynamical properties in cubic SrTiO3 with first-principles anharmonic force constants. Phys. Rev. B 92, 054301 (2015).
Tadano, T. & Tsuneyuki, S. Ab initio prediction of structural phase-transition temperature of SrTiO3 from finite-temperature phonon calculation. J. Ceram. Soc. Jap. 127, 404–408 (2019).
Guennou, M., Bouvier, P., Kreisel, J. & Machon, D. Pressure-temperature phase diagram of SrTiO3 up to 53 GPa. Phys. Rev. B 81, 054115 (2010).
Giannozzi, P. et al. Advanced capabilities for materials modelling with quantum ESPRESSO. J. Phys.: Condens. Matter 29, 465901 (2017).
Blöchl, P. E. Projector augmented-wave method. Phys. Rev. B 50, 17953 (1994).
Dal Corso, A. Pseudopotentials periodic table: from H to Pu. Comput. Mater. Sci. 95, 337–350 (2014).
Perdew, J. P., Burke, K. & Ernzerhof, M. Generalized gradient approximation made simple. Phys. Rev. Lett. 77, 3865–3868 (1996).
Ng, A. Y. Feature selection, L1 vs. L2 regularization, and rotational invariance. In Proc. twenty-first international conference on Machine learning, 78 (IEEE, 2004).
Simon, N., Friedman, J., Hastie, T. & Tibshirani, R. A sparse-group LASSO. J. Comput. Graph. Stat. 22, 231–245 (2013).
Tibshirani, R. Regression shrinkage and selection via the LASSO. J. R. Stat. Soc. Series B Stat. Methodol. 58, 267–288 (1996).
Hastie, T., Tibshirani, R. & Friedman, J. The Elements of Statistical Learning: Data Mining, Inference, and Prediction 2nd edn, Vol. 745 (Springer, 2009).
Nelson, L. J., Hart, G. L. W., Zhou, F. & Ozoliņš, V. Compressive sensing as a paradigm for building physics models. Phys. Rev. B 87, 035125 (2013).
Carrier, P., Justo, J. F. & Wentzcovitch, R. M. Quasiharmonic elastic constants corrected for deviatoric thermal stresses. Phys. Rev. B 78, 144302 (2008).
Giannozzi, P., de Gironcoli, S., Pavone, P. & Baroni, S. Ab initio calculation of phonon dispersions in semiconductors. Phys. Rev. B 43, 7231–7242 (1991).
Zhou, F., Sadigh, B., Åberg, D., Xia, Y. & Ozoliņš, V. Compressive sensing lattice dynamics. II. Efficient phonon calculations and long-range interactions. Phys. Rev. B 100, 184309 (2019).
Duffy, T. S., Hemley, R. J. & Mao, H.-k Equation of state and shear strength at multimegabar pressures: magnesium oxide to 227 GPa. Phys. Rev. Lett. 74, 1371–1374 (1995).
Belonoshko, A. B., Arapan, S., Martonak, R. & Rosengren, A. MgO phase diagram from first principles in a wide pressure-temperature range. Phys. Rev. B 81, 054110 (2010).
Dubrovinsky, L. S. & Saxena, S. K. Thermal expansion of periclase (MgO) and tungsten (W) to melting temperatures. Phys. Chem. Miner. 24, 547–550 (1997).
Anderson, O. L. & Zou, K. Thermodynamic functions and properties of MgO at high compression and high temperature. J. Phys. Chem. 19, 69–83 (1990).
Tange, Y., Nishihara, Y. & Tsuchiya, T. Unified analyses for P-V-T equation of state of MgO: A solution for pressure-scale problems in high p-t experiments. J. Geophys. Res. Solid Earth 114, B03208 (2009).
Otero-de-la-Roza, A. & Luaña, V. Equations of state and thermodynamics of solids using empirical corrections in the quasiharmonic approximation. Phys. Rev. B 84, 184103 (2011).
Oganov, A. R., Gillan, M. J. & Price, G. D. Ab initio lattice dynamics and structural stability of MgO. J. Chem. Phys. 118, 10174–10182 (2003).
Wang, Y., Liu, Z.-K., Chen, L.-Q., Burakovsky, L. & Ahuja, R. First-principles calculations on MgO: phonon theory versus mean-field potential approach. J. Appl. Phys. 100, 023533 (2006).
West, D. & Estreicher, S. K. First-principles calculations of vibrational lifetimes and decay channels: hydrogen-related modes in Si. Phys. Rev. Lett. 96, 115504 (2006).
Gonze, X., Charlier, J.-C., Allan, D. & Teter, M. Interatomic force constants from first principles: the case of α-quartz. Phys. Rev. B 50, 13035–13038 (1994).
Gonze, X. & Lee, C. Dynamical matrices, born effective charges, dielectric permittivity tensors, and interatomic force constants from density-functional perturbation theory. Phys. Rev. B 55, 10355–10368 (1997).
Wang, Y., Shang, S.-L., Fang, H., Liu, Z.-K. & Chen, L.-Q. First-principles calculations of lattice dynamics and thermal properties of polar solids. npj Comput. Mater. 2, 16006 (2016).
Hooton, D. J. A new treatment of anharmonicity in lattice thermodynamics: I. London Edinburgh Philos. Mag. & J. Sci. 46, 422–432 (1955).
Klein, M. L. & Horton, G. K. The rise of self-consistent phonon theory. J. Low Temp. Phys. 9, 151–166 (1972).
Monacelli, L. et al. The stochastic self-consistent harmonic approximation: calculating vibrational properties of materials with full quantum and anharmonic effects. J. Phys.: Condens. Matter 33, 363001 (2021).
Souvatzis, P., Eriksson, O., Katsnelson, M. I. & Rudin, S. P. Entropy driven stabilization of energetically unstable crystal structures explained from first principles theory. Phys. Rev. Lett. 100, 095901 (2008).
Souvatzis, P., Eriksson, O., Katsnelson, M. & Rudin, S. The self-consistent ab initio lattice dynamical method. Comput. Mater. Sci. 44, 888–894 (2009).
Hellman, O., Steneteg, P., Abrikosov, I. A. & Simak, S. I. Temperature dependent effective potential method for accurate free energy calculations of solids. Phys. Rev. B 87, 104111 (2013).
Angel, R. J., Miozzi, F. & Alvaro, M. Limits to the validity of thermal-pressure equations of state. Minerals 9, 562 (2019).
Flórez, M., Recio, J. M., Francisco, E., Blanco, M. A. & Pendás, A. M. First-principles study of the rocksalt–cesium chloride relative phase stability in alkali halides. Phys. Rev. B 66, 144112 (2002).
Fiquet, G., Richet, P. & Montagnac, G. High-temperature thermal expansion of lime, periclase, corundum and spinel. Phys. Chem. Miner. 27, 103–111 (1999).
Sinogeikin, S., Jackson, J., O’Neill, B., Palko, J. & Bass, J. Compact high-temperature cell for Brillouin scattering measurements. Rev. Sci. Instrum. 71, 201–206 (2000).
White, G. K. & Anderson, O. L. Grüneisen parameter of magnesium oxide. J. Appl. Phys. 37, 430–432 (2004).
Jacobsen, S. D. et al. Compression of single-crystal magnesium oxide to 118 GPa and a ruby pressure gauge for helium pressure media. Am. Mineral. 93, 1823–1828 (2008).
Perdew, J. P. et al. Restoring the density-gradient expansion for exchange in solids and surfaces. Phys. Rev. Lett. 100, 136406 (2008).
Becke, A. D. On the large-gradient behavior of the density functional exchange energy. J. Chem. Phys. 85, 7184 (1986).
Becke, A. D. & Johnson, E. R. Exchange-hole dipole moment and the dispersion interaction revisited. J. Chem. Phys. 127, 154108 (2007).
Johnson, E. R. The exchange-hole dipole moment dispersion model. In, Non-Covalent Interactions in Quantum Chemistry and Physics (Otero-de-la-Roza, A. & DiLabio, G. A.) 169–194 (Elsevier, 2017).
Otero-de-la-Roza, A. & Johnson, E. R. Van der Waals interactions in solids using the exchange-hole dipole moment. J. Chem. Phys. 136, 174109 (2012).
Yamanaka, T., Kittaka, K. & Nagai, T. B1-b2 transition in CaO and possibility of CaSiO3-perovskite decomposition under high pressure. J. Mineral. Petrol. Sci. 97, 144–152 (2002).
Acknowledgements
This study was funded by the Spanish MICIU/AEI/10.13039/501100011033 and European Union Next Generation EU/PRTR (TED2021-130874B-I00, PID2022-138063OB-I00, PID2021-122585-NB-C21-2, TED2021-129457B-I00, PGC2021-125518NB-I00, RED2022-134388-T, and CNS2023-144958) and by the Principado de Asturias (FICYT) and FEDER (AYUD/2021/51036). Computer resources/technical support from the Spanish Supercomputing Network (RES) are acknowledged: Xula/CIEMAT (RES-QHS-2023-1-0027), Lusitania/Cénits-COMPUTAEX (QHS-2022-3-0032), and MareNostrum5/BSC (RES-AECT-2024-2-0010). The funders played no role in study design, data collection, analysis and interpretation of data, or the writing of this manuscript.
Author information
Authors and Affiliations
Contributions
E.J.B., J.J.P., and A.O.R. conceived and designed the study, with help from A.L., A.M.M., J.M.R., P.N., and J.J.P.. E.J.B., A.L., F.I.R., and P.N. performed first-principles calculations, supervised by J.M.R. and A.M.M.. E.J.B. and A.O.R. designed the methodology with help from J.J.P. and P.N.. E.J.B., A.O.R., and J.J.P. performed the data analysis and investigation. J.J.P., A.M.M., and J.M.R. secured funding and resources for the study. E.J.B. and A.O.R. wrote the initial version manuscript. E.J.B., A.L., FIR, A.M.M., J.M.R., P.N., and J.P.P. discussed the results and contributed to the final version of the manuscript. All authors read and approved the final manuscript.
Corresponding author
Ethics declarations
Competing interests
The authors declare no competing interests.
Additional information
Publisher’s note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Supplementary information
Rights and permissions
Open Access This article is licensed under a Creative Commons Attribution-NonCommercial-NoDerivatives 4.0 International License, which permits any non-commercial use, sharing, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if you modified the licensed material. You do not have permission under this licence to share adapted material derived from this article or parts of it. The images or other third party material in this article are included in the article’s Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by-nc-nd/4.0/.
About this article
Cite this article
Blancas, E.J., Lobato, Á., Izquierdo-Ruiz, F. et al. Thermodynamics of solids including anharmonicity through quasiparticle theory. npj Comput Mater 10, 267 (2024). https://doi.org/10.1038/s41524-024-01447-8
Received:
Accepted:
Published:
DOI: https://doi.org/10.1038/s41524-024-01447-8