Towards nano-mechanical simulations of ceramics containing realistic defects via machine-learning potentials: the example of TiB₂

Abstract

Transition metal diboride (TMB₂) ceramics combine high hardness with excellent thermal and chemical stability, making them attractive for protective coating applications. TMB₂s commonly grow as largely off-stoichiometric—containing vacancies or other simple crystallographic defects—and a particularly essential question is how such defects alter the response to mechanical loads at the nanoscale. We exploit molecular dynamics (MD) equipped with here-trained machine-learning interatomic potential (MLIP) to reveal effects of point and planar defects during nano-mechanical tests of TiB₂, being a representative of the most common AlB₂-type TMB₂s. MLIP training consists of active learning on configurations from finite temperature ab initio molecular dynamics, including equilibrium and uniaxially loaded structures as well as various defective and/or extremely strained environments. Transferability to near-indenter-tip environments is achieved via on-the-fly training on extrapolative clusters from simple nanoindentation runs. Following MLIP's validation, we simulate room-temperature nanoindentation of TiB_2−x structures, where B sub-stoichiommetry is realized by disordered B vacancies, single- and double-planar defects previously revealed by electron microscopy. A somewhat non-intuitive prediction is that TiB_2−x structures can exhibit hardness comparable to the ideal TiB₂, challenging traditional assumptions about weakening effects of sub-stoichiometry. This behavior is ascribed to Ti-rich planar defects, contrarily to vacancies which, at the same chemistry, notably deteriorate mechanical properties. We hypothesize that similar effects may be expected for other group 5–6 transition metal diborides.

---

1 Introduction

Due to their high melting points, excellent hardness and distinguished chemical stability, transition metal diboride (TMB₂) thin films have been widely used in the field of cermat and wear-resistant materials.^1–4 TMB₂ s typically adopt layered hexagonal structures with the AlB₂-type (α; group IV–V TMs), WB₂-type (ω, group VI TMs), and ReB₂-type (γ, group VII TMs) phase prototype,⁵ where the phase preference may be tailored via synthesis conditions—especially for TMB₂ s exhibiting low energy barriers between phase polymorphs (e.g., NbB₂ and TaB₂ ^4,5)—or via shear-induced transformations (e.g., TaB₂ and WB₂ ⁶). TiB₂ has been the most widely researched diboride, with hardness of 30–50 GPa and elastic modulus of 450–600 GPa.^7–11 This significant spread of reported mechanical properties is seen not only in TiB₂ but in many other diborides,^12–14 reflecting various synthesis conditions which promote deviations from the ideal 1 : 2 TM-to-boron ratios.^15,16 However, experimentally resolving the effects of individual defect types on mechanical response is highly challenging, owing to the difficulty of isolating specific defects (e.g., boron vacancies and planar defects^12,17–20) and the limited control over synthesis conditions.²¹

Atomistic simulations, in principle, could shed light on the relationships between specific compositional defects and mechanical properties of TMB₂s, guiding the synthesis, micromechanical testing, and electron microscopy investigations. Density functional theory (DFT) calculations have provided basic insights, including predictions of 0 K elastic constants.²² Length scales accessible to DFT, however, do not allow for modeling more complex defects (being typically restricted to simple point defects) and loading scenarios. In particular, simulations of nanoindentation—the state-of-the-art tool for hardness measurements in TiB₂ and ceramic thin films in general—are unfeasible for DFT. A suitable avenue for such simulations is molecular dynamics (MD), with the main limitation being the availability of reliable interatomic potentials. Only few empirical interatomic potentials exist for transition metal nitride/boride based ceramics, such as embedded atom method (EAM) potential for Ti–Al–N,²³ MEAM (modified EAM) potential for Ti–N system with various stoichiometry,²⁴ and 2NN-MEAM (second nearest-neighbor MEAM) potential for Ti–B.²⁵ The last mentioned has been used to predict nano-hardness of defect-free TiB₂ as a function of temperature (T = 0–1500 K) via nanoindentation simulations yielding 70 GPa at 0 K and 65 GPa at 300 K (≈10–20% overestimated compared with experiment). The rapidly developing filed of machine learning interatomic potentials (MLIPs)—promising DFT-level of accuracy and computational efficiency of (semi-)empirical potentials—has boosted MD simulations which are now becoming feasible for more materials of practical relevance.^26–28 Recently, several MLIP training strategies suitable for nano-scale mechanical simulations of ceramics have been proposed,^29,30 providing a basis for our work. While showing promise for reliable predictions of deformation and fracture mechanisms as well as quantities such as theoretical tensile strength and fracture toughness,⁶ these strategies have only been applied to ideal defect-free structures and need to be further tested/extended for simulations including highly defective atomic environments.

According to experimental findings, mechanical properties of TiB₂ thin films vary significantly with B stoichiometry. For instance, the hardness of DC magnetron sputtered TiB_2−x films (with a B/Ti ratios 1.44–2.06) was reported to decrease from 43.9 GPa to 40.7 GPa, while the elastic modulus increased from 511 GPa to 541 GPa.^15,31 Despite a substantial lack of boron, the hardness and Young's modulus are less affected. The mechanisms driving this phenomenon require detailed exploration. Based on the Ti–B phase diagram, some studies have treated TiB₂ as a line compound,³² suggesting that under B-deficient conditions, the lack of B can be offset by formation of Ti-rich planar defects, locally alternating the crystalline structure and dislocation mobility. DFT calculated formation energy by Dahlqvist et al.²² suggested that B vacancies are energetically unfavorable, yielding ≈0.3 eV at⁻¹ differences from the ideal stoichiometric TiB₂ as B/Ti ratio vary from 1.5 to 2, whereas the formation of locally Ti-rich planar defects along the prismatic planes of TiB_2−x lattices is energetically more favorable than the formation of vacancies. These indications have been supported by transmission electron microscopy (TEM) investigations on various synthesized TiB_2−x thin films, showing that unpaired Ti atoms—resulting in B sub-stoichiometry—create planar defects on the prismatic plane of the TiB₂.^1,17

To elucidate how planar defects and randomly distributed boron vacancies impact the mechanical properties of boron-deficient TiB_2−x systems, in this work, we predict room-temperature mechanical properties of single-crystal TiB_2−x with various B contents realized via either point or planar defects, which were previously suggested to be largely present in magnetron-sputtered TiB₂ coatings.^17,18,33 Specifically, the elastic constants, hardness and indentation moduli are derived from nanoscale MD simulations of uniaxial tensile loading and nanoindentation. The underlying interatomic potential is machine-learned on various strained and/or defective ab initio bulk structures (active learning on configurations), up-fitted using highly extrapolative clusters representing environments near the indenter tip (on-the-fly training), and validated against ab initio calculations, and experimental data. The calculated mechanical properties as well as phenomena such as pop-in effects during nanoindentation align well with available experimental data. Further, the experimentally observed—and somewhat surprising—minor changes in hardness between the (nearly) ideal stoichiometric TiB₂ and largely B-deficient TiB_2−x, x ≈ 0.5, are attributed to Ti-rich planar defects. In a stark contrast, B vacancies cause much more significant hardness drop.

2 Methods

2.1 Ab initio datasets

Ab initio datasets for initial MLIP fitting, validation, and active learning were generated employing the Vienna ab initio Simulation Package (VASP^34–36), together with the projector augmented wave (PAW) method and the Perdew–Burke–Ernzerhof exchange–correlation functional revised for solids.^34,37 These included 0 K density functional theory (DFT) and finite-temperature Born-Oppenheimer Ab initio MD (AIMD) calculations³⁸ using the Machine-Learning-Force-Field (MLFF) feature embedded in VASP. The on-the-fly learning in each AIMD run was initiated from scratch (with default parameters) and served only to accelerate generation of atomic environments for own MLIP training detailed in the following subsection. Testing the accuracy of total energies (forces, stresses, structural parameters) during VASP MLFF-aided runs against equivalent fully ab initio simulations yielded plausible accuracy. For illustration, the Root-Mean-Square-Error (RMSE) of total energies, atomic forces, and stresses were ≈1 meV per atom, 0.06 eV Å⁻¹, 0.24 GPa, respectively. Regarding other parameters of our Ab initio calculations, the plane wave cut-off energy of 300 eV was applied together with the Γ point sampling of the reciprocal space.

