# First-Principles Calculations to Investigate the Third-Order Elastic Constants and Mechanical Properties of Mg, Be, Ti, Zn, Zr, and Cd

Authors: Xiaoqing Yang; Zhenya Meng; Hailin Cao
Journal: Advances in Materials Science and Engineering
 (2021)

---

## Abstract

We present theoretical studies for the third-order elastic constants of Mg, Be, Ti, Zn, Zr, and Cd with a hexagonal-close-packed (HCP) structure. The method of homogeneous deformation combined with first-principles total-energy calculations is employed. The deformation gradientFij is applied to the crystal lattice vectors ri, and the elastic strain energy can be obtained from the first-principles calculation. The second- and third-order elastic constants are extracted by a polynomial fit to the calculated energy-strain results. In order to assure the accuracy of our method, we calculated the complete set of the equilibrium lattice parameters and second-order elastic constants for Mg, Be, Ti, Zn, Zr, and Cd, and our results provide better agreement with the previous calculated and experimental values. Besides, we have calculated the pressure derivatives of SOECs related to third-order elastic constants, and high-pressure effects on elastic anisotropy, ductile-to-brittle criterion, and Vickers hardness are also investigated. The results show that the hardness model Hv=1.877k2G0.585 is more appropriate than Hv=2k2G0.585−3 for HCP metals under high pressure.


---

## Body

## 1. Introduction

In the theory of linear elasticity, infinitesimal deformation strains are assumed, and using finite deformation on materials is very significant for practical application. The second-order elastic constants (SOECs) are sufficient to describe the elastic stress-strain response and wave propagation in solids [1], whereas the third- and higher-order elastic constants (TOECs and HOECs) must be considered to characterize the strain properties of materials [2–5]. In general, the SOECs and TOECs describe the response of materials to the linear and nonlinear elasticity, respectively. For single crystals, TOECs reflect variation in acoustic velocities as a result of elastic strain [6, 7], so not only can TOECs describe mechanical phenomena when large amplitude stress are acted, but also it can describe other anharmonic properties including thermal expansion, temperature dependence of elastic properties, the interactions of phonon and phonon [6, 8]. In the past research, many experiments have been implemented to determine TOECs [3]. However, it is very difficult to obtain a complete set of TOECs from experimental methods for HCP metals, so it is eager to predict the TOECs of HCP metals by using theoretical methods. There are quite a few of the available theoretical methods for determining TOECs, including empirical force-constant models [9–13], molecular-dynamics simulations [14, 15], and first-principles total-energy methods [1, 16]. Nielson and Martin first applied to the methods of first-principles total-energy calculations determining TOECs of materials [17, 18]. Recently, the method from cubic symmetric crystals has been extended to arbitrary symmetric crystals by Zhao et al. [19].In the present study, we describe the method combined homogeneous deformation with first-principles total-energy calculations for determining TOECs of HCP metals (Mg, Be, Ti, Zn, Zr, and Cd). It is well known that the properties in HCP materials have been studied for a long time, especially for Mg. Mg has good ductility, better characteristics suppressing vibration and noise than aluminium, and excellent castability [20], and alloying magnesium can increase the strength-to-weight ratio due to their lightweight and high strength, making them important materials for structural applications [21–24]. Furthermore, zirconium and titanium have also been broadly applied, particularly in structural applications, for example, in the automotive and aerospace industries due to their corrosion resistance, lightweight, high strength, and so on [25]. This paper systematically studies the elastic characterization for HCP metals including Mg, Be, Ti, Zn, Zr, and Cd to better understand the strain behavior by using first-principles calculations. In order to test the feasibility of our model, we have calculated the equilibrium lattice parameters and SOECs, which are compared with previous theoretical and experimental results, and it is shown that our results are in perfect agreement with acceptable data. This paper is organized as follows: Section 2 gives a general overview of the nonlinear elasticity theory. Section 3 describes employed methodology and results discussion for TOECs obtained from first-principles calculations. Section 4 deals with the determination of the pressure dependent elastic constants related to TOECs, high-pressure effects on elastic anisotropy, ductile-to-brittle criterion, and hardness of Mg, Be, Ti, Zn, Zr, and Cd. Finally, we present our conclusions in Section 5.


## 2. Nonlinear Elasticity Theory

In this work, we will recall some basic facts from the nonlinear theory of elasticity [5, 26–30]. Let ai be the initial coordinates of some crystal element. After applying a finite deformation to a crystal, the initial crystal element will move to the position xi. After introducing the deformation gradient Fij,(1)Fij=∂xi∂aj,where i and j (=1, 2, 3) represent the Cartesian coordinates. Then, we may define the Lagrangian strain tensor η:(2)ηij=12∑pFpiFpj−δij,which is a convenient estimation of deformation for an elastic body, and the Lagrangian strains tensor does not contain information referring to rigid rotation of the crystal element based on its symmetry.In the nonlinear elasticity theory, the Lagrangian strain tensor is related to the internal energy, and the free energy is related to the Lagrangian strain tensor. The elastic constants can be obtained by expanding the internal energy as a Taylor series in terms of the Lagrangian strain tensor at constant entropy by Bruuger [29]:(3)UV,η,S=UV,0,S+V2!∑ijklCijklSηijηkl+V3!∑ijklmnCijklmnSηijηklηmn+…,where V is the volume of the system, S is the entropy, and UV,0,S is the corresponding internal energy of the ground state. An expansion UV,η,S in terms of symmetric Lagrangian strains tensor is appropriate because the internal energy is unaltered under rigid rotation [2]. The SOECs, TOECs, and HOECs are defined as the second-, third-, and higher-order derivatives of equation (3) with regard to the Lagrangian strain, respectively. Thus, the isentropic elastic constants can be written as(4)CijklS=1V∂2U∂ηij∂ηkl|η=0SOEC,CijklmnS=1V∂3U∂ηij∂ηkl∂ηmn|η=0TOEC.And the isothermal elastic constants are(5)FV,η,T=UV,0,T+V2!∑ijklCijklTηijηkl+V3!∑ijklmnCijklmnTηijηklηmn+…,which is obtained by expanding the Helmholtz free energy as a Taylor series of η at constant temperature [2, 6]. Since our first-principles calculations are performed at T=0K, F=U−TS=U, then CS=CT. We will not distinguish between isentropic elastic constants CS and isothermal elastic constants CT in this work.Only six of each set of the nine Lagrangian strain tensors are independent due to its symmetry, and it is convenient to introduce the Voigt notation (11⟶ 1, 22 ⟶ 2, 33 ⟶ 3, 23 ⟶ 4, 31 ⟶ 5, 12 ⟶ 6), for the Lagrangian strain tensors, η11⟶η1,η22⟶η2,η33⟶η3,η23⟶η4/2,η31⟶η5/2,η12⟶η6/2 . So, the Lagrangian strain tensors η can be written as follows:(6)η=η1η62η52η62η2η42η52η42η3.Equations (3) and (5) now can be written as(7)ΔUV,ηV=12!∑i,j=1,6Cijηiηj+13!∑i,j,k=1,6Cijkηiηjηk+…,where ΔUV,η/V=UV,η−UV,0/V as a function of applied finite strain is the elastic energy per unit volume. As shown above, the second and third orders of elastic constants can be obtained via using the method of homogeneous deformation by applying various simple deformations to the crystal, and the internal energy for the deformed crystal is calculated. Then, SOECs and TOECs are extracted by fitting equation (8) to the first-principles calculated energy-strain results.Besides, the Lagrangian stresst is defined as the first-order derivative of the internal energy with regard to the Lagrangian strain tensor η:(8)tij=1V∂U∂ηij,which can also be used to evaluate the elastic constants. In the following section, we will obtain the SOECs and TOECs by applying to the system ten simple deformation strains in the first-principles calculations.


