Straightforward computation of high-pressure elastic constants using Hooke’s law: A prototype of metal Ru

According to Hooke’s law, within the linear elastic regime of a crystal the relation between the stresses $\sigma_{i}$ and its induced strains $\varepsilon_{j}$ is,

$$\sigma_{i}=\sum_{j=1}^{6}C_{ij}\varepsilon_{j}$$

(1)

where the $C_{ij}$ are the elastic constant of the crystal.

We can obtain all the elastic constants of a material from Eq.(1), by applying the strain $\varepsilon_{j}$ and calculating the corresponding stresses, The deformation matrix applied is

$$\boldsymbol{D}=\boldsymbol{I}+\boldsymbol{\varepsilon},$$

(2)

where $\boldsymbol{I}$ is the $3\times 3$ unit matrix, and $\boldsymbol{\varepsilon}$ is the strain-matrix in Voigt notation. In the 3D case, the strain matrix is

$$\boldsymbol{\varepsilon}=\left[\begin{array}[]{l}\varepsilon_{1}\quad\frac{\varepsilon_{6}}{2}\quad\frac{\varepsilon_{5}}{2}\\ \frac{\varepsilon_{6}}{2}\quad\varepsilon_{2}\quad\frac{\varepsilon_{4}}{2}\\ \frac{\varepsilon_{5}}{2}\quad\frac{\varepsilon_{4}}{2}\quad\varepsilon_{3}\end{array}\right].$$

(3)

The deformed crystal lattice vector is

$$\boldsymbol{A}^{\prime}=\boldsymbol{A}\cdot\boldsymbol{D}$$

(4)

where $\boldsymbol{A}$ is the crystal lattice vector without deformation.

From the elastic constants, the elastic moduli can be readily derived. The Voigt and Reuss elastic moduli for different crystal systems are calculated according to Ref. Wu et al. (2007). According to the Voigt–Reuss–Hill approximations Hill , the arithmetic average of Voigt and Reuss bounds is

$$B=B_{\operatorname{VRH}}=\frac{B_{V}+B_{R}}{2},$$

(5)

and

$$G=G_{\operatorname{VRH}}=\frac{G_{V}+G_{R}}{2}$$

(6)

respectively.

Young’s modulus $E$ is calculated by

$$E=\frac{9BG}{3B+G},$$

(7)

and Poisson’s ratio is

$$\nu=\frac{3B-2G}{2(3B+G)}.$$

(8)

The elastic anisotropy is a direct reflect of spatial anisotropy of chemical bonding. Chung and Buessem defined an elastic anisotropy index Chung and Buessem ,

$$A^{C}=\frac{G_{V}-G_{R}}{G_{V}+G_{R}}.$$

(9)

Ranganathan and Ostoja-Starzewski proposed a more general anisotropy index called the universal elastic anisotropy index Ranganathan and Ostoja-Starzewski (2008),

$$A^{U}=5\frac{G_{V}}{G_{R}}-\frac{B_{V}}{C_{R}}-6\geqslant 0.$$

(10)

$A^{U}=0$ corresponds to the locally isotropic crystals  Ranganathan and Ostoja-Starzewski (2008).

The crystal Ru has a hexagonal closely-packed (hcp) structure up to 600 GPa. Before the calculations of elastic constants, the crystal structure of Ru was first optimized within the framework of DFT at each pressure. Then the stress components were calculated accurately after atomic positions were relaxed under specific deformations applied according to our recently proposed “the optimized high efficiency strain-matrix sets (OHESS)”Liu (2020a). The relaxations stopped after the forces acted on each atom were less than 0.02 eV/Å. The projector augmented wave (PAW) method Blöchl (1994) as implemented in the Vienna ab initio simulation package (VASP) Kresse and Joubert (1999); Kresse and Furthmüller (1996) was adopted in all the structural optimizations and stress computations in this work. The Perdew, Becke and Ernzerhof (PBE) Perdew et al. (1997) generalized gradient approximation (GGA) was used for the exchange-correlation functional. The energy cutoff values were set to ensure energy to be converged to $10^{-6}$ eV.

The calculated lattice parameters accords very well with both experimental and others’ theoretical results, as shown in Table 2. Also listed in Table 2 are our calculated elastic constants of Ru at 0 GPa. We see that our elastic constants are in reasonably good with experimental data and others’ calculated results. Especially, we note that our stress-strain results are in very good agreement with the results from the energy-strain method by Lugovskoy et al.Lugovskoy et al. (2014). It is not surprising because they are the zero-pressure elastic constants, not considering pressure effects in both the stress-strain and energy-strain methods.

