Evaluation of thermodynamic equations of state across chemistry and structure in the materials project

Thermodynamic equations of state (EOS) for crystalline solids describe material behaviors under changes in pressure, volume, entropy and temperature, making them fundamental to scientific research in a wide range of fields including geophysics, energy storage and development of novel materials. Despite over a century of theoretical development and experimental testing of energy–volume (E–V) EOS for solids, there is still a lack of consensus with regard to which equation is indeed optimal, as well as to what metric is most appropriate for making this judgment. In this study, several metrics were used to evaluate quality of fit for 8 different EOS across 87 elements and over 100 compounds which appear in the literature. Our findings do not indicate a clear “best” EOS, but we identify three which consistently perform well relative to the rest of the set. Furthermore, we find that for the aggregate data set, the RMSrD is not strongly correlated with the nature of the compound, e.g., whether it is a metal, insulator, or semiconductor, nor the bulk modulus for any of the EOS, indicating that a single equation can be used across a broad range of classes of materials. A systematic comparison between the performances of several thermodynamic equations of state revealed the superiority of three equations. Equations of state are widely used to describe materials properties based on variables like temperature, pressure, volume, etc. Now, a team from University of California Berkeley and the Lawrence Berkeley National Lab aim to determine the most suitable one for various conditions. The authors used DFT calculations to model the properties of hundreds of elemental, binary and ternary crystalline solids and subsequently fit them with the most commonly-used equations of state. The Birch, Tait and Vinet equations showed the lowest deviation from calculated points, while fitting reasonably well experimental data; this holistic approach underlines that there is not one equation of state to fit all cases.


INTRODUCTION
Thermodynamic equations of state (EOS) for crystalline solids describe material behaviors under changes in pressure, volume, entropy and temperature; making them fundamental to scientific research in a wide range of fields, including geophysics, energy storage and development of novel materials. [1][2][3] Despite over a century of theoretical development and experimental testing of EOS for solids, 4,5 there is still a lack of consensus on the most appropriate EOS under various conditions or even the metric to evaluate appropriateness. Previous attempts 6,7 have focused on comparison with the experimental data, which limits the range of the accessible data. By contrast, computational studies can span chemical and structural space in a consistent and methodical manner, and are often able to probe extreme conditions beyond the range of current experimental techniques. Density-functional theory (DFT) is gaining grounds as a first principles methodology that can accurately and efficiently search material space for complex functionality. In well-characterized systems, DFT calculations have been shown to be precise and reliable. 8 In particular, Lajaeghere et al. recently benchmarked results for DFT-calculated thermodynamic properties of the elements, and found overall satisfactory agreement between prediction and experiment (excepting cases where magnetic, correlative, or relativistic effects are significant; or in cases where van der Waals forces are dominant). 9 While direct calculations of systems under pressure are possible in most DFT codes, the equations of state provide a simple mathematical handle with minimum computational cost. These considerations pave the way for high-throughput studies that probe extreme conditions using DFT generated EOS.
There are three key questions to answer in order to utilize highthroughput DFT-based EOS. Which EOS best encapsulates the fundamental physics presented by DFT and are there rules or metrics that define EOS applicability? Are the optimal DFT-based EOS consistent with experimental values, and are therefore likely to reflect our physical reality? What are the limitations of these EOS?

