Computational prediction of ideal strength for a material

The ideal strength is crucial for predicting material behavior under extreme conditions, which can provide insights into material limits, guide design and engineer for enhanced performance and durability. In this work, we present a method within an allows for the estimation of tensile, shear, and indentation strengths in any crystallographic direction or plane. We have examined the strain-stress relationships of several well-known structures and compared our findings with previous work, demonstrating the effectiveness of our approach. Moreover, we performed extensive investigations into the indentation strength of hexagonal WC, $\beta$-SiC, and MgAl₂O₄. The current study uncovers the modes of structural deformation and the underlying atomistic mechanisms. The insights gained from this study have significant implications for the further exploration and design of superhard materials.

Keywords: First-principles calculations; Ideal strength; Tensile strength; Shear strength; Indentation strength.

The ideal elasticity strength is the minimum stress to plastically deform an infinite dislocation-free crystal [1], which gives an upper bound on the strength of a perfect crystal. The ideal strength could be classified as tensile, shear, and indentation strength according to the type of external force posted on the material. It is an important parameter to evaluate the shear and indentation strengths are believed to be closely related to the Vickers hardness of an ideal crystal. Therefore, it is of technical and scientific interest to investigate the ideal elasticity strength of a crystal. In industrial practice, stress-strain relationships are commonly employed to determine ideal strength [2], where various strains are applied to the crystal and the corresponding stress is measured. Theoretical simulation of the ideal strength under complex loading conditions by using stress-strain relationship can also provide fundamental insights into the evolution of chemical bonding and physical property changes of a materials, which is helpful for elucidating and predicting anomalous behaviors of material at extreme conditions [2, 3, 4, 5].

Great efforts had been dedicated to the accurate calculation of the ideal strength since the mid of 1980s[1]. Initially, some experimentally proposed semi-empirical models were used. It is commonly believed that the tensile and shear strength of materials are directly proportional to their Vickers hardness. However, the results are always inaccurate. Recent advances in computational physics made it possible to ab-initially calculate the ideal strength, and the ideal elasticity strength of tremendous materials has been precisely calculated within the framework of electronic structure[1, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29]. Indentation is a method used in experiments to measure the hardness and strength of materials. A more complex calculation of indentation strength has been performed in the framework of DFT [16, 27, 28]. Energy and force can be calculated faster and more accurately using the DFT approach. Therefore, this calculation can be applied in more scenarios. The calculations of the elastic strength have been extensively applied to an increasing number of materials, solving a wide range of engineering and physical problems [1, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29]. Li et al. investigated boron-rich tungsten borides, and demonstrated novel ultrahard properties with stress-strain behaviors revealing diverse and anomalous patterns. The profound influence of boron concentration on bonding configurations and deformation modes in WB_n (n = 2, 3, 4) results in distinct stress responses and unexpected variations in indentation strength, highlighting the unique deformation mechanisms and the need to explore unconventional structure-property relations in ultrahard materials.[21]Lu et al. designed as new-generation superhard materials, transition-metal light-element compounds have been hindered by indentation strain softening, limiting their intrinsic hardness below 40 GPa. Moreover, the calculations revealed that hP4-WN and hP6-WN₂ tungsten nitrides exhibit extraordinary strain stiffening, leading to enhanced indentation strengths exceeding 40 GPa, opening possibilities for nontraditional superhard materials.[24]

In this work, we provided a comprehensive exposition of the underlying principles of our newly proposed methodology on calculating the tensile, shear and indentation strengths based on Vienna Ab initio Simulation Package (VASP) code, accompanied by comparative analysis with other approaches to establish its robustness. We selected hexagonal WC, $\beta$-SiC, and MgAl₂O₄ for an in-depth investigation of ideal strength. Additionally, we examined the indentation strength of hexagonal WC, $\beta$-SiC, and MgAl₂O₄.This comprehensive exploration has afforded us a deeper comprehension of their mechanical characteristics.

