Thermo-Elasticity of Materials from Quasi-Harmonic Calculations

Abstract

An effective algorithm for the quasi-harmonic calculation of thermo-elastic stiffness constants of materials is discussed and implemented into the Crystal program for quantum-mechanical simulations of extended systems. Two different approaches of increasing complexity and accuracy are presented. The first one is a quasi-static approximation where the thermal dependence of elastic constants is assumed to be due only to the thermal expansion of the system. The second one is fully quasi-harmonic, takes into account thermal expansion, and explicitly computes Helmholtz free energy derivatives with respect to strain. The conversion of isothermal into adiabatic thermo-elastic constants is also addressed. The algorithm is formally presented and applied to the description of the thermo-elastic response of the forsterite mineral.

1. Introduction

In the last several decades, computational modeling has rapidly gained in popularity as a powerful complementary tool to the experimental characterization of materials. Taking a step forward, computational materials design is an emerging field of research where simulations guide the process of discovering new materials or engineering specific materials’ properties [1]. In this respect, the density functional theory (DFT) arguably represents the method-of-choice because of an ideal balance between its high accuracy (due to its rigorous quantum-mechanical derivation) and modest computational cost, particularly so if compared to other quantum-mechanical methodologies of comparable accuracy [2].

The state-of-the-art of computational materials design has recently been reviewed, the progresses made in the last decade illustrated and future challenges outlined [3]. Among others, recent progresses in the field of computational materials design include: (i) the ability of some software programs for quantum-mechanical simulations to efficiently exploit the parallel computing capabilities of high-performance computing (HPC) centers, thus allowing for the study of larger (more realistic) structural models [4,5,6,7,8,9,10]; (ii) the possibility of effectively including weak, dispersive interactions into the DFT formalism, thus allowing for the reliable simulation of properties of weakly-bound systems (such as molecular crystals) [11,12,13,14,15,16].

Many of the challenges current research on computational materials design is facing are related to getting a more accurate description of the electronic structure of the system (spin–orbit coupling, electron–phonon coupling, electronic transport properties, etc.) and of the conditions on which materials function (anharmonic thermal effects, phonon–phonon coupling, interaction with electric and magnetic fields, etc.) [1,2,3]. In particular, the effective inclusion of thermal effects (fully neglected by static DFT calculations) is much needed in order to widen the application domain and enhance the predictiveness of the computational approach.

In this context, let us briefly review some of the progresses made in the last few years in the development of the Crystal program for quantum-mechanical simulations of the condensed matter [17,18]. A massively-parallel version of the code has been developed, which was recently shown to efficiently scale up to 32,000 cores and which has been used to study several large systems (the largest one containing about 14,000 atoms/cell) [19]. London-type dispersive interactions are included with Grimme’s a posteriori correction (namely, DFT-D3) [20,21]. A fully-automated algorithm has been implemented to compute several quasi-harmonic thermal properties of materials (including thermal expansion, thermal dependence of the bulk modulus, and constant-pressure thermodynamic properties) [22,23,24,25].

Within the quasi-harmonic approximation, the thermal dependence of the bulk modulus can be determined at a relatively low computational cost by computing the lattice dynamics of the system at expanded and compressed space group symmetry-preserving configurations. However, the bulk modulus is just an average mechanical feature of the system. A more complete description of the anisotropic mechanical response of a material requires the fourth-rank elastic tensor to be determined. The evaluation of thermo-elastic constants (i.e., the thermal dependence of the elements of the elastic tensor) is a much more complicated task from a computational point of view as it requires phonons to be computed at several space group symmetry-breaking strained configurations. Nevertheless, the quasi-harmonic approximation [26,27] provides a formal framework for the calculation of thermo-elastic constants [28,29,30,31].

In this paper, we present and discuss an algorithm for the quasi-harmonic calculation of thermo-elastic constants of materials, which we have implemented in the Crystal program. Thermo-elastic constants are computed in two steps: (i) the equilibrium structure of the system at a given temperature T is determined from the quasi-harmonic description of thermal expansion; (ii) second derivatives of the Helmholtz free energy with respect to strain are computed numerically. Besides a fully quasi-harmonic scheme, a simplified quasi-static approach is also implemented and discussed, where the thermal dependence of elastic constants is assumed to be due entirely to the thermal expansion of the lattice. The quasi-static approach is characterized by a much lower computational cost than the fully quasi-harmonic one but only takes into account a fraction of the thermal effects on the mechanical response. The conversion of isothermal thermo-elastic constants into adiabatic ones is also discussed within the framework of the quasi-harmonic approximation.

We apply our quasi-harmonic algorithm to the description of the thermo-elastic response of forsterite, $\alpha$-Mg₂SiO₄: an end-member of the olivine solid solution series that is one of the most abundant silicates in the upper mantle of the Earth. The single-crystal thermo-elasticity of forsterite has been accurately determined experimentally at several temperatures from 300 K to 1700 K so that it represents an ideal system to validate and discuss our methodology [32].

2. Theoretical Aspects

2.1. The Athermal Elastic Tensor

We treat the crystal as a homogeneous anisotropic medium and we consider it in a reference state in which the three fundamental direct lattice vectors are ${\mathbf{a}}_1^0$, ${\mathbf{a}}_2^0$ and ${\mathbf{a}}_3^0$ so that the corresponding cell volume is $V_0=$ ${\mathbf{a}}_1·({\mathbf{a}}_2\times {\mathbf{a}}_3)$. In standard density functional theory (DFT) calculations, the reference state is the fully-optimized structure that minimizes the static internal energy of the system E at the absolute zero with no zero-point motion corrections.

