Atomistic-to-Macroscopic Mechanical Modeling of (HfNbTaZr)B₂ High-Entropy Ceramics via a Developed Interatomic Potential

1. Introduction

In recent years, the concept of entropy stabilization has guided the development of various multi-component materials. By combining several elements in equal or similar atomic ratios, the increased configurational entropy helps to form stable solid solutions [1,2,3]. The resulting features (e.g., lattice distortion, sluggish diffusion and the cocktail effect) generally improve the physical and mechanical performance of the material [4]. Following this principle, high-entropy alloys (HEAs) have been first designed and investigated [2,3]. Compared to traditional alloys, HEAs commonly display better strength-ductility synergy [5], radiation resistance [6], and good cryogenic [5] and catalytic properties [7]. However, the use of HEAs in high-temperature environments is sometimes limited by high-temperature softening and oxidation issues [8,9].

To meet the material requirements for more extreme environments (e.g., higher temperature), researchers applied the entropy-stabilization approach to nonmetallic materials, developing high-entropy ceramics (HECs) [10,11]. These materials are generally synthesized from mixtures of basic oxides [11], nitrides [10], carbides [4,10], or diborides [12]. In HECs, the random distribution of metal cations [6,13] and the shielding effect of the anion network [1,4,11,14] contribute to good phase stability. Particularly, the high-entropy diborides (i.e., HEBs with the AlB2 lattice structure) are found to maintain their hardness, strength, oxidation resistance and certain degree of ductility at ultrahigh temperatures. Because of these thermal and mechanical features, HEBs provide a practical alternative to conventional ultrahigh temperature ceramics and refractory alloys, showing potential for use in demanding fields like aerospace and nuclear engineering [1,15,16].

Research on HEBs is currently at an early stage. Because these materials are relatively difficult to synthesize, many studies have focused on improving preparation methods to achieve dense and uniform structures [17,18,19]. While basic material properties (e.g., thermal conductivity [20], elastic constants [16,21], and oxidation resistance [22,23,24]) have been tested and reported, the relationship between their nanoscale structures and macroscopic properties has not been fully investigated. Detailed atomic level factors, including chemical short-range order (CSRO), elemental diffusion, and defect behavior, are still not well understood. Without such nanoscale information, predicting material performance can be difficult. Current predictions often rely on simple empirical methods [25,26,27], such as the rule of mixing, entropy forming ability, or valence electron concentration. However, the predicted values can differ from experimental results, and materials with similar empirical parameters may show different actual properties [28]. This inconsistency makes it harder to design HEBs for specific applications.

Molecular dynamics (MD) simulations provide a useful way to study the atomic level mechanisms. To perform reliable MD simulations, accurate interatomic potentials are needed to describe the interactions among different elements. Recently, machine learning (ML) models have been used to develop these potentials [29,30,31,32]. While ML potentials usually match their training data well, their data-driven nature can limit their accuracy when applied to situations outside their training set. As a result, they may be less capable of simulating dynamic or non-equilibrium processes, such as large deformation, crack growth, dislocation movement, or radiation damage. By comparison, classical empirical or semi-empirical potentials are based on physical models [8,9,14,33], which generally helps them perform more stably in new or extreme conditions outside the training domain. The main difficulty in developing classical potentials for HEBs is the complex bonding environment, which involves a mix of metallic, ionic, and covalent bonds [14,15]. Describing these mixed interactions requires mathematical formulas that account for bond angles and bond orders. These formulas contain a large number of internal parameters, making the fitting process too complicated for conventional optimization methods [9,14,33].

To enable MD simulations for the HEB materials, this study develops an analytical bond-order potential (ABOP) for Hf/Nb/Zr/Ta/B exemplified system. To handle the complex parameter fitting process, we use a modified distributed breeder genetic algorithm (DBGA) based on ab initio reference data. The remainder of this paper is organized as follows. Section 2 introduces the ABOP model and explains the parameter optimization procedure. Section 3 evaluates the accuracy of the developed potential by comparing it with ab initio reference data. In Section 4, we apply the potential to several typical dynamic simulations to test its practical performance. Overall, this work provides a basic tool for studying the nanoscale mechanical behaviors and the associated atomic mechanisms of high-entropy ceramics.