Mathematically, the stress tensors are the second-order symmetric tensor with nine components (three positive stress componentsand six shear stress components), in which only six of them are independent. They also can be represented by a 3×3 matrix, where each matrix element corresponds to a stress component. To calculate the stress-strain relationship, it should strain the crystal in a series of incremental simple tensile or shear, and then relax the lattice as a function of the strain to obtain the total energy E_tot. According to DFT, E_tot can be obtained by

where $\Psi$_i represents wavefunction, V_ext is external potential, E_H($\rho$) is the Hartree energy, E_xc($\rho$) is exchange and correlation energy, and U_II is the ion-ion interaction. $\rho$(r) is the density of the electron gas, which can be calculated by

Based on Hellmann-Feynman force theorem, each stress component is the derivative of the energy with respect to the strain

where H represents the Hamiltonian, and $\varepsilon$_mn represents strain.

According to the above method, VASP can simply determine the six independent stress tensors $\tau$_xx, $\tau$_yy, $\tau$_zz, $\tau$_xy, $\tau$_xz, and $\tau$_yz by a structural optimizations, if a strain is induced in a well optimized crystal lattice by modifying the lattice parameters, which are represented by a lattice matrix in VASP code. For instance, if calculating the tensile stress along [100] direction in fcc lattice, one just needs to add the tensile strain along the [100] direction (or x-axis of crystal lattice) by changing the corresponding lattice matrix element, and then optimize the structure by fixing the F_xx component of the optimization matrix F_mn (m, n = x, y, z) in VASP code without optimization. Then we can obtain the tensile stress of $\tau$_xx which corresponds to the studied tensile stress $\tau$_[100]. Considering that VASP can straightforwardly compute the stress tensor of $\tau$_xx, $\tau$_yy, $\tau$_zz, if the tensile strain was added to the coordinate axis of crystal lattice. Therefore, we rotate the studied tensile direction [uvw] to one of coordinate axis (here, we choose x-axis). In this way, after structural optimization by setting F_xx to zero, we can obtained the tensile stress of $\tau$_xx which corresponds to the studied tensile stress $\tau$_[uvw]. By gradually applying a series of continuous strain along the tensile direction and ultimately deriving the corresponding tensile stress, we can obtain the tensile stress-strain relationship.

It is obvious that during the tensile strength calculations, one important step is to rotate the studied tensile direction to the x-axis. Here, take the orthorhombic lattice as an example to introduce how to rotate the studied tensile direction to x-axis. Using Cartesian coordinate system, the actual stretching direction is OA as shown in Fig. 1(a). It can be achieved by rotating -$\alpha$ degree (here counterclockwise rotation is set to be positive direction) around the z-axis (see Fig. 1b), then rotate $\beta$ degree around the y-axis.

Shear strength is more complicated than tensile strength. Shear stress calculations need to rotate the studied shear direction [uvw] to z-axis and reorient the normal vector of the shear plane to align with x-axis. Consequently, by maintaining F_xz as zero during structural optimization, we can calculate the shear stress $\tau$_xz, which is corresponding to the studied shear stress $\tau$_[uvw]. Here, to give a more vivid illustration on the rotation process of shear plane and shear direction, orthorhombic lattice was also taken as an example. As it is can be seen in Fig. 2(a), the sheared plane is characterized by its normal vector OA and the shear direction OB lies within this plane. Firstly, the normal vector OA was rotated -$\alpha$ degree around z-axis [see Fig. 2(b)], and then $\beta$ degree around y-axis [see Fig. 2(c)]. Subsequently, rotate OA $\gamma$ degree around the x-axis[Fig. 2(d)]. Finally, the normal vector is align with the x-axis, and the shear stress direction aligns with the z-axis.

The indentation strength is a more complex case, where in actual compression the material receives not only tensile stress perpendicular to the plane, but also in-plane shear stress. Experiments have shown that there is a close relationship between nano-scale indentation and the geometry of the indenter [16, 22, 23, 29]. When using the first-principles method to calculate compression, we assume that the material is uniformly deformed within a small volume. And the deformation is related to the shape of the indenter. Such a uniformly deformed lattice is sufficient in the study of the breakage of interactions on the atomic microscopic scale. In the first-principles calculations, we can simplify the situation to a uniform deformed infinite and perfect lattice under a uniform stress field[30]. We assume that a sharp indenter squeezes the material by the way shown in the Fig. 3(a). The tensile stress $\tau$_xx and shear stress $\tau$_xz received by the material under a sharp indenter always maintain the following relationship: $\tau$_xz =$\tau$_xxtan$\theta$, where $\theta$ is determined by the shape of the indenter, as shown in Fig. 3(b).