The elasticity of the system is given by the stress–strain relation. Both stress and pure strain (i.e., without rotational components) are described by symmetric second-rank tensors. The symmetric pure strain tensor $\mathit{\eta }$ is dimensionless and defines a strained lattice configuration where the deformed lattice vectors ${\mathbf{a}}_{\alpha }$ are expressed in terms of the undeformed reference ones ${\mathbf{a}}_{\alpha }^0$ (with $\alpha$ = $1,2,3$) as follows:

$$a_{\alpha k}=\sum _{l=1}^3({\delta }_{kl}+{\eta }_{kl})a_{\alpha l}^0,$$

(1)

where $k,l$ = $x,y,z$ are Cartesian indices, $a_{\alpha k}$ are Cartesian components of lattice vector ${\mathbf{a}}_{\alpha }$, and ${\delta }_{kl}$ is the Kronecker delta. The elements of the stress tensor $\mathit{\sigma }$ are defined as energy density derivatives with respect to strain components and, through inversion of Equation (1), can be expressed in terms of energy gradients with respect to Cartesian components of lattice parameters:

$${\sigma }_{ij}=\frac{1}{V_0}{\left.\frac{\partial E}{\partial {\eta }_{ij}}\right|}_{\mathit{\eta }=\mathbf{0}}=\frac{1}{V_0}\sum _{\alpha =1}^3{a}_{\alpha j}^0{\left.\frac{\partial E}{\partial a_{\alpha i}}\right|}_{{\mathbf{a}}_{\alpha }={\mathbf{a}}_{\alpha }^0},$$

(2)

where $i,j=x,y,z$ are Cartesian indices. The stress tensor can thus be evaluated from energy gradients with respect to the cell parameters [33]. The stress–strain relation for an anisotropic crystal is given by the generalized Hooke’s law:

$${\sigma }_{ij}=\sum _{k,l}{C}_{ijkl}{\eta }_{kl},$$

(3)

which relates the second-rank stress tensor to the second-rank pure strain tensor through the elastic stiffness constants $C_{ijkl}$ that define the fourth-rank elastic tensor $\mathcal{C}$. Therefore, the elastic constants are defined in terms of derivatives of stress components with respect to strain components:

$$C_{ijkl}={\left.\frac{\partial {\sigma }_{ij}}{\partial {\eta }_{kl}}\right|}_{\mathit{\eta }=\mathbf{0}}.$$

(4)

By casting the first equality of Equation (2) into Equation (4), we get the following expression for the second-order elastic constants in terms of second energy density derivatives with respect to pairs of strain components:

$$C_{ijkl}=\frac{1}{V_0}{\left.\frac{{\partial }^{2}E}{\partial {\eta }_{ij}\partial {\eta }_{kl}}\right|}_{\mathit{\eta }=\mathbf{0}}.$$

(5)

As the pure strain tensor $\mathit{\eta }$ is symmetric, there are just six independent components to be considered. It is convenient to adopt Voigt’s notation, which maps the two Cartesian index notation $\left(ij\right)$ or $\left(kl\right)$ into a one index notation $\left(v\right)$, where $v=1\equiv xx$, $v=2\equiv yy$, $v=3\equiv zz$, $v=4\equiv yz$, $v=5\equiv xz$, and $v=6\equiv xy$ [34]. The elastic stiffness tensor can thus be expressed as a 6 × 6 matrix:

$$C_{vu}=\frac{1}{V_0}{\left.\frac{{\partial }^{2}E}{\partial {\eta }_v\partial {\eta }_u}\right|}_{\mathit{\eta }=\mathbf{0}},$$

(6)

where $v,u=1,\cdots ,6$ are Voigt indices. From Equation (6), the equivalence $C_{vu}\equiv C_{uv}$ is evident, which makes the 6 × 6 elastic stiffness matrix symmetric and reduces the number of its independent components from 36 to 21. Most quantum-mechanical solid-state programs can not compute the second energy derivatives in Equation (6) analytically, which thus need to be evaluated numerically either as first derivatives of analytical forces or as second derivatives of the energy. Let us briefly discuss these two approaches below.

2.1.1. Calculation from Analytical Forces

If first energy derivatives with respect to strain, $\partial E/\partial {\eta }_v$, can be computed analytically, all the elastic stiffness constants $C_{vu}$ can be determined by applying to the lattice only the six independent strain components as described below. The lattice is deformed according to the u-th independent component, ${\eta }_u$, of the strain and the analytical energy gradients $\partial E/\partial {\eta }_v$ computed at the strained lattice configurations (for $v=1,\cdots ,6$) so that the u-th column of the elastic tensor is evaluated from numerical finite differences (via a double-sided formula):

$$C_{vu}=\frac{1}{V_0}\left[{\left.\frac{\partial E}{\partial {\eta }_v}\right|}_{{\eta }_u=\delta }-{\left.\frac{\partial E}{\partial {\eta }_v}\right|}_{{\eta }_u=-\delta }\right]/2\delta ,$$

(7)

where $\delta$ represents the amplitude of the applied strain. In the Crystal program, the energy gradients with respect to strain are computed analytically [35,36] and the strain amplitude has been optimized to the default value of 1% [37,38]. The nuclear-relaxation contribution to the elastic tensor can be evaluated by either relaxing atomic positions at strained configurations or computing the internal-strain Hessian tensor [39,40]. Thus, if analytical energy gradients can be exploited, only the fundamental independent strain components need to be applied to the lattice. The symmetry of the lattice further reduces the number of components to be applied [34].