Equations of state
The equations investigated in this study are described in detail below, and summarized in Table 1. Birch (Euler and Lagrange) Both forms of the Birch equation were derived by Francis Birch in 1947 for crystalline solids of cubic symmetry. 10 The derivation was based on Francis Murnaghan's extensive tensor formalism for analyzing finite strains. 11 The first Birch equation uses the Eulerian strain metric; the second form is based on the Lagrangian metric (roughly the inverse of the Eulerian). As noted by both Birch and Murnaghan, the Eulerian metric is generally considered to be a more apt description of elastic behavior for non-infinitesimal strains, since it treats the final (rather than initial) coordinates as the independent variable. The energy-volume relation arising from the Eulerian viewpoint is: where ν V Vo , and E o and V o are the calculated equilibrium energy and volume, respectively, at zero pressure and absolute zero temperature. As for all subsequent formulations, the values of bulk modulus (K) and its pressure derivative (K′) at ν = 1 may be expressed in terms of the fitting parameters B and C: Birch's equation from the Lagrangian viewpoint is: where: Mie-Gruneisen The Mie-Gruneisen EOS is based on the well-established empirical form of a general interatomic potential: 4,5 Setting m = −1 gives a reasonable representation of the Coulombic interaction between two atoms; 12 under the restriction that E = E o at r = r o reduces the four parameters in Eq. (7) to two to give: where: Murnaghan Interestingly, the Murnaghan equation is based not on Francis Murnaghan's prodigious 1937 paper, but rather on a short communication from 1944 in which he proposed a relatively simple EOS based on a linear variation of the bulk modulus with respect to pressure. 13 This hypothesis leads to the equation: where: Pack-Evans-James Pack et al. posited an exponential variation of pressure with respect to changes in volume (based on quantum mechanics), and adjusted their expression to ensure proper limiting behavior as pressure approaches zero or infinity: 14 where: Poirier-Tarantola Poirier and Tarantola proposed that pressure is best expressed in terms of a logarithmic strain metric (as opposed to the Eulerian or Lagrangian metrics of Birch's equations). 15 The resulting energy-volume relationship is: where: Tait The Tait equation was developed nearly 150 years ago based entirely on the empirical observations of Peter Tait, who was investigating the compressibility of seawater, and modified a few decades later by Gustav Tammann from a linear average to a C − 2 10 Mie-Gruneisen ν ÀC À1 C þ ν À 1 B C + 1 13 Pack-Evans-James Evaluation of thermodynamic equations of state across chemistry and K Latimer et al.

2
npj Computational Materials (2018) 40 Published in partnership with the Shanghai Institute of Ceramics of the Chinese Academy of Sciences differential form, 16 leading to the relation: where:

Vinet
The Vinet equation was shown by Stacey to be equivalent in formulation to the Rydberg potential: 7,17 This leads to the following energy-volume equation: 18 where: The Jellium EOS, 19 a third-order polynomial fit to energies, has historically shown a high quality of fit when utilized with the DFT data. 8 The authors for this EOS motivate their equation (which is essentially a polynomial expansion in powers of x V V0 À1=3 ) by arguing that each power of x corresponds to a contribution to the energy of a crystal structure as a function of its cell volume. However, it does not appear that these justifications are entirely physical. For example, the x 2 term in the Jellium model is aimed to represent the kinetic energy of an electron gas surrounding the atoms in the crystal. However, if the total energy varies roughly in (as is the case with NaCl, for instance), this coefficient is almost certainly negative if one examines the Taylor expansion of E as a function of x. Hence, as the Jellium EOS may in some cases revert to a non-physical polynomial function, we did not include it in this study.
Optimal equation of state Metrics and correlations were derived from the equation of state and data compiled in the Materials Project. The relative root mean square deviation (RMSrD) was used as the primary metric for quality of fit, defined as: 20 The RMSrD values for the fits of all EOS across the investigated material set are shown for the elements (Fig. 1) (see Supplementary Fig. S1 in SI for data on the compound set). Noble gases have been excluded from the elemental set due to drastic deviations from ideal behavior for solids: in the most extreme case, oscillations in the E-V curve of Kr lead to an RMSrD of 37%, orders of magnitude larger than the mean values in Fig. 2. Other studies have shown that it is necessary to incorporate the effects of van der Waals forces in the calculations to adequately capture the energy curve behavior of noble gases, 21 however, these forces were not included in our simulations. There is an interesting periodic trend in the plot of Fig. 1, which indicates a correlation between error in EOS fit and number of valence electrons for a given element. We conjecture that this is due to the inability of the simplified EOS models to capture the physics of systems with greater numbers of interacting bodies (the implicit assumption being that DFT is able to do so more accurately, although previous comparison to experimental data has demonstrated that this  method is likely still prone to some error. 9 ) It is clear from this plot that the Lagrangian formulation of the Birch EOS (Eq. (5)) fares significantly worse than the other EOS under consideration and for this reason, it has been removed from subsequent plots in order not to obscure other data. This property has been remarked by previous authors, 6,10 however, we do not believe any physical interpretation has previously been offered. We therefore point out that in approximately harmonic systems, where the E-V curve is dominated by a quadratic strain term, the Eulerian metric more closely resembles the true E-V relation for most solids (Fig. 3). This observation might be taken into consideration for future theoretical work.
There is a large degree of variability in RMSrD from material to material, captured by the violin plot in Fig. 2 which displays the distribution for the combined set of elemental and compound materials. The average relative deviations are less than 1% regardless of material or theoretical equation used to fit the E-V data for most materials investigated, which suggests that for general purposes of approximation, the particular choice of EOS is not crucial. Because the mean RMSrD values for some of the equations were quite close to each other, we also compared EOS performance on the basis of a few other error metrics: maximum relative deviation (MrD), the greatest magnitude of fractional deviation for a given fitted E-V curve from DFT-calculated points; statistical uncertainty for each of the fit parameters (E 0 , V 0 , K, K′) as calculated from the least-squares. These analyses do not indicate a clearly superior EOS (see for example, Supplementary Fig. S2 in the SI), however, do suggest that the Birch (Eulerian), Tait, and Vinet equations tend to show lower deviation from calculated data points and less variability of fitting parameters than the other equations investigated.
The RMSrD values did not exhibit strong correlation with the bulk modulus or the band gap in the elemental and compound sets (see Supplementary Tables S1 and S2 and Supplementary Figs. S3 and S4 in SI). This suggests that neither the metallic/ insulating nature nor softness of the material are predictive of error in the fitted EOS. The effect of structure within a fixed chemical system was also investigated for several materials, examples of which are shown in Figs. 4 and 5 for various polymorphs of Al2O3 and Ga, respectively. Al2O3 shows essentially no variation of quality of fit for any equation across a range of crystal systems. This is especially notable in light of the fact that both forms of the Birch equation were derived for media of cubic symmetry, 10 with the rest assuming an isotropic medium. In contrast, Ga does exhibit noticeable variation in RMSrD between different polymorphs. Similar trends were observed for other chemical systems, with stiff metal oxides exhibiting low variability relative to softer metals (see Supplementary Figures S5 and S6 in SI). It is true that structures with different space groups may still have similar physical properties; this typically occurs when one structure is a slight distortion of another with higher symmetry. We examined the supergroup/subgroup relations of the given systems, and did not find this to be the determining factor in the similarity or difference in error trends. For example, the only relation of this type in the alumina polymorph set (besides the trivial P1, which is a subgroup of every space group) is P21/c and C2/m (numbers 14 and 12, respectively); whereas for gallium, Cmce (number 64) is a subgroup of Cmcm (number 63) and I4/ mmm (number 139). These findings suggest that variations in structure are not necessarily predictive of variations in error for fitted EOS.