2. Computational Details

2.1. The Analytical Bond-Order Potential (ABOP)

The fidelity of atomistic simulations depends on an interatomic potential that accurately reflects the underlying bonding nature of the material. Given the pronounced directional covalent bonding in high-entropy ceramics such as (HfNbTaZr)B₂, an angular-dependent formulation that captures local bond-order effects is essential. Here, we implement the analytical bond-order potential (ABOP) framework [14,34], defining the total potential energy as:

$$V = {\sum }_{i,j}{V}_{\mathit{ij}} = {\displaystyle \frac{1}{2}}\sum _{i,j}{f}_c\left(r_{\mathit{ij}}\right)\left[V_R\left(r_{\mathit{ij}}\right) - b_{\mathit{ij}}{V}_A\left(r_{\mathit{ij}}\right)\right],$$

(1)

where $r_{\mathit{ij}}$ denotes the interatomic distance, $f_c\left(r_{\mathit{ij}}\right)$ is a smooth cutoff function, $V_R\left(r_{\mathit{ij}}\right)$ and $V_A\left(r_{\mathit{ij}}\right)$ represent the pairwise repulsive and attractive contributions, respectively. The many-body environment is implicitly encoded in the bond-order term b_ij, which incorporates local coordination numbers and bond-angle dependencies. This mathematical architecture enables ABOP to model complex, hybridized metallic-covalent bonding within a unified framework.

While ABOP has been successfully parameterized for various intermetallic, ionic, and covalent systems [34,35], abrupt cutoff truncations frequently induce unphysical artificial ductility and overestimate theoretical strength during tensile or fracture simulations [14]. To resolve this artifact, the cutoff parameters in fc are dynamically refined during parameter optimization. Additionally, to accurately capture high-energy collision cascades and close-contact interactions, the short-range repulsive domain is smoothly connected to the Ziegler–Biersack–Littmark (ZBL) potential following the strategy reported previously [14].

2.2. Training Data Generation

To parameterize the ABOP, a comprehensive reference dataset containing elemental metals, intermetallic compounds, and diborides is constructed via ab initio calculations. The representative configurations and corresponding target properties are summarized in Table 1. Notably, atomic virial stress tensors are incorporated as an independent fitting objective, which is a strategy demonstrated to substantially improve the robustness and mechanical transferability of interatomic force fields [36]. Benefiting from the physically motivated functional form, the ABOP requires a significantly more compact dataset than typical ML potentials [9,14].

All density functional theory (DFT) calculations are performed using the Vienna Ab Initio Simulation Package (VASP) [37] within the projector augmented-wave (PAW) framework. Electronic exchange-correlation effects are treated using the Perdew–Burke–Ernzerhof (PBE) generalized gradient approximation (GGA) [38]. The plane-wave cutoff energy is set to 520 eV, and ionic relaxations are converged until the Hellmann–Feynman forces are below 10^-4 eV/Å. The Brillouin zone is sampled using Γ-centered Monkhorst–Pack k-point grids generated via VASPKIT with a k-mesh resolution of $0.02\times 2\pi {\AA }^{-1}$ [39]. Initial structural prototypes are retrieved from the Open Quantum Materials Database (OQMD) [40].

2.3. Global Optimization

With the ab initio database established, parameterizing the quinary Hf/Nb/Zr/Ta/B ABOP requires solving a 254-dimensional optimization problem. The optimal parameters are determined by minimizing a normalized objective loss function (i.e., ∆) formulated as:

$$\Delta ={\sum }_i\frac{{\omega }_i(E_{\mathrm{ABOP},i}-E_{\mathrm{ab} \mathrm{initio},i}{)}^2}{E_{\mathrm{ab} \mathrm{initio},i}^2}+{\sum }_j{\sum }_{\alpha }\frac{{\omega }_j(S_{\mathrm{ABOP},j,\alpha }-S_{\mathrm{ab} \mathrm{initio},j,\alpha }{)}^2}{S_{\mathrm{ab} \mathrm{initio},j,\alpha }^2}+{\sum }_k{\sum }_{\beta }\frac{{\omega }_k(M_{\mathrm{ABOP},k,\beta }-M_{\mathrm{ab} \mathrm{initio},k,\beta }{)}^2}{M_{\mathrm{ab} \mathrm{initio},k,\beta }^2}$$

(2)

where E, S, and M denote energetic, structural (e.g., lattice vectors, atomic coordinates), and mechanical quantities (e.g., elastic constants, virial stress components), with indices (i, j and k) covering the respective sizes of dataset. The directional indices α=1~3 and Voigt stress components β=1~6 resolve spatial and tensorial dimensions. Although weight factors ω critically govern both convergence rate and predictive power [9,14,33], their assignment is often challenging due to disparate physical dimensions and orders of magnitudes. By employing a relative error normalization in Eq. 2, magnitude disparities across different properties are reduced. Therefore, the default weights are set to unity (ω = 1), except for cohesive energies which receive an increased weighting of ω = 10 to emphasize phase stability in dynamic simulations.

Searching such an ultrahigh-dimensional, non-convex parameter landscape makes gradient-based approaches (e.g., conjugate-gradient or steepest-descent methods) susceptible to convergence to spurious local minima [38,41]. Therefore, we employ a derivative-free global optimizer based on a modified distributed breeder genetic algorithm (DBGA) [9,14,33], which accelerates convergence by emulating artificial selection principles. As depicted in Figure 1, the workflow initializes a population of NP = 250 candidate parameter sets, sampled via a quasi-random Sobol sequence [42]. This sampling method guarantees high parameter space coverage and initial genetic diversity. Thereafter, the generated population is partitioned into five parallel subpopulations (so-called breeder pools) for localized evolutionary operations, including tournament selection, extended intermediate recombination, and a two-stage hybrid (i.e., uniform and breeder) mutation [9,14]. This dual-mutation strategy balances broad re-initialization with fine local perturbation, thus preventing premature stagnation, particularly in late-stage runs. Finally, the independent breeder pools are merged, and the top NP individuals ranked by fitness (defined as 1/∆) are retained to initialize the next generation. This process iterates until the predefined convergence criteria are met. More details of the global optimization procedure can be found in Refs. [9,14,33].

2.4. Molecular Dynamics (MD) Simulations

All atomistic simulations are conducted within the LAMMPS package [43] using a 1.0-fs timestep. Nosé–Hoover thermostat and barostat are utilized to control the temperature and pressure of the simulation system. Post-processing and visualization are carried out via OVITO software [44]. Specifically, local crystalline lattice and phase transformations are monitored using Polyhedral Template Matching (PTM) [45].

3. Results

3.1. Performance in Metallic Systems

The predictive performance of the developed ABOP for elemental Hf, Nb, Zr, and Ta is summarized in Figure 2 and Figure 3, compared with data from ab initio calculations. Apparently, the ABOP correctly reproduce structural, energetic and mechanical properties of the four metallic materials. Notably, the formation energies of defects (i.e., surface, vacancies and interstitials) are not included in training process, but the developed potential still captures these properties in well agreement with ab initio results (Figure 2c-e). Such a cross-validation strategy effectively avoids overfitting, thereby elevating transferability of the ABOP dynamic simulations. In addition, the equations of states (EOS) predicted by the ABOP also consistent well with ab initio data. The potential accurately tracks the energy profiles across various crystalline polymorphs under wide hydrostatic strains (±50%, Figure 3), thus precluding spurious local minima near equilibrium and preserving structural robustness during dynamic simulations [9,14,33].