2.1.2. Calculation from the Energy

If the elastic constants in Equation (6) are to be evaluated only from the energy of different lattice configurations (i.e., from numerical second derivatives), one can not compute all of them by only deforming the lattice according to the six independent components of the strain ${\eta }_v$. More general strains are required that can be expressed as combinations of the fundamental ones. Let us introduce the following Voigt’s vector notation for the independent strain components:

$${\underline{\eta }}_1=\left(100000\right){\underline{\eta }}_2=\left(010000\right){\underline{\eta }}_3=\left(001000\right),$$

$${\underline{\eta }}_4=\left(000100\right){\underline{\eta }}_5=\left(000010\right){\underline{\eta }}_6=\left(000001\right).$$

These vectors clearly constitute a complete orthonormal basis. In order to ensure a unitary norm for all the independent components, their 3 × 3 Cartesian representations are:

$$\left(\begin{array}{ccc}1& 0& 0\\ 0& 0& 0\\ 0& 0& 0\end{array}\right)\left(\begin{array}{ccc}0& 0& 0\\ 0& 1& 0\\ 0& 0& 0\end{array}\right)\left(\begin{array}{ccc}0& 0& 0\\ 0& 0& 0\\ 0& 0& 1\end{array}\right)$$

$$\left(\begin{array}{ccc}0& 0& 0\\ 0& 0& \frac{1}{2}\\ 0& \frac{1}{2}& 0\end{array}\right)\left(\begin{array}{ccc}0& 0& \frac{1}{2}\\ 0& 0& 0\\ \frac{1}{2}& 0& 0\end{array}\right)\left(\begin{array}{ccc}0& \frac{1}{2}& 0\\ \frac{1}{2}& 0& 0\\ 0& 0& 0\end{array}\right)$$

Any general strain $\underline{\eta }$ can thus be defined as a linear combination of the fundamental strains above:

$$\underline{\eta }=\sum _v{k}_v{\underline{\eta }}_v,$$

(8)

where $k_v$ are the coefficients of the linear combination defining the strain shape. The amplitude of the strain is defined by an additional parameter $\delta$. Let us introduce an example to make this notation clear. We consider a strain $\underline{\eta }={\underline{\eta }}_1+2{\underline{\eta }}_2+{\underline{\eta }}_4$. The corresponding 3 × 3 tensor $\mathit{\eta }$ is:

$$\underline{\eta }\to \left(\begin{array}{ccc}1& 0& 0\\ 0& 2& \frac{1}{2}\\ 0& \frac{1}{2}& 0\end{array}\right)\mathrm{and}\delta \underline{\eta }\to \left(\begin{array}{ccc}\delta & 0& 0\\ 0& 2\delta & \frac{\delta }{2}\\ 0& \frac{\delta }{2}& 0\end{array}\right).$$

(9)

Once a given strain shape $\underline{\eta }$ is defined, the lattice is deformed accordingly for different values of the strain amplitude $\delta$ and the energy of the system computed at each distorted lattice configuration: $E(\delta ;\underline{\eta })$. The energy is thus a function of $\delta$ and parametrically depends on the strain shape $\underline{\eta }$. The energy computed for different values of $\delta$ can then be fitted to a polynomial function (of second-order or higher) and the corresponding fitting parameters $c_0,c_1,c_2,\cdots$ determined:

$$E(\delta ;\underline{\eta })=c_0+c_1\delta +c_2{\delta }^2+c_3{\delta }^3+\cdots$$

(10)

The second energy density derivative with respect to the strain amplitude of the expression above, at the equilibrium configuration, is obtained in terms of the $c_2$ fitting parameter and corresponds to a linear combination of elastic stiffness constants as follows:

$$\frac{1}{V_0}{\left.\frac{{\partial }^{2}E(\delta ;\underline{\eta })}{\partial {\delta }^2}\right|}_{\delta =0}=\frac{2c_2}{V_0}=\sum _v{k}_v^2{C}_{vv}+2\sum _v\sum _{u>v}{k}_v{k}_u{C}_{vu}.$$

As an example, for the particular strain shape $\underline{\eta }$ in Equation (9), it gives:

$$\frac{1}{V_0}{\left.\frac{{\partial }^{2}E(\delta ;\underline{\eta })}{\partial {\delta }^2}\right|}_{\delta =0}=C_{11}+4C_{22}+C_{44}+4C_{12}+2C_{14}+4C_{24}.$$

In order to determine the value of all the elastic stiffness constants $C_{vu}$, different strain shapes have to be applied, each providing the value of a linear combination of them, until a linear system of equations can be solved to get all the individual values. Let us note that each diagonal elastic stiffness constant $C_{vv}$ can be obtained as such by applying the corresponding fundamental strain ${\underline{\eta }}_v$. Let us further note that the symmetry of the lattice can significantly reduce the number of strain shapes to be actually applied. For instance, for a cubic lattice, only three strain shapes are necessary: ${\underline{\eta }}_1$ that gives $C_{11}$, ${\underline{\eta }}_4$ that gives $C_{44}$, and ${\underline{\eta }}_1+{\underline{\eta }}_2$ that gives $2C_{11}+2C_{12}$ to determine the three independent elastic stiffness constants $C_{11}$, $C_{12}$ and $C_{44}$.

2.2. The Elastic Tensor at a Finite Temperature

