Discovering Superhard B-N-O Compounds by Iterative Machine Learning and Evolutionary Structure Predictions.

We searched for new superhard B-N-O compounds with an iterative machine learning (ML) procedure, where ML models are trained using sample crystal structures from an evolutionary algorithm. We first used cohesive energy to evaluate the thermodynamic stability of varying B x N y O z compositions and then gradually focused on compositional regions with high cohesive energy and high hardness. The results converged quickly after a few iterations. Our resulting ML models show that B x+2N x O3 compounds with x ≥ 3 (like B5N3O3, B6N4O3, etc.) are potentially superhard and thermodynamically favorable. Our meta-GGA density functional theory calculations indicate that these materials are also wide bandgap (≥4.4 eV) insulators, with the valence band maximum related to the p-orbitals of nitrogen atoms near vacant sites. This study demonstrates that an iterative method combining ML and ab initio simulations provides a powerful tool for discovering novel materials.

Introduction

With the continuous increase in the global demand for superhard materials, searching for new compounds with outstanding hardness and stability is becoming an important research topic.[1] Materials of superhardness (with hardness value H ≥ 40 GPa) can be classified into two major groups. The first includes transition metal (TM) ceramics, especially the borides and carbides.[2−4] The second concerns light elements B, C, N, and O,[5−8] which form short and strong covalent bonds. Various superhard materials can be produced by mixing two or three of these light elements. Materials in the second group have the advantages of being abundant and low cost, but they may require special synthesis conditions. Solozhenko et al. synthesized diamond-like BC5 under high pressure and high temperature (HPHT) and reported a corresponding hardness H = 71 GPa.[9] They also found that BC5 has a relatively high fracture toughness and exceptional thermal stability up to 1900 K. Baker et al. used microwave plasma chemical vapor deposition to synthesize boron-incorporated diamond at low temperature and low pressure. They found a hardness as high as 62 GPa in their cubic phase sample with 7.7 at% boron content.[10]

Superhard B–C–N compounds also have been studied in the literature. Cubic BC2N was reported with a superhardness between 62 and 76 GPa,[11−13] which is harder than cubic boron nitride (c-BN, H ∼ 50–70 GPa).[14,15] Experimentally, BC2N remains stable up to 1800 K, which demonstrates its superior thermal stability compared to diamond (H ∼ 100 GPa). Other ternary compounds such as BCN and BC9N in the cubic phase also have been synthesized under extreme conditions.[16] Recently, Chen et al. have used machine learning and evolutionary searches to discover superhard B–C–N structures. Their newly predicted BC10N has an ultrahigh hardness of ∼87 GPa with a relatively low formation energy.[17]

B–N–O compounds also have attracted considerable attention due to their relevance in several research areas, such as chemical adsorption,[18] fluorescent dots,[19] water splitting,[20] wide bandgap insulators,[21] and electrochemical applications.[22] Although B–N–O compounds have been synthesized in various forms, including amorphous, porous, layered, thin-film, and nanoparticle, their mechanical properties remain largely unexplored. More recently, Bhat et al. have used hexagonal BN and B2O3 as initial materials to synthesize B6N4O3 under HPHT conditions.[23] They also conducted computational studies and indicated that the most stable structures contain ordered vacancies in a zinc-blende structure. Such structures have a high bulk modulus of 300 GPa, which implies that B–N–O compounds may be superhard as well.[23]

Motivated by the study of Bhat et al., here we employ first-principles calculation and machine learning (ML) simulation to search for new superhard B–N–O compounds. Data-driven approaches have proven to be powerful in materials discovery,[32−34] and several ML models have been applied to find superhard compounds.[17,35−38] The first important step in building ML models is sample data acquisition. However, although there are several online computational materials databases,[39−43] only limited information exists on B–N–O. The lack of relevant training data may lead to inaccurate prediction of universal ML models, when they are applied to the not-well-explored B–N–O compounds. For example, our previous universal ML model[17] would predict a bulk modulus of ∼180 GPa for B6N4O3, which largely deviates from the reported value of ∼300 GPa from first-principles calculation.[23] Therefore, we need a different ML scheme for predicting B–N–O systems.