The model of various TiB_2−x structures was based on a 720-atom (ideal α-TiB₂, space group P6/mmm, Material ID:mp-1145)³⁹ supercell, equilibrated at 300 K during simulations with the NpT (Isothermal–isobaric ensemble) ensemble (3 ps), followed by a simulation with the NVT (Canonical ensemble) ensemble (2.5 ps, using time-averaged lattice parameters from NpT).

To represent various off-stoichiometric environments, TiB_2−x structures with (i) randomly-distributed B vacancies, and (ii) uniformly distributed Ti-rich planar defects were simulated (2-step NpT → NVT equilibration at 300 K, similar as for TiB₂). For MLIP training, we included 50% and 25% vacancy concentrations in the point defect structures. Ti-rich planar defects were constructed based on previous experimental observations.^17,18 Specifically, these defects were created by removing boron atoms located between selected single or double (101̄0) prismatic Ti planes. These defect planes were oriented perpendicular to the [101̄0] crystallographic direction and uniformly distributed across the models. To generate environments representing various mechanical loading conditions, surfaces, and voids, tensile strains (2–10%; 30% “shock” tension inducing cleavage) and compressive strains (2% to 10%) were applied to selected TiB_2−x structures.

2.2 Machine-learning interatomic potential (MLIP) training

The mlip-3 package⁴⁰—implementing the moment tensor potential (MTP) formalism—was used to train our MLIP, targeted to nanoscale mechanical simulations of (defective) AlB₂-type TiB_2−x structures. The weights of energy and forces were set to 1 and 0.01, respectively, while stresses were not fitted due to the presence of atomic clusters in the training set. Our training procedure (Fig. 1) consisted of four main steps detailed below.

Fig. 1 Schematic of our MLIP training routine. (1) Generation of initial training, learning, and validation set by finite-temperature ab initio MD calculations (point and planar defects, various strained structures). (2) Iterative active learning (on AIMD configurations) and validation using the concept of extrapolation grade, γ. (3) Preliminary nanoindentation runs to extract locally highly extrapolative environments (γ > γ_thr) near the indenter tip, to be labeled by ab initio energies and forces in additional DFT calculations; (4) further MLIP re-training using atomic clusters from (3). Note that (3) and (4) represent on-the-fly training on environments near the indenter tip. In the schematic of atomic configurations gray and green spheres represent Ti and B atoms, respectively.

1. Generation of AIMD dataset (Fig. 1, Step 1). AIMD simulations (detailed in Sec. 2.1) were performed to generate a pool of ab initio configurations to be randomly split into the initial training set (TS₀, 5%), the validation set (VS, 20%), and the learning set (LS, the rest ≈75%).

2. Active learning on AIMD configurations (Fig. 1, Step 2). The initial MLIP₀ was trained on TS₀ using an untrained 16-level MTP. Consequently, the MaxVol algorithm^41,42 was employed to identify extrapolative configurations in the LS to be added to the TS based on their extrapolation grade (γ). In simple words, γ is the “uncertainty indicator” quantifying the degree to which an atomic configuration deviates from local environments in the training set. Higher γ indicates higher extrapolation into unfamiliar configuration space where energy (interatomic force, stress) predictions may no longer be reliable. From the definition (eqn (9) in ref. 42), γ ≤ 1 and γ > 1 imply MLIP's interpolation and extrapolation, respectively, and γ_thr = 2 has been previously used as accurate extrapolation.^29,42,43 In this stage, a threshold γ_thr = 2 is set for active learning. The final MLIP produced in Step 2 exhibited γ ≤ 1 for the TS, γ < 2 for the LS but also the—independent, unseen—VS. Training and validation errors for energies and forces, as quantified by Root-Mean-Square-Error (RMSE), were below 6 meV per atom and 0.25 eV Å⁻¹.

3. Preliminary nanoindentation simulations to identify extrapolative environments near the indenter tip (Fig. 1, Step 3). Using the MLIP from Step 2, we simulate nanoindentation tests (for computational details see Sec. 2.2.1) and output the per-atom γ value at each step. These allow identifying locally highly extrapolative environments (γ ≫ 10 in a spherical atomic region with a 9 Å radius). The clusters—embedded in a box with a 10 Å vacuum buffer—were labeled by ab initio energies and forces in additional DFT calculations (see Sec. 2.1).

4. MLIP refinement using extrapolative atomic clusters (Fig. 1, Step 4). Atomic clusters generated above were added to the TS and the MLIP was re-trained until no highly extrapolative environments (γ > 10) were present during MD nanoindentation test for both the ideal stoichiometric TiB₂ as well as B-substoichiometric structures (on-the-fly training). MLIP's transferability to various TiB_2−x structural variants (some of which were directly trained on already in Step 2) is further discussed in the Results section.

2.2.1 Molecular dynamics (MD) simulations. MD simulations were performed using the Large-scale Atomic/Molecular Massively Parallel Simulator (LAMMPS)^44,45 interfaced with the mlip-3 package (requiring custom compilation) to enable the use of MTP-type machine-learning force fields. The force field was trained by us, as detailed in the previous section. All simulations used the extrapolation control (γ, explained in the previous section) to monitor the reliability of the underlying potential. Specifically, γ = 10 was set for terminating the simulation, but this critical value has never been reached during mechanical tests presented in the Results section. The timestep was always set to 1 fs. The Open Visualization Tool software (OVITO⁴⁶) served to visualize and analyze the results.
TiB_2−x models. The ideal AlB₂-type α-TiB₂ was modeled in a 120 000-atom supercell (≈11 × 11 × 11 nm³). B-substoichiometric structures were considered via (i) randomly distributed B vacancies, or (ii) uniformly distributed planar defects, specifically, Ti-rich planar defects observed in experiment^17,18 and visualized in Fig. 2. All structures were fully equilibrated at 300 K during a 10 ps NpT simulation.

Fig. 2 Illustration of our B-deficient TiB_2−x models. (a) The ideal stoichiometric TiB₂, (b) TiB_2−x, x = {0.2; 0.3; 0.4; 0.5}, with randomly distributed B vacancies also referred to as Vac; (c) TiB_2−x, x = {0.67}, with single Ti-rich planar defect also referred to as SP; (d) TiB_2−x, x = {0.67}, with double Ti-rich planar defect also referred to as DP. The visualization is at the (101̄0) plane (top row) and the basal (0001) plane (bottom row), respectively, generated using the program VESTA.⁴⁷

Mechanical tests with fully periodic boundary conditions. For selected systems (TiB₂, TiB_1.8, TiB_1.7, TiB_1.6, and TiB_1.5), room-temperature uniaxial tensile and shear tests were simulated by applying strain along the [0001], [101̄0], and [1̄21̄0] direction or along the (0001)[101̄0] slip system, respectively. The strain rate was 0.5 m s⁻¹ and the simulation was continued until fracture (tensile test) or lattice slip (shear test). Tensile tests considered Poisson's contraction (via the NpT ensemble orthogonal to the loading direction), while shear test were conducted at a constant volume. Stress/strain data from the elastic region (strain ≤4%) served to derive the second-order elastic constants, C_ij. Based on C_ijs, we evaluated elastic moduli (e.g., the bulk, Young's, and shear modulus) following standard formulas.⁴⁸
Nanoindentation simulations for MLIP training (detailed in Sec. 2.2). To identify environments representing highly strained regions near the indenter tip in an efficient manner, nanoindentation tests were performed in medium-sized supercells (≈6000 atoms, 3.5 × 3.5 × 3.5 nm³) at a high loading rate (≈200 m s⁻¹), considering a spherical indenter with a 1 nm radius. If present, highly extrapolative environments (with γ ≫ 10) were extracted from each simulation step and used to re-train our MLIP.
Nanoindentation simulations for quantification of mechanical properties. The model (Fig. 3) included a TiB_2−x supercell (≈11 × 11 × 11 nm³) together with a virtual repulsive indenter. Periodic boundary conditions (PBCs) were considered along the x‖[1̄21̄0] and y‖[101̄0] directions, while z‖[0001] was non-periodic: with lower and upper boundary being fixed and shrink-wrapped, respectively, to ensure that all atoms remain within the simulation box.

