A deep neural network potential model for transition metal diborides

Dataset collection

In this study, we adopt the machine learning potential model proposed by Zhang et al., which is known as the Deep Potential (DP) model^[17]. The DP model represents the potential energy surface using deep neural networks and consists of two sets of neural networks. The first set maps local atomic environments to symmetry-preserving descriptors, while the second set maps these descriptors to atomic energies^[17]. We use the revised version of the DP model, called the DPA-1 model^[28], which incorporates an element embedding mapping to improve generalizability and suitability for developing medium-scale interatomic potentials.

Training a machine learning potential requires a high-accuracy DFT dataset, containing atomic configurations and their corresponding energies, forces, and virial stresses. We collect the DFT dataset using the DP GENerator (DP-GEN) software, which generates the dataset under a concurrent learning scheme^[29]. This scheme enriches the candidate DFT dataset iteratively in three steps: training, exploration, and labeling [Figure 1A]. Prior to the DP-GEN procedure, an initial dataset must be created. For this purpose, we produced 300 structures featuring random compositions and performed single-point first-principles calculations to generate the initial dataset. The approach for generating structures is consistent with the methodology employed in the second step of DP-GEN. The detailed DP-GEN processes are:

Figure 1. (A) Illustration of the DFT dataset generation process using a concurrent learning scheme in the DP-GEN software. The lower part of the image demonstrates the coverage of the configurational space by the current model, with accurate, candidate, and inaccurate regions distinguished by different colors; (B) Comparison of energy predictions between the DP model and DFT calculations; (C) Error distribution in predicted energy; (D) Comparison of force predictions between the DP model and DFT calculations; (E) Error distribution in predicted force. DFT: Density functional theory; DP-GEN: Deep Potential GENerator; DP: Deep Potential.

(1) During the training step, we train four DP models with different initialization parameters and activation functions based on existing data.

(2) In the exploration step, we explore the configurational space through MD simulations with DP models from the training step. We generate 100 diboride supercell structures with random compositions at the beginning of each iteration. A supercell with 72 atoms is created for each structure. The axes of the supercell are oriented along [420], [030], and [002] directions of TMB₂ (The crystal structure of TMB₂ is described in Supplementary Text 1 and illustrated in Supplementary Figure 1). Then, the metallic elements of the supercell are randomly occupied by the candidate elements. For example, the element is randomly generated from Ti, Zr, Hf, V, Nb, Ta, Cr, Mo, and W for each metallic site. Both bulk and surface structures are explored. For surface structures, the surface plane is randomly chosen from the x (normal to [420] direction), y (normal to [030] direction), and z (normal to [002] direction) planes of the supercell. For bulk structure sampling, we adopt the isothermal-isobaric (NPT) ensemble, with pressures ranging from 0 to 20 GPa (with a step of 5 GPa) and temperatures ranging from 300 to 6,000 K (with a step of 300 K). For surface structure sampling, we adopt the isothermal-isovolumic (NVT) ensemble, with temperatures ranging from 300 to 2,100 K.

(3) During the labeling step, we select candidate structures and label them using a DFT-based method. We randomly select the candidate structures from the MD trajectories based on “model deviation” (ε), which is the maximum standard deviation of forces predicted by the four DP models. We set two criteria for “model deviation”: ε_low and ε_high. When ε < ε_low, all models predict similar values for the configuration, and adding the configuration to the dataset may only marginally improve the model. When ε > ε_high, the models have strong conflicts on the configuration, and there is a risk that the configuration is non-physical due to the inaccuracy of the current DP model. Adding the configuration to the dataset may disrupt the model. Therefore, we choose configurations with “model deviation” between ε_low and ε_high as candidates. The ε_high is set to be 2 eV/Å. The ε_low changes based on the exploration temperature (T), which are 0.3 eV/Å (T < 2,000 K), 0.5 eV/Å (2,000 K < T < 4,000 K), and 0.6 eV/Å (4,000 K < T < 6,000 K). We then label these candidate configurations using the Vienna ab initio simulation package (VASP)^[30], adding them to the dataset. The calculations adopt the projector augmented wave (PAW)^[31] approach to describe ion-electron interactions, and the exchange-correction is described using the Perdew-Burke-Ernzerhof (PBE) generalized gradient approximation (GGA)^[32]. We set the cutoff energy of the plane wave basis to 900 eV and adopt the Monkhorst-Pack^[33]k-point mesh with a separation of 0.15 Å^-1 in the Brillouin zone. Self-consistent field iteration stops when the difference in total energy of consecutive iterations is less than 10^-6 eV.