## 3. Determination of the Elastic Constants

### 3.1. Computational Methodology

In the work presented here, we have determined the second-order elastic constants and third-order elastic constants for Mg, Be, Ti, Zn, Zr, and Cd combined with the first-principles calculations based on density functional theory (DFT). The deformation gradientFij is applied to the crystal lattice vectors ri to obtain the unit cell for the strained crystal, and the deformed lattice vectors can be determined:(9)ri′=∑jFijaj.To carry out the different deformation modes in our work, the deformation gradient matrixFij needs to be obtained. Fij is related to the strain η and is determined by inverting equation (2):(10)Fij=δij+ηij−12∑kηikηkj+12∑klηikηklηlj+….In general, the deformation gradientFij is not single for a given strain η, while the various possible solutions differ from one another by a rigid rotation. The lack of a one-to-one relationship between the deformation gradient matrix Fij and the strain η is insignificant due to the invariability of the calculated total energy under rigid deformation [19]. For a given specific deformation, to obtain the minimized energy for the strained lattice, relaxation of the crystal internal coordinates for the distorted unit cell was performed.For hexagonal structure, there are five independent SOECsC11,C12,C13,C14,C44 and ten TOECs C111, C222, C333, C112,C113,C123, C133,C144,C155,C344. To determine these elastic constants, we introduced Lagrangian strain tensors η for hexagonal crystals that result in an energy expansion (equation (3)), including a small quantity of second-order elastic constants and third-order elastic constants. So, we can confirm the nonzero components of each Lagrangian strain tensor ηij according to a single parameter ξ. Moreover, the selection of the different deformation modes results in different strains used in our work, which are written as ηα, α  =  A, B, …, J and presented in Table 1. Inserting these strains into equation (7), the elastic energy per unit volume which resulted in this selection of the deformation mode can be expressed as an expansion in the strain parameter ξ:(11)V−1Uξ−U0=12P2ξ2+16P3ξ3+Oξ4.Table 1
The coefficientsP2 and P3 in equation (12) of corresponding Lagrangian strains are shown as combinations of second-order and third-order elastic constants for the hexagonal crystal.
Strain typeP2P3ηA=ξ,0,0,0,0,0C11C111ηB=ξ,ξ,0,0,0,02C11+2C124C111+6C112−2C222ηC=ξ,ξ,ξ,0,0,02C11+2C12+4C13+C334C111+6C112+6C113+6C123+6C133−2C222+C333ηD=0,ξ,0,0,0,0C11C222ηE=0,0,ξ,0,0,0C33C333ηF=0,ξ,,ξ,0,0,0C11+2C13+C333C113+3C133+C222+C333ηG=0,0,ξ,0,0,2ξ2C11−2C12+C336C113−6C123+C333ηH=0,ξ,0,0,2ξ,0C11+4C4412C144+C222ηI=0,0,ξ,0,2ξ,0C33+4C44C333+12C344ηJ=ξ,ξ,0,0,2ξ,02C11+2C12+4C444C111+6C112+12C144+12C155−2C222CoefficientsP2 and P3 represent combinations of second-order elastic constants and third-order elastic constants of the crystal, respectively. Furthermore, coefficients P2 and P3 for the specific strain tensors η are listed in Table 1. We can obtain these coefficients by fitting equation (10) to plots of energy per unit volume versus strain ξ. In every case for η, to obtain accurate TOECs, ξ is varied between −0.08 and 0.08 with Step 0.008. For every deformation mode, the coordinates of atoms optimized and stress tensors and internal energy are calculated based on density functional theory (DFT).In this work, we implement first-principles total-energy calculations on the basis of the density functional theory (DFT), which is embodied in the Vienna ab simulation package (VASP) [31–33]. All calculations on Mg, Be, Ti, Zn, Zr, and Cd are carried out by using the Perdew–Burke–Ernzerhof (PBE) [34, 35] exchange-correlation functional for the generalized-gradient-approximation (GGA). The projector augmented wave (PAW) method has been carried out in the VASP package [36]. The Monkhorst-Pack special k-point scheme [37] represents reciprocal space with 21  ×21×  21 grid meshes is enough for obtaining accurate results. Meanwhile, the plane-wave cutoff was set to 600 eV for all the calculated HCP metals. The convergence of energy and force are set to 1.0×10−6 eV and 1.0  ×10−5 eV/Å, respectively.Taking Mg as an example, Figure1(a) shows the relationship between the internal energy convergence and the k-point grid size. The internal energy is well converged after 15×15×15 mesh size. The test of the cutoff energy affects the internal energy convergence displayed in Figure 1(b), and E\spmathrmcutoff\spmathrmMg  = 350 eV is sufficient for Mg. The equilibrium lattice parameters a and c/a for Mg can be determined by minimizing the stress on the unit cell and the Hellmann–Feynman force on the atoms. Figure 1(c) shows the internal energy calculations for Mg corresponding to the equilibrium lattice parameters a and c/a. It shows that our values of calculation agree well with the previous experimental and calculated values (see Table 2).Figure 1
In the calculated internal energy of Mg as a function of k-points grid size, the cutoff energy, and lattice constants, respectively. (a) The analogous dependence on the density of k-points mesh (energy cutoff of 600 eV is used for all points). (b) The dependence of internal energy on the cutoff energy (Monk-Pack sampling21×21×21 is applied to all points). The energy converges well when the k-point mesh size and the cutoff energy are beyond 15×15×15 and 300 eV, respectively. In our calculation, we choose the Monkhorst-Pack sampling 21×21×21 and the energy cutoff of 600 eV for Mg. (c) The equilibrium lattice parameters for Mg a and c/a are 3.193 Åand 1.627 Å, which is determined from the corresponding minimum value of the internal energy\enleadertwodots.
(a)(b)(c)Table 2
Calculated values for the lattice parametersa, in Å) and elastic constants (Cij, in GPa), and the available experimental and theoretical data are given.
MgBeTiZnZrCda3.1932.2622.9252.6483.2343.0323.21a,3.19b2.28a,2.2856c2.90a,2.9506c2.64a,2.65b3.20a,3.2312c3.03b,2.98dc/a1.6271.5781.5801.9101.6001.9101.62a,1.61b1.56a,1.5677c1.55a,1.5857c1.81a,1.93b1.56a,1.5931c1.82b,1.89dC1163.96390.06175.26167.71142.4692.1860b,59.3e292e,292.3f160e,176.1g163b,165e144e,155.4g70b,116eC3367.27367.95191.2561.79164.9840.0475b,61.5e349e,246.0f181e,190.5g67b,61.8e166e,172.5g57b,50.9eC4417.02160.4636.3526.3022.2517.3822b,16.4e163e,162.5f46.5e,50.8g46b,39.6e33.4e,36.3g24b,19.6eC1228.3023.9396.1437.8671.6644.9634b,25.7e24e,26.7f90e,86.9g53b,31.1e74e,67.2g60b,43eC1320.5113.82381.2148.4663.1931.8921b,21.4e6e,14f66e,64.6g36b,50e67e,64.6g40b,41eaReference [50] calculated from first-principles calculations. bReference [51] calculated from first-principles calculations. cReference [52] taken from Barrett and Massalski. dReference [53] obtianed from experiments. eReference [54] obtained from the data compiled by Hearmon. fReference [55] taken from Simmon and Wang. gReference [56] obtianed from an ultrasonic wave interference technique.