Benchmarking to Experiment
Given the overall quality of performance of the Birch (Eulerian), Tait, and Vinet equations, these were selected further benchmarking analysis. To verify the ability of DFT-based EOS to predict    Supplementary  Figs. S8-S11 of SI. There is limited experimental data on secondorder elastic constants (used to calculate K′), but a few systems are compared in Fig. 8. Agreement between experimental and calculated values of elemental and compound bulk moduli is good, with regression slope and R 2 values for both sets very close to 1 (Fig. 6). Previous studies have indicated that zero-point energy can lead to greater discrepancy between 0 K calculations and experiment for lighter elements, 22 however, the error in bulk modulus values was not correlated to atomic mass for the elemental data set in this study (see Supplementary Figure S7). A few elements show variation from experimental data of greater than 50 GPa, labeled on the figure. Of these, Cr, Np, and Pa have been studied elsewhere by both experiment and calculation, and have shown similar systematic deviations as reported here. The discrepancy for Cr arises from magnetic ordering, 23 which was not considered in this study. Previous experimental studies have shown that the bulk modulus of this element increases considerably at lower temperatures (nearly 20% from room temperature to 77 K). 24 The variations for Np and Pa have been attributed to bonding in 5f orbitals. 25,26 Previous studies using exact and overlapping muffin tin orbitals have been able to reproduce experimental elastic properties of actinides, in particular plutonium. 27 This all-electron DFT method, while effective, is computationally expensive and difficult to generalize, making it unsuitable for EOS generation in high-throughput schemes.
Comparison to experimental values of K′ is more difficult, since they are small dimensionless quantities (typically in the range of 1-10), and the combination of experimental uncertainties and inherent variability due to different methods of fitting data lead to a wide range of cited literature values even for a single material. Figure 8 shows a few representative examples. Since the bulk modulus pressure derivative can be used to derive other material properties, such as behavior at finite temperatures, 7 it is useful to have an internally consistent method for deriving K′ (although it  should be noted that, since K′ is the highest-order fit parameter for the EOS investigated here, it is subject to the largest statistical uncertainty; see Fig. 9e). In light of their demonstrated accuracy in predicting K values, we therefore propose that numericallyderived EOS are valid avenues for analyzing and predicting thermodynamic properties of solid crystalline materials.
Bulk moduli calculated from the EOS were also compared against bulk moduli calculated from the elastic tensor using the formation detailed in elastic tensor workflow for Materials Project. 28 For the most part, values are identical, with some elements showing noticeably lower values from the elastic tensor workflow (Fig. 7). This is likely due to higher order terms the E-V curves, or asymmetries of the unit cell (e.g., in the cases of Tb and Dy). The discrepancy in carbon is due to the nature of the ground state structure, which is layered (graphite). This leads to a highly anisotropic elastic tensor, and the Voigt-Reuss-Hill averaging scheme is skewed by the low shear modulus of this structure. 28 On the other hand, the shear modulus is not accounted for by the isotropic volumetric deformations employed in this study resulting in a far less significant effect of the lack of inter-layer cohesive forces on the calculated bulk modulus.
Since many applications for thermodynamic EOS involve materials subjected to extreme pressures (e.g., in geophysical research, 29 ) it is natural to inquire whether there are pressure limits beyond which the DFT-derived EOS are no longer valid. A rough measure of these limits is given by comparing the maximal value found in the literature for pressure exerted on a material to that reached by the DFT simulation. Our findings, listed in Supplementary Tables S5 and S6 in SI, indicate that in many cases (e.g., AgCl), the latter of these limits is greater than or equal to the former, so that no extrapolation of the EOS is necessary for comparison with experimental results. However, in some cases, extrapolation would still be necessary: for example, for BeO, a stiff ceramic material, the maximum calculated volumetric compression of 15.7% corresponds to a calculated pressure of 127.5 GPa, whereas pressures of up to 175 GPa have been achieved in laboratory settings 30 (see Supplementary Table S5 in SI). Because computations are capable of pressurizing a material well beyond the validity of the pseudopotential approximation, extrapolation of the EOS may in some cases give a more physically valid (and computationally cheaper) prediction than simply performing higher-pressure calculations.
Many of the materials exhibit phase transformations far below the pressure range of the DFT curves, in which case the question is no longer whether a single DFT-derived EOS is valid at high pressures, but whether we can use DFT-based EOS of multiple polymorphs to predict these transitions. This analysis calls for simulations at finite temperatures, and has not been included in the present study.
In summary, the energy-volume curves of over 200 elemental, binary, and ternary crystalline solids were calculated using DFT, and the resulting data were fit to several theoretical EOS commonly cited in the literature. Although the quality of fit, as measured by RMSrD, did not vary greatly across the investigated set of equations, the Birch (Eulerian), Tait, and Vinet equations were found to give the best overall quality of fit to calculated energy-volume curves as compared to the other equations under examination. Band gap and bulk modulus were not found to be useful indicators of EOS quality of fit for the aggregate set of materials, although we observed that quality of fit does vary with structure for polymorphs within certain chemical systems. Bulk modulus values derived from the calculated EOS were benchmarked against experimental data, and in general display good agreement.