Model training

After collecting the dataset, the DP model was trained using the DeePMD-kit software^[34,35], specifically employing the DPA-1 model. The architecture of the DP model was established as follows:

(1) The descriptor net maps atomic environments to symmetry-preserving descriptors, utilizing a network type set as “se_atten_v2”. The “se_atten_v2” type in DeePMD-kit software designation refers to the DPA-1 model architecture, as detailed in the reference^[28]. The element embedding net consists of a single-layer neural network with eight nodes. The descriptor net is composed of three layers of neural networks, each containing 25, 50, and 100 nodes, respectively. The projection dimension is set to 16, and the cutoff is set to 7.0 Å with a smooth function imposed from 6.0 Å. The activation functions used are tanh.

(2) The fitting net maps the descriptors to atomic energies, including three layers of neural networks with each having 128 nodes. A ResNet (residual network) architecture is adopted in the fitting net, and the activation functions are also tanh.

The model was trained in two steps with the “Adam” optimizer. During the first step, the learning rate, which scales the loss that feeds to the optimizer, decayed from 1.0 × 10^-3 to 1.0 × 10^-8, and the pre-factors of energy, force, and virial in the loss function changed from 0.1 to 1, from 100 to 1, and from 0.1 to 1, respectively. Subsequently, the model was retrained, inheriting the model parameters from the first training round. In this round, the learning rate decayed from 1.0 × 10^-4 to 1.0 × 10^-8, and the pre-factors of energy, force, and virial in the loss function were set to 10, 1, and 1, respectively, remaining constant. In each training round, the total training step amounted to 8 million. For further details regarding the parameters, one may refer to ref^[34,35] and the open-source software. The DeePMD method has been successfully integrated into the open-source MD simulation software LAMMPS^[36]. Subsequent simulations were performed using LAMMPS with the trained DP model.

Prediction on fundamental properties

In addition to the comparable predictions of energy and force between the DP model and DFT calculations, it is crucial to assess the reliability of the DP model in predicting material properties. Table 1 presents a comparison of the lattice parameters (a and c), cohesive energy (E), and elastic constants (C₁₁, C₁₂, C₁₃, C₃₃, and C₄₄) of TMB₂s (TM = Ti, Zr, Hf, V, Nb, Ta, Cr, Mo, W) as predicted by the DP model and the DFT calculations, where a good agreement can be observed. Although the model demonstrates high accuracy, it inadequately fits certain properties, such as WB₂’s C₃₃ depicted in Table 1. Fitting properties of unstable materials typically poses a greater challenge. The stability of TMB₂s was assessed by examining their phonon dispersion curves [Supplementary Figure 2], where imaginary modes were detected for CrB₂, MoB₂, and WB₂. The equation of state characterizes the variation of energy with respect to volume changes, which can be determined by relaxing the structures at different volumes. Figure 2 displays a comparison of the equations of state for TMB₂ (TM = Ti, Zr, Hf, V, Nb, Ta, Cr, Mo, W) as predicted by the DP model and the DFT calculations, all of which exhibit good agreements.

Figure 2. Comparison of equations of state of TMB₂ (TM = Ti, Zr, Hf, V, Nb, Ta, Cr, Mo, W) predicted by the DP model (dash line) with the DFT calculations (points). TMB₂: Transition metal diboride; DP: Deep Potential; DFT: density functional theory.

Melting point prediction