The EOS curves for some binary intermetallics with various structural phases (e.g., $B2$, $L{1}_2$, $D{0}_3$, $D{0}_{19}$, etc.) are presented in Figure 4. The developed potential accurately describes both the ground-state energies and the energetic responses under finite volumetric dilation and compression. Furthermore, to assess the ABOP’s predictive capacity beyond the training domain, calculations are extended to ternary compounds (Figure 5). The resulting EOS profiles match the ab initio reference data with high accuracy, which indicates that this potential successfully learns the underlying many-body cross-interactions. This demonstrates its applicability in compositional variations across complex alloying spaces.

Apart from static target properties, evaluating the ability of the potential to preserve structural stability during finite-temperature dynamical simulations is essential. As reported in previous studies [46,47], complex potentials with flexible mathematic forms and high parameter counts (e.g., ML-based and bond-order formalisms) often suffer unphysical structural breakdown in dynamic simulations, even though they accurately describe static DFT reference data. To examine this dynamic stability, chemical short-range order (CSRO) in quaternary HfNbZrTa HEA is studied at room temperature. The bulk HEA maintains its BCC lattice structure throughout the Monte-Carlo relaxation at 300 K. An apparent local clustering of Hf-Zr and Nb-Ta pairs is evident from both the optimized atomic configuration (Figure 6a) and the Warren-Cowley (WC) order parameters (Figure 6b). The WC parameters can be calculated using [9,14]

$${\alpha }_{\mathit{ij}}^n=\left(p_{\mathit{ij}}^n - C_j\right) / \left({\delta }_{\mathit{ij}} - C_j\right),$$

(3)

where n represents the n-th nearest-neighbor coordination shell of species i, p_ij is the probability of encountering a j atom around an i atom, C_j is the nominal concentration and δ_ij is the Kronecker delta function. Negative values of the WC parameters indicate a clustering tendency for the given pair. The clustering between Hf-Zr and Nb-Ta agrees with their respective ground-state crystal phases (i.e., HCP for Hf/Zr and BCC for Nb/Ta). Moreover, the CSRO pattern is also consistent with previous experiments [48,49]. Findings in Figure 6a and Figure 6b indicate that the ABOP correctly learns interactions among metals at finite temperatures.

In addition, atomic deposition simulations are also performed at room temperature (Figure 6c). The ABOP accurately captures the vapor deposition process, wherein vaporized atoms spontaneously crystallize into the equilibrium BCC lattice on the substrate. This further demonstrates the robustness of the developed potentials for modeling HEAs under far-from-equilibrium extreme conditions, such as radiation cascades and shock loading [8,33]. Overall, the results in Figure 6 confirm the reliability and broad transferability of the ABOP formulation.

3.2. Performance in Ceramic Systems

The basic properties of four unary diborides HfB₂, NbB₂, ZrB₂ and TaB₂ (hexagonal lattice, space group: P6/mmm) calculated using ABOP and ab initio methods are displayed in Figure 7. The cohesive energies, lattice and elastic constants of these diborides are correctly reproduced by ABOP. Notably, given the particular lattice structures of the AB₂-type diborides, they exhibit two different (0001) surfaces terminated by metals and boron atoms, respectively. Thus, a reliable potential must correctly differentiate these two surfaces and give proper predictions on the corresponding surface formation energies. This is highly relevant to the reliability of simulations involving surfaces and interfaces, such as nanopillar compression, nanoindentation, grain boundaries and surface ion implantations. Figure 7c and Figure 7d demonstrate that the ABOP can predict the metal- and boron-terminated (0001) surface energies in high consistence with ab initio data.

In addition, the EOS curves of unary diborides calculated using ABOP in comparison with ab initio data are shown in Figure 8a. The MD results are consistent with ab initio training data within a wide range of hydrostatic strain. Moreover, the EOS curves of six binary diborides are displayed in Figure 8b, which are not included in training process. The ABOP predictions are consistent with ab initio data. This cross-validation result suggests that the developed potential is able to capture underlying physics about elemental interactions within ceramic phases.

