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

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. ,e method of homogeneous deformation combined with first-principles total-energy calculations is employed. ,e deformation gradient Fij is applied to the crystal lattice vectors ri, and the elastic strain energy can be obtained from the firstprinciples calculation. ,e secondand third-order elastic constants are extracted by a polynomial fit to the calculated energystrain 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. ,e results show that the hardness model Hv � 1.877(k2G) 0.585 is more appropriate than Hv � 2(k2G) 0.585 − 3 for HCP metals under high pressure.


Introduction
In the theory of linear elasticity, infinitesimal deformation strains are assumed, and using finite deformation on materials is very significant for practical application. e 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][3][4][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. ere are quite a few of the available theoretical methods for determining TOECs, including empirical force-constant models [9][10][11][12][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][22][23][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]. is 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. is 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.

Nonlinear Elasticity Theory
In this work, we will recall some basic facts from the nonlinear theory of elasticity [5,[26][27][28][29][30]. Let a i 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 x i . After introducing the deformation gradient F ij , where i and j (�1, 2, 3) represent the Cartesian coordinates. en, we may define the Lagrangian strain tensor η: 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. e 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]: where V is the volume of the system, S is the entropy, and U(V, 0, S) is the corresponding internal energy of the ground state. An expansion U(V, η, S) in terms of symmetric Lagrangian strains tensor is appropriate because the internal energy is unaltered under rigid rotation [2]. e SOECs, TOECs, and HOECs are defined as the second-, third-, and higher-order derivatives of equation (3) with regard to the Lagrangian strain, respectively. us, the isentropic elastic constants can be written as And the isothermal elastic constants are 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 C S � C T . We will not distinguish between isentropic elastic constants C S and isothermal elastic constants C T 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: Equations (3) and (5) now can be written as 2 Advances in Materials Science and Engineering where ΔU(V, η)/V � [U(V, η) − U(V, 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. en, SOECs and TOECs are extracted by fitting equation (8) to the first-principles calculated energy-strain results.
Besides, the Lagrangian stress t is defined as the firstorder derivative of the internal energy with regard to the Lagrangian strain tensor η: 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.

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). e deformation gradient F ij is applied to the crystal lattice vectors r i to obtain the unit cell for the strained crystal, and the deformed lattice vectors can be determined: To carry out the different deformation modes in our work, the deformation gradient matrix F ij needs to be obtained. F ij is related to the strain η and is determined by inverting equation (2): In general, the deformation gradient F ij is not single for a given strain η, while the various possible solutions differ from one another by a rigid rotation. e lack of a one-toone relationship between the deformation gradient matrix F ij 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 SOECs (C 11 , C 12 , C 13 , C 14 , C 44 ) and ten TOECs (C 111 , C 222 , C 333 , C 112 , C 113 , C 123 , C 133 , C 144 , C 155 , C 344 . 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 ξ: Coefficients P 2 and P 3 represent combinations of second-order elastic constants and third-order elastic constants of the crystal, respectively. Furthermore, coefficients P 2 and P 3 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][32][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). e projector augmented wave (PAW) method has been carried out in the VASP package [36]. e 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. e 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, Figure 1(a) shows the relationship between the internal energy convergence and the k-point grid size. e internal energy is well converged after 15 × 15 × 15 mesh size. e test of the cutoff energy affects the internal energy convergence displayed in Figure 1(b), and E \spmathrmMg \spmathrmcutoff � 350 eV is sufficient for Mg. e 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).

Results and Discussion.
We list our calculated values of second-order elastic constants (SOECs) C ij and third-order elastic constants (TOECs) C ijk 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. e C ij values are slightly different for different fitted curves. Taking Mg as an example, C 11 from coefficients in f A (ξ) and f D (ξ) Table 1: e coefficients P 2 and P 3 in equation (12) of corresponding Lagrangian strains are shown as combinations of second-order and third-order elastic constants for the hexagonal crystal.   Table 2. It is shown that our results of SOECs agree well with the available data from other theoretical and experimental values. e energies-strain (E − ξ) curves for Mg, containing the results of the first-principles calculations and the fitted polynomials, are plotted in Figure 2. e 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 C ijk are negative (see Table 3), while it is the opposite for SOECs (see Table 2).

Strain type
e determination of C ijk 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. e 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. e maximum strain ξ max is an important parameter in these calculations as well. In the present study, for the second-order term, the fitted coefficient P 2 was almost independent of the range of fitting, while the coefficient P 3 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. e results show that linear elasticity is not appropriate when the maximum strain ξ max larger than 30% and the third-order effects must be considered. erefore, the Lagrangian strain tensor η can be expressed as η � χ − 1/2χ 2 , where χ is the linear strain tensor.     the nonlinear elastic properties employing the concept of effective elastic constants C ij (P). For many applications, it is sufficient to consider only the terms linear in the external hydrostatic pressure. e five effective SOECs have been obtained from the finite strain theory [1] on the basis of SOECs and TOECs [40]: 13 (p) � C 13 + η C 113 + C 123 + ξ C 13 + C 133 , C 33 (p) � C 33 + η 4C 13 − 2C 33 + 2C 133 + ξ 5C 33 + C 333 , C 44 (p) � C 44 + η 1/2 C 11 + C 12 + C 13 + C 144 + C 155 + ξ 1/2 C 13 + C 33 + C 44 + C 344 . (12) e Lagrangian parameters η (along the basal plane) and ξ (along the unique axis) obtained in terms of hydrostatic pressure are

The Effective Elastic Constants under
e calculated hydrostatic pressure derivatives of SOEC C ij ′ for Mg, Be, Ti, Zn, Zr, and Cd in terms of our prediction for SOECs and TOECs are shown in Table 4. Here, C ij ′ represents dC ij ′ /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. e first-principles calculated results for the total elastic energy combined with the strainenergy relation will enable us to determine C ij (P) and C ij ′ . erefore, we believe that the method used in our paper is correct.

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 equation (C ijkl n j n k − Mδ il )u i � 0. M is the modulus of propagation and is written as M � ρυ 2 . ρ is density, υ is the velocity, the fourth rank tensor C ijkl is used to describe elastic constants, and n is the propagation [41]. e acoustic anisotropy is defined as [42] , where i is the index for the three types of elastic waves (one longitudinal and the two polarizations of the shear wave) [42] and n x is the extremal propagation direction except [100]. e anisotropy of the compressional wave (P) can be obtained by solving the Christoffel equation for a hexagonal lattice as Δ p � C 33 /C 11 . For the shear waves, the anisotropies of the wave polarized perpendicular to the basal plane (S 1 ) are Δ s 1 � (C 11 + C 33 − 2C 13 )/(4C 44 ) and the ones polarized in the basal plane (S 2 ) are Δ s 2 � 2C 44 /(C 11 − C 12 ). For S 1 , the extremum occurs at an angle of 45°from the c axis in the a − c plane, and for S 2 and P waves, it occurs along the c axis. e calculated elastic anisotropy factors, Δ p , Δ s 1 , and Δ s 2 for Mg, Be, Ti, Zn, Zr, and Cd as a function of the applied pressure are shown in Figure 3. e 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 , Δ s 1 and Δ s 2 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 Δ s 1 and Δ s 2 varies gently generally in the wide range of pressure. However, Δ s 1 of Cd and Δ s 2 for Mg increase rapidly, particularly over 3.5 GPa, with the increase in pressure in  In order to investigate the elastic anisotropy of HCP metals systematically, the universal anisotropy A U has been investigated as well. e universal anisotropy A U introduced by Shivakumar et al. can be used to measure the anisotropy of hexagonal, trigonal, and monoclinic [43]. For an elastically isotropic solid, Δ p � Δ s 1 � Δ s 2 � 1, while A U is zero.
0, which can be obtained from [43]. G V and G R are the Voigt [44] shear modulus and Reuss [45] shear modulus, and B V and B R are the Voigt bulk modulus and Reuss bulk modulus, respectively. For the hexagonal structure, B V � (1/9)[2(C 11 + C 12 ) + 4(C 13 + 12 . e calculated universal anisotropy A U 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 A U for Be, Ti, and Zn varies gently in the pressure range of 0-5 GPa. A U for Cd and Zr becomes small with the increase of the pressure, while the universal anisotropy A U of our calculation for Mg increases rapidly as the pressure increases.

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 ratio G/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 C 11 , C 12 , C 13 , C 33 , and C 44 by using the relations G � (G V + G R )/2 and B � (B V + B R )/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. e 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.
Vickers hardness H v 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], H v � 2(k 2 G) 0.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 (H v ) 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 H v � 1.877(k 2 G) 0.585 [47], as shown in Figure 5(b). is 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. e hardness H v 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 H v ≈ 48.1 (at zero pressure) is the biggest among six kinds of HCP metals.

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 firstprinciples calculations. Our calculated values of SOECs and TOECs for Mg, Be, Ti, Zn, Zr, and Cd are listed in Tables 2   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. e pressure derivatives of SOECs C ij ′ are calculated by using the obtained TOECs. Besides, highpressure effects on elastic anisotropy, ductile-to-brittle criterion, and Vickers hardness are also investigated. It is easy to find that the hardness model H v � 1.877(k 2 G) 0.585 is more appropriate than H v � 2(k 2 G) 0.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.

Data Availability
e data used to support the findings of this study are included in the article.

Conflicts of Interest
e authors declare that they have no conflicts of interest.