The family of TMB₂s is well-known for their high melting points, with some exhibiting melting points higher than 3,000 °C and being classified as UHTCs^[2]. Solid solutions are crucial for tailoring the properties of TMB₂s, such as removing harmful impurities, enhancing sintering, and improving high-temperature strength^[4,6,37]. In recent years, the high-entropy concept was introduced to ceramics, leading to the synthesis of new high-entropy diborides^[11,13]. HE-TMB₂s are believed to possess high melting points, high hardness, and other properties, making them promising candidate materials for extreme environments. However, the melting points of these high-entropy diborides are unknown, and the effects of solid solutions on their melting points are unclear.

In this section, the DP model is used to predict the melting points of diborides and to investigate the effects of solid solutions on melting points. TMB₂s (TM = Ti, Zr, Hf, V, Nb, Ta, Cr, Mo, W) were simulated, along with high-entropy compounds Zr_0.25Hf_0.25Nb_0.25Ta_0.25B₂ (HEMB₂_4), Ti_0.2Zr_0.2Hf_0.2Nb_0.2Ta_0.2B₂ (HEMB₂_Ti), V_0.2Zr_0.2Hf_0.2Nb_0.2Ta_0.2B₂ (HEMB₂_V), and W_0.2Zr_0.2Hf_0.2Nb_0.2Ta_0.2B₂ (HEMB₂_W). The solid-liquid coexistence method was employed to simulate the melting points of these compounds. The simulation box was aligned along the [010], [001], and [210] directions, with lengths of 8 × 8 × 10 (X × Y × Z) along each direction. The solid-liquid interfaces were parallel to the XY plane. The solid, solid-liquid coexistence, and liquid states are illustrated in Figure 3A.

Figure 3. (A) The snapshots of solid, solid-liquid coexistence, and liquid states of HEMB₂_Ti during MD simulations; (B) Illustration of how the melting point value is determined based on V~T data; (C) Comparison of experimentally measured and calculated melting points of TMB₂s, as well as the calculated melting points of HE-TMB₂s, and the estimated values obtained by averaging the melting points of their constituent compounds. The composition of HEMB₂_4, HEMB₂_Ti, HEMB₂_V and HEMB₂_W are Zr_0.25Hf_0.25Nb_0.25Ta_0.25B₂, Ti_0.2Zr_0.2Hf_0.2Nb_0.2Ta_0.2B₂, V_0.2Zr_0.2Hf_0.2Nb_0.2Ta_0.2B₂, and W_0.2Zr_0.2Hf_0.2Nb_0.2Ta_0.2B₂, respectively. Melting point values are listed in Supplementary Table 1. MD: Molecular dynamics; HE-TMB₂s: high-entropy transition metal diborides.

During MD simulations, the timestep was set to 1 fs, and the NPT ensemble was adopted with P = 0 Pa and pressure control applied independently in each direction. The damping parameters were set to 500 fs for pressure and 100 fs for temperature. The system was first equilibrated at a solid temperature (T_solid) close to the melting point for 50 ps. Then, half of the system was heated to a high temperature (T_liquid) and melted, while the other half remained frozen for 50 ps. This process resulted in a solid-liquid coexistence state. Finally, the entire system was set to another temperature (T) and equilibrated for 1 ns. T_solid and T_liquid were fixed for each compound, and the system with a series of T values was analyzed. For different Ts, the final state could be solid, solid-liquid coexistent, or liquid [Figure 3A]. The volume (V) vs. T data was collected from the simulations and fitted with a sigmoid-like function:

where a, b, c, T₀, and V₀ are fitting parameters. The linear expansion effect is accounted for by cT + V₀, while the sigmoid function $$ \frac{b}{1+\exp \left(-a\left(T-T_{0}\right)\right)} $$ represents the volume jump during melting, and T₀ is the melting point. Figure 3B shows the V~T data of HEMB₂_Ti and how well the equation fits the data. The midpoint of the jump is defined as the melting point.

Figure 3C displays the predicted melting points for all mono-compounds, experimentally measured data for TiB₂, ZrB₂, HfB₂, NbB₂, and TaB₂ by Ruby et al., and estimated melting points for HEMB₂s by averaging the melting points of their constituent compounds^[38]. The results demonstrate that the melting points predicted by our DP model closely match experimental measurements, indicating the model’s accuracy. For high-entropy compounds, their melting points closely align with the average values of their constituent compounds. This finding suggests that the mean value is a suitable initial estimate for the melting point of HE-TMB₂s. The melting points of unstable compounds MoB₂ and WB₂ are also plotted in Figure 3C, which cannot be measured experimentally. With these values for mono-compounds and the mean value approach, the melting point of any composition can be estimated, which can be applied in the future for a better design of high-entropy diborides.