The inclusion of thermal effects on computed properties of materials requires the description of the lattice dynamics. The commonly adopted harmonic approximation (HA) for the description of lattice vibrations—as combined with the laws of statistical thermodynamics through the definition of a canonical partition function—is rather effective in the evaluation of constant-volume thermal properties of solids. However, the HA is unable to describe volume-dependent thermal properties so that thermal expansion would be null, constant-pressure thermodynamic functions would coincide with constant-volume ones, and the elastic response of the system would be independent of temperature [41,42]. In order to go beyond the above-mentioned limitations of the HA, the so-called quasi-harmonic approximation (QHA) can be invoked, which, due to its formal simplicity and relatively low computational cost, represents an effective method of choice for describing many volume-dependent thermal features of solids [22,23,26,27,42,43].

In order to compute isothermal elastic stiffness constants at a finite temperature T, Equation (6) needs to be generalized as follows:

$$C_{vu}^{\mathrm{T}}\left(T\right)=\frac{1}{V\left(T\right)}{\left.\frac{{\partial }^{2}F}{\partial {\eta }_v\partial {\eta }_u}\right|}_{\mathit{\eta }=\mathbf{0}},$$

(11)

where the reference lattice configuration is now the equilibrium one at temperature T, whose lattice parameters are ${\mathbf{a}}_{\alpha }\left(T\right)$, with $\alpha =1,2,3$, and whose corresponding volume is $V\left(T\right)$, and where the second derivatives with respect to strain are now expressed in terms of the Helmholtz free energy F of the system and no longer in terms of the static internal energy E. The evaluation of Equation (11) can thus be performed in two steps:

• The determination of the equilibrium structure of the system at temperature T;
• The calculation of the second free-energy derivatives with respect to strain.

We discuss these two steps in Section 2.2.1 and Section 2.2.2, respectively. Finally, we stress that the superscript T in Equation (11) is needed to distinguish isothermal $C_{vu}^{\mathrm{T}}\left(T\right)$ from adiabatic $C_{vu}^{\mathrm{S}}\left(T\right)$ elastic constants (see Section 2.2.3 below for more details on this).

2.2.1. The Equilibrium Structure at a Finite Temperature

The QHA introduces the explicit volume dependence of phonon frequencies by retaining the simple harmonic expression for the Helmholtz free energy of the system [26,27]:

$$F^{\mathrm{QHA}}(V;T)=U\left(V\right)+k_{B}T\sum _{\mathbf{k}p}\left[ln\left(1-e^{-\frac{\hslash {\omega }_{\mathbf{k}p}\left(V\right)}{k_{B}T}}\right)\right],$$

(12)

where $k_B$ is Boltzmann’s constant, $U\left(V\right)=E\left(V\right)+E_{\mathrm{zp}}\left(V\right)$ is the zero-temperature internal energy of the crystal, which includes the zero-point energy $E_{\mathrm{zp}}\left(V\right)={\sum }_{\mathbf{k}p}\hslash {\omega }_{\mathbf{k}p}\left(V\right)/2$, and where the sum runs over phonons, whose vibration frequencies are ${\omega }_{\mathbf{k}p}$. The Helmholtz free energy in Equation (12) can thus be evaluated from harmonic phonon frequencies computed at several volumes. The equilibrium volume $V\left(T\right)$ of the system at a given temperature T is then obtained by minimizing the isothermal Helmholtz free energy of Equation (12) with respect to volume by keeping T as a fixed parameter. From the volume–temperature relation $V\left(T\right)$, the volumetric thermal expansion coefficient ${\alpha }_V\left(T\right)$ can be defined as:

$${\alpha }_V\left(T\right)=\frac{1}{V\left(T\right)}{\left.\frac{\partial V\left(T\right)}{\partial T}\right|}_P,$$

(13)

where the derivative is taken at constant pressure. The isothermal bulk modulus of the system, $K^{\mathrm{T}}\left(T\right)$, can also be obtained as an isothermal second derivative of Equation (12) with respect to the volume:

$$K^{\mathrm{T}}\left(T\right)=V\left(T\right){\left.\frac{{\partial }^2{F}^{\mathrm{OHA}}(V;T)}{\partial V^2}\right|}_T.$$

(14)

A fully-automated algorithm for computing the thermal expansion of the system (and derived properties) within the QHA has recently been implemented in the Crystal program [18,22,23], which relies on computing and fitting harmonic vibration frequencies at different compressed and expanded volumes (four volumes, including the equilibrium one, being enough for most systems) after having performed volume-constrained structural relaxations. This implementation of the QHA has been successfully applied to study several thermal properties (structural, thermodynamic, electronic) of various systems with different chemical features in the last couple of years: fully-covalent diamond [22], ionic MgO and CaO [23], LiF [44], mixed ionic/covalent corundum $\alpha$-Al₂O₃ [24], forsterite $\alpha$-Mg₂SiO₄ [25], calcium stannate CaSnO₃ [45], and the molecular crystals of urea, purine and carbamazepine [43,46,47].

2.2.2. Free Energy Derivatives with Respect to Strain