### 3.2. Results and Discussion

We list our calculated values of second-order elastic constants (SOECs)Cij and third-order elastic constants (TOECs) Cijk for Mg, Be, Ti, Zn, Zr, and Cd in Tables 2 and 3, respectively. We also provide values of calculation for SOECs and compare them with previous results to ensure the correctness of our calculation. The Cij values are slightly different for different fitted curves. Taking Mg as an example, C11 from coefficients in fAξ and fDξ is 64.02 GPa and 63.90 GPa, respectively; the C44 values are 17.15 GPa, 17.03 Gpa, and 16.88 GPa from coefficients in fHξ, fIξ, and fJξ, respectively. In such cases, the average of all calculated values is presented in Table 2. It is shown that our results of SOECs agree well with the available data from other theoretical and experimental values. The energies-strain E−ξ curves for Mg, containing the results of the first-principles calculations and the fitted polynomials, are plotted in Figure 2. The hollow circle denotes results of DFT results, solid lines represent the curves obtained from nonlinear elasticity theory, and dashed lines indicate the curves obtained from linear elasticity, respectively. For the curves obtained from nonlinear elasticity theory, E-ξ with positive strains are always smaller than those with negative strains, so the majority of values of TOECs Cijk are negative (see Table 3), while it is the opposite for SOECs (see Table 2).Table 3
Predictions for the third-order elastic constants (TOECs) of Mg, Be, Ti, Zn, Zr, and Cd. The unit of all data is GPa.
MgBeTiZnZrCdC111−879.14−3865.41−1591.00−3007.54−856.89−1978.36C222−801.04−2931.07−1457.22−2988.01−835.38−2340.52C333−670.98−2932.89−1635.42−528.17−1303.69−414.87C112−801.58250.623309.761−289.74−224.99−368.93C113−74.08−21.59−750.07−125.04−220.521−134.95C123−53.82−176.83−1112.21−81.77−556.40−42.00C133−193.3181.44419.96−265.92−94.83−185.52C144−28.16−256.838.22−33.8828.6114.33C155−36.03111.16−345.29−101.2822.94−131.33C344−160.09−930.03−109.11−139.51−77.4656.53Figure 2
The strain-energy relations for Mg. The hollow circle denotes the results of DFT results; solid lines represent the curves obtained from nonlinear elasticity theory; dashed lines indicate the curves obtained from linear elasticity. (a–j) Lagrangian strainsηA,ηB,ηC,ηD,ηE,ηF,ηG,ηH,ηI, and ηJ respectively.
(a)(b)(c)(d)(e)(f)(g)(h)(i)(j)The determination ofCijk is sensitive to errors in the calculation methods and the parameters and needs superconvergence of parameters controlling the precision of computations. In all of our calculations, the energy cutoff of 600 eV and the Monkhorst-Pack sampling 21 ×21× 21 are very reliable for Mg, Be, Ti, Zn, Zr, and Cd on the basis of our tests. The projector augmented wave (PAW) method chosen to solve Kohn–Sham equations does not affect the computation markedly [38] because the calculations of the static and dynamical properties for a wide range of solids within PAW have been performed properly [39]. Besides, we carry out GGA-PBE exchange-correlation functional and have proven its effectiveness for calculating the TOECs. As shown in Table 2, our results show that the equilibrium lattice parameters and SOECs perfectly agree with the theoretical values.Next, our concern is to investigate for which range of strainsξ the third-order effects dominant the properties of solids. The maximum strain ξmax is an important parameter in these calculations as well. In the present study, for the second-order term, the fitted coefficient P2 was almost independent of the range of fitting, while the coefficient P3 was more sensitive to ξmax. To illustrate this feature, we also show the curves of the nonlinear elasticity comparison with the linear elasticity for Mg in Figure 2. The results show that linear elasticity is not appropriate when the maximum strain ξmax larger than 30% and the third-order effects must be considered. Therefore, the Lagrangian strain tensor η can be expressed as η=χ−1/2χ2, where χ is the linear strain tensor.