According to most previous experimental investigations [1,4,11], the HEB usually exhibits outstanding stability and microscale homogeneity due to the boron anion sublattice, which screens metal cations from one another and effectively promotes configurational disorder. However, owing to the differences in atomic radii and electronegativities among the constituent transition metals, local CSRO should inherently develop at the atomistic level. Figure 9a and Figure 9b display the ABOP-predicted CSRO within the HEB, which suggest Hf-Nb and Zr-Ta pairing. Such a clustering tendency among cations, influenced by the interlocked anion sublattice, is different from that within the HfNbZrTa HEA (Figure 6). To verify the formation of cation CSRO, we evaluated the solution energies corresponding to various cation substitutions in the four unary diborides. The energetic profiles in Figure 9c reasonably explain results in Figure 9b. Nb, Zr and Ta show lowest solution energy in HfB₂, TaB₂ and ZrB₂. Hf shows negative solution energies in NbB₂, ZrB₂ and TaB₂, which is consistent with the negative WC parameters of Hf-Nb, Hf-Zr and Hf-Ta. Because CSRO is known to exert a major influence on macroscopic mechanical performance in complex concentrated systems [29,46,50], the developed ABOP provides a powerful computational tool to uncover these fine atomistic features and clarify their influence on mechanical responses.

Besides, Figure 9d shows that the ABOP-predicted lattice constant a = 3.17 Å at 300K (c = 3.42 Å), corresponding to a density of 8.77 g/cm³ reasonably consistent with experimentally measured values [12,51]. The variation of lattice constant with temperature described by ABOP is also consistent with ab initio results. Meanwhile, the ABOP-predicted temperature-dependence of elastic constants of HEB is shown in Figure 9e. C₁₁, C₁₂ and C₁₃ slightly decreases with increasing temperature, while C₃₃ and C₄₄ are almost invariant as temperature rises. The room temperature elastic constants predicted by ABOP and ab initio method are in reasonable agreement. Furthermore, the elastic constants follow the criteria of mechanical stability: C₁₂ > 0, C₄₄ > 0, C₁₁ > C₁₂, and (C₁₁ + C₁₂)C₃₃ – 2C₁₃² > 0. Overall, results in Figure 9 validate the ABOP’s applicability in dynamic simulations of complex ceramic systems at finite temperatures.

4. Discussions

As presented in the previous section, the developed ABOP is able to predict the properties of elements, intermetallic compounds and diborides in Hf/Nb/Zr/Ta/B system. Additionally, it is also proven capable of maintaining structural stability in simulations under finite temperatures. The robustness and reliability of the ABOP have been manifested. In this section, we further validate the broad applicability of the ABOP using several dynamic simulation scenarios.

4.1. Compressive Response of Nanopillars

The mechanical behavior of (HfNbTaZr)B₂ is further assessed under uniaxial compression using the developed ABOP; the simulation model is displayed in Figure 10a. The temperature is maintained at 300 K during the whole loading process; a typical MD strain rate of 10⁹/s is used [46,52,53]. To isolate the influence of CSRO, comparative simulations are conducted on both a random solid solution model and a MC-optimized model. As illustrated by the stress-strain responses in Figure 10b, both models display overall brittle failure behavior, consistent with prior assumptions [14,16,19,20,26]. The MC and Random models exhibit similar strength and flow stress. This indicates that CSRO exerts a limited effect on the mechanical properties of the HEB, which differs from the apparent strength-ductility synergy induced by CSRO in HEAs [46,54].