In this paper, we develop an iterative procedure involving crystal structure prediction (CSP), density functional theory (DFT) calculation, and machine learning (ML) simulation, in order to discover new ternary superhard B–N–O materials. The details of CSP, DFT, and ML simulations are discussed in the Computational Methods section. Below, we will first address our iterative calculation process, which involves four major steps as summarized in Figure : (I) Generating crystal structures, (II) computing physical properties from first-principles, (III) building ML models, and (IV) predicting properties using the ML models.

Figure 1

Schematic iterative process of machine learning and structure prediction for superhard B–N–O compounds. I. Generating B–N–O crystal structures from evolutionary algorithm implemented in USPEX.[24−26] II. Computing physical properties from first-principles density functional theory software VASP.[27,28] The target properties include cohesive energy, density, and mechanical properties. III. Building machine learning (ML) models with SCIKIT-LEARN.[29] The ML features from Meredig et al.[30] are generated by MATMINER.[31] IV. Predicting B–N–O physical properties with random forests ML models. Promising compositions with high stability (high cohesive energy) and high hardness are selected for the next round of iterative calculation.

In step (I), we use an evolutionary algorithm for CSP to create new B–N–O structures with varying chemical compositions. In step (II), we use DFT to compute the cohesive energy, volumetric density, and elastic tensor, for crystal structures found in the previous step. The elastic tensor enables the evaluation of a structure’s mechanical properties and elastic stability. The physical properties of stable structures from step (II) provide the sample data for training ML models. In step (III), we build separate ML models for the cohesive energy, density, and hardness. The cohesive energy helps to determine a structure’s thermodynamic stability. Finally, in step (IV), we construct triangular plots to map out the properties listed above for arbitrary B–N–O composition, using ML models from the previous step. Based on the ML predictions, we then select promising chemical compositions with high thermodynamic stability and high hardness as candidates for the next iteration of calculation. We find that our ML models can be gradually improved with this procedure. We stop at the fourth iteration, as the results are quickly converged. As shown below, using this iterative procedure to combine CSP/DFT calculations and ML simulation, we are able to discover several new stable and superhard B–N–O structures.

Computational Methods

Crystal Structure Prediction (CSP)

We use CSP to generate B–N–O sample structures, which later serve as the input data for training machine learning (ML) models. The purpose of CSP is to find stable and/or metastable structures of a compound given only its chemical formula.[44−46] Here, we use the highly powerful and efficient evolutionary algorithm implemented in the USPEX package.[24−26] During different generations of evolutionary optimization, new structures are created by heredity (50%), mutation (30%), and random (20%) operators. The enthalpy computed by density functional theory (DFT) is used as the fitness.

We consider varying compositions of ternary B–N–O compounds, with a unit-cell size between 9 and 19 atoms. We ensure that all of the selected atomic compositions contain an even number of total electrons. This typically leads to a more stable, insulating phase. Since our goal is to discover superhard compounds, we apply an external pressure P = 15 GPa during the structure searches. Adding a small but finite pressure will help avoid low-hardness layered graphite-like structures. For each chemical composition, we perform two separate USPEX calculations and search over 1200 structures. The optimized structures from both USPEX searches are further fully relaxed in DFT calculations without any external pressure. All DFT results and figures shown below are from calculations performed without any external pressure. In the end, for a given chemical formula, we have two sample structures for training ML models.

We note that USPEX also provides a variable-composition mode, which can sweep the whole B–N–O compositional space (constrained by the number of atoms in the unit cell) to find the most promising stable compositions, which could then be explored in more detail by separate fixed-composition calculations. It is expected that this method will lead to similar results as our iterative procedure.

Density Functional Theory (DFT) Calculation

