Computational discovery of hard and superhard materials

Superhard materials are a class of materials with unique mechanical properties that are widely used in many industrial applications, like mining, defense, and space industries. A material can be called superhard if its Vickers hardness is greater than 40 GPa.^1–3 The hardest material known to date is diamond, with Vickers hardness in the range of 60–120 GPa.^4,5 Some metastable carbon allotropes also have a high Vickers hardness,^6–9 but none of them surpass diamond. Among the noncarbon superhard materials, cubic boron nitride displays the high Vickers hardness of ∼60 GPa.^1,2,10,11 Materials with unique mechanical properties include some borides, carbides, and nitrides of transition metals, such as chromium,^12–14 rhenium,^15–17 molybdenum,^18,19 tungsten,^20–27 etc. Some of these carbides (WC) and nitrides (TiN) are already widely used in manufacturing and mining, e.g., in drilling equipment.

Diamond (as well as other sp³ allotropes of carbon) and cubic BN are high-pressure phases. In this work, we pay more attention to hard and superhard transition metal borides, carbides, and nitrides, many of which are stable at ambient conditions and easy to synthesize. We also cover recent advances in the prediction of new materials that may find a variety of applications. Crystal structures of some hard and superhard materials are shown in Fig. 1.

Titanium nitride's Vickers hardness of 22 GPa,²⁹ though not extraordinary, makes it suitable for various applications. Tungsten carbide (WC) has a higher Vickers hardness of ∼30 GPa,^30,31 enabling its use in “hard alloys” (WC-based composites) that are actively used in drilling, metal- and woodworking (in drawing tools, drill bits, cutters), etc. The composite of tungsten carbide and cobalt with WC:Co = 9:1 ratio is used for making metal cutting and drilling tools since 1930s. This hard alloy's manufacturing process is relatively simple, making it possible to easily scale its production: A fine powder of tungsten carbide or other refractory carbide and a fine powder of cobalt or nickel metal binder are mixed, pressed into molds, and sintered at a temperature close to the melting point of the binder metal, producing a very dense and hard material.

However, many transition metal borides and carbides pose challenges. From the experimental point of view, it is extremely difficult to determine the correct composition and crystal structure of transition metal borides and carbides by X-ray diffraction methods. This leads to errors in determining the mechanical properties of synthesized samples of an unclarified compound. Examples include the ongoing controversy around the exact composition and crystal structure of the highest tungsten boride and the superhardness of rhenium diboride, reported to have a hardness of 48 GPa¹⁷ which was later refuted,³² in agreement with our theoretical calculations (see below).

Crystal structure prediction methods develop rapidly^33–35 and can already be used to discover novel superhard materials. To facilitate this, simple and reliable ways of computing hardness and fracture toughness are necessary.

Several developed empirical models,^36–39 which require as an input the crystal structure and chemical composition or the elastic properties, can be used to calculate the hardness of different materials. One of the models proposed by Chen et al.³⁸ is based on the assumption that the indentation size is correlated with the shear modulus of the material (G), while the width of a formed imprint is proportional to the bulk modulus (B),³⁸ 

Analysis of experimental data on many materials³⁸ has led to the following empirical formula for Vickers hardness (H_V):

where k is the Pugh ratio (k = G/B), G is the shear modulus, and B is the bulk modulus; H_v, B, and G are expressed in GPa. Calculations of Vickers hardness for a number of materials using Chen's model agree well with experimental data:³⁸ the calculated value for diamond is 98 (experimental value ∼96^4,40), for TiN—22.6 (20.5²⁹), and for c-BN—56.9 (∼55^1,41). However, this reliable model requires calculations of the elastic constants.

The Lyakhov-Oganov model³⁹ is more convenient for high-throughput searches, numerically stable, usually reliable and can be successfully used for fast Pareto-screening for superhard materials. This model requires only the crystal structure and chemical composition; the hardness is computed from individual bond hardnesses.