Grain boundary segregation

Grain boundary engineering through solid solution segregation has emerged as a promising strategy to improve the high-temperature strength of diboride and carbide UHTCs. This approach has been validated by both experimental evidence^[6,37,39-41] and theoretical simulations^[25,26], and has proven effective for both conventional and high-entropy UHTCs^[40,41].

In this section, we demonstrate the capabilities of the DP model in simulating grain boundary segregation and its enhancement of grain boundaries. We constructed 14 symmetric tilt grain boundaries, similar to our previous work^[26]. The construction method is detailed in Supplementary Text 2, Supplementary Figure 3 and Supplementary Table 2. A grain boundary plane can be defined by two in-plane vectors, which are contained in its name. One is the rotation axis of the tilt grain boundary. The other is a combination of the other two basis vectors. For example, the 010_21 grain boundary has a rotation axis of [010] and another in-plane direction of 2[001]+$$ [\overline{21}{0}] $$, while the 210_12 grain boundary has a rotation axis of [210] and another in-plane direction of [001]+2[010]. For more information of the constructed grain boundaries [Supplementary Table 2]. The relaxed atomic structure of each grain boundary in ZrB₂ is shown in Supplementary Figure 4. Hybrid Monte Carlo/molecular dynamics (MC/MD) simulations were employed to examine grain boundary segregations. The simulations were performed under the NPT ensemble with P = 0 Pa, T = 2,000 K, and a time step of 1 fs. The damping parameters were set to 500 fs for pressure and 100 fs for temperature. Any pair of different metallic elements was swapped every five MD steps, with the acceptance probability of swapping following the Metropolis criterion: min(1, exp(-ΔE/k_BT)), where ΔE is the energy change due to swapping and k_B is Boltzmann’s constant. The total MC/MD simulation duration was 500,000 steps. A virtual tensile simulation was adopted to evaluate grain boundary strength, aligning the tensile direction normal to the grain boundary plane. The last snapshot of the MC/MD simulation was used to calculate grain boundary strength. The tensile simulations were conducted under the NPT ensemble at 2,000 K and an in-plane pressure of 0 Pa, with no constraints on stress along the tensile direction. The time step and damping parameters were consistent with those used in the segregation simulations, and the strain rate was set to 1 × 10⁹ s^-1. During tensile simulation, stress variation is approximately 1 GPa in the tensile direction. Consequently, when the strength difference is less than 1 GPa, distinguishing which grain boundary is stronger becomes inconclusive due to the margin of error.