In addition to the similar stress-strain responses, the nanostructural evolution processes of the Random and MC models are also analogous. Here, only the atomic structure of Random model at varying deformation levels is shown in Figure 10c. Yielding of the nanopillar is triggered via formation of several local shear bands (Figure 10c2), which promptly releases internal strain energy and causes a significant stress drop. These shear bands are oriented at ~45^o to the pillar axis, corresponding to resolved maximum shear stress. During the post-yielding stage, deformation proceeds through multiplication and propagation of shear bands as shown in Figure 10c3. Shear strain localizes and accumulates within shear bands, which finally incubate cracks, as illustrated in Figure 10c4 and Figure 10c5. Figure 10 displays a typical manner of brittle failure; the transformation from shear bands to cracks have been reported in many other brittle ceramics [14].

4.2. Nanoindentation

To facilitate a direct comparison with macroscopic mechanical properties of the HEB, we simulated a nanoindentation test on the (0001) surface (Figure 11a). A 60×60×25 supercell of single crystalline HEB is generated containing 270,000 atoms. The indentation is performed using a virtual spherical indenter with radius of 2.5 nm, which moves at a velocity of 20m/s. The repulsive force between the indenter and nearby atoms can be expressed as $F\left(r) = K\theta \right(R - r\left)\right(R - r{)}^2$, where K is a specified force constant which equals 3.3eV/ų [14]. From the simulated force-displacement curve (Figure 11b), the Young’s modulus (E_ABOP) and hardness (H_ABOP) are extracted as 435 GPa and 27.7 GPa, respectively. In addition, E_ABOP and H_ABOP can be verified using theoretic methods. Based on the previously calculated elastic constants (Figure 9e), the mechanical properties can be evaluated analytically:

$$B_{\mathrm{V}} = {\displaystyle \frac{2}{9}}\left(C_{11} + C_{12} + 2C_{13} + {\displaystyle \frac{1}{2}}{C}_{33}\right)$$

(4)

$$B_{\mathrm{R}}={\displaystyle \frac{\left(C_{11}+C_{12}\right)C_{33} - 2C_{13}^2}{C_{11}+C_{12}+2C_{33 }- 4C_{13}}}$$

(5)

$$G_{\mathrm{v}}={\displaystyle \frac{1}{15}}\left(2C_{11}+C_{33} - C_{12} - 2C_{13}\right)+{\displaystyle \frac{1}{5}}(2C_{44}+C_{66})$$

(6)

$$C^2=\left(C_{11}+C_{22}\right)C_{33} - 2C_{13}^2$$

(7)

$$G_{\mathrm{R}}={\displaystyle \frac{5C^2{C}_{44}{C}_{66}}{2\left[3B_V{C}_{44}{C}_{66}+{ C}^2\left(C_{44}+C_{66}\right)\right]}}$$

(8)

$$B_{\mathrm{H}}={\displaystyle \frac{1}{2}}\left(B_{\mathrm{V} }+B_{\mathrm{R}}\right)$$

(9)

$$G_{\mathrm{H}}={\displaystyle \frac{1}{2}}\left(G_{\mathrm{V}}+G_{\mathrm{R}}\right)$$

(10)

$$E={\displaystyle \frac{9B_{\mathrm{H}}{G}_{\mathrm{H}}}{3B_{\mathrm{H}}+G_{\mathrm{H}}}}$$

(11)

$$\nu ={\displaystyle \frac{3B_{\mathrm{H} }- 2G_{\mathrm{H}}}{2(3B_{\mathrm{H}}+G_{\mathrm{H}})}}$$

(12)

$$H=2{\left({\displaystyle \frac{G_{\mathrm{H}}^3}{B_{\mathrm{H}}^2}}\right)}^{0.585} - 3$$

(13)

$$K_{\mathrm{IC}}=V_0^{1/6}{G}_{\mathrm{H}}{\left({\displaystyle \frac{B_{\mathrm{H}}}{G_{\mathrm{H}}}}\right)}^{1/2},G_{\mathrm{IC}}=K_{\mathrm{IC}}^2\left({\displaystyle \frac{1-{\nu }^2}{E}}\right).$$

(14)