METHODS
The calculation of EOS is automated using self-documenting workflows compiled in the atomate 31 code base that couples pymatgen 32 for materials analysis, custodian 33 for just-in-time debugging of DFT codes, and Fireworks 34 for workflow management. The EOS workflow begins with a structure optimization and subsequently calculates the energy of isotropic deformations, including ionic relaxation with volumetric strain ranging from from −33 to 33% (−10 to 10% linear strain) of the optimized structure. Density-functional-theory (DFT) calculations were performed as necessary using the projector augmented wave (PAW) method 35,36 as implemented in the Vienna Ab Initio Simulation Package (VASP) 37,38 within the Perdew-Burke-Ernzerhof (PBE) Generalized Gradient Approximation (GGA) formulation of the exchange-correlation functional. 39 A cutoff for the plane waves of 520 eV is used and a uniform k-point density of approximately 7,000 per reciprocal atom is employed. Convergence tests were performed for selected materials across a range of bulk modulus values, with KPPRA values between 1,000 and 11,000; plots are provided in the SI (Supplementary Figure S12). In addition, standard Materials Project Hubbard U corrections are used for a number of transition metal oxides, 40 as documented and implemented in the pymatgen VASP input sets. 32 We note that the computational and convergence parameters were chosen consistently with the settings used in the Materials Project 41 to enable direct comparisons with the large set of available MP data.

Data availability
Numerical EOS for the data set investigated in this study have been included as a Supplementary json File (see explanation of file structure in SI). The equivalent data are also available for free to the public, via the Materials Project online database.