In order to compare the stress-strain method with the energy-strain method (Ref. Lugovskoy et al. (2014)), we accurately tuned our pressures to correspond to those calculated in Ref. Lugovskoy et al. (2014). The maximum difference in pressure is not larger than 0.3 GPa, as shown in Table 3. The corresponding mean atomic volumes are also in good agreement with each other. The maximum of mean atomic volume differences are not larger than 0.04 ų, below 500 GPa. The high-pressure elastic constants are listed in Table 3, from which we can numerically compare the elastic constants from the two different methods, the stress-strain and energy-strain methods. The pressure derivative of elastic constants $\mathrm{d}C_{ij}/\mathrm{d}P$ are shown in Table 2. It is surprisingly good that our stress-strain values agree with the energy-strain values that were calculated with complex pressure corrections by Lugovskoy et al. Lugovskoy et al. (2014). This indicates that, although the stress-strain method is rather straightforward to calculate the high-pressure elastic constants, it can reach the same accuracy of the energy-strain method with complicated pressure corrections.

To clearly compare the high-pressure elastic constants from the two methods, we presented all the elastic constants versus pressure data in Fig. 1. Except for the slight differences in $C_{11}$ and $C_{12}$ beyond 400 GPa, other elastic constants ($C_{13},C_{33}$, and $C_{44}$) are in excellent agreement with each other for the two methods, in the whole pressure range from 0 up to $\sim$600 GPa. We also compare totally all the high-pressure elastic constants from the two methods in Fig. 2, from which we again find the satisfactory agreements between the stress-strain and energy-strain methods. In order to statistically evaluate the agreement of our results with the previous ones, we calculated the coefficient of determination $R^{2}$ score R (2) of present compared to previous results. The evaluated $R^{2}=0.9982$, implying an excellent agreement quantitatively.

From our high-pressure elastic constants, we calculated the bulk modulus $B$, shear modulus $G$, and Young’s modulus $E$ according to Eqs. (5), (6), and (8), respectively, and show them in Table. 4. These three elastic moduli all increase monotonously with pressure, implying the rigidity of Ru increases with pressure. We also calculated the Poisson’s ratio of Ru and found that it increases from 0.2431 to 0.3367 at pressure ranging from 0 to $\sim$600 GPa.

In order to understand the variation of metallic bonding in Ru with pressure, we analyzed its elastic anisotropy properties. First, according to Eqs. (9) and (10), we calculated the elastic anisotropy indexes which can overall reflect the bonding characteristics of materials. As shown in Table 4, both $A^{U}$ and $A^{C}$ have very small variations at pressures from 0 up to $\sim$600 GPa. This indicates that the bonding properties keeps almost invariant in the hcp-Ru.

Second, we plotted the spatial variations of the Young’s modulus, shear modulus, and Poisson’s ratio at different pressures in Fig. 3. It can be clearly seen that the spatial elastic anisotropy of Ru almost remains unchanged with pressure. This accords well with the conclusion drawn from the elastic anisotropy indexes analysis.

The phase velocity $v$ and polarization of the three waves along a fixed propagation direction defined by the unit vector $n_{i}$ are given by Cristoffel equation,

$$(C_{ijkl}n_{j}n_{k}-\rho v^{2}\delta_{ij})u_{i}=0,$$

(11)

where $C_{ijkl}$ is the fourth-rank tensor description of the elastic constants, $n$ is the propagation direction, and $u$ the polarization vector.

The longitudinal wave velocity is

$$v_{P}=\sqrt{\frac{B+4G/3}{\rho}},$$

(12)

and the body wave velocity is

$$v_{B}=\sqrt{\frac{B}{\rho}}.$$

(13)

The shear wave velocity is calculated by

$$v_{S}=\sqrt{\frac{G}{\rho}}.$$

(14)

From $v_{S}$ and $v_{P}$, we can obtain the mean wave velocity via

$$v_{m}=\left[\frac{1}{3}\left(\frac{2}{v_{S}^{3}}+\frac{1}{v_{P}^{3}}\right)\right]^{-1/3}.$$

(15)

From the mean wave velocity data, we can calculate the Debye temperature from

$$\Theta_{D}=\frac{h}{k}\left[\frac{3n}{4\pi}\left(\frac{N_{A}\rho}{M}\right)\right]^{1/3}v_{m}$$

(16)

where $h$ is Planck’s constant, $k$ the Boltzmann’s constant, $N_{A}$ the Avogadro’s number, $n$ the number of atoms in the unit cell, $M$ the weight of the unit cell, and $\rho$ the density.

We show the calculated sound velocities in Fig. 4, from which we see the monotonous increases in all these sound velocities. This indirectly reflects the strengthening of the metallic bondings in Ru with increasing pressure.

Debye temperature $\Theta_{D}$ corresponds to the temperature of a crystal’s highest normal mode of lattice vibration, and it closely correlates the elastic properties with the thermodynamic properties Luo and Wang (2008). The calculated Debye temperatures are listed in Table 5. The increasing of $\Theta_{D}$ with pressure again reflects the enhancing of the metallic bonding in Ru.