Accurate thermoelastic tensor and acoustic velocities of NaCl

Despite the importance of thermoelastic properties of minerals in geology and geophysics, their measurement at high pressures and temperatures are still challenging. Thus, ab initio calculations are an essential tool for predicting these properties at extreme conditions. Owing to the approximate description of the exchange-correlation energy, approximations used in calculations of vibrational effects, and numerical/methodological approximations, these methods produce systematic deviations. Hybrid schemes combining experimental data and theoretical results have emerged as a way to reconcile available information and offer more reliable predictions at experimentally inaccessible thermodynamics conditions. Here we introduce a method to improve the calculated thermoelastic tensor by using highly accurate thermal equation of state (EoS). The corrective scheme is general, applicable to crystalline solids with any symmetry, and can produce accurate results at conditions where experimental data may not exist. We apply it to rock-salt-type NaCl, a material whose structural properties have been challenging to describe accurately by standard ab initio methods and whose acoustic/seismic properties are important for the gas and oil industry.

Despite the importance of thermoelastic properties of minerals in geology and geophysics, their measurement at high pressures and temperatures are still challenging.Thus, ab initio calculations are an essential tool for predicting these properties at extreme conditions.Owing to the approximate description of the exchangecorrelation energy, approximations used in calculations of vibrational effects, and numerical/methodological approximations, these methods produce systematic deviations.Hybrid schemes combining experimental data and theoretical results have emerged as a way to reconcile available information and offer more reliable predictions at experimentally inaccessible thermodynamics conditions.Here we introduce a method to improve the calculated thermoelastic tensor by using highly accurate thermal equation of state (EoS).The corrective scheme is general, applicable to crystalline solids with any symmetry, and can produce accurate results at conditions where experimental data may not exist.We apply it to rock-salt-type NaCl, a material whose structural properties have been challenging to describe accurately by standard ab initio methods and whose acoustic/seismic properties are important for the gas and oil industry.C 2015 Author(s).All article content, except where otherwise noted, is licensed under a Creative Commons Attribution 3.0 Unported License.
[http://dx.doi.org/10.1063/1.4938550]2][3][4] However, interpretation of seismic velocity models requires knowledge of equations of state (EoS) and thermoelastic properties of rock-forming minerals. 5,6As such, accurate estimates of thermoelastic tensors of minerals are fundamental in geophysics.They are also important for the gas/oil industry.Several gas fields have been discovered with the use of seismic surveys. 7n the context of gas/oil exploration, rock-salt NaCl has been attracting great interest since the discovery of large reserves beneath thick layers of salt on continental shelves, the so called pre-salt fields.In this letter, we report accurate thermoelastic properties of rock-salt NaCl (B1 structure) at relevant pressures (up to 0.5 GPa) and temperatures (up to 400 K) and beyond, which should be useful for surveying pre-salt fields.First, we perform ab initio free energy calculations based on the quasi-harmonic approximation (QHA) 8,9 using state-of-the-art density-functional theory (DFT) based codes (Quantum ESPRESSO 10 and VASP 11 ).Then, we use a recently proposed method to reconcile ab initio results and experimental data for thermal EoS. 12 Here, this method is further extended to reconcile ab initio results and experimental data for elastic coefficients.Experimental techniques to measure elastic properties at high pressures and temperatures (PT) are very well established.Ultrasonic interferometry [13][14][15] and Brillouin scattering [16][17][18][19] have provided fundamental information for interpretation of wave-velocities in Earth's interior.These properties can be accurately determined at low pressures and temperatures.However, experimental data become scarcer at combined high pressures and temperatures, including for NaCl.][22] The limiting factors in these calculations are: a) systematic deviations originating in the choice of exchange-correlation functional in DFT and b) limitations associated with the QHA -absence of intrinsic anharmonic effects.It is therefore desirable to combine experimental data and ab initio results in a sensible manner to produce optimally predicted thermal elastic properties.
Ab initio calculations were perfomed within the QHA, with three different functionals: LDA, 23,24 GGA-PBE, 25 and HSE06. 26We used the projector augmented wave (PAW) method 27 to describe the electronic wave function with a Brillouin zone sampling on a displaced 8x8x8 k-points mesh.Dynamical matrices at each sampled volume were calculated by the finite displacement method, using a 4x4x4 supercell for LDA, and GGA-PBE, and 2x2x2 for HSE06.Force constants were interpolated on a 12x12x12 mesh to produce phonon density of states (VDoS) to be used in conjunction with QHA free energies calculations: where q is the phonon wave vector, m is the normal mode, and E is the static internal energy at volume V .The second and third terms of the r.h.s are the zero point and the thermal contributions to the free energy, respectively.Rock-salt type NaCl has cubic structure with three independent elastic coefficients.At a given volume, negative and positive strains of 1% were applied and static elastic coefficients were calculated using the relation: Isothermal elastic coefficients can be separated into static and temperature dependent (vibrational) parts: The vibrational contribution to C vib i j k l was calculated using the semi-analytical approach by Wu and Wentzcovitch. 22This method avoids the calculation of VDoS for strained configurations, by assuming an isotropic distribution of phonon frequencies.9][30] With the isothermal elastic coefficients, the adiabatic ones can be obtained as: where C V is the heat capacity at constant volume and S the entropy at P and T. 21 It has been well documented 5,6 that the popular exchange correlation functionals LDA and PBE-GGA in general give static volumes smaller and larger than the experiments, respectively.These systematic deviations impact directly on the elastic properties.For instance, deviations in the bulk modulus correlate inversely with deviations in equilibrium volume. 31LDA with vibrational effects accounted by the QHA offers most accurate EoS for weakly correlated and strongly bonded materials.2][33] Some of these corrections to the theoretical EoS have adjustable parameters that minimize the difference between corrected ab initio results and experimental data. 12,32Other corrective schemes consistently rescale EoS variables to bring calculated and measured EoS parameters into close agreement. 31,33Here we extend one of these corrective schemes 12 to the elastic coefficients.Differences between ab initio results and experimental data are pressure dependent, therefore, the 300 K corrected volume is written as: V cor (P, 300) = V DFT (P, 300) + ∆V (P, 300) where ∆V (P, 300) is a numerical or analytical function modeling the difference between ab initio results and experimental data.Assuming that the error in thermal expansion is small implies that ∆V (P, 300) ∆V (P, 0), then: V cor (P, 0) = V DFT (P, 0) + ∆V (P, 300).
Therefore, corrected static volumes are obtained by subtracting the extra volume difference due to zero point motion: V cor st (P) = V DFT st (P) + ∆V (P, 300) − ∆V z p (P) which can be written as: f (P) is then used to correct the static energy that enters in QHA calculations.This corrective procedure offers highly accurate hybrid EoS 12 in the range of validity of the QHA and here we extend it to elastic coefficients.These calculations start by computing static coefficients (first term on the r.h.s of Eq.( 3)) using Eq. ( 2).Isothermal elastic coefficients can also be written as: where G(P,T) is the Gibbs free energy: Thus, the static elastic coefficients are: Conversely, the static compliance tensor is: Since static volumes are calculated via the Gibbs free energy: the corrected volume must be the derivative of a "corrected" Gibbs free energy: as in Eq. ( 8).Integrating Eq. ( 14) to obtain G cor st and using Eq. ( 12), we arrive at the corrected compliance tensor: V cor st (P) from which the elastic coefficients can be calculated: Similarly to corrected thermal EoS, as long as the thermal expansivity is well described by the uncorrected ab initio calculation and the temperature is within the range of validity of the QHA, there is no need to correct the thermal contribution to the elastic coefficient.
Calculated and measured 34,35 300 K EoS of NaCl are compared in Figure 1.LDA and GGA-PBE volumes are smaller and larger, respectively, than experimental measurements.HSE06 improves the EoS, but still produces a large deviation from experiments.Figure 2 shows the thermal expansion coefficients.GGA-PBE and HSE06 overestimate the thermal expansion, but LDA describes it well up to 600 K. 5,6,21 This is essentially the LDA limit of QHA for this material, where anharmonic effects start causing deviations of corrected results from experimental data. 12Corrective FIG. 4. a) Pressure dependence of corrected LDA elastic coefficients using at 300 K, b) and their temperature dependence at 0 GPa.][38][39] procedures for the EoS depend on a good description of the thermal expansivity.Thus, here we perform corrections only on LDA results.Calculated adiabatic elastic coefficients are shown in Voigt notation (C 11 = C 1111 , C 12 = C 1122 and C 44 = C 2323 ) in Figure 3 compared to experimental data from Kinoshita et.al., 36 Voronov et.al., 37 Whitfield et.al. 38 and Yamamoto et.al. 39 As LDA underestimates the experimental volume, C 11 is overestimated, similarly to pressure. 31The opposite occurs with GGA-PBE and HSE06.HSE06 agrees best with experiments, but still has large errors.C 44 is less volume-sensitive and agrees well with experimental data.
We used the corrected EoS 12,40 with this corrective scheme to produce the elastic coefficients shown in Figure 4.The pressure dependence of C 11 , C 12 and C 44 at 300 K is in excellent agreement with single crystal data by Kinoshita et.al. 36 and Voronov et.al. 37 Temperature dependent results at 0 GPa are also improved for C 11 , but it appears that even at 300 K the data by Yamamoto et.al. 39 differs from those by Kinoshita et.al. 36 and Voronov et.al. 37 As shown in Figure 2, above about ≈ 500 K, anharmonic effects cause deviations from experiments and further anharmonic corrections 12,32 or fully anharmonic calculations 41,42 are needed.With these elastic coefficients, bulk and shear moduli were calculated using Voigt-Reuss-Hill (VRH) averages [43][44][45] and are shown in Figure 5(a).Excellent agreement with data by Kinoshita et.al. 36 and Voronov et.al. 37 is expected given the good agreement between elastic coefficients.The corrected bulk modulus is not significantly changed, compared to the uncorrected one.Experimental bulk modulus obtained from the EoS by Decker 46 and from the hugoniot by Fritz et.al. 35 deviate somewhat from those calculated using corrected elastic coefficients.Some differences are also observed with respect to those obtained from measured velocities in polycrystalline samples. 47Those by Voronov were obtained using elastic coefficients measured on single crystals. 37The shaded areas in both figures are the upper and lower bonds of the Voig-Reuss-Hill averaging scheme.As expected, experimental data are within calculated bounds.
In summary, the best agreement between corrected and measured elastic coefficients is obtained by comparing results with single crystal elasticity data, such as those by Kinoshita et.al., 36 Voronov et.al. 37 and Whifield et.al. 38 The C 11 mean squared error is significantly reduced from 144 to 5, comparing to Voronov's data. 37Data obtained on polycrystaline NaCl 47 do not agree as well, nevertheless, they are within the Voigt-Reuss-Hill bounds.Such aggregates should be more isotropic than single crystals.The anisotropy of single crystal NaCl increases with pressure, since the deviation from the isotropic condition, C 11 − C 12 = 2C 44 increases with pressure.This increased anisotropy can also be seen in Fig. 5(a), where Reuss, G R , and Voigt, G V , bounds 43 are shown (K R = K V ).For an isotropic system, G R = G V , a condition held only at 0 GPa.With increasing pressure, G R and G V depart from each other, a clear signature of the pressure enhanced anisotropy.The difference between corrected results and data on aggregates with increasing pressure, as shown in Fig. 5(b), is therefore likely to be caused by the single crystal anisotropy.Further work is needed to address this issue.
We have extended a scheme to combine/reconcile experimental data and ab initio results for thermal EoS 12 to thermoelastic coefficients.Using this approach, we were able to predict quite accurately the single crystal elastic tensor and, within uncertainties, acoustic velocities in poly-crystalline rock-salt NaCl at relevant thermodynamics conditions.The proposed corrective scheme is general, applies to other crystalline systems, and can extend data with predictive accuracy to thermodynamics conditions difficult to be reached by experiments, as long as anharmonicity does not compromise the accuracy of QHA-based predictions.

FIG. 5 .
FIG. 5. (a)Bulk and shear moduli and (b) acoustic velocities compared to experimental data.Elastic moduli by Frankel were obtained via polycrystalline velocity measurements47 and Decker's EoS.46Those by Voronov were obtained using elastic coefficients measured on single crystals.37The shaded areas in both figures are the upper and lower bonds of the Voig-Reuss-Hill averaging scheme.As expected, experimental data are within calculated bounds.