## 3.1. Computational Methodology

In the work presented here, we have determined the second-order elastic constants and third-order elastic constants for Mg, Be, Ti, Zn, Zr, and Cd combined with the first-principles calculations based on density functional theory (DFT). The deformation gradientFij is applied to the crystal lattice vectors ri to obtain the unit cell for the strained crystal, and the deformed lattice vectors can be determined:(9)ri′=∑jFijaj.To carry out the different deformation modes in our work, the deformation gradient matrixFij needs to be obtained. Fij is related to the strain η and is determined by inverting equation (2):(10)Fij=δij+ηij−12∑kηikηkj+12∑klηikηklηlj+….In general, the deformation gradientFij is not single for a given strain η, while the various possible solutions differ from one another by a rigid rotation. The lack of a one-to-one relationship between the deformation gradient matrix Fij and the strain η is insignificant due to the invariability of the calculated total energy under rigid deformation [19]. For a given specific deformation, to obtain the minimized energy for the strained lattice, relaxation of the crystal internal coordinates for the distorted unit cell was performed.For hexagonal structure, there are five independent SOECsC11,C12,C13,C14,C44 and ten TOECs C111, C222, C333, C112,C113,C123, C133,C144,C155,C344. To determine these elastic constants, we introduced Lagrangian strain tensors η for hexagonal crystals that result in an energy expansion (equation (3)), including a small quantity of second-order elastic constants and third-order elastic constants. So, we can confirm the nonzero components of each Lagrangian strain tensor ηij according to a single parameter ξ. Moreover, the selection of the different deformation modes results in different strains used in our work, which are written as ηα, α  =  A, B, …, J and presented in Table 1. Inserting these strains into equation (7), the elastic energy per unit volume which resulted in this selection of the deformation mode can be expressed as an expansion in the strain parameter ξ:(11)V−1Uξ−U0=12P2ξ2+16P3ξ3+Oξ4.Table 1
The coefficientsP2 and P3 in equation (12) of corresponding Lagrangian strains are shown as combinations of second-order and third-order elastic constants for the hexagonal crystal.
Strain typeP2P3ηA=ξ,0,0,0,0,0C11C111ηB=ξ,ξ,0,0,0,02C11+2C124C111+6C112−2C222ηC=ξ,ξ,ξ,0,0,02C11+2C12+4C13+C334C111+6C112+6C113+6C123+6C133−2C222+C333ηD=0,ξ,0,0,0,0C11C222ηE=0,0,ξ,0,0,0C33C333ηF=0,ξ,,ξ,0,0,0C11+2C13+C333C113+3C133+C222+C333ηG=0,0,ξ,0,0,2ξ2C11−2C12+C336C113−6C123+C333ηH=0,ξ,0,0,2ξ,0C11+4C4412C144+C222ηI=0,0,ξ,0,2ξ,0C33+4C44C333+12C344ηJ=ξ,ξ,0,0,2ξ,02C11+2C12+4C444C111+6C112+12C144+12C155−2C222CoefficientsP2 and P3 represent combinations of second-order elastic constants and third-order elastic constants of the crystal, respectively. Furthermore, coefficients P2 and P3 for the specific strain tensors η are listed in Table 1. We can obtain these coefficients by fitting equation (10) to plots of energy per unit volume versus strain ξ. In every case for η, to obtain accurate TOECs, ξ is varied between −0.08 and 0.08 with Step 0.008. For every deformation mode, the coordinates of atoms optimized and stress tensors and internal energy are calculated based on density functional theory (DFT).In this work, we implement first-principles total-energy calculations on the basis of the density functional theory (DFT), which is embodied in the Vienna ab simulation package (VASP) [31–33]. All calculations on Mg, Be, Ti, Zn, Zr, and Cd are carried out by using the Perdew–Burke–Ernzerhof (PBE) [34, 35] exchange-correlation functional for the generalized-gradient-approximation (GGA). The projector augmented wave (PAW) method has been carried out in the VASP package [36]. The Monkhorst-Pack special k-point scheme [37] represents reciprocal space with 21  ×21×  21 grid meshes is enough for obtaining accurate results. Meanwhile, the plane-wave cutoff was set to 600 eV for all the calculated HCP metals. The convergence of energy and force are set to 1.0×10−6 eV and 1.0  ×10−5 eV/Å, respectively.Taking Mg as an example, Figure1(a) shows the relationship between the internal energy convergence and the k-point grid size. The internal energy is well converged after 15×15×15 mesh size. The test of the cutoff energy affects the internal energy convergence displayed in Figure 1(b), and E\spmathrmcutoff\spmathrmMg  = 350 eV is sufficient for Mg. The equilibrium lattice parameters a and c/a for Mg can be determined by minimizing the stress on the unit cell and the Hellmann–Feynman force on the atoms. Figure 1(c) shows the internal energy calculations for Mg corresponding to the equilibrium lattice parameters a and c/a. It shows that our values of calculation agree well with the previous experimental and calculated values (see Table 2).Figure 1
In the calculated internal energy of Mg as a function of k-points grid size, the cutoff energy, and lattice constants, respectively. (a) The analogous dependence on the density of k-points mesh (energy cutoff of 600 eV is used for all points). (b) The dependence of internal energy on the cutoff energy (Monk-Pack sampling21×21×21 is applied to all points). The energy converges well when the k-point mesh size and the cutoff energy are beyond 15×15×15 and 300 eV, respectively. In our calculation, we choose the Monkhorst-Pack sampling 21×21×21 and the energy cutoff of 600 eV for Mg. (c) The equilibrium lattice parameters for Mg a and c/a are 3.193 Åand 1.627 Å, which is determined from the corresponding minimum value of the internal energy\enleadertwodots.
(a)(b)(c)Table 2
Calculated values for the lattice parametersa, in Å) and elastic constants (Cij, in GPa), and the available experimental and theoretical data are given.
MgBeTiZnZrCda3.1932.2622.9252.6483.2343.0323.21a,3.19b2.28a,2.2856c2.90a,2.9506c2.64a,2.65b3.20a,3.2312c3.03b,2.98dc/a1.6271.5781.5801.9101.6001.9101.62a,1.61b1.56a,1.5677c1.55a,1.5857c1.81a,1.93b1.56a,1.5931c1.82b,1.89dC1163.96390.06175.26167.71142.4692.1860b,59.3e292e,292.3f160e,176.1g163b,165e144e,155.4g70b,116eC3367.27367.95191.2561.79164.9840.0475b,61.5e349e,246.0f181e,190.5g67b,61.8e166e,172.5g57b,50.9eC4417.02160.4636.3526.3022.2517.3822b,16.4e163e,162.5f46.5e,50.8g46b,39.6e33.4e,36.3g24b,19.6eC1228.3023.9396.1437.8671.6644.9634b,25.7e24e,26.7f90e,86.9g53b,31.1e74e,67.2g60b,43eC1320.5113.82381.2148.4663.1931.8921b,21.4e6e,14f66e,64.6g36b,50e67e,64.6g40b,41eaReference [50] calculated from first-principles calculations. bReference [51] calculated from first-principles calculations. cReference [52] taken from Barrett and Massalski. dReference [53] obtianed from experiments. eReference [54] obtained from the data compiled by Hearmon. fReference [55] taken from Simmon and Wang. gReference [56] obtianed from an ultrasonic wave interference technique.