Fig. 3 Schematic of our nanoindentation model. Atoms in the simulation box were grouped into (i) a 7 Å thick rigid layer with frozen atoms to prevent movement along the z axis, (ii) a 7 Å thick thermostat layer in which the temperatures fixed at 300 K using temperature scaling thermostat, and (iii) a Newtonian layer.

Prior to nanoindentation, the simulation box was equilibrated at the target temperature (T = 300 K for most simulations) for 10 ps using NpT ensemble, followed by a 2.5 ps simulation with in NVT ensemble with Nosé–Hoover thermostat. Nanoindentation was simulated in a displacement-controlled mode, by moving a virtual indenter (3.5 nm radius, force constant K = 10 eV Å⁻³) along the [0001] axis for 140 ps at loading speed of 15 m s⁻¹, until reaching ≈20 Å indentation depth. The hardness was evaluated from the hardness result averaged at the indentation depth of (15 ± 0.5) Å, corresponding to approximately 15% of the film thickness, in order to minimize the influence of the infinitive hard substrate on the hardness measurement. After holding the indenter for 2 ps, unloading was simulated for another 100 ps.

The hardness, H, was evaluated following Oliver–Pharr method,^21,49 i.e., from the relationship between the indentation load, P, and the projected contact area, A, during indentation process,

where R is the (spherical) indenter radius and h_c is the current projected contact depth. The reduced elastic modulus, E_r, was derived from the unload phase aswhere S = dP/dh is the stiffness (extracted by fitting a linear function using the linear elastic interval of load–displacement data during the unload) and β = 1 is the correction factor. Elastic deformation is generally observed in both indented material and indenter tip during the indentation process. The relationship between the reduced elastic modulus (E_r), the elastic modulus of the indented material (E_s), and that of the indenter tip (E_i) is given by:where ν is the Poisson's ratio of the indented material and ν_i the indenter. For an ideal virtual infinitely rigid indenter, E_i = +∞, thus, eqn (3) simplifies as⁵⁰

In this work, Poisson's ratios for the investigated material systems were obtained from elastic constants, C_ijs, calculations and are listed in Table 2.

3 Results and discussion

3.1 Validation of the trained MLIP

To carry out nanoscale MD nanoindentation tests for various defective TiB_2−x structures, we trained machine learning interatomic potential (MLIP) based on moment tensor potentials (MTP) formalism, building on previous training routine developed for tensile and shear simulations of single-crystal ceramics.^29,43 Enabling accurate description of complex, highly strained regions near the indenter tip, this routine was essentially upgraded by on-the-fly training on extrapolative atomic clusters from simplified nanoindentation models using a pre-trained MLIP (for details see the Methods and Fig. 1). Prior to employing the trained MLIP, we show its capability to reproduce material properties consistently with 0 K and finite-temperature ab initio results, as well as available experimental data.

Validation against an independent ab initio dataset—that contains ≈10⁴ snapshots from finite-temperature AIMD simulations of TiB_2−x featuring a wide range of point and planar defect types/concentrations/distributions as well as these structures under tensile, shear, or compressive loads—supports the credibility of our MLIP (Fig. 4). Fig. 4a depicts the histogram of total energies in the reference set, while the corresponding parity plot is given in Fig. 4b. The AIMD vs. MLIP-predicted energies fall at or near the x = y line (Fig. 4b), showing that the model closely reproduces ab initio values, as quantified also by statistical errors RMSE = 6 meV and R² = 0.999. Additional comparison with the ab initio MD dataset (300 and 1200 K data) from our previous work²⁹ yields RMSE and R² values of 13 meV and 0.999, respectively.

Fig. 4 Validation of the here-trained MLIP against a reference finite-temperature ab initio dataset (≈10⁴ configurations). (a) The total energy distribution in the reference set containing various defective and/or highly-strained TiB_2−x configurations, as described in the Methods. The B off-stoichiometry is realized via randomly distributed vacancies or planar defects (Vac, SP, DP–see Fig. 2). (b) Linear correlation between the corresponding AIMD and MLIP-predicted energies (blue data points) together with values predicted for another reference dataset from ref. 29 containing snapshots from finite-temperature (300 and 1200 K) AIMD simulations of TiB₂ subject to uniaxial tensile and shear loading until fracture or lattice slip.

Furthermore, our MLIP accurately reproduces stress/strain curves and failure mechanisms derived from room-temperature tensile and shear tests for defect-free TiB₂ in ref. 29, in particular, yielding accuracy in stresses of <1 GPa where slightly higher deviations from AIMD data may occur near the fracture point. To assess the MLIP's transferability beyond DFT accessible length scales and point/planar defect concentrations and configurations not directly trained on, we simulate room-temperature equilibration of TiB₂ as well as newly generated TiB_2−x structures (e.g., TiB_1.6, TiB_1.8), employing supercells with ≈11 000 atoms (i.e., beyond DFT accessible length scales). Comparison of the calculated lattice parameters (Table 1) with reference ab initio or experimental values shows excellent agreement. Specifically, the differences are 0.02 Å (0.6%) and 0.006 eV at⁻¹ (0.07%) in terms of lattice parameters and energies, respectively. Importantly, the MLIP's extrapolation grade during all equilibrations remains within the reliable extrapolation regime (γ < 10 (ref. 29 and 51)).

Table 1 Validation of the here-trained MLIP in terms of structural parameters of various TiB_2−x structures. Comparison of time-averaged lattice parameters (a and c in Å) obtained by MLIP-MD and equivalent ab initio MD simulations at temperature T = 300 K, or experimental values. Defects are categorized as vacancies (Vac) and planar defects (SP, DP), as schematically depicted in Fig. 2

Structure	Defect type	Method	T	a	c
	—	ab initio - MD	300	3.037	3.218
		MLIP-MD	300	3.040	3.236
		DFT^29	0	3.027	3.213
TiB2		DFT^52	0	3.028	3.226
		Experiment^53	300	3.030	3.230
		Experiment^54	300	3.030	3.240
TiB1.8	Vac	ab initio – MD	300	3.048	3.183
		MLIP-MD	300	3.049	3.197
TiB1.7	Vac	ab initio – MD	300	3.055	3.155
		MLIP-MD	300	3.054	3.179
TiB1.6	Vac	ab initio – MD	300	3.058	3.139
		MLIP-MD	300	3.056	3.154
TiB1.5	Vac	ab initio – MD	300	3.059	3.137
		MLIP-MD	300	3.061	3.132
TiB1.33	SP	ab initio – MD	300	3.046	3.259
		MLIP-MD	300	3.056	3.261
TiB1.33	DP	ab initio – MD	300	3.034	3.236
		MLIP-MD	300	3.055	3.251

---

Using the here-trained MLIP, we have additionally calculated 0 K formation energies, E_f, of representative non-stoichiometric TiB_2−x structures (Fig. 2 and Fig. 5). The results indicate that all defects cost energy (ΔE_f > 0 considering TiB₂ as a reference state) and TiB_2−x variants with disordered vacancies are significantly more expensive than their chemically equivalent counterparts containing single- or double-planar defects (SP, DP). These predictions are consistent with previous 0 K DFT indications²² and experimental (electron microscopy) findings.^17,18 The defects considered here are those frequently mentioned in experimental investigations of diborides. While there are other ways to realize B substoichiometry—such as Ti substitutional defects, Ti interstitials, or planar defects at different planes—none of those has been evidenced in experiment is therefore excluded from our models.

Fig. 5 Formation energy difference (ΔE_f) between various B-substoichiometric TiB_2−x (x = 0.2, 0.3, 0.4, 0.5, 0.66) structures calculated at 0 K using the here-trained MLIP. Defects are categorized as vacancies (green square), planar defects (SP: blue cross, DP: yellow cross), as schematically depicted in Fig. 2. Note that the MLIP-predicted E_f values differ from own reference DFT calculations by less than 0.01 eV at⁻¹, underpinning the credibility of our force field.