Based on this relationship, we only need to calculate the shear stress by using the method of calculating shear stress proposed in this work. Note that since we consider the crystal is an infinitely large perfect lattice with uniform deformation, we must keep the angle of the deformation constant in both directions. This can be achieved by modifying the optimization program in VASP code. For a blunt indenter, as long as we know that the direction of the tangent at each point of contact surface between the indenter and the material, we can use the approach of a sharp indenter to calculate the indentation stress at each point of the blunt indenter, as shown in Fig. 4[30].

To thoroughly validate the accuracy and reliability of our approach, we conducted an extensive analysis involving diamond, Si, hexagonal AlN, and $\alpha$-Al₂O₃. The tensile, shear and indentation strengths were calculated separately by setting the relative strain at each step to 0.01 times the elements of the lattice matrix. We used the VASP package with the Perdew-Burke-Ernzerhof (PBE) version of the generalized gradient approximation (GGA) for the exchange-correlation functional and a plane-wave basis set [31, 32, 33, 34]. For the integration in the Brillouin zone, we set the k-points to 6 × 6 × 6 for all the calculated material. The cutoff energy is set to 600 eV, which is sufficient to make the total energy converge to an accuracy of less than 0.001 eV/atom. The paramount objective of our investigation is to establish the veracity of our proposed methodology. In pursuit of this goal, we performed a meticulous comparison of our computed results against those obtained using established methods. This validation process was meticulously conducted under conditions consistent with prior research, ensuring a fair and direct comparison.

The crystallographic configurations of the chosen materials are presented in Table S1, while the outcomes of the tensile, shear, and indentation are elaborated in Tables 1, respectively. The significance of these tables lies in their role as a tangible demonstration of the effectiveness of our method across various loading scenarios. In the case of diamond, our calculations yielded the weakest tensile and shear directions consistent with previous research. The weakest tensile direction is [111], and the weakest shear direction is [111] <$11\overline{2}$>. Additionally, our Vickers hardness value is 95.7 GPa, which closely aligns with the prior result of 97.6 GPa, showcasing minimal deviation. This conforms to the characteristics of diamond as a superhard material. For cases Si, the weakest tensile and shear directions are [111] and [111] <$11\overline{2}$>, respectively. And its indentation strength is 6.1 GPa. Additionally, the weakest tensile direction for hexagonal AlN is [$\overline{1}2\overline{1}0$], and the weakest shear direction is [$\overline{1}2\overline{1}0$] <$10\overline{1}0$>. The indentation strength is 11.7 GPa. Moreover, the weakest tensile and shear directions are [$10\overline{1}0$] and [$10\overline{1}0$] <0001> for $\alpha$-Al₂O₃. The indentation strength is 28.9 GPa. Our computed results also closely match those from the previous work. The specific stress-strain curves for tension, shear, and indentation of these four materials are presented in Figures S1-4, respectively.

We subsequently placed particular emphasis on an in-depth investigation of the ideal strengths of three distinct materials: WC, $\beta$-SiC,and MgAl₂O₄. We would like to calculate the Vickers hardness, so the indenter we consider to be a Vickers indenter with an angle of 136°between the opposing surfaces.

Hexagonal tungsten carbide (WC) is crucial for cutting-edge industries like manufacturing, aerospace, and electronics, thanks to its exceptional properties. Its unique combination of hardness, strength, and thermal stability significantly enhances overall performance and efficiency. We extensively undertook an investigation into the ideal tensile strength, ideal shear strength, and indentation strength of hexagonal tungsten carbide (WC). According to academic standards, Fig.5(a) illustrates several distinctive crystallographic orientations of the WC crystal, which predominantly guided our calculations in this study. We conducted tensile calculations to explore these properties along the [0001], [$10\overline{1}0$] and [$\overline{1}2\overline{1}0]$ crystallographic directions. The tensile strengths were found to be 101.3 GPa, 72.6 GPa, and 63.5 GPa in the [0001], [$10\overline{1}0$] and [$\overline{1}2\overline{1}0]$ orientation, respectively. Fig.5(b) visually presents the stress-strain curves for the investigated directions, highlighting [$\overline{1}2\overline{1}0]$ as the least resistant tensile orientation within WC.