## 3.2. Results and Discussion

We list our calculated values of second-order elastic constants (SOECs)Cij and third-order elastic constants (TOECs) Cijk for Mg, Be, Ti, Zn, Zr, and Cd in Tables 2 and 3, respectively. We also provide values of calculation for SOECs and compare them with previous results to ensure the correctness of our calculation. The Cij values are slightly different for different fitted curves. Taking Mg as an example, C11 from coefficients in fAξ and fDξ is 64.02 GPa and 63.90 GPa, respectively; the C44 values are 17.15 GPa, 17.03 Gpa, and 16.88 GPa from coefficients in fHξ, fIξ, and fJξ, respectively. In such cases, the average of all calculated values is presented in Table 2. It is shown that our results of SOECs agree well with the available data from other theoretical and experimental values. The energies-strain E−ξ curves for Mg, containing the results of the first-principles calculations and the fitted polynomials, are plotted in Figure 2. The hollow circle denotes results of DFT results, solid lines represent the curves obtained from nonlinear elasticity theory, and dashed lines indicate the curves obtained from linear elasticity, respectively. For the curves obtained from nonlinear elasticity theory, E-ξ with positive strains are always smaller than those with negative strains, so the majority of values of TOECs Cijk are negative (see Table 3), while it is the opposite for SOECs (see Table 2).Table 3
Predictions for the third-order elastic constants (TOECs) of Mg, Be, Ti, Zn, Zr, and Cd. The unit of all data is GPa.
MgBeTiZnZrCdC111−879.14−3865.41−1591.00−3007.54−856.89−1978.36C222−801.04−2931.07−1457.22−2988.01−835.38−2340.52C333−670.98−2932.89−1635.42−528.17−1303.69−414.87C112−801.58250.623309.761−289.74−224.99−368.93C113−74.08−21.59−750.07−125.04−220.521−134.95C123−53.82−176.83−1112.21−81.77−556.40−42.00C133−193.3181.44419.96−265.92−94.83−185.52C144−28.16−256.838.22−33.8828.6114.33C155−36.03111.16−345.29−101.2822.94−131.33C344−160.09−930.03−109.11−139.51−77.4656.53Figure 2
The strain-energy relations for Mg. The hollow circle denotes the results of DFT results; solid lines represent the curves obtained from nonlinear elasticity theory; dashed lines indicate the curves obtained from linear elasticity. (a–j) Lagrangian strainsηA,ηB,ηC,ηD,ηE,ηF,ηG,ηH,ηI, and ηJ respectively.
(a)(b)(c)(d)(e)(f)(g)(h)(i)(j)The determination ofCijk is sensitive to errors in the calculation methods and the parameters and needs superconvergence of parameters controlling the precision of computations. In all of our calculations, the energy cutoff of 600 eV and the Monkhorst-Pack sampling 21 ×21× 21 are very reliable for Mg, Be, Ti, Zn, Zr, and Cd on the basis of our tests. The projector augmented wave (PAW) method chosen to solve Kohn–Sham equations does not affect the computation markedly [38] because the calculations of the static and dynamical properties for a wide range of solids within PAW have been performed properly [39]. Besides, we carry out GGA-PBE exchange-correlation functional and have proven its effectiveness for calculating the TOECs. As shown in Table 2, our results show that the equilibrium lattice parameters and SOECs perfectly agree with the theoretical values.Next, our concern is to investigate for which range of strainsξ the third-order effects dominant the properties of solids. The maximum strain ξmax is an important parameter in these calculations as well. In the present study, for the second-order term, the fitted coefficient P2 was almost independent of the range of fitting, while the coefficient P3 was more sensitive to ξmax. To illustrate this feature, we also show the curves of the nonlinear elasticity comparison with the linear elasticity for Mg in Figure 2. The results show that linear elasticity is not appropriate when the maximum strain ξmax larger than 30% and the third-order effects must be considered. Therefore, the Lagrangian strain tensor η can be expressed as η=χ−1/2χ2, where χ is the linear strain tensor.


## 4. The Effective Elastic Constants under Pressure and Mechanical Properties for Mg, Be, Ti, Zn, Zr, and Cd

### 4.1. The Effective Elastic Constants