In what follows we present MD simulations—particularly nanoindentation tests for various off-stoichiometric TiB_2−x systems—that are beyond length scales accessible to ab initio calculations. Their reliability is inferred from (i) the previously presented MLIP's validation against ab initio data, and (ii) MLIP's extrapolation grades, γ (see the Methods), monitored during some of these simulations. Importantly, γ values indicate reliable extrapolation (γ < 10) even under far-from-equilibrium conditions characterized by high stresses, local distortions, Ti/B intermixing, and crack nucleation beneath the indenter tip. MLIP's reliability at the nanoscale is further supported by comparing the calculated material's properties, such as elastic constants and hardness, with reference DFT or experimental data (see following sections).

3.2 Elastic constants, tensile and shear strength of defective TiB_2−x

Employing the here-trained MLIP, we simulate tensile and shear deformation of various defective TiB_2−x structures, considering low-index crystallographic directions and slip systems. The simulations are carried out at 300 K, for supercells with ≈110 000–150 000 atoms, under fully periodic boundary conditions. Prior to deformation, all structures are fully relaxed, exhibiting excellent thermal stability. In particular, no B diffusion is observed even for largely off-stoichiometric models (e.g. TiB_1.6), mirroring the intrinsic structural stability of the α phase which imposes high energy barrier for vacancy migration and diffusion. According to the Arrhenius relationship, this results in diffusion at 300 K unlikely, making macroscopic observation virtually impossible.^55–57

The linearly-elastic part of stress/strain curves (Fig. S1 in SI) serves to derive second-order elastic constants, C_ij, and related mechanical properties (Table 2, Fig. 6). Furthermore, global stress maxima in the stress/strain data are used as estimates of the theoretical tensile and shear strength (Table S1 in SI) following ref. 58 and 59.

Fig. 6 Room-temperature directional Young's moduli, polycrystalline Young's modulus, bulk modulus and shear modulus of the ideal TiB₂ as well as of various B-substoichiometric TiB_2−x structures derived using the here-trained MLIP. (a) in-plane Young's modulus E_[1̄01̄0] = E_[1̄21̄0]; (b) out-of-plane Young's modulus E_[0001]; (c) bulk modulus, B; (d) Young's modulus, E; and (e) shear modulus, G were obtained from the elastic constants, C_ij. The B deficiency in TiB_2−x was realized by uniformly distributed single- or double-planar defects (SP: blue data points, DP: orange data points), or disordered vacancies (Vac: green data points), as schematically depicted in Fig. 2. The dashed horizontal line marks the ideal TiB₂. For the underlying elastic constants and numerical values of the elastic moduli, see Table 2.

Table 2 Room-temperature mechanical properties of the ideal TiB₂ as well as of various B-substoichiometric TiB_2−x structures derived from uniaxial tensile and shear tests at 300 K using the here-trained MLIP. Second-order elastic constants, C_ij (in GPa), were obtained from the linearly-elastic part of the stress/strain curves, and were subsequently used to calculate the polycrystalline bulk, shear, and Young's modulus (B, G, E in GPa) and the Poisson's ratio (ν). The B deficiency in TiB_2−x was realized by uniformly distributed single- or double-planar defects (SP, DP), or disordered vacancies (Vac), as schematically depicted in Fig. 2

Models	Defect type	T	C11	C12	C13	C33	C44	B	G	E	ν
TiB1.5	SP	300	461	86	100	336	156	206	165	392	0.183
	DP		466	78	98	371	183	205	179	417	0.161
	Vac		319	73	85	291	113	157	115	278	0.206
TiB1.6	SP	300	478	83	94	381	217	208	195	447	0.142
	DP		478	74	89	377	204	203	193	440	0.140
	Vac		357	70	78	319	143	165	139	326	0.171
TiB1.7	SP	300	503	86	103	397	224	220	204	467	0.146
	DP		497	79	101	392	205	216	196	452	0.151
	Vac		395	79	92	334	162	184	153	361	0.174
TiB1.8	SP	300	510	84	92	401	236	216	212	480	0.131
	DP		506	67	90	397	225	211	210	473	0.126
	Vac		434	76	83	366	191	191	177	406	0.145
TiB2	—	300	542	89	105	421	246	233	222	505	0.138
TiB2 ^48	—	0	670	64	101	473	267	263	269	601	0.119
TiB2 ^29	—	300	588	85	98	409	261	236	246	554	0.113
TiB2 ^60	—	300	654	56	98	454	263	250	262	583	0.111

---

Fig. 6 underscore the strong influence of the B content on TiB₂'s mechanical properties already in the linearly-elastic regime. Specifically, the in-plane Young's modulus, E_[1̄21̄0] (Fig. 6a), decreases from 507 GPa to 426 GPa (16% decrease) when going from TiB₂ to TiB_1.5, where the latter contains planar defects (SP, DP). Due to the prismatic plane isotropy, differences between E_[1̄21̄0] and E_[101̄0] are negligible. The out-of-plane Young's modulus, E_[0001] (Fig. 6b), drops from 386 GPa (TiB₂) down to 330 GPa (TiB_1.5, 15% decrease) with planar defects. While the B-sub-stoichiometry-mediated decrease in directional Young's moduli is slightly more pronounced for single planar defects (SP), the differences from double planar defects (DP) remain rather small (<10 GPa). Compared with planar defects, the directional Young's modulus of TiB_2−x containing disordered vacancies (Vac) is significantly lower. Trends in the polycrystalline B, G, and E (Fig. 6c–e) follow these observed for the directional Young's modulus, in particular, decrease as the B sub-lattice becomes less populated. This decrease is the least pronounced for the bulk modulus.

The polycrystalline Poisson's ratio (ν, shown in Table 2) increases from 0.14 (TiB₂) up to 0.16–0.2 (TiB_1.5 with planar defects or vacancies), indicating an improved ability to redistribute stresses via plastic deformation, as discussed later in relation to the results obtained from nanoindentation simulation. A Poisson's ratio (ν) of 0.25 has been previously assumed when deriving nanoindentation hardness in experimental studies.⁶¹ This may, however, not accurately reflect the actual TiB₂‘s stoichiometry, causing discrepancies in the literature. When deriving Young's modulus from our later presented nanoindentation tests, we therefore use Poisson's ratios representing the correct stoichiometry and defect type (i.e., SP, DP, or Vac). Using a constant ν = 0.25, the values would change by up to 5%.

Furthermore, Fig. 7 depicts trends in the ideal strengths, σ_hkil, of TiB_2−x, derived from the maximum stress withstood during our room-temperature uniaxial tensile simulations. Irrespective of the defect type (SP, DP, Vac), our predictions reveal σ_[101̄0] < σ_[0001] < σ_[1̄21̄0], suggesting that the basal plane isotropy does not extend beyond linearly-elastic regime. For the corresponding fracture strains, ε^F_[hkil], we predict ε^F_[101̄0] < ε^F_[1̄21̄0] < ε^F_[0001]. Although theoretical tensile strengths generally diminish with decreasing the B content, the trends seem more conspicuous than those derived from elastic constants. In particular, TiB_1.7 with single planar defects exhibits the second highest strengths—after the ideal stoichiometric TiB₂—followed by TiB_1.7 with double planar defects, which, however is notably below, by >20%. Increasing the fraction of double planar defects, tensile strength of TiB_1.5 remains essentially unchanged compared with TiB_1.7, while increasing the fraction of single planar defects, TiB_1.5, yields the overall lowest strength.

Fig. 7 Theoretical tensile strength (and the corresponding fracture strain) of the ideal TiB₂ as well as of various B-substoichiometric TiB_2−x structures derived using the here-trained MLIP. (a), (b) and (c) depict the tensile strength (σ_max; triangle data points, left y-axis) and fracture strain (σ_max; square data points, right y-axis) along the [0001], [101̄0] and [1̄21̄0] crystallographic directions, respectively. The B deficiency in TiB_2−x was realized by uniformly distributed single- or double-planar defects (SP, DP), or disordered vacancies (Vac), as schematically depicted in Fig. 2. The σ_max values of SP, DP, Vac are shown in blue, orange and green lines, respectively. The ε_max of SP, DP and Vac are shown in cyan, magenta and brown lines, respectively. The structure snapshots (color code based on the volumetric strain) show TiB_1.8 with SP near the fracture point, illustrating the fracture pattern. These general fracture patterns are similar for other stoichiometries and defect types.

