Hybrid ab-initio/experimental high temperature equations of state: Application to the NaCl pressure scale

Accurate equations of state (EoS) are important to describe materials properties in several scientific domains, from low temperature physics to astrophysics.^1–3 For instance, geophysics thrives on knowledge of equations of state of minerals and melts, and the oil/gas industry needs accurate equations of state to perform phase-behavior calculations in hydrocarbons.⁴ Nowadays, it is possible to measure equations of state parameters accurately, but the pressure and temperature range needed for many applications is still difficult to reach experimentally. An outstanding example is the burgeoning field of planetary modeling, which has gained a tremendous impulse with the discovery of exoplanets.^5,6 Today, ab initio calculations contribute an essential tool to predict phase behavior and equations of state at extreme conditions.^7,8 The essential methods involved are based on density functional theory (DFT)^9,10 and quasiharmonic approximation (QHA)¹¹ combined with ab initio phonon calculations. Molecular dynamics and thermodynamics integration also play an important role but are less practical for solid state calculations because of limited phase space sampling. There are two main problems with the DFT/QHA approach: uncertainties associated with DFT and anharmonicity. Static energies are calculated solving the one-electron Kohn-Sham equation¹⁰ self-consistently and then using the total charge density to obtain the total energy. However, the absence of an exact description of exchange correlation potential and energies produces systematic errors when compared with experimental data. There are numerous approximations for the exchange correlation functional, e.g., the popular local density approximation (LDA), the generalized gradient approximation (GGA), meta-functionals, hybrid functionals, and extended functionals (DFT+U, DFT+vdW), and theoretical equations of state are very sensitive to the choice of functional.¹² The second issue concerns the QHA, which works exceedingly well in its range of validity, but, as temperature rises, anharmonic effects become noticeable and the QHA becomes inadequate. Resolution of both issues is likely to take time and the need to generate practical and useful databases of equations of state and thermodynamics and thermoelastic properties has motivated analytical approaches to combine experimental data and computational results.^13–15 Experimental data at low pressures and high temperatures are usually very accurate, while results of calculations at high pressures are much less dependent on the choice of exchange correlation functional. Also, anharmonicity becomes less important at low temperatures, except, obviously, in the vicinity of dynamical instabilities.

By now, several analytical approaches to combine the best data of both worlds, i.e., low pressure and high temperature experimental data and high pressure and high temperature ab initio results, have been proposed. Wu et al.¹³ used an exponential function to correct MgO's theoretical 300 K compression curve followed by an approximate correction to anharmonic effects. Kunc and Syassen¹⁴ observed that experimental and theoretical (static GGA and/or LDA) compression curves plotted using scaled volume and pressure agree considerably better. Based on this observation, Otero-de-la-Roza and Luana¹⁶ developed an empirical correction to the static energies and Zhang et al.¹⁵ extended this approach to thermoelasticity. Except for the first work (Wu et al.¹³), the limitations of DFT and QHA have not been explicitly decoupled and addressed separately.

Here, we consider these issues and explore analytical ways of combining experimental data and computational results in modeling NaCl's equation of state.^13–15 Theoretical calculations of alkali halides in general have not been as successful as those for other ionic compounds in the same structure with stronger bonds, e.g., MgO.¹⁷ In the early 1980's, the recognition that alkali elements have large cores lead to development of the partial core correction in pseudopotential calculations.¹⁸ As it will be illustrated here, LDA, GGA, and hydrid calculations¹⁹ using projector augmented wave (PAW)^20,21 potentials have not solved the problem yet, as shown in NaF, LiCl, and LiF calculations.^22,23 Another motivation for studying NaCl comes from the discovery of huge oil/gas reservoirs below a thick layer of salt at the southeast coast of Brazil, coast of Africa, Mexico, etc.,^24,25 the so-called pre-salt oil fields. Seismic mapping of these oil fields and salt regions require good knowledge of its high temperature equation of state. Finally, accurate equations of state of NaCl are important in high pressure physics, since it is a popular pressure calibrants.^26–32 Improving methodologies to calculate and predict high temperature equations of state is an important step in advancing high pressure technology at conditions still challenging to experiments. High pressure mineral physics still struggles with conflicting measurements of pressure using different standards,⁵ and more accurate ways to determine pressure are still highly desirable, especially at high temperatures.

To validate theoretical/experimental hybrid approaches to generate reliable equations of state, we investigate the performance of four types of corrections applied to 300 K results and extended to other relatively low temperatures by correcting the common static energy first. Such corrected DFT results are followed by a semi-empirical correction to anharmonicity³³ at higher temperatures.

Ab initio calculations were performed with Quantum ESPRESSO³⁴ and VASP³⁵ computational packages using three different functionals for the exchange correlation energy, LDA,³⁶ GGA-PBE,³⁷ and HSE06.¹⁹ Electronic wave functions were described using the PAW^20,21 method. The plane wave cutoff used was 60 Ry. Brillouin zone sampling for electronic states was done on a displaced 6 × 6 × 6 k-mesh; dynamical matrices were obtained in a $Γ$-point centered 4 × 4 × 4 q-mesh using density functional perturbation theory³⁸ and then interpolated in a 12 × 12 × 12 q-mesh to produce the vibrational density of states. Thermodynamic properties were calculated using the QHA, with the Helmholtz free energy given by

