Enhancing efficiency and scope of first-principles quasiharmonic approximation methods through the calculation of third-order elastic constants

A novel computational strategy is presented to calculate from first principles the coefficient of thermal expansion and the elastic constants of a material over meaningful intervals of temperature and pressure. This strategy combines a novel implementation of the quasiharmonic approximation to calculate the isothermal-isochoric linear and nonlinear elastic constants of a material, with elementary equations of nonlinear continuum mechanics. Our implementation of the quasiharmonic approximation relies on finite deformations, the use of nonprimitive supercells to describe a material, a recently proposed technique to calculate generalized mode Grüneisen parameters, and the numerical differentiation of the stress tensor to calculate both second- and third-order elastic constants. The combination of this method with nonlinear continuum mechanics is shown to yield accurate predictions of lattice parameters and linear elastic constants of a material over finite intervals of temperature and pressure, at the cost of calculating isothermal second- and third-order elastic constants for a single reference state. Here, the validity and limits of our novel methods are assessed by carrying out calculations of MgO based on classical interatomic potentials. To demonstrate potential, our methods are then used in conjunction with a density functional theory approach to calculate thermal expansion and elastic properties of silicon, lithium hydrate, and graphite.

Figure 1

Schematic illustration of the key ingredients of conventional (left) and present (right) implementations of QHA to calculate thermoelastic parameters of a material. Both implementations rely on generating a set of deformed configurations (light red) of a reference state (light gray); $\mu$ is the Lagrangian strain, and for convenience, tensor components are indexed using the Voigt notation. In the conventional approach, SOECs are obtained via second-order differentiation of the Helmholtz free energy ($A$), whereas the present approach relies on the calculation of the generalized mode Grüneisen parameters (${\gamma }_i^{\left(\alpha \right)}$) to obtain SOECs and TOECs via first- and second-order differentiation of the second Piola-Kirchhoff stress tensor ($P$), respectively.

Figure 2

Lattice parameter (top) and bulk modulus (bottom) of MgO vs $T$ at zero pressure calculated by using a conventional QHA approach (red circles) and the methods presented in this work (black solid line), i.e., explicit QHA calculation of thermal stress, isochoric-isothermal SOECs and TOECs for a reference state yielding a zero static pressure, followed by the extrapolation procedure described in Sec. 4. MgO is described through the use cubic supercells contains 512 (red circles) and 64 (black solid line) atoms. Inset, rock-salt structure of MgO.

Figure 3

Helmholtz free energy (blue discs) and pressure (red discs) vs the lattice parameter of MgO at 1000 K. Symbols show results calculated as described in the text by using an energy scheme based on classical interatomic potentials [50], whereas solid black curves are the result of a polynomial interpolation.

Figure 4

Isothermal SOECs of MgO at zero pressure vs $T$ calculated by using the present methods. Red circles show results obtained by calculating SOECs explicitly through the use of our implementation of QHA; at each $T$, we used our QHA method to calculate the isothermal-isochoric SOECs of MgO with lattice parameter $a\left(T\right)$ (Fig. 2). The black solid lines show values of SOECs obtained by using the extrapolation procedure relying on nonlinear continuum mechanics discussed in Sec. 4, using thermal stress, SOECs and TOECs calculated explicitly via QHA for a reference state yielding a zero static pressure.

Figure 5

Percent variation of the lattice parameter (top) and bulk modulus (bottom) of MgO at 300 K vs pressure calculated by using a conventional QHA approach (red circles) and the methods presented in this work (solid lines), i.e., explicit calculation by using our QHA method of thermal stress, and isochoric-isothermal linear and nonlinear elastic constants for a reference state yielding a zero static pressure, followed by the extrapolation procedures relying on nonlinear continuum mechanics described in Sec. 4. The black solid line shows results obtained by combining Eqs. (8) and (11), whereas the thin blue solid line shows results obtained by including in Eq. (8) also the leading third-order terms in strain [Eq. (29)], whose coefficients are the isothermal-isochoric fourth-order elastic constants $C_{\alpha \beta \beta \beta }^{\left(4\right)}$, with $\alpha ,\beta =1,...,6$.

Figure 6

Percent variation of the lattice parameter (top), linear thermal expansion coefficient (middle), and bulk modulus (bottom) of Si at zero pressure and increasing $T$. Solid lines show results obtained from DFT by employing supercells containing 64 (gray) and 216 (black) atoms, and by carrying out QHA calculations of thermal stress, and isothermal-isochoric SOECs and TOECs for a reference state yielding a static pressure of $-$0.34 GPa, followed by the extrapolation technique discussed in Sec. 4. Red circles show results obtained by using a conventional QHA approach [42]. Experimental data for the linear thermal expansion coefficient are from Refs. [74] (blue discs), and [75] (light blue discs), and [27] (blue circles). Inset, the diamond structure of Si.

Figure 7

Percent variation of the lattice constant of LiH vs $T$ at zero pressure, and (inset, top left) vs pressure at 300 K, calculated by using both a conventional QHA approach (solid red circles) and our methods. Gray and black solid lines show results obtained by using supercells containing 64 and 216 atoms, respectively. Blue discs show selected experimental data extracted from Refs. [65, 76]. Inset, bottom right, the rock-salt structure of LiH.

Figure 8

(Top) Percent variation of the in-plane ($a$) and out-of-plane ($c$) lattice parameters of graphite vs $T$. Bottom panel, independent SOECs of graphite vs $T$; values are shifted as indicated in figure. These results are obtained from DFT by using our QHA method to calculate thermal stress, SOECs and TOECs for a reference state, followed by the use of the extrapolation techniques relying on nonlinear continuum mechanics discussed in Sec. 4. Thick gray curves show results obtained by considering a reference state yielding a zero static pressure at 0 K, whereas the solid green lines show results obtained by considering a reference state yielding a zero total pressure at 300 K. DFT calculations are carried out by using a hexagonal supercell containing 100 atoms, consisting of a $5\times 5\times 1$ array of primitive unit cells. Inset, top left corner, an image showing the layered structure of graphite.

Figure 9

Linear thermal expansion coefficients parallel (top) and perpendicular (bottom) to the $c$ axis of graphite. Black solid lines show results obtained by combining QHA calculations and nonlinear continuum mechanics, whereas blue circles show experimental data [77].