From the methodology outlined in Section 2.2.1, the equilibrium structure of the system at temperature T is determined, whose volume is $V\left(T\right)$. In order to evaluate the isothermal elastic constants in Equation (11), the second derivatives of the Helmholtz free energy given in Equation (12) with respect to strain have to be computed. Such derivatives can not be computed analytically but rather need to be evaluated numerically from the free energy of strained lattice configurations. A scheme similar to that described in Section 2.1.2 can be adopted, where a strain shape $\underline{\eta }$ is selected and the lattice distorted accordingly for different values of the strain amplitude $\delta$. At each strained configuration, the harmonic vibration frequencies have to be computed in order to get the corresponding free energy: $F(\delta ;\underline{\eta },T)$. The free energy is thus a function of $\delta$ and parametrically depends on the strain shape $\underline{\eta }$ and on the temperature T. The free energy computed for different values of $\delta$ can be fitted to a polynomial function and the corresponding fitting parameters $c_0,c_1,c_2,\cdots$ determined:

$$F(\delta ;\underline{\eta },T)=c_0+c_1\delta +c_2{\delta }^2+c_3{\delta }^3+\cdots$$

(15)

The second free energy density derivative with respect to the strain amplitude of the expression above, at the equilibrium configuration at temperature T, is given by $2c_2/V\left(T\right)$ and corresponds to a linear combination of isothermal elastic stiffness constants as follows:

$$\frac{1}{V\left(T\right)}{\left.\frac{{\partial }^{2}F(\delta ;\underline{\eta },T)}{\partial {\delta }^2}\right|}_{\delta =0}=\sum _v{k}_v^2{C}_{vv}^{\mathrm{T}}\left(T\right)+2\sum _v\sum _{u>v}{k}_v{k}_u{C}_{vu}^{\mathrm{T}}\left(T\right).$$

(16)

As discussed in Section 2.1.2, in order to get the full set of elastic constants, several strain shapes have to be applied, each providing a linear combination of them. In the present implementation, the evaluation of the derivatives in Equation (16) with respect to the strain amplitude $\delta$ is fully-automated and requires a single run of the program per each temperature T and per each strain shape $\underline{\eta }$.

2.2.3. Adiabatic versus Isothermal Elastic Moduli

From an experimental point of view, elastic constants at finite temperature are generally measured in adiabatic rather than isothermal conditions, for instance through Brillouin scattering experiments. Adiabatic elastic moduli are always larger in value than their isothermal counterparts. In this section, we are going to illustrate how the isothermal bulk modulus and elastic constants can be transformed into adiabatic ones by means of quasi-harmonic calculations. At any temperature T, the adiabatic bulk modulus is given in terms of the isothermal bulk modulus as follows:

$$K^{\mathrm{S}}\left(T\right)=K^{\mathrm{T}}\left(T\right)\left[1+\frac{T{\alpha }_V^2\left(T\right)\overline{V}\left(T\right)K^{\mathrm{T}}\left(T\right)}{{\overline{C}}_V\left(T\right)}\right],$$

(17)

where $K^{\mathrm{T}}\left(T\right)$ is the isothermal bulk modulus defined in Equation (14), ${\alpha }_V\left(T\right)$ is the volumetric thermal expansion coefficient introduced in Equation (13), ${\overline{C}}_V\left(T\right)$ is the molar constant-volume specific heat of the system, and $\overline{V}\left(T\right)$ is the molar volume of the system. All the quantities entering Equation (17) can be computed with the HA and QHA. Adiabatic elastic stiffness constants are related to isothermal ones through the following relation [28]:

$$C_{vu}^{\mathrm{S}}\left(T\right)=C_{vu}^{\mathrm{T}}\left(T\right)+\frac{T\overline{V}\left(T\right){\lambda }_v\left(T\right){\lambda }_u\left(T\right)}{{\overline{C}}_V\left(T\right)},$$

(18)

where the ${\lambda }_v\left(T\right)$ are defined as the isothermal derivatives of the entropy S with respect to strain or, equivalently, as the derivative of the stress with respect to temperature at constant strain, and, by exploiting the chain rule of differentiation:

$${\lambda }_v\left(T\right)={\left.\frac{\partial S}{\partial {\eta }_v}\right|}_T=-{\left.\frac{\partial {\sigma }_v}{\partial T}\right|}_{\eta }=-\sum _u{\left.\frac{\partial {\eta }_u}{\partial T}\right|}_{\sigma }{\left.\frac{\partial {\sigma }_v}{\partial {\eta }_u}\right|}_T=-\sum _u{\alpha }_u\left(T\right)C_{vu}^{\mathrm{T}}\left(T\right),$$

(19)

where $\mathit{\alpha }$ is the thermal expansion tensor [48]. For cubic, hexagonal, trigonal, tetragonal and orthorhombic crystals (i.e., those lattices where angles among fundamental lattice vectors are fixed by symmetry), we have:

$$\begin{array}{cc}{\alpha }_u\left(T\right)=\frac{1}{|{\mathbf{a}}_u\left(T\right)|}{\left.\frac{\partial |{\mathbf{a}}_u\left(T\right)|}{\partial T}\right|}_{P=0},\hfill & \mathrm{for}u=1,2,3,\\ {\alpha }_u\left(T\right)=0,\hfill & \mathrm{for}u=4,5,6.\end{array}$$

More complicated expressions for the elements of the thermal expansion tensor are obtained for lower-symmetry systems. In particular, the expression for the most general triclinic case can be found elsewhere [49,50].

2.2.4. The Quasi-Static Approximation to Thermo-Elastic Moduli