In the case of materials under larger hydrostatic pressure, it is of value to describe the nonlinear elastic properties employing the concept of effective elastic constantsCijP. For many applications, it is sufficient to consider only the terms linear in the external hydrostatic pressure. The five effective SOECs have been obtained from the finite strain theory [1] on the basis of SOECs and TOECs [40]:(12)C11p=C11+η4C11+2C12+C111+C112+ξ−C11+2C13+C113,C12p=C12+ηC111+2C112−C222+2C12+ξ−C12+C123,C13p=C13+ηC113+C123+ξC13+C133,C33p=C33+η4C13−2C33+2C133+ξ5C33+C333,C44p=C44+η1/2C11+C12+C13+C144+C155+ξ1/2C13+C33+C44+C344.The Lagrangian parametersη (along the basal plane) and ξ (along the unique axis) obtained in terms of hydrostatic pressure are(13)η=C11−C33pC11+C12C33−2C132,(14)ξ=2C13−C11−C12pC11+C12C33−2C132.The calculated hydrostatic pressure derivatives of SOECCij′ for Mg, Be, Ti, Zn, Zr, and Cd in terms of our prediction for SOECs and TOECs are shown in Table 4. Here, Cij′ represents dCij′/dp. In [40], the above method has been employed to confirm the pressure dependence of the SOECs. First, applying the hydrostatic pressure to a crystal, then the pressure-dependent elastic constants that can be obtained from the crystal have been additionally deformed. The first-principles calculated results for the total elastic energy combined with the strain-energy relation will enable us to determine CijP and Cij′. Therefore, we believe that the method used in our paper is correct.Table 4
Predictions for the pressure derivatives of second-order elastic constants for Mg, Be, Ti, Zn, Zr, and Cd.
MgBeTiZnZrCdC11′12.217.513.256.912.137.92C12′14.941.531.472.503.242.52C13′1.410.343.883.302.993.68C33′5.624.23-0.323.742.535.32C44′0.271.350.160.76−0.84−2.29


### 4.2. Elastic Anisotropy

In this section, elastic anisotropy of Mg, Be, Ti, Zn, Zr, and Cd has been investigated. It is well known that the acoustic velocities are related to the elastic constants by the Christoffel equationCijklnjnk−Mδilui=0. M is the modulus of propagation and is written as M=ρυ2. ρ is density, υ is the velocity, the fourth rank tensor Cijkl is used to describe elastic constants, and n is the propagation [41]. The acoustic anisotropy is defined as [42] Δi=Minx/Mi100, where i is the index for the three types of elastic waves (one longitudinal and the two polarizations of the shear wave) [42] and nx is the extremal propagation direction except [100]. The anisotropy of the compressional wave P can be obtained by solving the Christoffel equation for a hexagonal lattice as Δp=C33/C11. For the shear waves, the anisotropies of the wave polarized perpendicular to the basal plane S1 are Δs1=C11+C33−2C13/4C44 and the ones polarized in the basal plane S2 are Δs2=2C44/C11−C12. For S1, the extremum occurs at an angle of 45° from the c axis in the a−c plane, and for S2 and P waves, it occurs along the c axis. The calculated elastic anisotropy factors, Δp, Δs1, and Δs2 for Mg, Be, Ti, Zn, Zr, and Cd as a function of the applied pressure are shown in Figure 3. The five-pointed star, triangle-right, circle, triangle-up, diamond, and square represent Mg, Be, Ti, Zn, Zr, and Cd, respectively. We find that changing tendencies of the elastic anisotropy Δp, Δs1 and Δs2 in the pressure range 0–5 GPa is similar for all six metals investigated here. It is easy to see that they become small with the increase in pressure. From Figure 3(a), the anisotropy of the compression wave Δp for Zr, Cd, and Zn increases gradually with the increase of the pressure, while the variation trends of Δp for Be, Ti, and Mg are the opposite. Besides, Δp for all six metals from big to small, arrange Be>Zr>Ti>Mg>Cd>Zn. It is seen that the anisotropy of the shear waves Δs1 and Δs2 varies gently generally in the wide range of pressure. However, Δs1 of Cd and Δs2 for Mg increase rapidly, particularly over 3.5 GPa, with the increase in pressure in Figures 3(b) and 3(c). The different change tendencies are determined by the anharmonicity of the acoustic vibrations.Figure 3
Elastic anisotropy factors varied with the applied pressures. (a–d) describeΔp, Δs1, Δs2, and AU respectively. The five-pointed star, triangle-right, circle, triangle-up, diamond and square represent Mg, Be, Ti, Zn, Zr and Cd respectively.
(a)(b)(c)(d)In order to investigate the elastic anisotropy of HCP metals systematically, the universal anisotropyAU has been investigated as well. The universal anisotropy AU introduced by Shivakumar et al. can be used to measure the anisotropy of hexagonal, trigonal, and monoclinic [43]. For an elastically isotropic solid, Δp=Δs1=Δs2=1, while AU is zero. AU=5GV/GR+5BV/BR−6≥0, which can be obtained from [43]. GV and GR are the Voigt [44] shear modulus and Reuss [45] shear modulus, and BV and BR are the Voigt bulk modulus and Reuss bulk modulus, respectively. For the hexagonal structure, BV=1/92C11+C12+4C13+C33,GV=1/30M+12C44+12C66,BR=C2/M,GR=5/2C2C44C66/3BVC44C66+C2C44+C66M=C11+C12+2C33−4C13,C2=C11+C12C33−2C122. The calculated universal anisotropy AU for Mg, Be, Ti, Zn, Zr, and Cd as a function of the applied pressure is shown in Figure 3(d). It is seen that AU for Be, Ti, and Zn varies gently in the pressure range of 0–5 GPa. AU for Cd and Zr becomes small with the increase of the pressure, while the universal anisotropy AU of our calculation for Mg increases rapidly as the pressure increases.


### 4.3. Ductile-to-Brittle Criterion and Vickers Hardness