We next undertook an in-depth examination involving shear calculations of properties along the [$\overline{1}2\overline{1}0]$<0001>, [$10\overline{1}0$]<0001>, and [$\overline{1}2\overline{1}0$] <$10\overline{1}0$> crystallographic orientations. The shear strength along the [$\overline{1}2\overline{1}0]$<0001> direction emerged at an impressive 41.6 GPa, while the [$10\overline{1}0$]<0001> orientation showcased a noteworthy tensile strength of 33.9 GPa. Additionally, the [$\overline{1}2\overline{1}0$] <$10\overline{1}0$> orientation unveiled a distinctive tensile strength of 18.7 GPa. These significant findings are shown in Fig.5(c), elegantly capturing the essence of stress-strain profiles for each respective orientation. Remarkably, [$\overline{1}2\overline{1}0$] <$10\overline{1}0$> represents the weakest shear direction of WC. Therefore, we selected the weakest shear direction as the orientation for its indentation.

A comprehensive analysis of the indentation process was conducted and visualized in Fig.5(d). Intriguingly, as the strain surpassed 0.27 times the crystal length in the aforementioned direction, the chemical bonds between tungsten and carbon atoms suffered rupture due to increased atomic separation. At this point, the indentation strength of WC is 36.8GPa. Its stress-strain curve is displayed in Fig.5 (e). This phenomenon is attributed to the exceeding distance between the atoms, leading to bond fracture. Our study provides valuable insights into the mechanical response of hexagonal WC under various loading conditions and crystallographic orientations, shedding light on the intricate interplay between atomic bonds and mechanical deformation. The calculation results indicate that WC is a relatively hard material with the potential to become a superhard material.

At equilibrium condition, $\beta$-SiC adopts the cubic structure in F-43m symmetry [Fig. 6(a)], where each Si atom is connected to four C atoms, forming a strong sp³ covalent bond unit with Si-C bond lengths of 1.89 Å. We first set out to determine the easy cleavage plane of $\beta$-SiC by simulating its ideal tensile strength along several high-symmetry directions. The calculated results shown in Fig. 6(b) clearly indicate that the weakest peak tensile stress occurs in the [111] direction in $\beta$-SiC, giving rise the (111) easy cleavage planes. The result is 45 GPa in <111> direction for $\beta$-SiC and the Si-C bonds soften considerably under the tensile strain long before the bond breaking. Since a material’s hardness is primarily determined by its ability to resist shear deformation, we examine shear stress responses of $\beta$-SiC in its major low-index planes to establish key benchmarks. We show in Fig.6(c) the stress responses under pure and indentation shear strains along the (111) [11$\overset{\overline{}}{2}$] direction. The two sets of shear stresses reach peak values 29.7 and 19.7 GPa, respectively. At the indentation strain ($\varepsilon$ = 0.36), the stress has only decreased from the peak value. The peak pure shear stress occurs at the strain $\varepsilon$ = 0.22 but the interlayer sp³ bonding persists up to $\varepsilon$ = 0.32 where significant buckling still exists and the Si-C bond length increased to 2.076 Å. The corresponding structural snapshot shows significant buckling in the (111) planes, indicating a strong sp³ bonding character. The structural failure occurs at $\varepsilon$ = 0.38, where the interlayer bonds in the (111) [11$\overset{\overline{}}{2}$] direction break up and the layers become essentially flat, signaling a dominant sp² bonding character. Here we notice that all the Si-C bonds responsible for the structural failure also break simultaneously under shear strains.The calculation results show that $\beta$-SiC possesses favorable mechanical properties.