It is important to compare the accuracy of both empirical models with the available theoretical and experimental data. The comparison of the results obtained by using these models of hardness with the experimental data that we calculated for various predicted stable and metastable materials is shown in Fig. 2. Here, the largest studied class of materials is transition metal borides (TMBs). The two models are consistent with each other and with the reference experimental data for TMBs. Theoretical hardnesses are more often underestimated than overestimated, compared with the experimental data (Fig. 2). One of the exceptions is WB₄, which may be related to the incorrectly determined structure and/or composition of the synthesized material. For transition metal polynitrides with quasimolecular N_n groups (Cr-N and Mo-N systems, the blue region in Fig. 2), there is a significant difference between these models, showing that Chen's model³⁸ is more accurate. This is because the Lyakhov-Oganov model includes very strong N–N bonds in the calculation of hardness, while only metal-nitrogen bonds are broken during hardness tests. Expectedly, the Lyakhov-Oganov model overestimates the hardness of such polynitrides.

For a very large number of materials, including diamond, c-BN, SiC, borides, hydrides, etc., a much better agreement can be observed between both theoretical models and the experimental data. Among experimental data, an unusual situation can be noted for γ-B28; recent experimental measurement of low Vickers hardness of γ-B₂₈⁵⁷ (∼30 GPa) looks controversial: there are two earlier measurements of the hardness of γ-B₂₈,^28,58 indicating that it is superhard (H_v = 50–58 GPa). It is possible that the low value of hardness for γ-B₂₈ is due to systematic underestimation: in the same set of experiments, Vickers hardness of cubic BN was also underestimated (47 GPa vs ∼60 GPa from Refs. 1, 2, 10, and 11). Similar values (∼50 GPa) of hardness of γ-B₂₈ result from physically motivated empirical models; such models (which can be either analytical or based on machine learning) are justified, and one can see how well this model describes the hardness of very different materials (Fig. 2).

For industrial applications, fracture toughness, along with hardness, also plays a key role. Unfortunately, the hardest materials are usually brittle (covalent crystals), while the materials with the highest fracture toughness are metals whose hardness ranges from low to medium.

Fracture toughness can be calculated using the empirical model from Ref. 59,

where α is the enhancement factor accounting for the degree of metallicity, V is the volume per atom, and G and B are the shear and bulk moduli, respectively.⁵⁹ For insulators, semiconductors, transition metal carbides, nitrides, borides, and hydrides, α = 1.⁵⁹ The calculated values of fracture toughness of diamond, WC, TiN, and c-BN are close to those measured experimentally, e.g., 6.33 MPa m^0.5 (experiment: 4–7 MPa m^0.5)^60–62 for diamond, 5.37 MPa m^0.5 (experiment: 5–8 MPa m^0.5)^30,63 for WC, 3.9 MPa m^0.5 (experiment: 3.4–5 MPa m^0.5)⁶³ for TiN, and 5.41 MPa m^0.5 (experiment: 2–5 MPa m^0.5)^41,62 for c-BN.

Chromium-based materials reveal attractive mechanical properties. In the Cr-C system, only three stable carbides, Cr₂₃C₆, Cr₃C₂, and Cr₇C₃,^64–68 are known from experiments. Two metastable chromium carbides, CrC and Cr₃C, have also been synthesized.^69–72 In the Cr-B system, six different phases, Cr₂B, Cr₅B₃, CrB, Cr₃B₄, CrB₂, and CrB₄, are known from experiments,^{12,50,73–76} and their mechanical properties were examined theoretically.^13,14,50,77 A recent theoretical study of chromium carbides¹⁴ using the evolutionary algorithm USPEX^78–80 has led to the prediction of new stable phase Pmn2-Cr₂C, in addition to already known phases.^65,81 Pmn2-Cr₂C phase is anticipated to have the highest shear modulus of 292 GPa among all chromium carbides and the Vickers hardness of 27 GPa.¹⁴ The highest Vickers hardness among chromium-based materials is achieved in CrB₄ and was reported to be in the range of 29–44 GPa,^14,49 while our calculations give 48 GPa.