where the first term is the static internal energy given by DFT, the second is the zero point motion energy, and the last is the contribution from thermal excitations. We used a third order finite strain Birch–Murnaghan equation of state

to fit results by three different functionals. In this equation, V₀, K₀, and $K′$⁠, are zero pressure equilibrium volume, bulk modulus, and pressure derivative of bulk modulus, respectively.

As shown in Fig. 1, 300 K ab initio compression curves differ considerably from experiments. The equilibrium volume obtained with LDA is about 7% smaller than the experimental one, while GGA-PBE overestimates it in about 6%, as shown in Table I. This is a very large difference and will produce a large error in pressure if used as a pressure calibrant. With the increase in pressure, the difference between LDA and PBE predictions and experimental data decreases. At 25 Gpa, the difference is 3.6% and 1.7% for LDA and PBE, respectively. These differences are still large compared to required pressure scale accuracy. Although HSE06 improves the results, the equilibrium volume still differs by 3.5% and at 25 GPa the error is about 1.0%. It is also worth noting that the calculated and experimental $K′$ differ from each other and from the standard value of 4. This is an important factor when considering corrections to DFT.

In general, the QHA gives very accurate results at low temperatures, so that the difference between experiments and calculations are due to DFT errors at 300 K.¹⁷ This can be verified by comparing calculated and measured thermal expansion coefficients.^13,41,42 Figure 2 shows the thermal expansivity, α, calculated with these functionals. Although PBE and HSE06 overestimate this quantity, LDA describes it well up to 600 K. Therefore, we will exemplify the corrective procedure for LDA results only. As will be indicated, large differences in thermal expansivity require further elaboration of the corrective procedure. This is not a hindrance, in principle, but are more laborious.

Our approach consists in correcting first the calculated 300 K compression curve and then working backwards, i.e., subtracting vibrational effects, finding the corrected static compression curve, and correcting static energies. Then, starting from the latter, we compute high temperature compression curves, which should be in agreement with high temperature measurements if the DFT errors have been properly removed. In principle, if the procedure is consistent, it should be possible to correct results obtained with any functional, as long as errors do not push calculations into anharmonic regions, e.g., into very large volumes due to thermal expansion.

The results presented in Figs. 1 and 2 show that some sort of correction is essential to get accuracy comparable to experiments. Several analytical schemes have been proposed to correct ab initio results for DFT errors. Here, we investigate four different corrections. Wu et al.¹³ proposed a simple correction based on the fact that, at high pressures, the EoS is better described by the standard DFT functionals

where $ΔV0$ is the difference between calculated and measured volumes at 0 GPa and P_c is a parameter to be determined by comparing calculated and measured equations of state. The P_c that minimized DFT errors for NaCl is 16.9. This correction worked very well for MgO,¹³ but DFT errors for NaCl are relatively larger. Equation (3) derives from the following assumption:

Here, we generalize further this correction and introduce a second free parameter

The resulting expression takes then the form

where $ΔV0$ has the same meaning. Optimal value for P_c and A is 12.53 and 1.152, respectively.

Another type of correction was proposed by Kunc and Syassen.¹⁴ This scheme is based on the observation that the second order finite strain equation of state, the Birch–Murnaghan equation,

for scaled quantities,

produces a unique curve, if a second order finite strain equation of state is suitable for the material in question. In Eq. (8), β stands for LDA, GGA, etc., or experiment, $V0β$ and $K0β$ are the respective zero pressure equilibrium volume and bulk modulus. Equation (7) is usually a suitable equation of state if the pressure range of measured volumes is small compared to the bulk modulus or if $K′=4$⁠. The well-known correlation between underestimated equilibrium volume and overestimated bulk modulus, or vice-versa, observed in virtually all calculations, is therefore a consequence of Eq. (7). This insight was useful to design an analytical correction for calculated results. As indicated in Eq. (7),

Then,

This scaling motivated the Kunc and Syassen (KS) correction;

Such correction based on the second order Birch–Murnaghan equation (Eq. (7)), or third order with $K′=4$⁠, will not be adequate if $K′exp$ differs from $K′DFT$⁠, which is usually the case when a third order Birch–Murnaghan equation of state (Eq. (2)) is necessary, i.e., when the compression spans a pressure range comparable to the K₀. The KS correction can then be generalized by using relationships derived from Eq. (2) for scaled quantities, including scaled $K′$⁠, or simply put, replacing $K′DFT$ by $K′exp$⁠, producing to the following correction:

where V is calculated solving the equation of state (Eq. (2)) at scaled pressures.

These are the four types of volume corrections to 300 K DFT results we will investigate in this paper, and we refer to them as WW (Eq. (3)), MWW (Eq. (6)), KS, and KSr (for Kunc & Syassen revised). However, we need corrections at all temperatures. To achieve this, we need to start from correct static compression curves. The current corrective scheme assumes that the error at 300 K and 0 K are the same, since thermal volume expansion is about 2.5% at 300 K and is well described by calculations, in this case LDA (see Fig. 2)

