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.

Most information about Earth’s interior is obtained by analyzing the internal propagation of seismic waves.1–4 However, interpretation of seismic velocity models requires knowledge of equations of state (EoS) and thermoelastic properties of rock-forming minerals.5,6 As 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.7 In 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 ESPRESSO10 and VASP11). 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 interferometry13–15 and Brillouin scattering16–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. In contrast, ab initio calculations of thermal elastic coefficients can reach much higher pressures and temperatures.5,6,20–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.26 We used the projector augmented wave (PAW) method27 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:

F V , T = E ( V ) + 1 2 q , m ħ ω q , m ( V ) + k B T q , m ln 1 exp ħ ω q , m ( V ) k B T
(1)

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:

σ i j P δ i j = C i j k l s t ϵ kl .
(2)

Isothermal elastic coefficients can be separated into static and temperature dependent (vibrational) parts:

C i j k l T ( P , T ) = C i j k l s t ( P ) + C i j k l v i b ( P , T ) .
(3)

The vibrational contribution to C i j k l v i b was calculated using the semi-analytical approach by Wu and Wentzcovitch.22 This method avoids the calculation of VDoS for strained configurations, by assuming an isotropic distribution of phonon frequencies. This method provided accurate high temperature/pressure elasticity of several minerals.22,28–30 With the isothermal elastic coefficients, the adiabatic ones can be obtained as:

C i j k l S ( P , T ) = C i j k l T ( P , T ) + T V C V S ( P , T ) e i j S ( P , T ) e kl
(4)

where CV is the heat capacity at constant volume and S the entropy at P and T.21 It has been well documented5,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.31 LDA with vibrational effects accounted by the QHA offers most accurate EoS for weakly correlated and strongly bonded materials. Several corrections have been designed to reconcile experiments and calculations and to accurately predict high temperature EoS.12,31–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,32 Other corrective schemes consistently rescale EoS variables to bring calculated and measured EoS parameters into close agreement.31,33 Here we extend one of these corrective schemes12 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 D F T ( P , 300 ) + Δ V ( P , 300 )
(5)

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 D F T ( P , 0 ) + Δ V ( P , 300 ) .
(6)

Therefore, corrected static volumes are obtained by subtracting the extra volume difference due to zero point motion:

V st cor ( P ) = V st D F T ( P ) + Δ V ( P , 300 ) Δ V z p ( P )
(7)

which can be written as:

V st cor ( P ) = V st D F T ( P ) + f ( P ) .
(8)

f(P) is then used to correct the static energy that enters in QHA calculations. This corrective procedure offers highly accurate hybrid EoS12 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:

C i j k l T ( P , T ) = 1 V 2 G ( P , T ) e i j e kl P
(9)

where G(P, T) is the Gibbs free energy:

G ( P , T ) = F + P V = U 0 ( V ) + U z p ( V ) + U t h ( V , T ) + P V st + P Δ V z p + P Δ V t h = G st ( P ) + G v i b ( P , T ) .
(10)

Thus, the static elastic coefficients are:

C i j k l s t ( P ) = 1 V 2 G st ( P ) e i j e kl P
(11)

Conversely, the static compliance tensor is:

M i j k l s t = 1 V 2 G st ( P ) σ i j σ kl P .
(12)

Since static volumes are calculated via the Gibbs free energy:

G st ( P ) P = V st ( P )
(13)

the corrected volume must be the derivative of a “corrected” Gibbs free energy:

G st cor ( P ) P = V st cor ( P ) = V st D F T ( P ) + f ( P )
(14)

as in Eq. (8). Integrating Eq. (14) to obtain G st cor and using Eq. (12), we arrive at the corrected compliance tensor:

M i j k l cor st ( P ) = V st ( P ) V st cor ( P ) M i j k l s t + δ i j δ kl 9 V st cor ( P ) f ( P ) P
(15)

from which the elastic coefficients can be calculated:

C i j k l cor st = M i j k l cor st 1
(16)

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 measured34,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.12 Corrective 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 (C11 = C1111, C12 = C1122 and C44 = C2323) 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, C11 is overestimated, similarly to pressure.31 The opposite occurs with GGA-PBE and HSE06. HSE06 agrees best with experiments, but still has large errors. C44 is less volume-sensitive and agrees well with experimental data.

FIG. 1.

Calculated and measured NaCl EoS. Volumes are underestimated by LDA and overestimated by GGA-PBE. Experimental data are from Boehler and Kennedy34 and Fritz et al.35 

FIG. 1.

Calculated and measured NaCl EoS. Volumes are underestimated by LDA and overestimated by GGA-PBE. Experimental data are from Boehler and Kennedy34 and Fritz et al.35 

Close modal
FIG. 2.

Calculated thermal expansivity compared to experimental data from Enck and Dommel48 and Rao et al.49 

FIG. 2.

Calculated thermal expansivity compared to experimental data from Enck and Dommel48 and Rao et al.49 

Close modal
FIG. 3.

a) Pressure dependence of calculated elastic coefficients at 300 K. Experimental data from Kinoshita et. al.,36 Voronov et. al.37 and Whitfield et. al.38 b) Temperature dependence of calculated elastic coefficients at 0 GPa. Experiments from Yamamoto et. al.39 

FIG. 3.

a) Pressure dependence of calculated elastic coefficients at 300 K. Experimental data from Kinoshita et. al.,36 Voronov et. al.37 and Whitfield et. al.38 b) Temperature dependence of calculated elastic coefficients at 0 GPa. Experiments from Yamamoto et. al.39 

Close modal

We used the corrected EoS12,40 with this corrective scheme to produce the elastic coefficients shown in Figure 4. The pressure dependence of C11, C12 and C44 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 C11, 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 corrections12,32 or fully anharmonic calculations41,42 are needed. With these elastic coefficients, bulk and shear moduli were calculated using Voigt-Reuss-Hill (VRH) averages43–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 Decker46 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.47 Similar agreement and deviations are observed between calculated and measured acoustic velocities.

FIG. 4.

a) Pressure dependence of corrected LDA elastic coefficients using at 300 K, b) and their temperature dependence at 0 GPa. Experiments from.36–39 

FIG. 4.

a) Pressure dependence of corrected LDA elastic coefficients using at 300 K, b) and their temperature dependence at 0 GPa. Experiments from.36–39 

Close modal
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.46 Those by Voronov were obtained using elastic coefficients measured on single crystals.37 The 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.

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.46 Those by Voronov were obtained using elastic coefficients measured on single crystals.37 The 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.

Close modal

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 C11 mean squared error is significantly reduced from 144 to 5, comparing to Voronov’s data.37 Data obtained on polycrystaline NaCl47 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, C11C12 = 2C44 increases with pressure. This increased anisotropy can also be seen in Fig. 5(a), where Reuss, GR, and Voigt, GV, bounds43 are shown (KR = KV). For an isotropic system, GR = GV, a condition held only at 0 GPa. With increasing pressure, GR and GV 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 EoS12 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.

We thank Raphael Ranna from Brazilian Petroleum Association (ANP), Rio de Janeiro, and Prof. Miguel A. Mane from Department of Applied Geology of Rio de Janeiro State University for suggesting calculations in NaCl. We thank Prof. Lucy Assali for helpful comments on the manuscript. MLM thanks support by CAPES from Brazil through Science without Borders program process number BEX 14456/13-3. RMW thanks support by NSF grant EAR-1341862. We acknowledge computer allocations at the Minnesota Supercomputing Institute and at the Blue Waters grant through the Great Lakes Consortium.