The results in Fig. 7 also suggest that both types of planar defects deteriorate tensile strength parallel to the Ti-rich defect region. Analysis of [101̄0] tensile tests reveals cleavage along the entire Ti-rich plane (structure snapshots in Fig. 7). Possible explanation is that the Ti–Ti bond strength within the Ti-rich region, perpendicular to the planar defect, is weaker than the intact hexagonal structure with not only Ti–Ti bonds but also Ti–B bonds, in stoichiometric TiB₂. This differences in distribution of mechanical properties throughout the simulation box increases stress concentration at the Ti-rich plane, causing cleavage fracture along these planes.

Furthermore, considering the typical orientation of diboride thin films—with (0001) planes perpendicular to the growth direction⁶²—we simulate shearing along the (0001)[101̄0] slip system. These are also representative of local stress states during nanoindentation. Fig. 8 presents the theoretical shear strength (together with the corresponding slip strain) derived as the maximum stress withstood by the material during shear deformation. The ideal TiB₂ exhibits the highest shear strength but only an intermediate fracture strain compared with the B-deficient systems. Consistently with the shear modulus trends in Fig. 6, TiB_2−x with disordered vacancies show a linear shear strength increase from ≈10 to 20 GPa. The correlation of the corresponding slip strain and B content, however, is rather weak. For single-planar defects, the shear strength initially increases with the B content, reaching a plateau for TiB_1.7. The corresponding slip strain decreases significantly, as the B/Ti ratio decrease. For double-planar defects, the shear strength is the highest and remains stable despite stoichiometry changes, which can be attributed to a more favorable Ti–B bonding environment and improved p–d orbital hybridization, as shown in the SI (Fig. S3). The defective layers also help distribute the applied shear over a wider region, reducing local stress concentration. As a result, the reduction in slip strain is more moderate compared with the SP case.

Fig. 8 Room-temperature shear strength of TiB_2−x systems derived from (0001)[101̄0] shear tests using the here-trained MLIP. (a) Shear strength and (b) the corresponding slip strain. Results for the ideal TiB₂ and defective TiB_2−x structures—single-planar defects, double-planar defects, and randomly distributed B vacancies (schematically depicted in Fig. 2)—are marked in blue, orange and green, respectively. The underlying values are shown in Table S1 in SI.

Overall, simulations of tensile and shear tests indicate that the mechanical properties of TiB_2−x are controlled by both the B content and the underlying defect type. In some experimental observations, however, sub-stoichiometric TiB_2−x thin films—somewhat surprisingly—were reported to possess hardness as high as the stoichiometric TiB₂.^31,33 To further investigate this discrepancy, nanoindentation simulations are carried out in the following section.

3.3 Nanoindentation of defective TiB_2−x

To shed light on mechanisms driving the variation in hardness and Young's modulus of perfect TiB₂ and sub-stoichiometric TiB_2−x thin films, room-temperature MLIP-MD nanoindentation simulations for both ideal stoichiometric as well as B-sub-stoichiometric structures were performed. In line with most TiB₂ (and other transition metal diboride) films being preferentially [0001] oriented during the synthesis,^1,2,63 our model considers an indenter moving along the [0001] crystallographic axis.³ Detailed description of the simulation setup and derivation of the corresponding hardness and Young's modulus values is given in the Methods.

Fig. 9 presents the obtained load-depth or load–displacement (L–D) curves for the ideal TiB₂ as well as its defective TiB_2−x counterparts, represented by TiB_1.5, in which the B-sub-stoichiometry is realized by planar defects (SP, DP) or B vacancies (Vac). Irrespective of the B content, pressing the indenter into the (0001) surface initially induces an elastic response. As the indentation depth (h_c) increases, the reactive force on the indenter exhibits a small bump indicating the “pop-in” effect, hence hinting at the onset of plastic deformation via dislocation motion. These observations closely align with L–D curves obtained from lab-scale nanoindentation tests.^64,65 The pop-in event observed in the load–displacement curve is primarily due to the sudden activation of specific slip planes (mainly basal and pyramidal planes) beneath the indenter tip.⁶⁶ These “pop-in” effects seem to reflect the B vacancy distribution, in particular, are notable for structures containing planar defects—which can be viewed as highly ordered vacancy arrangement—but almost negligible for those with vacancies randomly distributed in the B network.

Fig. 9 Representative load–depth (L–D) curves derived from room-temperature nanoindentation simulations for TiB₂ as well as various B-substoichiometric TiB_2−x structures. These nanoscale simulations were carried out using the here-trained MLIP. The B deficiency in TiB_2−x was realized by uniformly distributed single- or double-planar defects (SP: orange data points, DP: blue data points), or disordered vacancies (Vac: green data points), as schematically depicted in Fig. 2.

The calculated L–D data was further used to determine the hardness, H, and Young's modulus, E (Fig. 10). With H = 65 GPa and E = 563 GPa, TiB₂ ranks as the second highest and the highest in terms of H and E, respectively, which is ≈10–20% overestimated with respect to experiment. Such discrepancy can be ascribed to the fact that real samples are never perfectly pure and ideally defect-free. Besides simple point and planar defects considered in our work, H(E) measurements are influenced by other macroscopic features like grain boundaries, micro-cracks as well as porosity arising from various deposition conditions. As previously shown^67–71 the yield strength, thus hardness of ceramics, is dependent on the grain size following the Hall–Petch relationship. Quantitative assessment of these effects for TiB₂, however, is outside of our scope.

Fig. 10 The hardness (a) and Young's modulus (b) derived from room-temperature nanoindentation simulations using the here-trained MLIP. The results for TiB₂ (the dashed horizontal line) are compared with various B-substoichiometric TiB_2−x structures. The SP and DP refer to single and double-planar defects, respectively, as schematically depicted in Fig. 2.

While both the predicted H and E decrease upon decreasing the B content, i.e., by introducing more defects, there is an apparent difference between effects of B vacancies and planar defects, where the former cause a more dramatic hardness and Young's modulus deterioration. For structures with B vacancies, H and E decrease steadily by 49% and 40%, respectively, when the Ti-to-B ratio changes as 1 : 2 → 1 : 1.5, i.e., when 25% of the B sub-lattice becomes unpopulated. The H and E decrease in TiB_2−x with planar defects occurs at a significantly slower rate. In the presence of single-planar defects, H = 59–62 GPa (TiB_1.5 → TiB_1.8), which is only a 10% to 5% decrease compared with the ideal TiB₂, thus, the effect of decreasing the B content massively—down to TiB_1.5 —causes only minimal H changes. For TiB_2−x with double-planar defects, H = 66 GPa of TiB_1.7 even slightly exceeds that of the defect-free TiB₂.

Our simulations (Fig. 10b) further indicate that the Young's modulus E = 563 GPa of TiB₂ decreases to 338 GPa (B vacancies) or 472 GPa (SP) or 422 GPa (DP) for the TiB_1.5 stoichiometry. The differences between SP and DP defects are generally minor, especially compared with the significant E decrease in vacancy-containing TiB_2−x structures. Note that B vacancies are more detrimental not only to TiB₂'s mechanical performance but also its energetic stability, as suggested by us as well as by previous DFT formation energy calculations.²² These simulations help to understand recent experimental findings³³ demonstrating that the hardness of TiB_2−x thin films deposited via Direct Current Magnetron Sputtering (DCMS) is maintained or slightly enhanced upon reducing the B/Ti ratio from 2.12 to 1.46, whereas the elastic modulus exhibits a noticeable decline. In particular, according to our predictions, the hardness remains nearly constant or undergoes a slight reduction as the B/Ti ratio decreases from 2 to 1.7 and eventually to 1.5. Conversely, the Young's modulus decreases significantly as the B/Ti ratio varies from 2 to 1.5, particularly in models with SP and DP defects.

