An efficient method for calculating thermoelastic properties
AAn efficient method for calculating thermoelastic properties
Zhongqing Wu Department of Chemical Engineering and Materials Science and Minnesota Supercomputing Institute, University of Minnesota, Minneapolis, MN 55455, USA First-principles quasi-harmonic calculations play a very important role in mineral physics because they can accurately predict the structure and thermodynamic properties of materials at pressure and temperature conditions that are still challenging for experiments. It also enables calculations of thermoelastic properties by obtaining the second-order derivatives of the free energies with respect to strain. However, these are exceedingly demanding computations requiring thousands of large jobs running on 10 processors each. Here we introduce a simpler approach that requires only calculations of static elastic constants and phonon density of states for unstrained configurations. This approach decreases the computational time by more than one order of magnitude. We show results on MgO and forsterite that are in very good agreement with previous first-principles results and experimental data. . Introduction Elasticity is a basic property of materials important in many processes such as brittle failure, flexure and the propagation of elastic waves. The study of the elasticity of the Earth’s materials is very important to geophysics. The primary source of information about the Earth’s interior is provided by seismological observation. Understanding these observations requires knowledge of the elastic properties of materials at high pressures and temperatures. Although experimental progress has been steady, measuring the elastic constant at the pressure and temperature conditions of the Earth’s mantle remains a considerable challenge. Recently, a complementary approach to experiments [Karki et al. et al. et al. et al. . Method
The isothermal elastic constant for the infinitesimal strain [Barron and Klein, 1965] is given by ijkl ij kl il jk ik jlij kl Fc PV e e δδ δδ δδ∂= + − −∂ ∂ (1) Here F is Helmholtz free energy, which is expressed in the quasi-harmonic approximation as
1( , ) ( ) ( ) 2 ln{1 exp[ ( ) / ]} qmq mB qm Bq m
F V T U V Vk T V k T ωω (cid:1)(cid:1) = ++ − − (cid:1) (cid:1) (2) where subscript q is the wave vector and m is index of normal modes. Usually, in order to obtain the derivative of the Helmholtz free energy with respect to all kinds of strains, we need to calculate the phonon DOS not only for a configuration without strain but also for different strain configurations, which number 15 for orthorhombic structures and 5 for cubic ones. All these calculations must be repeated for several volumes (generally >8) to get the pressure and temperature dependence of the elastic constant. Phonon calculation is time-consuming. For example, the computing time for phonon DOS are almost three orders of magnitude greater than that for the self-consistent calculation of MgO. Therefore, calculating the thermal elastic constant by the method described above will necessitate a large computation effort and may even be impracticable for materials with a large primitive cell, which is usually the case for minerals. Obviously, we can save more than one order of magnitude of computation time if we can calculate the thermal elastic constant just from the phonon DOS for configuration without strain. Phonon DOSs for configuration without strain at several volumes provides information about the volume dependence of phonon frequencies. If we can know how the frequency varies with infinitesimal strain based on its volume dependence, then we can obtain the thermal elastic onstant without performing the phonon calculations for strain configuration. Here we derived the relations between frequency dependence on volume and strain for a crystal with orthorhombic symmetry. The relation can also be applied to cubic and tetragonal crystals. The relations will be a little different for crystals with other symmetries; however, the basic idea should be still hold. The volume dependence of the frequencies is described by mode Grüneisen parameters as qm qmqm d dVV ω γω = − (3) If the volume change is only due to longitudinal strain, we also can show that the frequencies change with
1, 1 2, 2 3, 3 ( ) qm qm qm qmqm d e e e ω γ γ γω = − + + (4) where, we use the Voigt’s notation, e , e , and e is the longitudinal strain along the three crystal axes, x , x and x , for orthorhombic crystals, respectively. In order to avoid confusion, we call qm γ , qm γ , and qm γ strain Grüneisen parameters. We can either express the frequency change by the power of the volume change
12 [ ( ) 2]1 ( )( 21 ) 2 qm qmqm qm qm qmq qmm e e ee edV dVV Ve e e e e eeV V ω ωω γ ωγγ ω ω∂∂∂ ∂∆ = ⋅ + ⋅∂ ∂ + ++ += − + + − + +−+ (5) r express the frequency change by the power of the strain
2, 2 2 21, 1 2, 2 3, 3 1, 1 2, 2 3, 31, 2, 3, 1, 2, 32 2 21 2 3 1 2 2 31 2 3 2 3
12 1 1 1 ( )2 2 22 2 21 ( 2 qm qmqm i i ji i ji i jqm qm qm qm qm qm qmqm qm qm qm qmqm e e ee e ee e e e e ee e e e e e ee e e e e ω ωω ω γ γ γ γ γ γγ γ γ γ γ γω∂ ∂∆ = ⋅ + ⋅∂ ∂ ∂= − + + − − −∂ ∂ ∂ ∂ ∂ ∂− + + + + +∂ ∂ ∂ ∂ ∂ (cid:1) (cid:1) , 3 112 2 2 2 2 21, 1 2, 2 3, 3 1, 2, 1 2 2, 3, 2 3 3, 1, 3 1 )1 [( ) ( ) ( ) 2 2 2 ]2 qmqm qm qm qm qm qm qm qm qm qm e eee e e e e e e e e ω γ γ γ γ γ γ γ γ γ∂+ + + + + + (6) For any longitudinal strain, both equations should have the same frequency change. According to the fist-order term of two equations, we get one of relationships
31 21, 2, 3,1 2 3 1 2 3 1 2 3 ( ) ( ) ( ) qm qm qm qm ee ee e e e e e e e e γ γ γ γ+ + =+ + + + + + (7-1) This relation also ensures that the last term in both equations are equal each other. In order to keep the second term in both equations equal for any longitudinal strain, the corresponding coefficients in both equation should have relation as following , 2, , i qm qmji qm j qm qm
Ve V γ γγ γ γ∂ ∂∂ ∂= (7-2) Therefore we can introduce two azimuth parameters and qm qm θ φ , which we call the Grüneisen azimuth, and rewrite the equation (7-1) c o s s insin s inc o s qm qm qm qmqm qm qm qmqm qm qm e e eee e eee e ee γ γ θ φγ γ θ φγ γ φ+ += + += + += (8) ach mode m and wave vector q corresponded to a Grüneisen azimuth. All these Grüneisen azimuth will form a distribution function ( , ) f θφ . What we need in calculating elastic constant is the average value about strain Grüneisen parameters, namely something like , , ,, i qm q i qm i qmq m w f d dN γ γ γ θφ θφ φ θ φ= = (cid:1) (cid:2) (9) where, q w is the weight at wave vector q and are the number of normal modes. The distribution function ( , ) f θφ is unknown for most materials except one case, isotropic materials. For isotropic materials, ( , ) f θφ is a constant and should be equal to
1/ 4 π . The symmetry reduction seems not change the distribution function ( , ) f θφ significantly. For example, the tetragonal system, produced by applying a small strain in a cubic system, and its original cubic system should have almost the same distribution function ( , ) f θφ since the small strain does not change the strain Grüneisen parameter, which is reflected by equation (4). Therefore we assume that the uniform distribution function ( , ) f θφ is a good approximation for the crystal. Then we can calculate all kinds of average values of strain Grüneisen parameters, which are needed in calculating the thermal elastic constant. i qm qm qm qmi qm qm qmi qmi e e e dee e e dee e ee γ γ θ φπγ θ φπγ (cid:2)(cid:2)(cid:2)(cid:2) + += Ω+ + Ω+ + (cid:2) (10) Similarly, we have ( )1 ( ) if 5 ( )1 ( ) if 15 qmi ji qm j qm qmi j e e e i je ee e e i je e γγ γ γ (cid:3) + + = (cid:4)(cid:4) = (cid:5) + + (cid:4) ≠ (cid:4)(cid:6) (11) By using the relation 2 (Eq. (7-2)), we can get
21 2 3, 21 2 3 ( )1 if 5 ( )1 if 15 qmi ji qmj qmi j e e e V i je e Ve e e e V i je e V γγ γ (cid:3) ∂+ + = (cid:4) ∂∂ (cid:4) = (cid:5) ∂ ∂+ + (cid:4) ≠ (cid:4) ∂ (cid:6) (12) where i,j =1, 2, 3. Although we assume the isotropic distribution function, the average value about the strain Grüneisen parameters are anisotropy, which are reflected by the axis compressive ratio under pressure, namely e : e : e . We know the axis compressive ratio e : e : e at different pressures from first-principles static calculations. For cubic e = e = e For other crystals, they are not equal in general and change with volume but are known. The volume dependence of qm γ is known since we had calculated the frequencies at different pressures (or volume). Therefore, we can calculate all kinds of average values about strain Grüneisen parameters using the above equation. For orthorhombic crystals, there are nine elastic constants. Using Voigt’s notation, the isothermal elastic constant in equation (1) excluding the shear elastic constant can be rewritten as ij ij ij ijPHi j
F V Tc V T P V T c V c V TV e e δ (cid:7) (cid:8) ∂= + − = + (cid:9) (cid:10)(cid:9) (cid:10) ∂ ∂ (cid:11) (cid:12) (13) here i,j = 1, 2, 3. The first term
20 , 0 ( ) (1 ) ( ) ij i ji j
Uc V P Ve e δ∂= + −∂ ∂ is the elastic constant from the local-density approximation (LDA) static calculation, where ( ) P V is the pressure from the LDA static calculation. All phonon contributions to the elastic constant are included in the second term, which can be written as ( , ) ( ) ( , ) (1 ) ( , ) ijPH ij ijT ij T c V T c V c V T P V T ω δ = + + − (14) where ( , ) ( , ) ( ) T P V T P V T P V = − is the pressure resulting from the phonon contribution. Zero motion and temperature contributions to the elastic constant are described by the first and second terms of equation (14), respectively, ( ) 3( ) ( )2 2 qm i qmij i qm j qm i j i qmqm i j j
V Nc V V e e V e ω ω γγ γ δ γ∂ ∂= = − +∂ ∂ ∂ (cid:1) (cid:1) (cid:1) (15) ln(1 )( , ) e 1[ ( )](e 1) e 1 qmqmqm qm QBijT q m i jQ i qmB qm i qm j qm qm i qm j qm i qm ijQ Qq m j k T ec V T V e ek T Q QV e γγ γ γ γ γ δ − (cid:1)(cid:1) ∂ −= ∂ ∂ ∂= − ⋅ ⋅ + ⋅ ⋅ − +− − ∂ (16) where qm qm B Q k T ω= (cid:1) and 3N is the number of normal modes. e /(e 1) qm qm Q Qqm Q − and /(e 1) qm Qqm Q − are weakly frequency dependent especially at high temperature. Therefore, it is a good approximation to assume that e /(e 1) qm qm Q Qqm Q − and /(e 1) qm Qqm Q − do not correlate with , i qm γ . Then we have qmqm qm Q i qmBijT qm i qm j qm qm i qm j qm i qm ijQ Q j k TNc V T Q QV e γγ γ γ γ γ δ∂− ⋅ ⋅ + ⋅ ⋅ − +− − ∂ (cid:2) (17) he relation between the adiabatic and isothermal elastic constant is described [Davies 1974] by ( , ) ( , ) i jSij ij V
VTc V T c V T C λλ= + , (18) where i,j = 1, 2, 3 and V C is heat capacity at constant volume, and Q Qqm qmB BQ Qi qm i qm qm i qmqm qmq mi
S V T k e k N eQ QV e V e V e λ γ γ∂= = =∂ − − (cid:1) (19)
Similarly, the shear elastic constant can also be expressed as ( , ) ( ) ( , ) kk kk kkPh c V T c V c V T = + , (20) where k = 4, 5, 6. Up to now, we still do not know how to calculate the second term of the equation, the phonon contribution to the shear elastic constant, because we only discuss the longitudinal strain. However, as Figure 1 shows, the shear strain in the x and x planes can be regarded as a combined strain with a compressive strain along the x (cid:1) axis and a tensile strain along the x (cid:1) axis. Energy change are be described by ijkl ij kl E c e eV ∆ = (21) If volume is conserve in strain. In X1 and X2 coordinate, it is
E c e e c e e c e e c e eV c e e ∆ = + + += (22) In X1’ and X2’ coordinate, it is '1'1'1' 1'1' 1'1' 1'1'2'2' 1'1' 2'2' 2'2'1'1' 2'2' 1'1' 2'2'2'2' 2'2' 2'2'1'1' 1' 1' 1'2' 1' 2' 2'1' 2' 1' 2'2' 2' 2'
E c e e c e e c e e c e eV c e e c e e c e e c e e ∆ = + + += + + + (23) The strain relation in two coordinates is
1' 6 2' 6 , e e e e = = − (24) With Eq. (21) –(24), we got
66 1'1' 2'2' 1'2' ( 2 ) / 4
Ph Ph Ph Ph c c c c = + − (25-1) The equation (25) builds the relation between the shear elastic constant and longitudinal and off-diagonal elastic constant, which we already know hot to calculate. The formula used to calculated ijPh c can be applied to ' ' i j Ph c , here i,j=1.3. The only difference is that i e qm γ should be replaced by ' i e ' qm γ The strain in X , X and X i e (i=1,3) transformed into the stain in X , X and X as fellows: ; ; ; 2 2 e e e ee e e e e e e + += = = = − (26) When we use the equations in section 2.2.1. to calculate the ' ' i j Ph c by replacing e , e and e with e (cid:1) , e (cid:1) and e (cid:1) , in principle we need new volume Grüneisen parameters , ' q m γ corresponding to strain ; ; ; 02 2 e e e ee e e e e + += = = = . , ' q m γ are different from , q m γ which correspond to strain ; ; ; 2 2 e e e ee e e e e e e + += = = = − because non zero e Fortunately, the difference between , ' q m γ and , q m γ should be small because (i) in contrast to the compressive strain, the shear strain e should be relatively small as shown in equation (26); (ii) the frequency changes caused by the shear strain are far smaller than the frequency changes caused by the compressive strain ence contribution to Grüneisen parameters from the shear strain are relatively small. Therefore we can still use , q m γ to calculate the shear elastic constant without sacrificing too much precision. c and c can be calculated similarly.
3. Results and discussion
Our results for the thermal elasticity in two major minerals, cubic MgO (periclase) and orthorhombic forsterite ( (cid:2) phase of Mg SiO ), are presented here. The calculations are performed using the density functional theory with the local-density approximation [Perdew and Zunger, 1981]. Calculation details are similar to those reported in previous works [Karki et al. et al. k -points ( q -points) for periclase and forsterite, respectively. Phonon frequencies were calculated using density functional perturbation theory [Baroni 2001]. As shown in Figure 2, our calculated elastic moduli of MgO at zero pressure are almost the same as our previous calculations [Kariki et al. above 2000 K and c above 1200 -K. Both results agree well with the experimental data [Isaak et al. , 1989a] for a very wide temperature range. At temperatures above ~1000 K, the c appear to deviate more sharply from the experimental trend. The deviations result from overestimating the thermal expansivity by using the quasi-harmonic approximation (QHA) at those temperatures [Kariki et al. , although our new method adopts the approximation. The possible reason is that our method avoids the direct calculation of the second derivative of the free energy with respect to infinitesimal strain, which is sensitive to fitting conditions such as hich equation of state is adopted and how many points are used in fitting. In contrast, our results are insensitive to fitting conditions. Furthermore the approximation should work very well for cubic crystals. The elastic moduli of forsterite at zero pressure also are consistent with the experimental data (see Fig. 3). At room temperature, the longitudinal elastic constants c , c and c are slightly smaller than the experimental data [Issak et al. , 1989b]. However, their temperature derivatives are in excellent agreement with the experimental results. At temperatures above ~1000 K, the longitudinal elastic constant of forsterite decreases with temperature more slowly than the experimental data. This similar temperature behavior also appears in the bulk modulus, K s, (see Fig. 3). The bulk moduli are calculated without using any of the approximations mentioned above, the deviation from experimental data should also be attributed to the anharmonic effect as in the case of MgO. However, MgO and forsterite show completely different anharmonic behavior. The longitudinal elastic constant of MgO decreases more quickly than the experimental data. On the contrary, the longitude elastic constant of forsterite decrease more slowly than the experimental data. This behavior is consistent with the thermodynamic observation of anharmonicity. The thermal expansivity is overestimated in MgO but underestimated in forsterite by the QHA [Li et al. V C of forsterite exceeds the Dulong and Petit limit at about 1400 K [Gillet 1991, Bouhifd 1996 , Anderson 1996)]. In contrast, the heat capacity V C of MgO is predicted to be 2.6% smaller than the Dulong and Petit limit even at T = 2000 K [Anderson and Zou, 1990]. Our calculations can also reproduce high-pressure experimental data well. Figure 4 shows the dependence of velocity of on pressure at various temperatures for forsterite. Our first principles static calculation is in agreement with previous results by da Silva et al. [1997]. Both alculations overestimate the compressive velocity. After including zero motion and the temperature contributions in the elastic constant, the velocity becomes consistent with the experiment observations [Zha et al. et al. et al.
4. Conclusion
We derived the relation between the strain Grüneisen parameters and volume Grüneisen parameters, from which we introduce a Grüneisen azimuths for each phonon mode. By assuming the Grüneisen azimuths have a homogenous distribution, we use the formulas to calculate all kinds of average values of the strain Grüneisen parameters just requiring the axis compressive ratio e : e : e and volume Grüneisen parameters qm γ at different pressures, which can be obtained from first-principles static calculations and phonon DOS of unstrain configurations. Then we can calculate the thermal elastic constant without requiring the phonon DOS of the strain configuration, which reduces the computing time by more than an order of magnitude compared with the previous method. This kind of saving in computing time is very useful in investigating the elasticity of minerals because the phonon DOS calculation is very demanding and usually becomes too unwealdy for minerals because of their large primitive cell size. The calculated elasticity on periclase and forsterite shows that our method includes the phonon contribution well both at zero pressure and at high pressure and works well both for cubic and orthorhombic crystals. cknowledgements References
Anderson, O.L. 1996 Anharmonicity of forsterite and the thermal pressure of insulators.
Geophys. Res. Lett . J. Phys. Chem. Ref. Data , 69–83 Baroni, S., de Gironcoli, S., Dal Corso, A. and Giannozzi, P. 2001 Phonons and related crystal properties from density-functional perturbation theory, Rev. Mod. Phys.
Proc. Phys. Soc.
Geophys. Res. Lett.
Geophys. Res. Lett . , 1963–1966 Davies, G.F. 1974 Effective elastic moduli under hydrostatic stress—I. quasi-harmonic theory, J. Phys. Chem. Solids , 1513–1520 Gillet, P, Richet, P, Guyot, F, Fiquet, G. 1991 High-temperature thermodynamic properties of forsterite. J. Geophys. Res.
Phys. Chem. Miner . , 704 saak, G., Anderson, O.L. and Goto, T. 1989b Elasticity of single-crystal forsterite measured to 1700 K, J. Geophys. Res. , 5895 Karki, B.B., Wentzcovitch, R.M., De Gironcoli, S. and Baroni, S. 1999 High-pressure lattice dynamics and thermoelasticity of MgO, Phys. Rev. B , 8793–8800 Li, B., Kung, J., Liebermann, R.C. 2004 Modern techniques in measuring elasticity of Earth materials at high pressure and high temperature using ultrasonic interferometry in conjunction with synchrotron X-radiation in multi-anvil apparatus. Phys. Earth Planet. Inter. , 559–574 Li, L., Wentzcovitch, R.M., Weidner, D.J. and da Silva, C. 2007 Vibrational and thermodynamic properties of forsterite at mantle conditions,
J. Geophys. Res. , B05206 Perdew, J.P., and Zunger, A. 1981 Self-interaction correction to density-functional approximations for many-electron systems,
Phys. Rev. B , postperovskite at D'' conditions Proc. Natl. Acad. Sci. USA (3), 543–546 Wentzcovitch, R.M., Karki, B.B., Cococcioni, M. and de Gironcoli, S. 2004 Thermoelastic properties of MgSiO -perovskite: Insights on the nature of the Earth's lower mantle, Phys. Rev. Lett.
J. Geophys. Res . , 17535–17545 igure Captions Figure 1.
The strain relations in different coordinate axis. The structure before and after strain are denoted by green and red lines, respectively. The structure experiences a shear strain as shown by the solid line if viewing from the x and x axis. This shear strain can also be viewed as the combined strain with a compressive strain along the x (cid:1) axis and a tensile strain along the x (cid:1) axis as shown by the dashed lines. Figure 2.
The temperature dependence of the elastic constant of MgO at zero pressure (solid line), compared to previous results [Karki et al. et al.
Figure 3.
The temperature dependence of the elastic modulus of forsterite at zero pressure (solid line), compared to experimental data [Issak et al.
Figure 4.
The pressure dependence of the velocity of forsterite at various temperatures, compared with experiment data from Issak et al. [1989b] (stars), Zha et al. [1996] (open circles) and Li et al. [2003] (solid circles). Fig .1 X X (cid:1) X X (cid:1) E l a s t i c c on s t an t ( G P a ) Temperature (K)
MgO C C C Fig. 2
500 1000 1500 2000 250040608010012050150200250300350400 c c E l a s t i c m odu l u s ( G P a ) c c c c Temperature (K) Ks c c c Fig. 3 V e l o c i t y ( k m / s ) Pressure (GPa) static result 300K 1070K 1073K 1100K 2000k Li (2004) Issak (1989) Zha(1996)static result 300K 1070K 1073K 1100K 2000k Li (2004) Issak (1989) Zha(1996)