Tungsten borides are extremely interesting phases, displaying mechanical properties suitable for industrial applications. One of the first reports on their synthesis was made by Kiessling⁸² in 1947. In 1963, the first phase diagram of the W-B system was drawn.⁸³ Since then, the W-B phase diagram was reinvestigated and clarified several times.^84–87 Five stable tungsten borides were reported in experiments: W₂B,^20,82 WB (α and β phases),^44,82 WB₂,⁸⁸ and the much debated WB₄.^20,21,23 Numerous theoretical studies of the stability and properties of these phases were published recently.^{15,22,25,26,89,90} Wide homogeneity regions of W-B phases, mentioned in theoretical and experimental works,^20,25,91 may at least partly be caused by an extensive polysomatism.⁹² This often leads to difficulties in the synthesis of stoichiometric phases and to inaccurate crystallographic descriptions of synthesized phases by experimental methods, especially in determining the exact positions of light boron atoms using X-ray diffraction.

Because of this, the phase originally claimed in 1947 as W₂B₅⁸² was correctly identified as W₂B₄ with the $P63/mmc$ space group only after 60 years.^91,93 Another example of contradiction is boron-rich phases observed for the first time by Chretien and Helcorsky,⁹⁴ where the composition was initially reported to be WB₄,⁹⁴ while other works presented it as W₂B₉,⁹⁵ WB₁₂,⁹⁶ and W_1−xB₃.⁸⁶ Various experimental techniques were used to determine and characterize the obtained samples.^95,97 This phase with unknown composition and crystal structure has made the W-B system very popular in recent studies, with Gu et al.²¹ suggesting that WB₄ may be superhard. Its hardness was measured experimentally and calculated theoretically. It was determined that its Vickers hardness is higher than 40 GPa^23,98 and can even reach 57.3 ± 1.9 GPa (under the load of 0.49 N) with the addition of chromium.⁹⁹ Later, Cheng et al.¹⁰⁰ showed that this material's composition is WB_3 + x, not WB₄. A large excess of boron can provide good crystallinity and stability of WB_3 + x.²⁰ 

A recent evolutionary computational search²⁶ for new stable superhard tungsten borides shows that in the boron-rich region of the W-B phase diagram, a new Pmmn-WB₅ phase is stable at ambient conditions. WB₄ was predicted to be thermodynamically stable only at pressures above 1 GPa.²⁷ The most intriguing properties of the newly predicted WB₅ are a combination of very high Vickers hardness, ∼45 GPa, and high fracture toughness, 4 MPa m^0.5. In addition, it is thermodynamically stable at ambient pressure and predicted to retain excellent mechanical properties even at 2000 K. According to these predictions, WB5 is a promising material for practical applications.

Modern technology requires new materials displaying unique combinations of mechanical, electronic, and other properties to replace traditional materials widely used in the industry. Recently, we developed a new method¹⁰¹ that allows one to predict optimal materials from among all possible combinations of all the elements across the Periodic Table. This nonempirical method combines our coevolutionary approach with a carefully restructured “Mendelevian” chemical space, energy filtering, and Pareto optimization¹⁰² of target properties and stability, to ensure that the predicted materials have attractive properties and a high chance to be synthesizable. Variation operators (i.e., the “moves” in the structural and compositional space) are of central importance in every evolutionary/coevolutionary algorithm and need to be designed carefully. For this reason, the variation operators of this method work in the space of two important elemental properties (electronegativity and atomic radius) and are designed in a way that newly sampled compounds are based on elements that are similar in their properties to elements present in already sampled compounds with good properties. This definition of chemical space originates from Goldschmidt's law^103,104 of crystal chemistry (which states that crystal structure is determined by stoichiometry, atomic size, electronegativity, and polarizability of atoms/ions) and was inspired by the Mendeleev number of Pettifor^105,106 (for more details on our new method, Mendelevian search, or MendS, see Ref. 101). This way, step by step, the algorithm learns about promising regions of the chemical space and focuses on them at the expense of unpromising regions. Using our MendS method, we carefully studied the chemical space constructed from 74 elements (all the elements, excluding the noble gases, rare-earth elements, and elements heavier than Pu) and searched for materials with optimal hardness (calculated by the Lyakhov-Oganov empirical model³⁹) and stability (computed using Maxwell's convex hull construction). As Chen's model³⁸ is more accurate, it was used to obtain the final theoretical hardness, which was then used for constructing the Ashby plot of the Vickers hardness vs fracture toughness (Fig. 3). This results in a map on which materials with an optimal combination of the Vickers hardness and fracture toughness can be easily found.