To develop basic level of understanding why disordered B vacancies and Ti-rich planar defects trigger a distinct mechanical response, Fig. 11 shows the volumetric shear strain distribution in the indented models in which the indenter is pressed fairly deep in the sample (indentation depth ≈2 nm). The analysis indicates qualitatively different stress redistribution due to the presence of planar defects compared with disordered B vacancies as well as the ideal TiB₂. In particular, TiB_2−x structure with B vacancies reveals highly concentrated shear strain below the indenter tip and fewer slip events, in contrast to the defect-free TiB₂ as well as TiB_2−x featuring Ti-rich planar defects. This may diminish the occurrence of “pop-in” effects in models with B vacancies. Unlike that, TiB_2−x with planar defects (both SP and DP) restricts lattice slip between several planar defects beneath the indentation zone. These indications are consistent with Thörnberg et al.,³¹ suggesting that Ti-rich planar defects act as barriers to dislocation glide, contributing to anisotropic deformation behavior and high hardness.

Fig. 11 Atomic shear strain distribution at the final stage of nanoindentation tests simulated with the here-trained MLIP. (a) Stoichiometric TiB₂, (b–d) off-stoichiometric TiB_1.8 with disordered vacancies (b), single planar defects (c), and double planar defects (d), as schematically depicted in Fig. 2. The color code is mapped using the atomic strain analysis function in OVITO.⁴⁶ The supercell slice has thickness of ≈1 nm, oriented perpendicular to the prismatic plane. The tripod in color (top right) indicates crystallographic directions.

4 Conclusion and outlooks

Employing MLIP-MD simulations, we predicted fundamentally different mechanical response of a paradigm hard coating material, TiB₂, depending on (i) the B content—within the experimentally reported B/Ti ratios of 2 to 1.5—and, (ii) point vs. planar defects realizing the same stoichiometry. The defect types were motivated by other DFT studies (vacancies) and experimental high-resolution electron microscopy investigations (Ti-rich single and double planar defects). Our machine-learning force field (MTP-type MLIP) for finite-temperature nanoindentation simulations was trained by active learning on both stoichiometric and off-stoichiometric TiB_2−x structures subject to tensile, shear, or compressive loads, as well as on-the-fly (OTF), on atomic clusters from simplified nanoindentation tests.

Main conclusions of this work are summarized below:

1. Our MLIP accurately reproduced finite-temperature ab initio energies of AlB₂-type TiB_2−x structures with point or planar defect types and contents not directly trained on (RMSE_max ≈ 6 meV at⁻¹). Further, it showed good agreement with ab initio and—where applicable—experimental lattice parameters (±0.02%) and second-order elastic constants (±10%), as well as reliable extrapolation grades during mechanical tests at both AIMD length scales (uniaxial tensile or shear simulations under fully periodic boundary conditions) and beyond (nanoindentation).

2. Room-temperature tensile tests for supercells with ≈140 000 atoms, uni-axially loaded along [0001], [1̄21̄0], or [101̄0] direction, ascribed the highest theoretical strength (54.5 GPa) to the stoichiometric TiB₂. In B-sub-stoichiometric structures, the strength gradually reduced by 57%, 53% and 28% in case of vacancies, single-planar defects and double-planar defects, respectively, as B content decreased from the ideal TiB₂ to TiB_1.5. Notably, the most significant reduction occurred in the direction normal to the planar defect, associated with brittle cleavage perpendicular to Ti-rich (101̄0) planes. For the shear test along the (0001)[101̄0] slip system, the perfect TiB₂ model exhibited the highest shear strength but only an intermediate slip strain compared with all other tested models.

3. Room-temperature nanoindentation tests—considering commonly reported [0001]-oriented samples—showed a significant hardness and modulus decrease for TiB_2−x with disordered B vacancies: by 48% and 37%, respectively, for TiB_1.5, as compared with the ideal TiB₂. Contrarily, the H decrease for TiB_2−x with Ti-rich planar defects was way less dramatic, about 8%, for TiB_1.5 with single-planar defects. TiB_1.7 with double-planar defects exhibited H = 66 GPa even slightly above H = 65 GPa of the ideal TiB₂. The predicted H values were generally 10%–20% above experimental data, as our models only account for one defect type. Further, structural analysis suggested that planar defects serve as barriers impeding the dislocation glide or slip beneath the indentation region, thereby contributing to the increased H or preventing a rapid H decrease.

Our study indicated that MLIPs transferable to complex stress states near indenter tips require on-the-fly training, additionally to training on bulk and surface structures from ab initio MD. The here-proposed training strategy can be straightforwardly applied to other AlB₂-type diborides and is likely suitable also for transition metal carbides and nitrides (typically crystallizing with different, cubic, crystal symmetry). Further, the wide spread of H data and the fact that largely B-substoichiometric compositions—though less energetically favorable compared with the fully occupied lattice—can even exceed H of the ideal TiB₂ was explained by different mechanical response of TiB_2−x with B vacancies and Ti-rich planar defects. We hypothesize that similar H trends and non-stoichiometry driven deformations mechanisms can be found in diborides energetically favoring the AlB₂ structure, particularly ZrB₂ and HfB₂ isoelectronic to TiB₂.

All in all, the unique capability of atomistic simulations to evaluate effects of individual defect types, for desired chemical compositions, under well defined and comparable loading conditions, was demonstrated. Challenges for follow-up research include modeling more complex (combinations of) defects, such as grain boundaries or amorphous B-rich tissue phase frequently reported in diborides. Following our MLIP training and MD simulation routine, a fundamentally interesting question is whether TM-rich planar defects have similar effects in other TMB₂s and can be used in rational structural design.

Author contributions

CD: conceptualization, methodology, investigation, data curation, visualization, writing – original draft. SL: investigation, validation, writing – review & editing. NK: conceptualization, methodology, writing – review & editing, supervision. PHM: resources, writing – review & editing, supervision.

Conflicts of interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Data availability

The supporting data has been provided as part of the supplementary information (SI). Supplementary information: Table S1, Fig. S1, S2 and S3. See DOI: https://doi.org/10.1039/d6nr00466k. The relevant MLIP and lattice structures are stored in the TU Wien Research Data repository. See DOI: https://doi.org/10.48436/3mat6-bp059.

Acknowledgements

NK and PHM acknowledge the support by the Austrian Science Fund, FWF, via projects 10.55776/PAT4425523 and 10.55776/RIC2714224. The computations handling were enabled by resources provided by the Vienna Scientific Cluster (VSC) in Austria, as well as by National Academic Infrastructure for Supercomputing in Sweden (NAISS) and the Swedish National Infrastructure for Computing (SNIC) at the National Supercomputer Center (NSC) partially funded by the Swedish Research Council through grant agreements no. 2022-06725. The authors acknowledge TU Wien Bibliothek for financial support through its Open Access Funding Programme. Further, the authors thank D.G. Sangiovanni from Linköping University and D. Holec from Montanuniversität Leoben for insightful discussions.

References