As we have discussed above, the explicit quasi-harmonic calculation of thermo-elastic constants, as in Equation (11) requires the evaluation of two quantities: (i) the thermal expansion of the system $V\left(T\right)$; and (ii) the second free energy derivatives with respect to strain. From a computational point of view, the latter is much more expensive than the former as it involves several phonon calculations at strained lattice configurations. A simplified way to obtain thermo-elastic constants has been suggested, which is based on the assumption that the thermal dependence of the elastic moduli is mostly due to thermal expansion [51,52,53]. Within this methodology, referred to as “quasi-static approximation” (QSA), the thermo-elastic constants are computed as follows:

$$C_{vu}^{\mathrm{T}}\left(T\right)=\frac{1}{V\left(T\right)}{\left.\frac{{\partial }^{2}E}{\partial {\eta }_v\partial {\eta }_u}\right|}_{\mathit{\eta }=\mathbf{0}},$$

(20)

where the thermal expansion of the system $V\left(T\right)$ still has to be determined from the quasi-harmonic approximation, but the energy derivatives with respect to strain are computed in terms of static internal energy E. This approximation drastically reduces the computational cost of thermo-elastic calculations and often allows to take into account a significant fraction of thermal effects on the elastic response of the system.

3. The Implemented Algorithm

In this section, we sketch the algorithm that we have devised to compute quasi-harmonic thermo-elastic constants of materials. In particular, we explicitly illustrate the sequence of calculations that are required by stressing what steps can be performed automatically. The algorithm that we propose is the following:

• A full structural relaxation of the system is performed (both atomic positions and lattice parameters are optimized). The static equilibrium structure, with volume $V_0$, is obtained.
• A space group symmetry-preserving QHA calculation is performed, which provides the thermal expansion of the system. A fully-automated algorithm is implemented in the Crystal program to perform this task [22,23], where four different volumes are explored (compressed and expanded with respect to $V_0$). For each volume, a volume-constrained, lattice symmetry-preserving structural relaxation is performed and phonon frequencies computed. By minimizing $F(V;T)$ at several temperatures, the $V\left(T\right)$ relation is determined. As a result, lattice parameters as a function of temperature are obtained.
• A value of temperature T is selected. Starting from the values of the lattice parameters at this temperature obtained at the end of the previous step, a volume-constrained, lattice symmetry-preserving structural relaxation is performed to get the equilibrium structure (also in terms of atomic positions) at the desired temperature.
• A given strain shape $\underline{\eta }$ is chosen, which will provide a linear combination of elastic stiffness constants according to Equation (16).
• The second free energy derivatives with respect to the strain are computed as discussed in Section 2.2.2. A fully-automated algorithm has been implemented in the Crystal program for this task. The starting point is represented by the optimized structure obtained at the end of step 3 above (i.e., the equilibrium structure at temperature T). The structure is deformed, in terms of the strain shape $\underline{\eta }$, into four strained configurations (two with positive and two with negative strain amplitude $\delta$). At each strained configuration, atomic positions are relaxed and phonon frequencies computed. The computed quasi-harmonic free energy as a function of strain amplitude is fitted to a second-order polynomial and the corresponding second-derivative determined.

4. Computational Parameters

All the calculations reported in this study have been performed with algorithms implemented in a developmental version of the Crystal17 program [18]. An all-electron atom-centered Gaussian-type-function (GTF) basis set is used for all atoms; silicon, oxygen and magnesium atoms are described by (8s) (6311$sp$) (1d), (8s) (411$sp$) (1d) and (8s) (511$sp$) (1d) contractions of primitive GTFs, respectively. The exponents (in bohr⁻² units) of the most diffuse $sp$ shells are 0.32 and 0.13 (Si), 0.59 and 0.25 (O), and 0.68 and 0.22 (Mg); the exponents of the d shells are 0.6 (Si), 0.5 (O), and 0.5 (Mg). The same basis set has already been successfully utilized in a couple of recent studies of vibrational and spectroscopic properties of forsterite [25,54].

Thresholds controlling the accuracy of Coulomb and exchange series are set to default values. Reciprocal space is sampled using a Monkhorst–Pack mesh with a shrinking factor of 6 for the primitive cell of forsterite, corresponding to 64 independent k-points in the irreducible portion of the Brillouin zone. A pruned grid with 1454 radial and 99 angular points is used to calculate the DFT exchange-correlation contribution through numerical integration of the electron density over the unit cell volume. The self-consistent-field (SCF) convergence on energy was set to a value of 10⁻¹⁰ hartree for all geometry optimizations and phonon frequency calculations. In the calculation of phonon frequencies, the Hessian matrix is computed numerically from finite differences of analytical forces computed at a set of displaced nuclear configurations: a step of 0.003 Å is used.

Harmonic thermodynamic properties, such as the constant-volume specific heat used in Equation (17), are computed by sampling phonon dispersion over 27 k-points (corresponding to a supercell containing 756 atoms). All quasi-harmonic properties are computed by working in terms of the primitive cell of the system (containing 4 formula units or equivalently 28 atoms), which already ensures convergence.

5. Results and Discussion

We apply the quasi-harmonic methodology outlined in Section 2 to compute the thermo-elastic response of the orthorhombic forsterite mineral, $\alpha$-Mg₂SiO₄, that is characterized by nine independent elastic constants: $C_{11}$, $C_{22}$, $C_{33}$, $C_{44}$, $C_{55}$, $C_{66}$, $C_{12}$, $C_{13}$, $C_{23}$. The first step of our algorithm, as sketched in Section 3, consists of a full structural relaxation of the system (i.e., both atomic positions and lattice parameters are optimized by minimizing the static internal energy E of the system). The second step of the algorithm is a lattice symmetry-preserving quasi-harmonic calculation, where harmonic phonon frequencies are evaluated at four volumes (one compressed, one equilibrium and two expanded) according to the fully-automated scheme previously suggested and implemented in the Crystal program by some of the authors [22,23]. From the calculation of phonon harmonic frequencies at different volumes, the Helmholtz free energy of the system as a function of temperature and volume can be evaluated as in Equation (12). By minimizing this function with respect to volume at each temperature, the $V\left(T\right)$ relation is obtained that provides the thermal expansion of the system.