Initially, we investigated W segregation in ZrB₂ grain boundaries and the grain boundary strengthening effect it caused. The average W concentration at the metallic sites is established at 5 at%, marginally exceeding the experimentally determined W levels in ZrB₂ bulk, which range between 2 at% and 4 at%^[8]. Following MC/MD simulations, W content in the bulk decreases owing to grain boundary segregation, aligning the concentration with experimental observations. Figure 4A displays the occupation probability of W at each metallic site for the 010_21 and 210_11 grain boundaries. Snapshots were extracted every 1,000 steps from steps 250,000 to 500,000 for statistical analysis. The occupation probability of a given site was calculated by N_W/N_frames, where N_frames is the total number of selected frames and N_W accounts for the number of times the site is occupied by W. Noticeable W segregation around grain boundaries can be observed in Figure 4A, with W forming a monolayer at the grain boundary plane in the 010_21 grain boundary and an off-center bilayer pattern in the 210_11 grain boundary. In our earlier studies^[26,42], we demonstrated that size effects predominantly govern segregation behavior, leading to W atoms favoring smaller-sized grain boundary sites. The segregation pattern is mainly governed by the strain state of each grain boundary^[26]. Additionally, the concentration profiles of W (R_W) were calculated from the selected frames. Slices parallel to the grain boundary plane with a thickness of 10 Å were cut, and the W concentration in the slice was adopted as the concentration at the midpoint of the slice. Figure 4B displays the W concentration profiles of all calculated grain boundaries, revealing prominent W segregation around these boundaries, with W concentrations around grain boundaries reaching as high as ~20 at% in metallic sites, which are much higher than W contents in the bulk. Figure 4C compares the ideal grain boundary strength between pure ZrB₂ and ZrB₂ with W segregation. For pristine ZrB₂ grain boundaries, the strength of the grain boundaries investigated in this study mainly depends on two factors. First, grain boundaries such as the 010_32, 210_13, and 210_23 types demonstrate low strength attributed to improperly self-accommodated structures, as evidenced by the strained B−B bonds near the grain boundaries. The relaxed structures in Supplementary Figure 4 show some bonds extending over 2 Å, signaling the weak bonding character of these grain boundaries and hence their reduced strength; Second, the strength of well-self-accommodated grain boundaries is influenced by their orientation relative to the boron net plane. Typically, a smaller angle between the grain boundary normal and the boron net plane correlates with higher strength, as exemplified by the robust 010_31 grain boundary. With W segregation, the results indicate that most of the grain boundaries are significantly enhanced (improved by ~30%), except 010_13, 010_31, and 210_21, where the ideal strengths decrease slightly within the error level of calculation. Due to the large size of Zr, some of the boron bonds around grain boundaries are stretched and weakened, leading to low strength. Replacing Zr by W around grain boundaries is equivalent to introducing additional compressive strains to grain boundaries, reducing boron bond lengths around grain boundaries, which then enhances grain boundaries^[26]. Especially for the 010_21 boundary, segregation of W will increase its strength substantially, which fractures shifting by one atomic layer into one of the crystals instead of the grain boundary plane^[26]. The enhancement induced by segregation is highly dependent on the atomic structure of a grain boundary. For instance, the 010_21 grain boundary exhibits the smallest in-plane periodicity [Figure 4A and Supplementary Figure 4], which allows near-total coverage by segregated W atoms at the boundary sites. This near total occupancy, if it indeed reinforces local strength, would result in considerable strengthening of the entire boundary. Conversely, a larger in-plane periodicity, as in the 010_13 grain boundary, results in tungsten atoms replacing only some sites^[26], leading to a lower degree of segregated coverage. Subsequently, the strength of this grain boundary is marginally influenced by segregation. The reduction in the ideal strength of the 010_31 boundary is consistent with our previous findings^[26], where the W segregation level in this boundary is minimal [Figure 4B]. After MC/MD simulation, we obtained an asymmetrical segregation pattern of W around the boundary. Such asymmetry destroys the well-accommodated structure, potentially weakening the grain boundary. For the 210_21 grain boundary, additional space permits increased W segregation. As indicated in Supplementary Figure 4, the metal elements shift to one side of the grain boundary in the 210_21 boundary, creating a large volume on the opposite side. The MC/MD methods employed here are incapable of filling this free volume, likely contributing to the slightly diminished grain boundary strength. Nevertheless, most of the grain boundaries can be enhanced due to grain boundary segregation, which may improve high temperature strength of ZrB₂, agreeing with experimental observations^[6,37,39].

Figure 4. (A) The occupation probability of W at each metallic site, with noticeable W segregation around grain boundaries. W forms a monolayer at the grain boundary plane in the 010_21 grain boundary, while it forms an off-center bilayer pattern in the 210_11 grain boundary. Gray atoms are B, and colored atoms are metals; (B) The concentration profile of W (R_W) at the metallic site surrounding various grain boundaries, displaying prominent W segregation around these boundaries. The unit of horizontal axis is Å; (C) A comparison of the ideal grain boundary strength between pure ZrB₂ and ZrB₂ with W segregation. Ideal grain boundary strength values are listed in Supplementary Table 3.