We perform DFT calculations using the VASP software,[27,28] which adopts a pseudopotential method and plane-wave basis sets. We use the projector augmented wave (PAW)[47,48] method and generalized gradient approximation (GGA) functional based on the Perdew–Burke–Ernzerhof (PBE) formalism.[49] The kinetic energy cutoff for wave function expansion is 520 eV, and the k-points are sampled by a Γ-centered Monkhorst–Pack mesh with a resolution ∼0.02 × 2π/Å. The convergence criteria of electronic self-consistency and structural relaxation are set to 10–6 eV/unit cell and 10–3 eV/Å, respectively.

For fully relaxed crystal structures, we further utilize the strain–stress method[50] implemented in VASP to compute the elastic constants C, which in turn help to evaluate the elastic stability and mechanical properties. In particular, the eigenvalues of the elastic tensor (for elastic stability), as well as the Voigt-Reuss-Hill (VRH) averaged[51−53] bulk modulus K and shear modulus G, are computed via the MechElastic python library.[54,55] The Vickers hardness H is also evaluated by using Chen’s empirical hardness model:[56]where k = K/G is the Pugh’s ratio. Other empirical hardness models[57,58] based on elastic moduli in general produce similar hardness predictions. Finally, phonon density of states (for dynamical stability) are computed using the PHONOPY package,[59] via the density functional perturbation theory scheme implemented in VASP. Convergence tests of all calculations are examined carefully.

Machine Learning (ML) Prediction

Since the ML predictive capability depends crucially on the training data, we first exclude elastically unstable sample structures, by examining whether the elastic tensor is positive-definite.[60] After preparing data from CSP and DFT calculations, we build three ML models, respectively, for predicting cohesive energy, density, and hardness.

We first use the python library PYMATGEN[61] to create structure objects using VASP’s POSCARs as inputs. The structure objects are in turn used to create the ML features or descriptors using the python library MATMINER.[31] For cohesive energy and density, we adopt the compositional features of Meredig et al.[30] A total number of 16 features derived from a given chemical composition are considered: atomic fractions of B, N, and O; mean atomic weight; mean column number; mean row number; range of atomic number; mean atomic number; range of atomic radius; mean atomic radius; range of electronegativity; mean electronegativity; average numbers of s and p electrons; fractions of s and p electrons. For the hardness model, we further consider the volumetric density as an additional feature, in order to improve the model performance.

Our ML models are based on random forests, implemented in the SCIKIT-LEARN library.[29] Compared to a single decision tree, ensemble trees can improve model prediction accuracy and meanwhile avoid high variance. When training the ML models, we use 100 estimators and 10-fold cross validation on 80% of the samples (the training-validation set). The maximum tree depth is restricted to 6 layers to further avoid overfitting.

Results and Discussion

Figure shows triangular plots for B–N–O compositions considered in different iterations of our evolutionary structure prediction (top panels) and the cohesive energy predicted by our machine learning (ML) random forests models (bottom panels). In the first iteration, we uniformly sample the B–N–O compositional space, with 12-atom unit cells for binary and ternary compounds. Corner points of the graphs correspond to elemental compounds. For example, pure B here represents α-B12.

Figure 2

[Top panels] Distribution of B–N–O compositions in evolutionary structure prediction and their elastic stabilities computed using density functional theory. Only elastically stable structures (green diamonds) are used for constructing machine learning models; elastically unstable structures (red squares) are excluded during data selection. [Bottom panels] Random forests prediction of cohesive energy. Based on the predicted cohesive energy and hardness (not shown), promising B–N–O compounds are selected for calculation in the next iteration. Ternary graphs are visualized by the PYTHON-TERNARY[62] library.

In the first iteration of data sampling, the ratio between elastically stable and unstable structures is roughly 1:1. Even though the sampling grid is coarse, it can be seen that ternary compounds composed of mostly nitrogen and oxygen tend to be elastically unstable. In contrast, the systems with high boron content are elastically more stable. To ensure high prediction accuracy, we only consider elastically stable structures to build our ML models. In the four iterations, the numbers of elastically stable structures involved in building ML models are 27, 78, 104, and 153, respectively.