The calculated values are listed in Table 2. The theoretical moduluse E (481.18 GPa) and hardness H (26.30 GPa) are comparable to the values measured from the simulation. The Poisson’s ratio of 0.215 falls between typical values for strong covalent bonds (~0.1) and purely ionic materials (> 0.25) [14,15], indicating a mixed bonding character in the HEB. Moreover, existing studies have reported the following ranges for the key mechanical properties [16,18,19,20,55]: bulk modulus B = 276-301 GPa, shear modulus G = 228-237 GPa, Yongs’ modulus E = 340-503 GPa, hardness H = 22-31.9 GPa, and fracture toughness K_IC = 3.72-4.37 MPa·m^1/2. These consistent mechanical parameters confirm that the developed ABOP reliably captures the underlying physical correlations.

Additionally, Figure 11c illustrates the atomic level structural evolution during indentation, in which atoms are colored according to local shear strain. Clearly, atoms with significant shear strain locate close to the indented surface region. No subsurface plastic zone is observed, indicating markedly brittle nature. This phenomenon differs from that of certain hard ceramics exhibiting nanoscale plasticity; for instance, in the (HfNbZrTa)C HEC, a dislocation network penetrates the bulk material to a depth nearly twice the radius of the indenter tip [14]).

4.3. Multiscale Cohesive Zone Modeling of Polycrystalline (HfNbZrTa)B₂

Previous sections have demonstrated the reliability and applicability of the developed ABOP. Nonetheless, there exists a long-standing issue for MD simulation, namely, the MD results are difficult to directly compare with macroscale experimental tests or continuum modeling. Here, we propose a simple but effective multiscale method to connect MD to finite element method (FEM); the workflow is shown in Figure 12. The hexagonal HEB supercell is converted to an orthogonal cell using ATOMSK (Figure 12a1), which searches for an equivalent orthogonal cell while preserves the periodicity of the configuration [56]. Numerous bicrystalline models are generated by combining two orthogonal cells at different tilt angles (Figure 12a2). The regions near the central grain boundaries (GBs) are divided into bins to record evolution of displacement and force during subsequent uniaxial tensile deformation, as detailed in Refs. [57,58,59]. The traction-separation (TS) data from different bicrystalline models are gathered, as displayed in Figure 12a3. An atomic scale bilinear TS law is obtained by fitted these simulation data, characterized by the cohesive strength T_max = 52.24 GPa, the softening-onset separation δ₀ = 1.51 Å, the critical separation δ_c = 8.04 Å, and the fracture energy G_c = 0.5(T_max · δ_c) = 21.01 J/m². Notably, the critical energy strain rate is close to the theoretic value in Table 2, indicating that it is reasonable to obtain the TS law from bicrystalline tensile simulations.

The atomic scale TS law extracted from MD data quantifies the cohesive strength and fracture energy of GBs. However, directly applying this atomic scale TS law to a macroscopic FEM calculation is impractical due to the large difference in scales between the atomistic process zone and standard FE meshes [59,60,61]. According to linear elastic fracture mechanics, the length of the process zone ahead of a crack tip l_cz is calculated as l_cz = E · G_c / T_max². With a bulk elastic modulus E = 481 GPa (Table 2), l_cz is about 3.7 nm. Resolving this zone by at least n = 2-3 interface elements requires an element length l_e ≤ l_cz / n ≈ 1.2 nm [59,60]; even for a two-dimensional 100 × 100 μm² polycrystal this would demand O(10¹⁰) elements, far beyond computational feasibility. Moreover, since the entire failure separation (δ_c ≈ 8 Å) is much smaller than a standard macroscopic load step, the elements will fail instantly within one step; the fracture energy is then released in an unstable burst.