MgAl₂O₄ holds significant importance across diverse fields due to its remarkable combination of properties. As a transparent ceramic, it finds applications in optical systems, lasers, and transparent armor. Its high-temperature stability, electrical insulating properties, and chemical resistance make it valuable in aerospace, electronics, and corrosion-resistant coatings. Moreover, spinel’s biocompatibility has fueled its adoption in medical implants. Its versatile attributes and extensive utility underscore spinel’s role as a multifunctional material with widespread impact.

Our investigation centered around the characterization of spinel MgAl₂O₄ in terms of its ideal tensile strength, ideal shear strength, and indentation strength. To adhere to academic conventions, we adhered to Fig.7(a), which visually illustrates distinct crystallographic orientations within the MgAl₂O₄ crystal. These orientations served as the foundation for our computations in this study. Tensile calculations were carried out along the [100], [110], and [111] crystallographic directions, yielding corresponding strengths of 89.2 GPa, 28.28 GPa, and 28.3 GPa, respectively. The stress-strain curves for these orientations are depicted in Fig.7(b), with the [111] orientation revealing the least resistance to tensile forces within MgAl₂O₄.

Our subsequent focus was on an exhaustive exploration of shear calculations along the $[100]$<001>, $[110$]<001>, and [111] <$11\overline{2}$> crystallographic orientations. Notably, the shear strength along the$[100]$<001> direction registered an impressive 29.7 GPa, while the $[110$]<001> orientation displayed a remarkable tensile strength of 31.0 GPa. Furthermore, the [111] <$11\overline{2}$> orientation exhibited a distinctive tensile strength of 14.2 GPa. Fig6 (c) elegantly encapsulates these pivotal findings, showcasing stress-strain profiles for each orientation. Interestingly, the [111] <$11\overline{2}$> orientation emerged as the weakest shear direction within MgAl₂O₄, thereby being chosen for indentation analysis.

Our comprehensive analysis of the indentation process is vividly presented in Fig. 7(d). Intriguingly, when the strain surpassed 0.46 times the crystal length in the aforementioned direction, the increased atomic separation led to the rupture of chemical bonds. At this juncture, the indentation strength of MgAl₂O₄ was determined to be 22.2 GPa, as depicted in the stress-strain curve presented in Fig.7(e). This phenomenon was attributed to the exceeding atomic distance, which induced bond fracture. Our study thus yields valuable insights into the mechanical behavior of hexagonal MgAl₂O₄ under varied loading conditions and crystallographic orientations, unraveling the intricate interplay between atomic bonds and mechanical deformation. Its hardness enables it to be applied in more diverse applications.

In summary, we introduced a novel methodology for computing the ideal strength of materials based on first-principles calculations. The approach enabled the determination of tensile, shear, and indentation strengths in arbitrary crystallographic directions, facilitating a comprehensive understanding of material response under different loading conditions. Through a meticulous comparison with prior works, we demonstrated the accuracy and reliability of our proposed method. Furthermore, our investigation delved into the specific cases of hexagonal WC, cubic SiC, and spinel MgAl₂O₄, shedding light on their ideal strengths and the structural transformations that occur under applied stress. Based on our computational results, the Vickers hardness of WC reaches 36.8 GPa, which is closed to 40 GPa. Therefore, it could not be classified as a superhard material. While WC has the potential to become a superhard material through various methods, such as nanocrystallization, high-pressure synthesis, and alloying. Meanwhile, cubic SiC and MgAl₂O₄ exhibit Vickers hardness within the range of 20-30 GPa, rendering them suitable for diverse applications in various fields. This detailed analysis not only enriched our understanding of these materials’ mechanical properties but also highlighted the versatility of our methodology in exploring a diverse range of crystalline structures. The findings presented in this study have broad implications for material design and engineering, offering a rigorous foundation for predicting and optimizing the mechanical behaviors of various materials. As we continue to refine and expand our computational techniques, the insights gained from this research will contribute to advancements across multiple disciplines, driving innovation in materials science and engineering.

This work is supported by the Natural Science Foundation of China (Grant No.12074138), and the Jilin Provincial Science and Technology Development Joint Fund Project (Grant No. YDZJ202201ZYTS581).