The ML prediction of cohesive energy in Figure also indicates that compounds with more boron atoms show higher cohesive energy (i.e., higher thermodynamic stability). The behaviors of elastic stability and thermodynamic stability are thereby consistent with each other. Therefore, starting from the second iteration, we do not consider compositions with simultaneous high nitrogen and high oxygen contents. Instead, we manually select B–N–O compositions based on ML prediction from the previous iteration, by gradually zooming in compositional space with both high predicted cohesive energy and high hardness. With increasing number of iterations, the ML results then have higher and higher accuracy and resolution. Importantly, we find that a region connecting BN and B2O3 (forming BNO3) is thermodynamically favorable with high cohesive energy. This result is consistent with the recently reported B6N4O3 compound by Bhat et al., where a mixture of hexagonal BN and B2O3 is used as the starting material for HPHT synthesis.[23]

Figure shows the predictions of our random forests ML models, respectively, for cohesive energy, density, and hardness, built from training samples in the fourth iteration of evolutionary structure searches. Here, 80% of the samples are used as the training-validation set to construct random forests models, and the remaining 20% are used as the test set for a final unbiased evaluation of model performance. Our ML model for the cohesive energy shows high prediction accuracy, with a Pearson correlation coefficient r = 0.986 between the ML values and actual DFT calculations. The model performances for predicting density and hardness are also comparable, with r scores of 0.953 and 0.975, respectively.

Figure 3

Prediction [top panels] and evaluation [bottom panels] of random forests machine learning (ML) models for (a) cohesive energy, (b) density, and (c) hardness. The Pearson correlation coefficient (r) between the ML predicted value and the density functional theory (DFT) calculation is utilized as the evaluation metric. The ternary graphs indicate that (i) BNO3 [from linear combinations of (BN) and B2O3] has higher cohesive energy, and (ii) hardness is strongly correlated with density.

In general, a material’s hardness is correlated with its volumetric density.[68] For example, cubic diamond as the hardest material also has the largest volumetric density at ambient condition. Therefore, to enhance the ML prediction accuracy, we also consider density as an additional feature when building the hardness model. We first use only compositional features to train a density ML model, which in turn is utilized to generate the density feature for the subsequent hardness ML model. We note that creating another model to predict missing or lacking feature data is a common practice in machine learning. In principle, error propagation could be a concern for such a workflow. In practice, it is not a concern here because of the high accuracy of our ML models. As shown in Figure b,c, the ML predicted density and hardness values indeed show a nearly linear relationship.

Our iterative calculation combining an evolutionary structure search and ML prediction shows that B–N–O compounds can exhibit superhardness in a compositional region around BN. With the aid of the ternary graphs shown in Figure , we can directly inspect the sample crystal structures generated from the evolutionary algorithm. Table lists various superhard B–N–O compounds discovered in our study. In particular, we find mainly two kinds of ternary compositions. In the first kind, the number of boron is equal to the total number of nitrogen and oxygen (e.g., B5N3O2). The second kind contains one more oxygen, where the number of boron is equal to the total number of nitrogen and oxygen minus 1 (e.g., B5N3O3). In principle, the two different kinds of B–N–O compositions can have different stable crystal structures.

Table 1

Physical Properties of Superhard B–N–O Compounds with Cohesive Energy >6.75 eV Discovered in This Studya