To solve these problems, the TS law is rescaled while keeping the total fracture energy G_c constant. We select the element size l_e = 1 μm, which gives a rescaled strength $T_{\max }^{\prime }$ = 1.83 GPa and a larger separation ${\delta }_c^{\prime }$ = 22.9 nm. The rescaled process zone $l_{\mathrm{cz}}^{\prime }$ = 3.0 μm (Figure 12b1). This rescaling method is reasonable because brittle fracture is mainly controlled by the energy release rate G_c [62]. Similar multiscale approaches, which transfer the energy-based fracture mechanism from the atomic level to the engineering level, have been proven reliable previously [59,60,61]. Thereafter, a 100 × 100 μm² representative volume element (RVE) containing 100 grains (mean grain size 11 μm consistent with experiments [16,25,51]) is generated as shown in Figure 12b2. Periodic boundary conditions are applied to the RVE edges. The model is discretized with a 1μm grid of bilinear plane-stress elements (CPS4R). To explicitly simulate intergranular fracture, zero-thickness cohesive elements (COH2D4) are inserted along all GBs. Consequently, the grains deform linearly and elastically; structural damage is governed by the cohesive interfaces.

The engineering stress-strain curve (Figure 12b3) initially shows a linear elastic behavior, reaching a peak stress of 2,683 MPa at a strain of 0.66%. This is followed by a sharp stress drop, where the load-carrying capacity falls below 20 MPa by 2% strain. This rapid failure indicates a brittle intergranular fracture, which is consistent with findings in the previous sections. Two aspects of this mechanical response are notable. First, the critical strain of 0.66% is slightly higher than theoretical elastic limit of the grains (~ 0.55%), which is caused by the early softening of GBs that adds compliance before fracture. Second, the peak stress of 2,683 MPa exceeds the assigned GB strength of 1.83 GPa. Because the load redistributes across the complex GB network during deformation, the overall strength differs from a simple average of individual GBs. The corresponding microstructural damage evolution is illustrated in Figure 13, including the Mises stress and GB damage across different deformation stages. During the initial loading stage (ε ≤ 0.49%), discrete damage begins at weaker GBs perpendicular to the tensile direction. As the stress approaches its peak (ε = 0.66%), these damaged GBs connect to form horizontal bands across several grains. The intact boundaries between these bands continue to carry the load, indicating that initial failure is driven by widespread microcracking rather than the propagation of a single dominant crack. Following the peak stress (ε ≥ 0.95%), the damage bands quickly merge into continuous cracks. As these cracks open along the GB network, the surrounding intact grains unload their stored elastic energy. The resulting fracture path is entirely intergranular, which is consistent with the fundamental assumption of the cohesive zone model.

5. Conclusions

In this study, an analytical bond-order potential is developed for the complex Hf/Nb/Zr/Ta/B system to investigate the nanoscale structural and mechanical behaviors of HEBs. The training process against an extensive ab initio dataset is achieved via a DBGA parameterization. The main conclusions are listed as follows:

1. The developed ABOP successfully reproduces the structural, energetic, and mechanical properties of unary metals, binary and ternary intermetallics, and diborides over a wide range of strains. Additionally, a cross-validation test demonstrates that this potential ensures structural stability during finite-temperature dynamic simulations and captures the CSRO pattern.

2. Using the developed ABOP, simulations of single crystalline HEB nanopillars under uniaxial compression revealed a brittle failure behavior. Yielding and subsequent fracture are dominated by the formation, multiplication, and cracking of local shear bands oriented at approximately 45° to the loading axis. Unlike in HEAs, CSRO exerts a limited influence on the mechanical response of the HEB.

3. Simulated nanoindentation on the (0001) surface gives a Young’s modulus of 435 GPa and a hardness of 27.7 GPa, consistent well with theoretical evaluations and existing experimental data. The structural evolution during indentation showed highly localized deformation without deep subsurface plastic zones, further denoting the intrinsic brittleness of the material.

4. A multiscale framework is established to bridge atomic scale MD simulations with macroscopic FEM modeling. By extracting an atomic TS law from bicrystalline MD models and rescaling it with conserved fracture energy G_c, the FEM-CZM simulations reproduce the brittle intergranular failure of polycrystalline HEB.