The multiobjective Pareto technique¹⁰² and energy filtering (discarding the structures with energy above the convex hull by more than 0.5 eV/atom) were used to ensure that the evolutionary algorithm generates hard phases that have low energy. In this calculation, 600 binary systems were studied in 20 MendS generations (which is only about one-fifth of all the systems that can be made of 74 elements). Impressively, all the hard binary materials reported in the literature and materials claimed to be potentially hard were found in our calculation. Diamond (and its polytype lonsdaleite) as the hardest, and boron allotropes as the second hardest elemental phases, were found correctly in our Mendelevian search. B_xC_y,^109–111 B_xO_y,^108–112 C_xN_y,^113,114 B_xN_y,^1,2,10,11 Cr_xC_y,^53,54,115 Cr_xB_y,^{13,48–50,77} Cr_xN_y,^53,54,115 Re_xB_y,^15–17,32 W_xB_y,^20–27 Si_xC_y,^{36,37,116,117} W_xC_y,^36,37,118 Al_xO_y,^36,117 Ti_xC_y,¹¹⁹ Si_xN_y,¹¹⁸ Ti_xN_y,¹¹⁸ Be_xO_y,¹¹⁸ Ru_xO_y,^119–121 Os_xO_y,¹²² Rh_xB_y,¹²³ Ir_xB_y,¹²³ Os_xB_y,^124–126 and Ru_xB_y^124–126 are some of the examples of binary systems reported to be hard in the literature and found by us in this single calculation. Moreover, our calculation revealed several binary systems with very promising hardness: Mo_xB_y,^18,19 B_xP_y,¹²⁷ V_xB_y,^128–130 Fe_xB_y,^131,132 Mn_xB_y,^133,134 and Mn_xH_y.

In the W-B system, the most promising material is WB₅, predicted²⁶ to be thermodynamically stable and to have high hardness and fracture toughness. The Mn-H system appeared in our list unexpectedly; it was never thought to contain hard phases. However, we showed that many phases in this system are stable or low-energy metastable, have high-symmetry structures, and hard indeed, having the Vickers hardness in the range of 20–30 GPa (Fig. 3). Figure 3, the Ashby plot of Vickers hardness and fracture toughness, is a “treasure map” of superhard materials, summarizing much of what we know and expect today and pointing to future directions of research. Concerning low-pressure phases, the best combination of hardness and fracture toughness is achieved in TiB₂, ZrB₂, VB, V₃B₄, VB₂, VB₁₂, CrB₄, WB₅, MnB₄, ReB₂, and WC.

We have discussed here several important examples of transition metal boride systems that were computationally predicted to be promising for practical applications. The new method, Mendelevian search, shows great predictive power in the search for hard and superhard materials. Modern computational techniques allow the prediction and relatively fast development of new materials (e.g., superhard materials) with enhanced properties, destined to replace traditional materials.

A.R.O. was supported by the Russian Science Foundation (Grant No. 19-72-30043). A.G.K. thanks the Russian Science Foundation (Grant No. 17-73-20038). This work was carried out using our Rurik supercomputer at MIPT and the Arkuda supercomputer of the Skolkovo Foundation.