A universal ductile-to-brittle criterion: to evaluate material ductility or brittleness, the classical criteria of pressure and of Pugh’s modulus ratioG/B were proposed by Pugh et al. [46]. G indicates shear modulus and B bulk modulus. According to Voigt [44] and Reuss [45], G and B can be calculated from elastic constants C11, C12, C13, C33, and C44 by using the relations G=GV+GR/2 and B=BV+BR/2. G/B=0.571 corresponds to critical Pugh’s modulus ratio defined by Pugh [46]. Brittleness and ductility for Mg, Be, Ti, Zn, Zr, and Cd as functions of pressure are presented in Figure 4. It is seen that Figure 4 is divided into three regions, I, II, and III, which represent ductility region, transition region, and brittleness region, respectively. The material is ductile when the ratio is below the critical value of 0.5, and the material is considered brittle when the ratio G/B is larger than 0.571. It is evident that Mg, Ti, Zr, and Cd belong to the ductility region, in which metallic bonding is stronger, and the majority of Zn is between ductility region and brittleness region, while Be in the brittleness region. Furthermore, Pugh’s modulus ratio G/B for Mg, Be, Ti, Zn, Zr, and Cd decreases with the increase in pressure. In other words, for Mg, Be, Ti, Zn, Zr, and Cd, larger pressure leads to higher ductility, namely, lower brittleness.Figure 4
Brittleness and ductility as functions of pressure. The material is ductile when the ratio is below the critical value of 0.5, and the material is considered brittle when the ratio G/B is larger than 0.571. The five-pointed star, triangle-right, circle, triangle-up, diamond, and square represent Mg, Be, Ti, Zn, Zr, and Cd, respectively.Vickers hardnessHv for HCP metals under pressure ranges 0–5 GPa. Hardness is defined as the resistance of a material to deformation. Recently, Chen et al. proposed a new hardness model [47, 48], Hv=2k2G0.585−3, where k=G/B. We can obtain that the hardness is connected not only with the bulk modulus but also with the shear modulus. Vickers hardness Hv as a function of pressure is presented in Figure 5(a). However, we find that the hardness values of Mg and Zr are negative. Generally speaking, the hardness values are positive, so we adopt another hardness model Hv=1.877k2G0.585 [47], as shown in Figure 5(b). This model is also obtained by Chen et al. via investigating the hardness of polycrystalline materials and bulk metallic glasses and analyzing these data. Comparing Figure 5(a) with Figure 5(b), it is easy to find that variation trend of hardness value as the change of pressure is approximately consistent. From Figure 5(b), we can obtain that the hardness values of all six metals are mainly linearly related to the pressure. The hardness Hv for Be and Cd varies by a small amount, while that for Mg, Ti, Zn, and Zr decreases as the pressure increases. Besides, the hardness value of Hv≈48.1 (at zero pressure) is the biggest among six kinds of HCP metals.Figure 5
Vickers hardnessHv as a function of pressure, in which (a, b) adopt various methods. The five-pointed star, triangle-right, circle, triangle-up, diamond, and square represent Mg, Be, Ti, Zn, Zr, and Cd, respectively.
(a)(b)


## 4.1. The Effective Elastic Constants

In the case of materials under larger hydrostatic pressure, it is of value to describe the nonlinear elastic properties employing the concept of effective elastic constantsCijP. For many applications, it is sufficient to consider only the terms linear in the external hydrostatic pressure. The five effective SOECs have been obtained from the finite strain theory [1] on the basis of SOECs and TOECs [40]:(12)C11p=C11+η4C11+2C12+C111+C112+ξ−C11+2C13+C113,C12p=C12+ηC111+2C112−C222+2C12+ξ−C12+C123,C13p=C13+ηC113+C123+ξC13+C133,C33p=C33+η4C13−2C33+2C133+ξ5C33+C333,C44p=C44+η1/2C11+C12+C13+C144+C155+ξ1/2C13+C33+C44+C344.The Lagrangian parametersη (along the basal plane) and ξ (along the unique axis) obtained in terms of hydrostatic pressure are(13)η=C11−C33pC11+C12C33−2C132,(14)ξ=2C13−C11−C12pC11+C12C33−2C132.The calculated hydrostatic pressure derivatives of SOECCij′ for Mg, Be, Ti, Zn, Zr, and Cd in terms of our prediction for SOECs and TOECs are shown in Table 4. Here, Cij′ represents dCij′/dp. In [40], the above method has been employed to confirm the pressure dependence of the SOECs. First, applying the hydrostatic pressure to a crystal, then the pressure-dependent elastic constants that can be obtained from the crystal have been additionally deformed. The first-principles calculated results for the total elastic energy combined with the strain-energy relation will enable us to determine CijP and Cij′. Therefore, we believe that the method used in our paper is correct.Table 4
Predictions for the pressure derivatives of second-order elastic constants for Mg, Be, Ti, Zn, Zr, and Cd.
MgBeTiZnZrCdC11′12.217.513.256.912.137.92C12′14.941.531.472.503.242.52C13′1.410.343.883.302.993.68C33′5.624.23-0.323.742.535.32C44′0.271.350.160.76−0.84−2.29


## 4.2. Elastic Anisotropy

In this section, elastic anisotropy of Mg, Be, Ti, Zn, Zr, and Cd has been investigated. It is well known that the acoustic velocities are related to the elastic constants by the Christoffel equationCijklnjnk−Mδilui=0. M is the modulus of propagation and is written as M=ρυ2. ρ is density, υ is the velocity, the fourth rank tensor Cijkl is used to describe elastic constants, and n is the propagation [41]. The acoustic anisotropy is defined as [42] Δi=Minx/Mi100, where i is the index for the three types of elastic waves (one longitudinal and the two polarizations of the shear wave) [42] and nx is the extremal propagation direction except [100]. The anisotropy of the compressional wave P can be obtained by solving the Christoffel equation for a hexagonal lattice as Δp=C33/C11. For the shear waves, the anisotropies of the wave polarized perpendicular to the basal plane S1 are Δs1=C11+C33−2C13/4C44 and the ones polarized in the basal plane S2 are Δs2=2C44/C11−C12. For S1, the extremum occurs at an angle of 45° from the c axis in the a−c plane, and for S2 and P waves, it occurs along the c axis. The calculated elastic anisotropy factors, Δp, Δs1, and Δs2 for Mg, Be, Ti, Zn, Zr, and Cd as a function of the applied pressure are shown in Figure 3. The five-pointed star, triangle-right, circle, triangle-up, diamond, and square represent Mg, Be, Ti, Zn, Zr, and Cd, respectively. We find that changing tendencies of the elastic anisotropy Δp, Δs1 and Δs2 in the pressure range 0–5 GPa is similar for all six metals investigated here. It is easy to see that they become small with the increase in pressure. From Figure 3(a), the anisotropy of the compression wave Δp for Zr, Cd, and Zn increases gradually with the increase of the pressure, while the variation trends of Δp for Be, Ti, and Mg are the opposite. Besides, Δp for all six metals from big to small, arrange Be>Zr>Ti>Mg>Cd>Zn. It is seen that the anisotropy of the shear waves Δs1 and Δs2 varies gently generally in the wide range of pressure. However, Δs1 of Cd and Δs2 for Mg increase rapidly, particularly over 3.5 GPa, with the increase in pressure in Figures 3(b) and 3(c). The different change tendencies are determined by the anharmonicity of the acoustic vibrations.Figure 3
Elastic anisotropy factors varied with the applied pressures. (a–d) describeΔp, Δs1, Δs2, and AU respectively. The five-pointed star, triangle-right, circle, triangle-up, diamond and square represent Mg, Be, Ti, Zn, Zr and Cd respectively.
(a)(b)(c)(d)In order to investigate the elastic anisotropy of HCP metals systematically, the universal anisotropyAU has been investigated as well. The universal anisotropy AU introduced by Shivakumar et al. can be used to measure the anisotropy of hexagonal, trigonal, and monoclinic [43]. For an elastically isotropic solid, Δp=Δs1=Δs2=1, while AU is zero. AU=5GV/GR+5BV/BR−6≥0, which can be obtained from [43]. GV and GR are the Voigt [44] shear modulus and Reuss [45] shear modulus, and BV and BR are the Voigt bulk modulus and Reuss bulk modulus, respectively. For the hexagonal structure, BV=1/92C11+C12+4C13+C33,GV=1/30M+12C44+12C66,BR=C2/M,GR=5/2C2C44C66/3BVC44C66+C2C44+C66M=C11+C12+2C33−4C13,C2=C11+C12C33−2C122. The calculated universal anisotropy AU for Mg, Be, Ti, Zn, Zr, and Cd as a function of the applied pressure is shown in Figure 3(d). It is seen that AU for Be, Ti, and Zn varies gently in the pressure range of 0–5 GPa. AU for Cd and Zr becomes small with the increase of the pressure, while the universal anisotropy AU of our calculation for Mg increases rapidly as the pressure increases.


