#### 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 gradient is applied to the crystal lattice vectors , 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 is more appropriate than for HCP metals under high pressure.

#### 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 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 . After introducing the deformation gradient ,where and (=1, 2, 3) represent the Cartesian coordinates. Then, 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. 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]:where is the volume of the system, is the entropy, and is the corresponding internal energy of the ground state. An expansion 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

And the isothermal elastic constants arewhich 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 , , then . We will not distinguish between isentropic elastic constants and isothermal elastic constants 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, . So, the Lagrangian strain tensors can be written as follows:

Equations (3) and (5) now can be written aswhere 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 stress is defined as the first-order 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.

#### 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 gradient is applied to the crystal lattice vectors 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 needs to be obtained. is related to the strain and is determined by inverting equation (2):

In general, the deformation gradient 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 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 and ten TOECs , , , , . 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 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 , = , , …, 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 and represent combinations of second-order elastic constants and third-order elastic constants of the crystal, respectively. Furthermore, coefficients and 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 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 eV and 1.0 eV/Å, respectively.

Taking Mg as an example, Figure 1(a) shows the relationship between the internal energy convergence and the -point grid size. The internal energy is well converged after mesh size. The test of the cutoff energy affects the internal energy convergence displayed in Figure 1(b), and = 350 eV is sufficient for Mg. The equilibrium lattice parameters and 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 and . It shows that our values of calculation agree well with the previous experimental and calculated values (see Table 2).

**(a)**

**(b)**

**(c)**

##### 3.2. Results and Discussion

We list our calculated values of second-order elastic constants (SOECs) and third-order elastic constants (TOECs) 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 values are slightly different for different fitted curves. Taking Mg as an example, from coefficients in and is 64.02 GPa and 63.90 GPa, respectively; the values are 17.15 GPa, 17.03 Gpa, and 16.88 GPa from coefficients in , , and , 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 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 are negative (see Table 3), while it is the opposite for SOECs (see Table 2).

**(a)**

**(b)**

**(c)**

**(d)**

**(e)**

**(f)**

**(g)**

**(h)**

**(i)**

**(j)**

The determination of 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 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 is an important parameter in these calculations as well. In the present study, for the second-order term, the fitted coefficient was almost independent of the range of fitting, while the coefficient was more sensitive to . 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 larger than 30% and the third-order effects must be considered. Therefore, the Lagrangian strain tensor can be expressed as , 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 constants . 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]:

The Lagrangian parameters (along the basal plane) and (along the unique axis) obtained in terms of hydrostatic pressure are

The calculated hydrostatic pressure derivatives of SOEC for Mg, Be, Ti, Zn, Zr, and Cd in terms of our prediction for SOECs and TOECs are shown in Table 4. Here, represents . 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 and . Therefore, we believe that the method used in our paper is correct.

##### 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 equation . is the modulus of propagation and is written as . is density, is the velocity, the fourth rank tensor is used to describe elastic constants, and is the propagation [41]. The acoustic anisotropy is defined as [42] , where is the index for the three types of elastic waves (one longitudinal and the two polarizations of the shear wave) [42] and is the extremal propagation direction except [100]. The anisotropy of the compressional wave can be obtained by solving the Christoffel equation for a hexagonal lattice as . For the shear waves, the anisotropies of the wave polarized perpendicular to the basal plane are and the ones polarized in the basal plane are . For , the extremum occurs at an angle of from the axis in the plane, and for and waves, it occurs along the axis. The calculated elastic anisotropy factors, , , and 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 , and 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 for Zr, Cd, and Zn increases gradually with the increase of the pressure, while the variation trends of for Be, Ti, and Mg are the opposite. Besides, for all six metals from big to small, arrange . It is seen that the anisotropy of the shear waves and varies gently generally in the wide range of pressure. However, of Cd and 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.

**(a)**

**(b)**

**(c)**

**(d)**

In order to investigate the elastic anisotropy of HCP metals systematically, the universal anisotropy has been investigated as well. The universal anisotropy introduced by Shivakumar et al. can be used to measure the anisotropy of hexagonal, trigonal, and monoclinic [43]. For an elastically isotropic solid, , while is zero. , which can be obtained from [43]. and are the Voigt [44] shear modulus and Reuss [45] shear modulus, and and are the Voigt bulk modulus and Reuss bulk modulus, respectively. For the hexagonal structure, . The calculated universal anisotropy 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 for Be, Ti, and Zn varies gently in the pressure range of 0–5 GPa. for Cd and Zr becomes small with the increase of the pressure, while the universal anisotropy 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 ratio were proposed by Pugh et al. [46]. indicates shear modulus and bulk modulus. According to Voigt [44] and Reuss [45], and can be calculated from elastic constants , , , , and by using the relations and . 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 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 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 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], , where . We can obtain that the hardness is connected not only with the bulk modulus but also with the shear modulus. Vickers hardness 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 [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 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 (at zero pressure) is the biggest among six kinds of HCP metals.

**(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 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. The pressure derivatives of SOECs 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 is more appropriate than 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

The data used to support the findings of this study are included in the article.

#### Conflicts of Interest

The authors declare that they have no conflicts of interest.

#### Acknowledgments

This research was funded by the Fundamental Research Funds for the Central Universities (Grant no. 2019CDXYTX0023); Chongqing Natural Science Foundation (Grant no. cstc2018jscx-msyb1002); and the Doctoral Startup Fund of Chongqing University of Posts and Telecommunications (Grant no. E012A2021224).