1.
K. E.
Bullen
,
Phys. Chem. Earth
1
,
68
(
1956
).
2.
K. E.
Bullen
and
R. A.
Haddon
,
Proc. Nat. Acad. Sci.
58
,
846
(
1967
).
3.
A. M.
Dziewonski
and
D. L.
Anderson
,
Phys. Earth Planet. Int.
25
,
297
(
1981
).
4.
D. L.
Anderson
and
A. M.
Dziewonski
,
Scientific American
251
,
60
(
1984
).
5.
R. M.
Wentzcovitch
,
Y. G.
Yu
, and
Z.
Wu
,
Reviews in Mineralogy and Geochemistry
71
,
59
(
2010
).
6.
R. M.
Wentzcovitch
,
Z.
Wu
, and
P.
Carrier.
,
Reviews in Mineralogy and Geochemistry
71
,
99
(
2010
).
7.
J.
Gluyas
and
R.
Swarbrick
,
Petroleum Geoscience
(
Blackwell Science Ltd
,
2004
)
8.
D.
Wallace
,
Thermodynamics of Crystals
(
John Willey & Sons Inc
,
New York
,
1972
).
9.
P.
Carrier
,
R. M.
Wentzcovitch
, and
J.
Tsuchiya
,
Phys. Rev. B
76
,
064116
(
2007
).
10.
P.
Giannozzi
,
S.
Baroni
,
N.
Bonini
,
M.
Calandra
,
R.
Ca
,
C.
Cavazzoni
,
D.
Ceresoli
,
G. L.
Chiarotti
,
M.
Cococcioni
,
I.
Dabo
,
A. D.
Corso
,
S.
Fabris
,
G.
Fratesi
,
S.
de Gironcoli
,
R.
Gebauer
,
U.
Gerstmann
,
C.
Gougoussis
,
A.
Kokalj
,
M.
Lazzeri
,
L.
Martin-Samos
,
N.
Marzari
,
F.
Mauri
,
R.
Mazzarello
,
S.
Paolini
,
A.
Pasquarello
,
L.
Paulatto
,
C.
Sbraccia
,
S.
Scandolo
,
G.
Sclauzero
,
A. P.
Seitsonen
,
A.
Smogunov
,
P.
Umari
, and
R. M.
Wentzcovitch
,
J.Phys.:Condens.Matter
21
,
395502
(
2009
).
11.
G.
Kresse
and
J.
Furthmuller
,
Phys. Rev. B
54
,
11169
(
1996
).
12.
M. L.
Marcondes
and
R. M.
Wentzcovitch
,
J. Appl. Phys.
117
,
215902
(
2015
).
13.
D. B.
Fraser
and
R. C.
LeCraw
,
Rev. Sci. Instrum.
35
,
1113
(
1964
).
14.
B.
Li
and
R. C.
Liebermann
,
Proc. Nat. Acad. Sci.
104
,
9145
(
2007
).
15.
Y.
Higo
,
T.
Inoue
,
B. S.
Li
,
T.
Irifune
, and
R. C.
Liebermann
,
Proc. Nat. Acad. Sci.
104
,
9145
(
2007
).
16.
D. J.
Weidner
,
K.
Swyler
, and
H. R.
Carleton
,
Geophys. Res. Lett.
2
,
189
(
1975
).
17.
J. D.
Bass
and
D. J.
Weidner
,
J. Geophys. Res. Solid Earth
89
,
4359
(
1984
).
18.
M.
Murakami
,
Y.
Ohishi
,
N.
Hirao
, and
K.
Hirose
,
Earth Planet. Sci. Lett.
277
,
123
(
2009
).
19.
S.
Speziale
,
H.
Marquardt
, and
T. S.
Duffy
,
Rev. Min. Geochem.
78
,
543
(
2014
).
20.
B. B.
Karki
,
R. M.
Wentzcovitch.
,
S.
de Gironcoli
, and
S.
Baroni
,
Science
286
,
1705
(
1999
).
21.
R. M.
Wentzcovitch
,
M. C. B. B.
Kark
, and
S.
de Gironcoli
,
Phys. Rev. Lett.
92
,
018501
(
2004
).
22.
Z.
Wu
and
R. M.
Wentzcovitch
,
Physical Review B - Condensed Matter and Materials Physics
83
,
1
(
2011
).
23.
J. P.
Perdew
and
A.
Zunger
,
Physical Review A
23
,
5048
(
1981
).
24.
D.
Ceperley
and
B. J.
Alder
,
Phys. Rev. Lett
45
,
566
(
1980
).
25.
J.
Perdew
,
K.
Burke
, and
M.
Ernzerhof
,
Physical Review Letters Phys. Rev. Lett. (USA)
77
,
3865
(
1996
).
26.
J.
Heyd
,
G. E.
Scuseria
, and
M.
Ernzerhof
,
The Journal of Chemical Physics
118
,
8207
(
2003
).
27.
P. E.
Blochl
,
Physical Review B
50
,
17953
(
1994
).
28.
M.
Nunez-Valdez
,
M. K.
Umemoto
, and
R.
Wentzcovitch
,
Appl. Phys. Lett.
101
,
170912
(
2012
).
29.
R.
Yang
and
Z.
Wu
,
Earth Planet. Sci. Lett.
404
,
14
(
2014
).
30.
G.
Shukla
,
Z.
Wu
,
H.
Hsu
,
A.
Floris
,
M.
Cococcioni
, and
R. M.
Wentzcovitch
,
Geophys. Res. Lett
42
,
1741
(
2015
).
31.
K.
Kunc
and
K.
Syassen
,
Physical Review B
81
,
134102
(
2010
).
32.
Z.
Wu
,
R. M.
Wentzcovitch
,
K.
Umemoto
,
B.
Li
,
K.
Hirose
, and
J.-C.
Zheng
,
Journal of Geophysical Research
113
,
B06204
(
2008
).
33.
Z.
Zhang
,
L.
Stixrude
, and
J.
Brodholt
,
Earth and Planetary Science Letters
379
,
1
(
2013
).
34.
R.
Boehler
and
G. C.
Kennedy
,
J. Phys. Chem. Solids
41
,
517
(
1980
).
35.
J. N.
Fritz
,
S. P.
Marsh
,
W. J.
Carter
, and
R. G.
McQueen
,
Accurate Characterization of the High Pressure Enviroment
326
,
201
(
1972
).
36.
H.
Kinoshita
,
N.
Hamaya
, and
H.
Fujisawa
,
J. Phys. Earth
27
,
337
(
1979
).
37.
F. F.
Voronov
,
E. V.
Chernyshera
, and
V. A.
Goncharova
, in
Proceedings of the Seventh Institute AIRACT Conference of High Pressure Science Technology
, edited by
B.
Vodar
and
A.
Morteau
(
Pergamon
,
1980
).
38.
C. H.
Whitfield
,
E. M.
Brody
, and
W. A.
Basset
,
Rev. Sci. Instrum.
47
,
942
(
1976
).
39.
S.
Yamamoto
,
I.
Ohno
, and
O. L.
Anderson
,
Journal of Physics and Chemistry of Solids
48
,
143
(
1987
).
41.
D.
Zhang
,
T.
Sun
, and
R. M.
Wentzcovitch
,
Phys. Rev. Lett.
24
,
058501
(
2014
).
42.
T.
Sun
,
D. B.
Zhang
, and
R. M.
Wentzcovitch
,
Phys. Rev. B
89
,
094109
(
2014
).
43.
R.
Hill
,
Proc. Phys. Soc. A
65
,
349
(
1952
).
44.
A.
Reuss
,
Z. Angew. Math. Mech.
9
,
49
(
1929
).
45.
W.
Voigt
,
Lehrbuch der Kristallphysic
(
Vieweg+Teubner Verlag
,
1928
)
46.
D. L.
Decker
,
Journal of Applied Physics
42
,
3239
(
1971
).
47.
J.
Frankel
,
F. J.
Rich
, and
C. G.
Homan
,
J. Geophys. Res
81
,
6357
(
1976
).
48.
F. D.
Enck
and
J. G.
Dommel
,
Journal of Applied Physics
36
,
839
(
1965
).
49.
A. S. M.
Rao
,
K.
Narender
,
G. K. K.
Rao
, and
N. G.
Krishna
,
Journal of Mod. Physics
4
,
208
(
2013
).