The first two steps of the algorithm, which we have just briefly described, have been performed by using four different functionals of the DFT, which belong to three different classes: the local-density approximation (LDA), two generalized-gradient approximation (GGA) functionals (PBE and PBEsol) and one hybrid functional (PBE0). The results in terms of the thermal expansion are presented in Figure 1, where the volume of forsterite is reported as a function of temperature up to 2000 K. Three independent sets of experimental data are reported from Bouhifd et al. [55] (filled circles), Ye et al. [56] (empty circles) and Hazen [57] (empty squares). The following considerations can be made: (i) as expected, different functionals provide different absolute values for the volume of forsterite, with the PBEsol and PBE0 functionals providing almost identical results in this case and the closest to the experimental values; (ii) all functionals provide a very satisfactory description of the thermal expansion of the system up to 2000 K. In particular, we want to stress that, despite predicting different absolute values for the equilibrium volume, different functionals provide an extremely consistent description of the thermal expansion of the system. This aspect has very interesting computational implications and is not specific of forsterite but is rather a common feature of all quasi-harmonic properties [22,23,24,45]. The computational implication is the following: while a more accurate (and costly) functional (like a hybrid one) is often needed to get a better description of the static absolute value of the equilibrium volume, one can safely work in terms of a less costly functional (like a GGA one) to get the thermal dependence of the volume.

From the second step of our algorithm, we can also compute the thermal dependence of the bulk modulus of the system by evaluating the second derivatives with respect to volume of the Helmholtz free energy, as in Equation (14). The bulk modulus evaluated in this way is the isothermal one, while the adiabatic one can be obtained from the correction as in Equation (17). The isothermal and adiabatic bulk modulus of forsterite are reported as a function of temperature up to 2000 K in the left and right panels of Figure 2, respectively. Experimental data are from Anderson et al. [58] (filled and empty black circles for isothermal and adiabatic moduli). The same considerations made above on the volume hold here on the bulk modulus: (i) different functionals provide different absolute values, with PBEsol and PBE0 resulting in very similar descriptions, almost matching the experimental values; (ii) the computed quasi-harmonic thermal dependence of the bulk modulus is basically insensitive to the particular functional of the DFT used. Let us stress that the adiabatic correction to the isothermal bulk modulus in Equation (17) requires a good description of both quasi-harmonic (such as $V\left(T\right)$, ${\alpha }_V\left(T\right)$, and $K^{\mathrm{T}}\left(T\right)$) and harmonic (such as $C_V\left(T\right)$) thermodynamic properties. We have computed quasi-harmonic properties by working in terms of the primitive cell (28 atoms/cell), which already ensures convergence, while we have used a supercell containing 756 atoms to converge harmonic thermodynamic properties, which show a much slower convergence with sampling of phonon dispersion. The results reported in the right panel of Figure 2 show the good description of all quasi-harmonic and harmonic thermodynamic properties.

We compute the thermal dependence of the elastic constants of forsterite by using the steps 3 to 5 of our quasi-harmonic algorithm, as described in Section 3. Eight different temperatures are considered from 50 K to 1450 K. In order to compute the nine independent elastic constants of forsterite, nine different shapes of deformations are considered. According to the formalism introduced in Section 2.1.2, from the deformation $\underline{\eta }={\underline{\eta }}_1$ we get the $C_{11}$ constant, from $\underline{\eta }={\underline{\eta }}_2$ we get $C_{22}$, from $\underline{\eta }={\underline{\eta }}_3$ we get $C_{33}$, from $\underline{\eta }={\underline{\eta }}_4$ we get $C_{44}$, from $\underline{\eta }={\underline{\eta }}_5$ we get $C_{55}$, from $\underline{\eta }={\underline{\eta }}_6$ we get $C_{66}$, from $\underline{\eta }={\underline{\eta }}_1+{\underline{\eta }}_2$ we get the linear combination $2C_{12}+C_{11}+C_{22}$, from $\underline{\eta }={\underline{\eta }}_1+{\underline{\eta }}_3$ we get the linear combination $2C_{13}+C_{11}+C_{33}$, and from $\underline{\eta }={\underline{\eta }}_2+{\underline{\eta }}_3$ we get the linear combination $2C_{23}+C_{22}+C_{33}$. At each temperature T, and for each of the nine deformation shapes above, the second free energy derivatives in Equation (16) with respect to the deformation amplitude $\delta$ have to be evaluated numerically: five strained configurations corresponding to five different strain amplitudes are explored (including the equilibrium configuration), where harmonic phonon frequencies are computed. Therefore, a total of $9\times 8\times 5=360$ harmonic frequency calculations have been performed to characterize the thermo-elastic response of the orthorhombic crystal of forsterite. It is clear that this methodology is rather computationally expensive. Later on, we are going to discuss how the computational cost of the evaluation of thermo-elasticity can be drastically reduced by adopting a quasi-static approximation rather than a fully quasi-harmonic one.