crystal	ρ	K	G	E	ν	AU	KIC	H	EgPBE/EgmBJ	Ecoh	Eform
B5N3O2	0.151	286	257	593	0.155	0.294	4.38	42	4.3/5.5	6.761
B5N3O3	0.150	279	255	586	0.150	0.088	4.09	43	6.3/7.6	6.935	83
B6N4O2	0.156	322	284	658	0.160	0.121	5.26	44	3.0/4.4	6.810
B6N4O3	0.155	292	268	616	0.149	0.311	4.33	45	4.5/5.7	6.900	120
B7N5O2	0.146	272	241	559	0.157	0.336	4.09	40	3.6/4.6	6.776
B7N5O3	0.158	306	283	649	0.147	0.276	4.54	47	4.1/5.3	6.904	117
B9N7O2	0.161	330	296	683	0.154	0.139	5.31	46	3.2/4.5	6.802
B9N7O3	0.160	318	301	688	0.140	0.245	4.42	50	3.9/5.1	6.926	97
c-BN	0.168	373	383	856	0.118	0.172	4.73	64	4.5/5.3	7.028	0
B2O3	0.101	35	33	75	0.127	2.347	0.14	12	6.3/8.9	7.008	0

Density ρ (atom/Å3), bulk modulus K (GPa), shear modulus G (GPa), Young’s modulus E (GPa), Poisson’s ratio ν, universal elastic anisotropy A, fracture toughness KIC (MPa·m1/2), hardness H (GPa), bandgap Eg (eV), cohesive energy Ecoh (eV/atom), and formation energy Eform (meV/atom). The fracture toughness is based on the empirical model by Mazhnik and Oganov.[58] The bandgaps are computed respectively with the standard Perdew–Burke–Ernzerhof (PBE)[49] functional and the Tran-Blaha modified Becke-Johnson (TB-mBJ)[63,64] exchange potential for improved bandgap estimation. For benchmark, the experimental hardness Hexp = 50–70 GPa for c-BN,[14] and Hexp = 1.5 GPa for B2O3.[65] The experimental bandgap Egexp=6.36 eV for c-BN,[66] and Egexp > 10 eV for B2O3.[67]

Using DFT calculations, we have studied various physical properties of B–N–O compounds listed in Table . In particular, their volumetric densities ρ (atom/Å3) are overall positively correlated with the hardness values, computed by Chen’s empirical hardness model based on the bulk modulus K and shear modulus G.[56] Moreover, by examining the electronic density of states (DOS), we find these new superhard B–N–O materials are insulators, with wide bandgaps ≥3.0 eV based on the PBE functional. It is known that GGA functionals like PBE will tend to underestimate the bandgaps. For a more accurate estimation, we further consider the Tran-Blaha modified Becke Johnson (TB-mBJ)[69,70] meta-GGA exchange potential. Our TB-mBJ results show that the superhard B–N–O compounds under study are wide-bandgap insulators with a gap size ≥4.4 eV.

We next evaluate the thermodynamic stabilities of the newly predicted superhard B–N–O compounds, by comparing their cohesive energy:With the total energies of isolated atoms as references, Ecoh is the energy gained by arranging atoms in a crystalline state. Therefore, the higher the cohesive energy, the higher the thermodynamic stability of a material. The second kind of compound BNO3 is found to exhibit a higher cohesive energy than the first kind. Moreover, the cohesive energies of c-BN and B2O3 are even higher than the predicted B–N–O compounds, which are thereby metastable at ambient condition, but may be stabilized under HPHT synthesis conditions.[23]

For BNO3, we also calculate the formation energy using the total energies of c-BN and B2O3 as references. While hexagonal boron nitride (h-BN) and B2O3 were utilized as the starting materials to synthesize B6N4O3,[23] c-BN and B2O3 would be the relevant energy minima of the convex hull, especially under HPHT synthesis conditions. In particular, the formation energy for these B–N–O compounds can be computed byBased on this formula, B5N3O3 exhibits the lowest formation energy among the B–N–O compounds under study. A corresponding superlattice structure computed by an evolutionary algorithm for B5N3O3 is shown in Figure a. The result is consistent with the report from Bhat et al. that models of ordered structure would agree better with their experiments.[23] With the previously reported success in the synthesis of B6N4O3,[23] it is likely that other metastable B–N–O structures from our predictions could be synthesized as well.

Figure 4