In addition to micro-alloying for creating conventional solid solutions, heavily alloying with multiple elements to produce medium- or high-entropy diborides is another effective method for tailoring material properties. Our previous work has reported grain boundary segregation in medium- and high-entropy carbides through simulations, as well as the grain boundary strengthening effect due to segregation^[25]. Here, the DP model is employed to investigate grain boundary segregation in a medium-entropy diboride Ti_1/3Zr_1/3Hf_1/3B₂ (ME-TMB₂). Figure 5A illustrates the occupation probability of Ti at each metallic site in the 010_21 and 210_11 grain boundaries, with results similar to W segregation in ZrB₂ grain boundaries. Ti forms a monolayer at the grain boundary plane in the 010_21 grain boundary, and adopts an off-center bilayer pattern in the 210_11 grain boundary with Ti depletion at the grain boundary plane. Figure 5B presents the concentration profile of Ti (R_Ti) at the metallic site surrounding various grain boundaries, indicating significant Ti segregation around most of these boundaries. The concentration increase of Ti is not as substantial as W in ZrB₂ grain boundaries, suggesting a lower segregation tendency of Ti. Figure 5C compares the grain boundary strength of pure ZrB₂ and Ti_1/3Zr_1/3Hf_1/3B₂ after MC/MD simulations, demonstrating that all the grain boundaries of Ti_1/3Zr_1/3Hf_1/3B₂ are stronger. This finding indicates that medium- or high-entropy approaches are also effective in designing diborides with robust high-temperature strength.

Figure 5. (A) The occupation probability of Ti at each metallic site, with noticeable Ti segregation around grain boundaries. Ti forms a monolayer at the grain boundary plane in the 010_21 grain boundary, while it forms an off-center bilayer pattern in the 210_11 grain boundary. Gray atoms are B, and colored atoms are metals; (B) The concentration profile of Ti (R_Ti) at the metallic site surrounding various grain boundaries, displaying prominent Ti segregation around most of these boundaries. The unit of the horizontal axis is Å; (C) A comparison of the ideal grain boundary strength between pure ZrB₂ and medium-entropy Ti_1/3Zr_1/3Hf_1/3B₂ (ME-TMB₂) after MC/MD simulation. Ideal grain boundary strength values are listed in Supplementary Table 3. MC/MD: Hybrid Monte Carlo/molecular dynamics.

To further enhance grain boundaries and improve high-temperature strength, we can incorporate micro-alloying into medium- or high-entropy diborides, such as by introducing a solid solution of W into Ti_1/3Zr_1/3Hf_1/3B₂. We added 5 at% of W to the metallic site of Ti_1/3Zr_1/3Hf_1/3B₂ and then simulated segregation. As illustrated in Figure 6A, substantial W segregations were obtained, similar to those in ZrB₂, with Ti depletion occurring around grain boundaries due to W’s higher segregation tendency. Figure 6B displays a comparison of the grain boundary strengths between Ti_1/3Zr_1/3Hf_1/3B₂ alloys with and without W segregation. It highlights that the grain boundaries 010_11, 010_21, 210_11, 210_23, and 210_32 experience a significant reinforcement from W micro-alloying. All of these grain boundaries have normals that form a small angle with the boron net plane. While certain other grain boundaries display slight improvements or weakenings owing to W micro-alloying, these changes fall within the limits of calculation error. In general, the findings demonstrate that, on average, W segregation contributes to the enhancement of grain boundary strengths. This result indicates that combining the medium/high-entropy design strategy with the micro-alloy strategy is effective in designing new high-performance materials.

Figure 6. (A) The concentration profile of W (R_W) and Ti (R_Ti) at the metallic site surrounding various grain boundaries, displaying prominent W segregation around these boundaries, which results in Ti depletion around these boundaries. The unit of the horizontal axis is Å; (B) A comparison of the ideal grain boundary strength between pure Ti_1/3Zr_1/3Hf_1/3B₂ (ME-TMB₂) and Ti_1/3Zr_1/3Hf_1/3B₂ with W segregation. Ideal grain boundary strength values are listed in Supplementary Table 3.