1. B. Paul, N. L. Okamoto, M. Kusakari, Z. Chen, K. Kishida, H. Inui and S. Otani, Plastic deformation of single crystals of CrB2, TiB2 and ZrB2 with the hexagonal AlB2 structure, Acta Mater., 2021, 211, 116857.
2. P. H. Mayrhofer, C. Mitterer and H. Clemens, Self-organized nanostructures in hard ceramic coatings, Adv. Eng. Mater., 2005, 7, 1071–1082.
3. M. Magnuson, L. Hultman and H. Högberg, Review of transition-metal diboride thin films, Vacuum, 2022, 196, 110567.
4. T. Fiantok, N. Koutná, D. G. Sangiovanni and M. Mikula, Ceramic transition metal diboride superlattices with improved ductility and fracture toughness screened by ab initio calculations, Sci. Rep., 2023, 13, 12835.
5. T. Leiner, N. Koutná, J. Janovec, M. Zelenỳ, P. H. Mayrhofer and D. Holec, On energetics of allotrope transformations in transition-metal diborides via plane-by-plane shearing, Vacuum, 2023, 215, 112329.
6. S. Lin, D. Holec, D. G. Sangiovanni, T. Leiner, L. Hultman, P. H. Mayrhofer and N. Koutná, Shear-activated phase transformations of diborides via machine-learning potential molecular dynamics, Acta Mater., 2025, 121606.
7. M. Mikula, B. Grančič, V. Buršíková, A. Csuba, M. Držík, Š. Kaveckỳ, A. Plecenik and P. Kúš, Mechanical properties of superhard TiB2 coatings prepared by DC magnetron sputtering, Vacuum, 2007, 82, 278–281.
8. N. Panich and Y. Sun, Mechanical properties of TiB₂ -based nanostructured coatings, Surf. Coat. Technol., 2005, 198, 14–19.
9. M. Berger, L. Karlsson, M. Larsson and S. Hogmark, Low stress TiB2 coatings with improved tribological properties, Thin Solid Films, 2001, 401, 179–186.
10. F. Lofaj, T. Moskalewicz, G. Cempura, M. Mikula, J. Dusza and A. Czyrska-Filemonowicz, Nanohardness and tribological properties of nc-TiB₂ coatings, J. Eur. Ceram. Soc., 2013, 33, 2347–2353.
11. M. Berger, M. Larsson and S. Hogmark, Evaluation of magnetron-sputtered TiB₂ intended for tribological applications, Surf. Coat. Technol., 2000, 124, 253–261.
12. C. Hu, S. Lin, M. Podsednik, S. Mráz, T. Wojcik, A. Limbeck, N. Koutná and P. H. Mayrhofer, Influence of co-sputtering AlB₂ to TaB₂ on stoichiometry of non-reactively sputtered boride thin films, Mater. Res. Lett., 2024, 12, 561–570.
13. A. Goncharov, S. Dub, A. Agulov and V. Petukhov, Structure, composition, and mechanical properties of thin films of transition metals diborides, J. Superhard Mater., 2015, 37, 422–428.
14. C. Fuger, R. Hahn, L. Zauner, T. Wojcik, M. Weiss, A. Limbeck, O. Hunold, P. Polcik and H. Riedl, Anisotropic super-hardness of hexagonal WB_{2 ± z} thin films, Mater. Res. Lett., 2022, 10, 70–77.
15. N. Hellgren, A. Sredenschek, A. Petruins, J. Palisaitis, F. F. Klimashin, M. A. Sortica, L. Hultman, P. O. Persson and J. Rosen, Synthesis and characterization of TiBx (1.2 × 2.8) thin films grown by DC magnetron co-sputtering from TiB₂ and Ti targets, Surf. Coat. Technol., 2022, 433, 128110.
16. C. Fuger, R. Hahn, A. Hirle, P. Kutrowatz, M. Weiss, A. Limbeck, O. Hunold, P. Polcik and H. Riedl, Revisiting the origins of super-hardness in TiB2 + z thin films–Impact of growth conditions and anisotropy, Surf. Coat. Technol., 2022, 446, 128806.
17. J. Palisaitis, M. Dahlqvist, A. Hall, J. Thörnberg, I. Persson, N. Nedfors, L. Hultman, J. Greene, I. Petrov, J. Rosen and O. Persson, Where is the unpaired transition metal in substoichiometric diboride line compounds?, Acta Mater., 2021, 204, 116510.
18. J. Palisaitis, M. Dahlqvist, L. Hultman, I. Petrov, J. Rosen and P. O. Persson, On the nature of planar defects in transition metal diboride line compounds, Materialia, 2022, 24, 101478.
19. V. Moraes, C. Fuger, V. Paneta, D. Primetzhofer, P. Polcik, H. Bolvardi, M. Arndt, H. Riedl and P. Mayrhofer, Substoichiometry and tantalum dependent thermal stability of α -structured W-Ta-B thin films, Scr. Mater., 2018, 155, 5–10.
20. S. Lin, Z. Chen, R. Janknecht, Z. Zhang, L. Hultman, P. H. Mayrhofer, N. Koutná and D. G. Sangiovanni, Machine-learning potentials predict orientation- and mode-dependent fracture in refractory diborides, Acta Mater., 2025, 301, 121568.
21. L. Xue, G. Feng and S. Liu, Molecular dynamics study of temperature effect on deformation behavior of m-plane 4H–SiC film by nanoindentation, Vacuum, 2022, 202, 111192.
22. M. Dahlqvist and J. Rosen, Impact of vacancies on structure, stability and properties of hexagonal transition metal diborides, MB₂ (M = Sc, Y, Ti, Zr, Hf, V, Nb, Ta, Cr, Mo, W, Mn, and Fe), Materialia, 2022, 26, 101629.
23. G. Almyras, D. G. Sangiovanni and K. Sarakinos, Semi-empirical force-field model for the Ti_{1 ( x} AlxN (0 × 1) system, Materials, 2019, 12, 215.
24. D. G. Sangiovanni, A. Kraych, M. Mrovec, J. Salamania, M. Odén, F. Tasnádi and I. A. Abrikosov, Descriptor for slip-induced crack blunting in refractory ceramics, Phys. Rev. Mater., 2023, 7, 103601.
25. S. Attarian and S. Xiao, Development of a 2NN-MEAM potential for TiB system and studies of the temperature dependence of the nanohardness of TiB₂, Comput. Mater. Sci., 2022, 201, 110875.
26. P. Rowe, G. Csányi, D. Alfè and A. Michaelides, Development of a machine learning potential for graphene, Phys. Rev. B, 2018, 97, 054303.
27. M. Qamar, M. Mrovec, Y. Lysogorskiy, A. Bochkarev and R. Drautz, Atomic cluster expansion for quantum-accurate large-scale simulations of carbon, J. Chem. Theory Comput., 2023, 19, 5151–5167.
28. Y. Lysogorskiy, A. Bochkarev, M. Mrovec and R. Drautz, Active learning strategies for atomic cluster expansion models, Phys. Rev. Mater., 2023, 7, 043801.
29. S. Lin, L. Casillas-Trujillo, F. Tasnádi, L. Hultman, P. H. Mayrhofer, D. G. Sangiovanni and N. Koutná, Machine-learning potentials for nanoscale simulations of tensile deformation and fracture in ceramics, NPJ Comput. Mater., 2024, 10, 67.
30. A. Shapeev, Moment Tensor Potentials: a class of systematically improvable interatomic potentials, Multiscale Model. Simul., 2016, 14, 1153–1173.
31. J. Thörnberg, J. Palisaitis, N. Hellgren, F. F. Klimashin, N. Ghafoor, I. Zhirkov, C. Azina, J.-L. Battaglia, A. Kusiak and M. A. Sortica, et al., Microstructure and materials properties of understoichiometric TiB_x thin films grown by HiPIMS, Surf. Coat. Technol., 2020, 404, 126537.
32. Y. V. Bogatov, V. Shcherbakov, D. Y. Kovalev and A. Sychev, Effect of Al on the Phase Composition, Microstructure, and Properties of SHS Composites Based on the Ti–B System, Inorg. Mater., 2022, 58, 1321–1327.
33. L. Enzlberger, A. Kirnbauer, T. Astecker, T. Wojcik, M. Podsednik, A. Limbeck, S. Kolozsvári and P. H. Mayrhofer, Tuning microstructure and mechanical properties in superstoichiometric to substoichiometric TiB₂ thin films: a DCMS and HiPIMS study, Surf. Coat. Technol., 2025, 132811.
34. G. Kresse and J. Hafner, Ab initio molecular dynamics for liquid metals, Phys. Rev. B: Condens. Matter Mater. Phys., 1993, 47, 558.
35. G. Kresse and J. Furthmüller, Efficiency of ab initio total energy calculations for metals and semiconductors using a plane-wave basis set, Comput. Mater. Sci., 1996, 6, 15–50.
36. G. Kresse and J. Furthmüller, Efficient iterative schemes for ab initio total-energy calculations using a plane-wave basis set, Phys. Rev. B: Condens. Matter Mater. Phys., 1996, 54, 11169.
37. R. Jinnouchi, F. Karsai and G. Kresse, On-the-fly machine learning force field generation: Application to melting points, Phys. Rev. B, 2019, 100, 014105.
38. G. Kresse and D. Joubert, From ultrasoft pseudopotentials to the projector augmented-wave method, Phys. Rev. B: Condens. Matter Mater. Phys., 1999, 59, 1758.
39. A. Jain, S. P. Ong, G. Hautier, W. Chen, W. D. Richards, S. Dacek, S. Cholia, D. Gunter, D. Skinner, G. Ceder and K. a. Persson, The Materials Project: A materials genome approach to accelerating materials innovation, APL Mater., 2013, 1, 011002.
40. I. S. Novikov, K. Gubaev, E. V. Podryabinkin and A. V. Shapeev, The MLIP package: moment tensor potentials with MPI and active learning, Mach. Learn. Sci. Technol., 2020, 2, 025002.
41. S. A. Goreinov, I. V. Oseledets, D. V. Savostyanov, E. E. Tyrtyshnikov and N. L. Zamarashkin, Matrix Methods: Theory, Algorithms And Applications: Dedicated to the Memory of Gene Golub, World Scientific, 2010, pp. 247–256.
42. E. V. Podryabinkin and A. V. Shapeev, Active learning of linearly parametrized interatomic potentials, Comput. Mater. Sci., 2017, 140, 171–180.
43. N. Koutná, S. Lin, L. Hultman, D. G. Sangiovanni and P. H. Mayrhofer, Machine-learning potentials for structurally and chemically complex MAB phases: strain hardening and ripplocation-mediated plasticity, Mater. Des., 2025, 114307.
44. A. P. Thompson, H. M. Aktulga, R. Berger, D. S. Bolintineanu, W. M. Brown, P. S. Crozier, P. J. In’t Veld, A. Kohlmeyer, S. G. Moore and T. D. Nguyen, et al., LAMMPS - a flexible simulation tool for particle-based materials modeling at the atomic, meso, and continuum scales, Comput. Phys. Commun., 2022, 271, 108171.
45. J. R. Gissinger, I. Nikiforov, Y. Afshar, B. Waters, M.-k. Choi, D. S. Karls, A. Stukowski, W. Im, H. Heinz and A. Kohlmeyer, et al., Type label framework for bonded force fields in LAMMPS, J. Phys. Chem. B, 2024, 128, 3282–3297.
46. A. Stukowski, Visualization and analysis of atomistic simulation data with OVITO-the Open Visualization Tool, Modelling Simul. Mater. Sci. Eng., 2010, 18, 015012.
47. K. Momma and F. Izumi, VESTA 3 for three-dimensional visualization of crystal, volumetric and morphology data, Appl. Crystallogr., 2011, 44, 1272–1276.
48. I. R. Shein and A. L. Ivanovskii, Elastic properties of mono-and polycrystalline hexagonal AlB2-like diborides of s, p and d metals from first-principles calculations, J. Phys.: Condens. Matter, 2008, 20, 415218.
49. W. C. Oliver and G. M. Pharr, Measurement of hardness and elastic modulus by instrumented indentation: Advances in understanding and refinements to methodology, J. Mater. Res., 2004, 19, 3–20.
50. W. Xia, J. Song, D. D. Hsu and S. Keten, Understanding the interfacial mechanical response of nanoscale polymer thin films via nanoindentation, Macromolecules, 2016, 49, 3810–3817.
51. A. Shapeev, Accurate representation of formation energies of crystalline alloys with many components, Comput. Mater. Sci., 2017, 139, 26–30.
52. G. Volonakis, L. Tsetseris and S. Logothetidis, Electronic and structural properties of TiB₂ : Bulk, surface, and nanoscale effects, Mater. Sci. Eng., B, 2011, 176, 484–489.
53. B. Post, F. W. Glaser and D. Moskowitz, Transition metal diborides, Acta Metall., 1954, 2, 20–25.
54. R. G. Munro, Material properties of titanium diboride, J. Res. Natl. Inst. Stand. Technol., 2000, 105, 709.
55. D.-l. Ouyang, T. Cheng and C. Xia, et al., Experiment and modeling of TiB2/TiB boride layer of Ti- 6Al- 2Zr- 1Mo- 1 V alloy, Trans. Nonferrous Met. Soc. China, 2021, 31, 3752–3761.
56. G. Kara and G. Purcek, Growth kinetics and mechanical characterization of boride layers formed on β -type Ti-45Nb alloy, Surf. Coat. Technol., 2018, 352, 201–212.
57. M. Keddam, S. Taktak and S. Tasgetiren, A diffusion model for the titanium borides on pure titanium, Surf. Eng., 2016, 32, 802–808.
58. D. Sangiovanni, F. Tasnadi, L. Johnson, M. Odén and I. Abrikosov, Strength, transformation toughening, and fracture dynamics of rocksalt-structure Ti1-x AlxN (0 × 0.75) alloys, Phys. Rev. Mater., 2020, 4, 033605.
59. A. Mei, H. Kindlund, E. Broitman, L. Hultman, I. Petrov, J. E. Greene and D. Sangiovanni, Adaptive hard and tough mechanical response in single-crystal B1 VN_x ceramics via control of anion vacancies, Acta Mater., 2020, 192, 78–88.
60. N. L. Okamoto, M. Kusakari, K. Tanaka, H. Inui and S. Otani, Anisotropic elastic constants and thermal expansivities in monocrystal CrB₂, TiB₂, and ZrB₂, Acta Mater., 2010, 58, 76–84.
61. G. N. Greaves, A. L. Greer, R. S. Lakes and T. Rouxel, Poisson's ratio and modern materials, Nat. Mater., 2011, 10, 823–837.
62. M. Berger, E. Coronel and E. Olsson, Microstructure of dc magnetron sputtered TiB₂ coatings, Surf. Coat. Technol., 2004, 185, 240–244.
63. M. Stüber, H. Riedl, T. Wojcik, S. Ulrich, H. Leiste and P. Mayrhofer, Microstructure of Al-containing magnetron sputtered TiB₂ thin films, Thin Solid Films, 2019, 688, 137361.
64. P. Rupa, P. Chakraborti and S. K. Mishra, Mechanical and deformation behaviour of titanium diboride thin films deposited by magnetron sputtering, Thin Solid Films, 2009, 517, 2912–2919.
65. O. Sobol, S. Dub, A. Pogrebnjak, R. Mygushchenko, A. Postelnyk, A. Zvyagolsky and G. Tolmachova, The effect of low titanium content on the phase composition, structure, and mechanical properties of magnetron sputtered WB₂ -TiB₂ films, Thin Solid Films, 2018, 662, 137–144.
66. A. Hasnaoui, P. Derlet and H. Van Swygenhoven, Interaction between dislocations and grain boundaries under an indenter–a molecular dynamics simulation, Acta Mater., 2004, 52, 2251–2258.
67. P. H. Mayrhofer, H. Clemens and C. Mitterer, Interfaces in nanostructured thin films and their influence on hardness, Int. J. Mater. Res., 2022, 96, 468–480.
68. W. Dai, T. Zhang, J. Yang, R. Sun and J. Xu, Morphological analysis of TiB₂ thin film prepared by rf magnetron sputtering, J. Vac. Sci. Technol., A, 2008, 26, 610–615.
69. H. Deng, J. Chen, R. Inturi and J. Barnard, Structure, mechanical and tribological properties of dc magnetron sputtered TiB₂ and TiB₂ (N) thin films, Surf. Coat. Technol., 1995, 76, 609–614.
70. I. Zyganitidis, N. Kalfagiannis and S. Logothetidis, Ultra sharp Berkovich indenter used for nanoindentation studies of TiB₂ thin films, Mater. Sci. Eng., B, 2009, 165, 198–201.
71. H. Jia, X. Liu, Z. Li, S. Sun and M. Li, The effect of grain size on the deformation mechanisms and mechanical properties of polycrystalline TiN: A molecular dynamics study, Comput. Mater. Sci., 2018, 143, 189–194.

---