The computed isothermal elastic constants are then transformed into adiabatic ones by means of Equation (18) in order to be compared to experimental values from Isaak et al. [32]. Figure S1 of the Supplementary Materials shows a comparison between isothermal and adiabatic computed values. Given that the fully quasi-harmonic description of the thermo-elasticity of forsterite is a rather computationally expensive task, and given that different functionals provide consistent trends for the thermal dependence of structural and mechanical properties of solids, we have decided to use the PBEsol functional for this task. This is because, from the analysis of the results shown in Figure 2, it provides a good description of the absolute value of the bulk modulus of the system, and, being a GGA functional, it has a lower cost than the hybrid PBE0 one.

The computed quasi-harmonic adiabatic elastic stiffness constants of forsterite are reported in Figure 3 (with lines) as a function of temperature, up to 1500 K, where they are compared with experimental data (circles). Very regular trends are obtained for all the nine independent elastic constants as a function of temperature. Given that each elastic constant required eight separate calculations at different volumes (corresponding to eight different temperatures), this is not trivial and shows the high numerical robustness of the implemented algorithm. The absolute values of all constants at low temperature is in good agreement with the experiments: for the $C_{11}$ at about 340 GPa, for $C_{33}$ at about 240 GPa, for $C_{22}$ at about 210 GPa, and for the bunch of other constants at about 70–80 GPa (see Figure S2 of the Supplementary Materials for a zoomed-in plot of these latter constants). This confirms the good description provided by the PBEsol functional already inferred from the analysis of both volume and bulk modulus in Figure 1 and Figure 2. The most important aspect of our thermo-elastic description is clearly the thermal dependence of the computed constants, which turns out to be in excellent agreement with the experimental data. Figure 3 confirms that a full quasi-harmonic approach provides an excellent description of the thermo-elasticity of crystalline materials.

The characterization of the thermal dependence of single-crystal elastic constants is clearly crucial if one wants to discuss the evolution with temperature of the mechanical anisotropy of a material. This is relevant to minerals of geophysical interest as well as to other classes of materials such as metal-organic frameworks, whose mechanical instability at room temperature often represents a limiting factor to their effective use [59]. Once all elastic constants are known, the spatial dependence of the Young modulus of a material can be computed [60]. In Figure 4, we report the spatial dependence of the Young modulus of forsterite as a function of temperature. The experimental evolution of the mechanical anisotropy is reported as well as the computed one. Forsterite has a very anisotropic mechanical response. In particular, the system has a larger rigidity along the x-axis than along y and z, which results in strong mechanical anisotropies in the $xy$ and $xz$ planes, while the mechanical anisotropy is way less pronounced in the $yz$ plane. From the analysis of the experimental data in Figure 4, we see that, as temperature increases, the mechanical response of the system changes in two respects: (i) the Young modulus shrinks, which corresponds to an overall softening of the system with temperature; and (ii) the mechanical anisotropy becomes more pronounced (clearly seen in the 2D projections of the Young modulus in the right panels). Both these effects are nicely reproduced by the full quasi-harmonic description of the thermo-elasticity, apart from a slight underestimation of the mechanical anisotropy in the $yz$ plane.

As we have discussed above, the full quasi-harmonic description of the thermo-elastic response of forsterite is rather computationally expensive as it required about 360 harmonic frequency calculations at different strained lattice configurations. In Section 2.2.4, we have introduced a simplified approach, called quasi-static approximation (QSA), to the calculation of thermo-elastic constants where the thermal dependence of the mechanical response is only due to the thermal expansion of the system. Within the QSA, isothermal elastic constants are computed according to Equation (20) by starting from the equilibrium structure at temperature T and by computing the second energy derivatives in terms of internal static energy E rather than of Helmholtz free energy F. This results in a dramatic reduction of the complexity of the algorithm as very few (just four, in our algorithm) phonon frequency calculations are needed to determine the $V\left(T\right)$ relation and none for the second energy derivatives. At the same time, the QSA is a further approximation with respect to the full QHA one that typically results in the underestimation of thermal effects. We have computed the nine thermo-elastic constants of forsterite as a function of temperature with the QSA (they are reported in Figure S3 of the Supplementary Materials). The bottom panels of Figure 4 show the thermal dependence of the Young modulus of forsterite as obtained from the QSA computed constants, as compared with experimental and fully quasi-harmonic values. It can be seen that the two main effects of temperature on the mechanical response of forsterite are both correctly described by the QSA from a qualitative point of view: the shrinking of the Young modulus and the increased anisotropy. From a quantitative point of view, the shrinking of the Young modulus with temperature (i.e., the thermal softening of the mechanical response) is underestimated, as can be clearly seen from the small separation of the different lines in the 2D projections, compared to the experimental ones.

6. Conclusions

Two computational approaches for the first-principles calculation of thermo-elastic stiffness constants of crystalline materials have been formally illustrated and implemented into the public Crystal program for quantum-mechanical simulations of the condensed matter. Both schemes are based on the quasi-harmonic description of the thermal expansion of the system. The first approach, referred to as the quasi-static approximation, assumes that the whole thermal dependence of the mechanical properties of a solid is due to its thermal expansion and therefore involves second energy derivatives with respect to the strain in terms of the purely static internal energy E. The second scheme is a full quasi-harmonic one where thermal expansion is taken into account and second energy derivatives with respect to strain are computed in terms of the Helmholtz free energy F. The latter approach is much more computationally expensive than the former but also ensures a better description of thermal effects on the mechanical response of the system. Both approaches have been applied to the description of the thermo-elasticity of the forsterite mineral.