Predicted crystal structures [top panels], electronic density of states (DOS) [middle panels], and phonon DOS [bottom panels] for (a) B5N3O3, (b) B6N4O3, (c) B7N5O3, and (d) B9N7O3. The results are obtained using the Perdew–Burke–Ernzerhof (PBE) functional. All four compounds are wide-bandgap insulators, with a peak at the valence band maximum originating from p-orbitals of nitrogen atoms near the vacant sites (see Figure ). The phonon spectra show only positive modes, indicating dynamical stability of all four compounds. The crystal structures are visualized by the VESTA software.[71]

Figure 5

B5N3O3 crystal structure: (a) side view and (b) top view. The unit cell is indicated by thin solid lines. The vacant sites around three oxygen atoms and one nitrogen atom are emphasized by dotted circles. The electron charge contour (cyan color) corresponds to the valence band maximum in the electronic density of states, and it is mainly contributed by nitrogen p-orbitals near the vacant sites.

Based on Table , we select a few more stable BNO3 compounds to further characterize their properties. Figure shows the crystal structures, electronic density of states (DOS), and phonon DOS, respectively, for B5N3O3, B6N4O3, B7N5O3, and B9N7O3. Among these compounds, only B5N3O3 assumes a wurtzite structure, and all the others have a zinc-blende structure. Compared to BN, BNO3 contains one less boron atom (due to B2O3), leading to vacant boron sites as indicated by the dotted circles in Figure .

As seen in the middle panels of Figure , the three zinc-blende structures have similar electronic profiles around the Fermi level, with a bandgap ≳4 eV. The electronic DOS of the wurtzite B5N3O3 is also similar, but it exhibits an even larger bandgap ≥6 eV. This result is unusual, because in binary BN structures, wutzite BN has a smaller bandgap than that of cubic BN, which is distinct from the ternary B–N–O compounds studied here. We also find that in the zinc-blende structures, the bandgap increases with increasing oxygen content. This result is reasonable, as B2O3 has an extremely large band gap ≥10 eV,[67] compared to that of 6 eV in c-BN. Overall, introducing oxygen atoms into binary BN compounds can substantially alter the electronic structures, which also could cause changes in thermal and chemical stabilities. Finally, the phonon DOS (bottom panels of Figure ) shows only positive phonon modes, ensuring that these B–N–O compounds are dynamically stable.

Finally, we note that the electronic DOS for all compounds in Figure shows a sharp peak at the valence band maximum. By studying the local electronic DOS, we find the sharp peak is mainly contributed by localized p-orbitals of nitrogen atoms next to the vacant boron sites. The electron charge distribution corresponding to the sharp peak in B5N3O3 is shown in cyan color in Figure .

Conclusion

We have developed an iterative machine learning procedure to discover new superhard B–N–O compounds, where the input training samples are generated from evolutionary structure searches and density functional theory calculations. Our combined machine learning and first-principles results revealed several stable and superhard B–N–O compounds of chemical compositions BNO3 (x ≥ 3), with hardness values ≳45 GPa. We also found that these newly predicted B–N–O systems are all wide bandgap insulators, with gap size ≥4 eV based on the Perdew–Burke–Ernzerhof GGA functional. Our additional meta-GGA calculations indicate that their actual bandgaps could be even larger. The electronic density of states for BNO3 all show a prominent peak around the valence band maximum, and it is related to localized p-orbitals of nitrogen atoms near vacant boron sites. Since B6N4O3 has already been reported in the literature, other BNO3 compounds in principle can be synthesized by high-pressure high-temperature techniques or by chemical vapor deposition methods. We expect these newly discovered superhard B–N–O materials to have a wide range of applications in extreme environments, and they may outperform diamond or cubic boron nitride in an oxidizing environment or in humid conditions at high temperature. In particular, superhard wide bandgap insulators can be more resilient to radiation and have a higher breakdown voltage, so they can be practical for reactor and high-voltage electronics. Higher thermal stability also enables industrial applications in cutting and machining of hard ferrous metals. Because of their superior mechanical properties, these materials may be superhard thermal conductors for heat management and transfer as well.