Therefore,

where the second term on the r.h.s. above is negligible. Thus, we start from the following corrected static volume:

where $δVzp(P)$ is the volume change caused by zero point motion. The corrected static compression curve yields $ΔPst(V)$ that is then used to correct the static energy

As indicated by Wu and Wentzcovitch,³³ after correcting for DFT errors, anharmonic effects⁴⁵ become more evident. To improve results at very high temperatures, we use a semi-empirical correction³³ to the QHA free energy. This method is based on the assumption that the Helmholtz free energy, Eq. (1), is still a good approximation for anharmonic phonons with “temperature renormalization frequencies”¹¹ given by

where Ω are T dependent frequencies at V and ω are static frequencies calculated at $V′$⁠. The relation between V and $V′$ is given by

where c is an empirical parameter to be determined, by minimizing the mean squared error (MSE) between the corrected V(0,T) and high temperature experimental measurements. This expression reflects the fact that anharmonic effects are more prominent at high temperatures and low pressures and decrease with increasing pressure. This relation implies that

from which anharmonic thermodynamic properties can be calculated.

Experimental and ab initio compression curves for all functionals investigated in this paper up to 3 GPa and 770 K and to 30 GPa at 300 K are shown in Fig. 1. Corresponding equation of state parameters are given in Table I. Results were fit with a third order Birch–Murnaghan equation of state, compared to experimental data from Boehler and Kennedy³⁹ and the Hugoniot from Fritz et al.⁴⁰ After applying all four corrections to these ab initio 300 K compression curves, they are essentially indistinguishable by the eye from the experimental ones. Figure 3 shows the mean square error between corrected and experimental compression curves for all cases.

It can be seen that the introduction of a second parameter (MWW) increases noticeably the accuracy of the correction compared to the correction based on just one free parameter (WW). The error in the KS correction is the largest compared to the others and can be significantly reduced using the third order Birch–Murnaghan equation of state with scaled $K′$⁠. These corrections are then transferred to the static energies following the procedure outlined in the previous session and the high temperature equations of state recalculated with the QHA. LDA results are presented in Fig. 4.

All corrections work extremely well in the pressure range considered, validating this procedure to construct equations of state with the accuracy of experimental data at low temperatures. As seen in Fig. 2, the LDA thermal expansion starts deviating from experiments at around 550 K. Beyond this temperature, the QHA starts to become inadequate. This behavior is evident in Fig. 4, which includes the mean squared errors for all corrections for several temperatures. It is clear that the errors increase with temperature for any of the corrections.

At this point, we add to the DFT corrected results shown in Fig. 4, the anharmonic corrections indicated in Eqs. (19) and (20). The resulting equation of state for MWW and KSr corrections are shown in Fig. 5, along with their mean squared errors. As can be seen, the larger errors at high temperatures are successfully removed.

The error in $ΔV$ between our results and Bohler's experimental data is 0.03% at 3.1 GPa and 300 K and 0.24% at 3.2 GPa and 770 K. This is a huge improvement compared to the −8% error for LDA results at ambient conditions before the correction. Since the experimental error from Bohler and Kennedy data is 0.7% in V, the procedures presented in this study were able to construct predictive equations of state with the same accuracy as experimental data. Table II compares the calculated pressures in this work with data from Boehler and Kennedy.³⁹ 

Figure 6 shows our equation of state up to ∼23 GPa at 500 K and 800 K compared to other semi-empirical high temperature equations. It indicates that this procedure could be used to construct highly accurate equations of state.

In addition, it is possible to obtain corrected thermodynamical properties with the anharmonic Helmholtz free energy. Figure 7 shows corrected thermal expansivity, α, Grüneiser parameter, γ, and constant pressure specific heat, C_P. Although corrected and uncorrected Grüneiser parameters fall within the experimental error bars, the corrected one seems to have a more correct temperature dependence. The specific heat is greatly improved. Uncorrected results deviate at small temperatures and diverge at high temperatures. These errors are nicely removed in the corrected curves. The thermal expansion shows a small deviation, but still has a good agreement with experiments.

In this paper, we test and validate a systematic procedures to reconcile experimental data on equations of sate with ab initio results, correcting the 300 K compression curve. The performance of previously proposed correction schemes was compared to newly proposed ones for NaCl, and significant improvement has been achieved with the new schemes. We further extend this procedure to remove DFT errors from the static energy, which is then used to provide relatively low temperature (up to 300 K) equations of state with experimental accuracy, i.e., <0.07 GPa in pressure. To account for errors at higher temperatures, a simple anharmonic correction successfully removed errors caused by the QHA, and extend the validity of our equation of state to temperatures up to near T_melting. This approach is intended to be applicable to weakly anharmonic systems, primarily.

We thank Gaurav Shukla and Professor Lucy Assali for helpful comments on the manuscript. M.L.M. thanks support by CAPES from Brazil through Science without Borders Program process No. BEX 14456/13-3. R.M.W. thanks support by NSF Grant No. EAR-1341862. We acknowledge computer allocation at the Minnesota Supercomputing Institute and a Blue Waters grant through the Great Lakes Consortium.