## 4.3. Ductile-to-Brittle Criterion and Vickers Hardness

A universal ductile-to-brittle criterion: to evaluate material ductility or brittleness, the classical criteria of pressure and of Pugh’s modulus ratioG/B were proposed by Pugh et al. [46]. G indicates shear modulus and B bulk modulus. According to Voigt [44] and Reuss [45], G and B can be calculated from elastic constants C11, C12, C13, C33, and C44 by using the relations G=GV+GR/2 and B=BV+BR/2. G/B=0.571 corresponds to critical Pugh’s modulus ratio defined by Pugh [46]. Brittleness and ductility for Mg, Be, Ti, Zn, Zr, and Cd as functions of pressure are presented in Figure 4. It is seen that Figure 4 is divided into three regions, I, II, and III, which represent ductility region, transition region, and brittleness region, respectively. The material is ductile when the ratio is below the critical value of 0.5, and the material is considered brittle when the ratio G/B is larger than 0.571. It is evident that Mg, Ti, Zr, and Cd belong to the ductility region, in which metallic bonding is stronger, and the majority of Zn is between ductility region and brittleness region, while Be in the brittleness region. Furthermore, Pugh’s modulus ratio G/B for Mg, Be, Ti, Zn, Zr, and Cd decreases with the increase in pressure. In other words, for Mg, Be, Ti, Zn, Zr, and Cd, larger pressure leads to higher ductility, namely, lower brittleness.Figure 4
Brittleness and ductility as functions of pressure. The material is ductile when the ratio is below the critical value of 0.5, and the material is considered brittle when the ratio G/B is larger than 0.571. The five-pointed star, triangle-right, circle, triangle-up, diamond, and square represent Mg, Be, Ti, Zn, Zr, and Cd, respectively.Vickers hardnessHv for HCP metals under pressure ranges 0–5 GPa. Hardness is defined as the resistance of a material to deformation. Recently, Chen et al. proposed a new hardness model [47, 48], Hv=2k2G0.585−3, where k=G/B. We can obtain that the hardness is connected not only with the bulk modulus but also with the shear modulus. Vickers hardness Hv as a function of pressure is presented in Figure 5(a). However, we find that the hardness values of Mg and Zr are negative. Generally speaking, the hardness values are positive, so we adopt another hardness model Hv=1.877k2G0.585 [47], as shown in Figure 5(b). This model is also obtained by Chen et al. via investigating the hardness of polycrystalline materials and bulk metallic glasses and analyzing these data. Comparing Figure 5(a) with Figure 5(b), it is easy to find that variation trend of hardness value as the change of pressure is approximately consistent. From Figure 5(b), we can obtain that the hardness values of all six metals are mainly linearly related to the pressure. The hardness Hv for Be and Cd varies by a small amount, while that for Mg, Ti, Zn, and Zr decreases as the pressure increases. Besides, the hardness value of Hv≈48.1 (at zero pressure) is the biggest among six kinds of HCP metals.Figure 5
Vickers hardnessHv as a function of pressure, in which (a, b) adopt various methods. The five-pointed star, triangle-right, circle, triangle-up, diamond, and square represent Mg, Be, Ti, Zn, Zr, and Cd, respectively.
(a)(b)


## 5. Conclusions

In this work, we have described a systematic scheme to calculate the SOECs and TOECs for Mg, Be, Ti, Zn, Zr, and Cd using the homogeneous deformation method and first-principles calculations. Our calculated values of SOECs and TOECs for Mg, Be, Ti, Zn, Zr, and Cd are listed in Tables2 and 3, respectively. To verify the reliability of the present method, we have compared our calculated results for the equilibrium lattice parameters and SOECs with previous theoretical and experimental results, and our results are in perfect agreement with the previous available calculated and experimental values. The pressure derivatives of SOECs Cij′ are calculated by using the obtained TOECs. Besides, high-pressure effects on elastic anisotropy, ductile-to-brittle criterion, and Vickers hardness are also investigated. It is easy to find that the hardness model Hv=1.877k2G0.585 is more appropriate than Hv=2k2G0.585−3 for Mg, Be, Ti, Zn, Zr, and Cd under the pressure range 0–5 GPa. A method of the nonlinear theory of elasticity for HCP metals must accelerate the improvement of experimental methods of measuring TOECs. We believe that our first-principle computation of TOECs can be a very useful tool in applying the material with the hexagonal structure to practical engineering.



---

*Source: 8726250-2021-12-03.xml*