When generating the training data for the DP model, no grain boundary structure is included. We examine the reliability of the DP model in predicting grain boundaries by evaluating grain boundary energies and segregation energies, with the results shown in Supplementary Text 3, Supplementary Table 4 and Supplementary Figure 5. Even though the DP model does not exhibit the same high accuracy as it does for other properties, its predictive capabilities on grain boundaries are reasonable, correctly capturing the elementary dependent grain boundary properties. In addition, our simulation results are in good agreement with previous simulations on ZrB₂ with a solid solution of W^[26], and they also correspond well with experimental results in ZrB₂ and Ti_1/3Zr_1/3Hf_1/3B₂^{[4,6,11,37,43]}. In experiments, the segregation of W into ZrB₂ grain boundaries has been confirmed by transmission electron microscopy^[39]. Furthermore, the high-temperature fracture surface transforms from being dominated by intergranular features to being dominated by transgranular features with the addition of WC^[6]. The grain boundaries also exhibit improved thermal corrosion resistance^[37], which indicates a reduction in grain boundary energy and an enhancement in grain boundary strength due to W segregation. As reported in ref^[11], similar to ZrB₂ with WC as a sintering additive, medium-entropy Ti_1/3Zr_1/3Hf_1/3B₂ also demonstrates good strength retention at elevated temperatures, and the addition of WC can further improve the strength of Ti_1/3Zr_1/3Hf_1/3B₂, which is consistent with our simulation results. In addition to the strengthening effect, the addition of WC can also help remove surface oxides, resulting in cleaner grain boundaries^[4]; pin grain boundaries, leading to an improved microstructure uniformity; and increase grain boundary oxidation resistance due to the low activation of W.

Discussion

In previous sections, we have demonstrated the accuracy and reliability of the DP model in predicting various properties, such as lattice parameters, elastic constants, equations of state, melting points, and grain boundary segregations. This indicates that our DP model can be effectively applied to TMB₂ mono-compounds or solid solutions, regardless of whether they are micro-alloyed, heavily alloyed, or even designed as high-entropy materials. Consequently, the model serves as an efficient and accurate tool for predicting material properties, exploring micro-mechanisms, and ultimately assisting in material design. This paves the way for the customization and design of these materials to meet the demands of various applications.

While the model shows good accuracy, a few properties have not been well fitted, such as the C₃₃ of WB₂, as shown in Table 1. In machine learning potentials, there is still a need to strike a balance between accuracy and efficiency. The DPA-1 model achieves a satisfactory balance between these two factors, making it a suitable choice for simulating systems with up to 10,000 atoms. The prediction errors in energy and force values for this model are 8.8 meV/atom and 387 meV/Å, respectively. If a more accurate model is required, one can opt for a more complex model, such as the DPA-2 model^[27] available in the DeePMD-kit software. This model can improve the prediction error by approximately 30%, with errors in energy and force being 6.0 meV/atom and 245 meV/Å, respectively. However, this increased accuracy comes at the cost of reduced efficiency, which is one to two orders of magnitude slower than the DPA-1 model. As a result, the DPA-2 model may only be suitable for simulations involving hundreds to thousands of atoms. For scenarios where simulations with millions of atoms are frequently conducted, a more efficient model is needed, such as the NVNMD model^[44] in the DeePMD-kit software. This model offers an efficiency that is one order of magnitude faster than the DPA-1 model. However, the trade-off is an increase in prediction errors for energy and force values, which rise to 13.2 meV/atom and 410 meV/Å, respectively.

In addition to the models, our highly precise DFT dataset also plays a significant role in the impending revolution of atomic-scale simulations. The potential of big machine learning models is set to bring about substantial changes to the field of material research. In order to develop reliable big machine learning potential models, it is crucial not only to create new artificial intelligence models but also to enrich valuable datasets for model training. While large high-throughput DFT datasets have already been developed under the support of the Materials Genome Initiative, such as the Materials Project^[45], AFLOW^[46], and OQMD^[47], there is still a need for many domain-specific datasets to cater to the varying demands of different fields. These datasets will not only aid in the development of big machine learning potential models but also assist in driving the revolution in materials research. By focusing on enriching our dataset resources and refining our artificial intelligent models, we can truly unlock the full potential of atomic-scale simulations and revolutionize the field of material research.