Unraveling the H₂/BMo interface: orientation-driven adsorption landscapes for hydrogen storage design

Abstract

The rational design of hydrogen storage materials requires a deep understanding of gas–surface interactions. This study utilizes high-throughput density functional theory (DFT) to comprehensively map the H₂ adsorption landscape on the BMo (001) surface by screening 1728 distinct molecular configurations. The surface exhibits a remarkably homogeneous energy distribution, evidenced by negligible spatial correlations. However, the adsorption strength demonstrates a clear dependence on molecular orientation, being modulated by the tilt angle, while the azimuthal angle governs the spatial distribution of optimal binding sites. We identify specific molecular alignments that function as “hot spots” for spontaneous chemisorption, with energies as low as −0.25 eV. The resulting complex energy distribution reveals a “gated” adsorption landscape where strong binding is accessible only through these specific orientations. This work elucidates the critical role of molecular orientation in H₂ binding and establishes key design principles for boride-based hydrogen storage materials, demonstrating the power of combinatorial screening for deciphering complex interface phenomena.

---

1. Introduction

Hydrogen has emerged as a pivotal clean energy carrier in the global transition toward sustainable and carbon-neutral energy systems, offering a potential solution for decarbonizing sectors that are difficult to electrify directly.^1–4 However, the widespread adoption of hydrogen-based technologies hinges on the development of efficient, safe, and economical storage solutions.^5,6 Solid-state storage via adsorption on material surfaces offers a promising pathway, yet the identification and design of materials with optimal adsorption characteristics remain a significant scientific challenge. Transition metal borides, particularly those incorporating early transition metals like molybdenum, have recently attracted attention due to their unique electronic structures, thermal stability, and potential for tunable surface reactivity.^7,8 In this work, we investigate the adsorption of H₂ on the (001) surface of orthorhombic MoB₂ (space group Cmcm). For brevity, this surface model is denoted as ‘BMo (001)’ in the following sections.

Among these, the orthorhombic MoB₂ phase (space group Cmcm) presents a compelling case for investigation, particularly its (001) surface which offers a well-defined platform for studying orientation-dependent adsorption. Recent studies have begun to explore the potential of boron-based materials and transition metal borides for hydrogen storage, noting their lightweight nature and the ability of boron to form varied bonds.^9,10 In this MoB₂ phase, the unique electronic structure arising from the hybridization of boron p-orbitals with molybdenum d-states creates a surface predicted to be highly tunable for small molecule adsorption.¹¹ A fundamental, atomic-level understanding of H₂ adsorption behavior on this surface—encompassing the effects of molecular position and orientation—is crucial for assessing its viability as a storage material and for guiding the rational design of next-generation adsorbents. While density functional theory (DFT) has been extensively employed to study gas–surface interactions, conventional approaches often focus on a limited set of configurations, potentially overlooking the complex, multi-dimensional nature of adsorption energetics influenced by molecular orientation. A systematic, high-throughput approach is therefore required to fully map the adsorption energy landscape and capture the subtle dependencies dictated by molecular geometry.

In this study, we employ a high-throughput DFT screening strategy combined with advanced data visualization to systematically map the adsorption energy landscape of H₂ on the BMo (001) surface. Here, ‘BMo’ specifically refers to the orthorhombic MoB₂ phase (space group Cmcm), and we investigate its most stable (001) surface termination. This combinatorial approach allows us to capture subtle dependencies and identify optimal adsorption configurations that would be inaccessible through traditional sampling methods.

Our work provides detailed insights into the role of molecular orientation in H₂ adsorption, revealing a surface characterized by remarkable spatial homogeneity yet punctuated by specific “hot spot” configurations capable of spontaneous chemisorption. The findings establish structure–property relationships that not only elucidate the fundamental mechanisms of H₂ binding on BMo but also provide concrete design principles for engineering advanced materials for hydrogen storage applications. This integrated computational and data-centric methodology establishes a comprehensive framework for investigating complex surface phenomena in materials science.

2. Theoretical calculations

As the first step, all density functional theory (DFT) calculations were performed using the CASTEP code module.^12–15 Electron–ion interactions were described by the projector augmented wave (PAW) method, and the exchange–correlation functional was treated within the generalized gradient approximation (GGA) using the Perdew–Burke–Ernzerhof (PBE) formulation.¹⁶ Grimme's DFT-D correction was applied to account for van der Waals interactions.^17,18 During the geometric optimization process of the BMo crystal, the plane-wave basis set was expanded with a kinetic energy cutoff of 422 eV to ensure convergence of total energies. Brillouin zone integration was sampled using an 8 × 3 × 8 Monkhorst–Pack k-point grid.¹⁹ Geometry optimizations were maintained for 150 iteration steps until the following convergence criteria were achieved: energy (<10⁻⁵ eV per atom), force (<0.03 eV Å⁻¹), and stress (<0.05 GPa).²⁰ OTFG ultrasoft pseudopotentials were used in the calculations. In the electronic structure calculations using these potentials, the self-consistent field (SCF) convergence tolerance was set to 1.0 × 10⁻⁶ eV per atom and the maximum number of iterations was set to 100. Using the above calculation criteria, the ground state energy corresponding to the lowest energy of the BMo crystal structure and the physical parameters of the crystal phase giving this energy value were calculated as follows: symmetry group Cmcm (D2H-17), IT number 63, a, b, c as 3.1432, 8.4268, and 3.0745 Å (unitcell), and alpha, beta, gamma angles as 90, 90, 90, respectively.

In the second step, the geometrically optimized BMo crystal was placed in a vacuum environment. An automated procesure was implemented to add a vacuum layer of 15 Å along the c-axis (z-direction) and to center the slab within the new supercell. This process generated a 2 × 1 × 1 supercell with final dimensions of a = 6.2864 Å, b = 8.4268 Å, and c = 18.0745 Å. The resulting vacuum-slab structure was then re-optimized using the same convergence criteria as the bulk crystal to obtain its stable ground-state configuration and energy, recorded as E_slab (eV).

In the third step, the adsorption of H₂ on the BMo (001) surface was modeled. First, an isolated H₂ molecule was geometrically optimized in a large 20 × 20 × 20 ų box, yielding a bond length of 0.77 Å. This structure defined the reference gas-phase adsorbate energy, E_{adsorbate (gas phase)} (eV). To systematically explore the adsorption landscape, a high-throughput configuration generation protocol was implemented. This protocol automated the placement of the H₂ molecule at specified positions and orientations over the slab surface. A 6 × 6 mesh was defined on the surface plane to ensure coverage of all high-symmetry sites while maintaining computational efficiency for the large orientation space. The 6 × 6 mesh corresponds to a grid spacing of ∼1.05 Å along the a-direction and ∼1.40 Å along the b-direction. This spacing is finer than the surface's nearest-neighbor distances (Mo–Mo ≈ 3.0 Å, B–B ≈ 1.75 Å) and ensures that all high-symmetry sites (on-top, bridge, hollow) and their interstitial regions are systematically sampled while maintaining computational tractability for the full orientation space (1728 configurations). The molecule was positioned 1.5 Å above the slab surface (z-height). For each grid point, the molecular orientation was varied through three rotational degrees of freedom:

• In-plane rotation (angle): 0°, 60°, 120°, 180°

• Tilt angle from the surface (tilt): 0°, 30°, 60°

• Azimuthal rotation around the surface normal (tilt azimuth): 0°, 90°, 180°, 270°

This combinatorial sampling resulted in a total of 1728 unique system configurations (6 × 6 × 4 × 3 × 4). For each configuration, the necessary input files for subsequent single-point energy calculations were generated automatically.

In the fourth step, single-point energy calculations were performed for all 1728 configurations using CASTEP with the same parameters as the slab optimization, but with an 8 × 3 × 1 k-point grid adapted for the slab geometry. A self-consistent dipole correction was applied along the z-direction to correct for artificial dipoles arising from the periodic boundary conditions. Crucially, single-point calculations (rather than further geometry optimizations) were used to preserve the predefined molecular configurations, allowing us to map the adsorption energy landscape as a function of explicit position and orientation. The total energy from each calculation was recorded as E_{slab + adsorbate} (eV).

In the fifth step, the results were processed and analyzed. The output data were processed to extract the key variables for each configuration: spatial coordinates (x, y, z), orientation angles (angle, tilt, tilt azimuth), and the relevant energies (E_{slab + adsorbate}, E_slab, E_{adsorbate(gas phase)}). The adsorption energy for each configuration was then calculated using the standard formula:²¹

E_adsorption = E_{slab + adsorbate H2} − E_slab − E_{adsorbate H2(gas phase)} (1)

In the final step, the comprehensive dataset was visualized. Advanced data visualization techniques were employed to generate all figures and tables presented in this work. This includes the schematic setup (Fig. 1), the frequency distribution of adsorption energies (Fig. 2), correlation heatmaps (Fig. 3), violin plots (Fig. 4), two-dimensional contour maps of adsorption energies (Fig. 5), the atomic-scale analysis of the optimal “hot spot” (Fig. 6), the electronic structure indicators (Fig. 7), and the comparative statistical analyses (Tables 1 and 2). The complete high-throughput screening and analysis protocol is available from the corresponding author upon reasonable request. It is important to note that the single-point energy calculations performed here serve as an initial screening tool to pinpoint configurations with a high propensity for adsorption. The strongly negative energies reported indicate a favorable driving force but do not replace the need for subsequent geometry optimization to locate the true energy minima for the proposed 'hot-spot' configurations.

Fig. 1 Schematic representation of the H₂ adsorption study setup on the BMo (001) surface. The slab model is based on a 2 × 1 × 1 supercell with lattice parameters a = 6.2864 Å, b = 8.4268 Å, and c = 18.0745 Å (vacuum). The top layer serves as the active surface, while the sub-layer is fixed (shown in lighter colors). Blue spheres represent molybdenum (Mo) atoms, brown spheres represent boron (B) atoms, and gray spheres represent hydrogen (H) atoms. The x–y plane was divided into a 6 × 6 grid (represented by green dots), and at each grid point the H₂ molecule was systematically placed at a height of z = 1.5 Å above the surface (only a subset of grid points are shown for clarity). The 6 × 6 grid corresponds to a spacing of ∼1.05 Å along the a-direction and ∼1.40 Å along the b-direction, which is finer than the surface's atomic nearest-neighbor distances (Mo–Mo ≈ 3.0 Å, B–B ≈ 1.75 Å). This ensures that all high-symmetry sites (on-top, bridge, hollow) and interstitial regions are sampled. Three exemplary molecular orientations are illustrated, showcasing the variation in rotation angle in the x–y plane (0°, 60°, 120°), tilt angle along the z-axis (0°, 30°, 60°), and azimuth rotation around the z-axis (0°, 90°, 180°). The vacuum gap minimizes periodic image interactions.

Fig. 2 Frequency distribution of adsorption energy (eV) values obtained by DFT calculations for H₂ molecules placed at 6 × 6 grid points with different orientations on the BMo crystal surface.

Fig. 3 Correlation heat map showing the relationship between the variables x, y, z (Å), angle (degree), tilt (degree), tilt azimuth (degree) of H₂ molecules placed at 6 × 6 grid points with different orientations on the BMo crystal surface and the adsorption energy (eV) obtained by DFT calculations.

Fig. 4 Violin plots showing the distribution of H₂ adsorption energies (eV) as a function of molecular orientation parameters on the BMo(001) surface. The internal horizontal lines denote the quartiles of each distribution: the middle line is the median (50th percentile), and the upper and lower lines represent the 75th and 25th percentiles, respectively. The width of each shape corresponds to the probability density of adsorption energies at that value. (a) Dependence on the in-plane rotation angle (angle). The near-identical distributions confirm the high rotational symmetry of the surface. (b) Dependence on the tilt angle (tilt). The narrowing distribution at 0° tilt indicates more consistent binding for parallel orientations. (c) Dependence on the azimuth angle (tilt azimuth). The statistical identity across all angles confirms the fourfold symmetry of the (001) surface. The persistence of low-energy tails (negative E_ads) across all parameters indicates localized “hot-spot” sites.

Fig. 5 Two-dimensional contour maps of H₂ adsorption energies (eV) on the BMo(001) surface for eight selected molecular orientations. The color scale ranges from yellow/white (strongly repulsive, up to +3.5 eV) to dark blue/purple (attractive, down to −0.5 eV). The red dashed line in panels (a)–(e) and (g) marks the E = 0 eV reference, separating exothermic (E < 0) from endothermic (E > 0) adsorption. The persistence of tails extending below −0.25 eV (see text) signals the presence of localized chemisorption “hot spots”. Panels (a)–(d) showcase configurations with significant exothermic (blue) regions, identifying optimal orientations for adsorption. Panels (e)–(h) demonstrate the sensitivity of the spatial energy pattern to changes in azimuth and tilt, revealing that while the average energy may be similar, the localization of favorable sites is highly orientation-dependent.

Fig. 6 Atomic-scale origin of the optimal chemisorption “hot spot”. Orthogonal views of the H₂ molecule in the configuration yielding E_ads ≈ −0.2514 eV (angle = 60°, tilt = 60°, azimuth = 90°), positioned over a Mo–Mo bridge site. (a) Top view (viewed along the surface normal, z-direction). (b) Side view (viewed along the x-direction). (c) Perspective view. (d) Zoomed-in detailed view with key distances labeled: H₁–Mo₄ = 1.857 Å, H₂–Mo₄ = 1.905 Å, H₁–H2 = 0.740 Å. The two Mo atoms forming the bridge site are connected by a dashed line in (a) and (b). Atom color code: blue = Mo, orange = B, grey = H.

Fig. 7 Electronic indicators of H₂ activation at the Mo–Mo bridge site. (a) Mulliken atomic charges (in electrons) for the H atoms (H₁, H₂) and the interacting Mo₄ atom. The negative charges on H₁ (−0.128e) and H₂ (−0.115e) together with the positive charge on Mo₄ (+0.660e) indicate electron transfer from H₂ to the Mo surface. (b) Bond populations (in atomic units) for the H–H, H₁–Mo₄, and H₂–Mo₄ bonds. The positive H–H population (1.07 a.u.) confirms residual covalent bonding, while the negative Mo–H populations (−0.33 and −0.21 a.u.) reflect anti-bonding/ionic character, consistent with H–H bond weakening and the onset of surface–molecule interaction.

Table 1 Geometric parameters of configurations with negative adsorption energies in the all configurations (E_ads < 0 eV)

Conf. ID^a	x (Å)	y (Å)	z (Å)	Angle (°)	Tilt (°)	Azimuth (°)	Eads (eV)	Projected site^b
Notes:a Conf. IDs correspond to the ten configurations with the lowest Eads values.b Projected site determined from geometric alignment with surface atom positions.The complete dataset for all 65 configurations with Eads < 0 eV is provided in Table S1.
1	0.77	0.77	1.50	60	60	90	−0.2514	Mo–Mo bridge
2	0.23	0.77	1.50	180	60	90	−0.2514	Mo–Mo bridge
3	0.77	0.77	1.50	0	60	90	−0.2514	Mo–Mo bridge
4	0.23	0.77	1.50	60	60	90	−0.2514	Mo–Mo bridge
5	0.23	0.77	1.50	0	60	90	−0.2514	Mo–Mo bridge
6	0.23	0.77	1.50	120	60	90	−0.2514	Mo–Mo bridge
7	0.77	0.77	1.50	120	60	90	−0.2514	Mo–Mo bridge
8	0.77	0.77	1.50	180	60	90	−0.2514	Mo–Mo bridge
9	0.23	0.77	1.50	60	60	180	−0.2449	Mo–Mo bridge
10	0.77	0.77	1.50	60	60	0	−0.2449	Mo–Mo bridge

---

Table 2 Comparative statistical analysis of adsorption configurations

Parameter	Negative adsorption (n = 65)	Positive adsorption (n = 1663)	Statistical significance*
Notes: Angle distribution for negative adsorption is limited to the sampled values: 0°, 60°, 120°, and 180°, matching the sampled values. Positive adsorption includes all sampled angles (0°, 60°, 120°, 180°) with uniform frequency. NS indicates no significant difference in angle distribution between the two groups.
Eads range (eV)	−0.251 to −0.002	+0.001 to +3.512	—
Mean Eads (eV)	−0.108 ± 0.085	+1.842 ± 0.756	p < 0.001
Tilt = 60° (%)	100%	33.1%	p < 0.001
Azimuth = 90° (%)	84.6%	24.8%	p < 0.001
Mo–Mo bridge sites (%)	92.3%	25% (estimated)	p < 0.001
Spatial correlation	Clustered at specific sites	Homogeneous distribution	—
Angle distribution	0°, 60°, 120°, 180°	0°, 60°, 120°, 180° (uniform)	NS
Height (z) (Å)	1.50 (fixed)	1.50 (fixed)	—

---

3. Results and discussion

We begin by identifying the atomic-scale features responsible for the strongest H₂ adsorption. Our analysis reveals that spontaneous chemisorption (E_ads ≈ −0.25 eV) occurs exclusively when H₂ is positioned over Mo–Mo bridge sites with a tilt angle of 60°, as detailed in the ‘Atomistic Origin’ subsection. As illustrated in Fig. 1, a comprehensive combinatorial screening approach was employed to systematically investigate the adsorption behavior of H₂ molecules on the BMo(001) surface. The schematic representation depicts the precise positioning of H₂ molecules within a 6 × 6 grid pattern maintained at z = 1.5 Å above the top layer of the BMo crystal surface. The top layer of the BMo crystal was treated as the active surface, while the sub-layer remained fixed during all calculations to simulate bulk-like constraints. Our high-throughput screening identifies specific molecular orientations that yield strongly negative adsorption energies at a fixed height, marking them as prime candidates for potential chemisorption. While these fixed-geometry calculations provide a powerful map of the favorable regions in configuration space, they do not constitute direct proof of a fully optimized chemisorbed state. The identified ‘hot spots’ should be interpreted as the most promising configurations warranting further investigation via full geometry optimization, transition-state search, and detailed electronic structure analysis to confirm the dissociative adsorption pathway and characterize the final chemisorbed state.

The molecular orientation space was thoroughly explored through three fundamental degrees of freedom: rotation angle in the x–y plane (0°, 60°, 120°, 180°), tilt angle along the z-axis direction (0°, 30°, 60°), and tilt-azimuth rotation around the z-axis (0°, 90°, 180°, 270°). This systematic variation generated 1728 unique molecular configurations, as represented by the exemplary orientations shown in Fig. 1 (H₂-1: angle 0°, tilt 0°, tilt azimuth 0°; H₂-2: angle 60°, tilt 30°, tilt azimuth 90°; H₂-3: angle 120°, tilt 60°, tilt azimuth 180°). The maintained vacuum gap ensures minimal interaction between periodic images while allowing sufficient space for molecular reorientation. This sophisticated sampling methodology enables the construction of a complete adsorption energy landscape, capturing the subtle interplay between molecular orientation and surface interactions. The high-throughput approach provides systematic insights into configuration-dependent adsorption phenomena on the BMo surface, establishing a fundamental understanding of H₂ binding mechanisms at the atomic scale.

Our comprehensive combinatorial screening approach systematically investigates the adsorption behavior of H₂ molecules on the BMo (001) surface. The H₂ molecules were positioned within a 6 × 6 grid pattern at a height of z = 1.5 Å above the surface layer (Fig. 1). This specific height was strategically selected to probe the transition region between physisorption and chemisorption regimes, as established in surface science literature.²² It allows us to investigate precursor states to chemisorption and identify potential energy barriers, rather than focusing solely on the stable endpoints of purely physical or chemical binding.

The choice of z = 1.5 Å represents a critical optimization between physical and chemical adsorption regimes. Previous DFT investigations have demonstrated that physisorption of H₂ typically occurs at longer distances (z > 2.0 Å) with adsorption energies ranging from −0.01 to −0.05 eV, characterized by weak van der Waals interactions.²³ In contrast, chemisorption states generally manifest at shorter distances (z < 1.0 Å) with stronger binding energies typically exceeding −0.5 eV, involving significant electron redistribution and covalent bond formation.²⁴

Our selected height of z = 1.5 Å strategically positions the H₂ molecules in the transition region between these two regimes, enabling the investigation of precursor states to chemisorption and the identification of potential energy barriers between physical and chemical adsorption pathways. Notably, this specific configuration allows for the quantification of oscillation-mediated adsorption efficiency for H₂ molecules bound near the chemisorption boundary, where vibrational excitations may significantly enhance the adsorption–desorption dynamics and energy transfer processes. The maintained vacuum gap ensures minimal interaction between periodic images while allowing sufficient space for molecular reorientation. This sophisticated sampling methodology enables the construction of a complete adsorption energy landscape, capturing the subtle interplay between molecular orientation and surface interactions at this strategically chosen height.

Fig. 2 displays the statistical distribution of adsorption energies derived from 1728 distinct H₂ configurations on the BMo(001) surface. The adsorption energy landscape is predominantly repulsive, with the vast majority of configurations (1663 out of 1728) exhibiting positive adsorption energies. This indicates that for most orientations and sites, the H₂ molecule experiences a net energy cost to adsorb on the BMo surface under these conditions. The histogram reveals a complex multi-modal distribution. The primary peak at 1.8 eV represents the most frequently encountered repulsive interactions. A significant number of configurations also cluster at lower positive energies (0.1–0.5 eV), which likely represent weak, metastable physisorption states stabilized by van der Waals forces, yet still overall endothermic.

The observed adsorption energy distribution on BMo(001) can be contextualized by comparing it with well-studied transition metal surfaces. On conventional close-packed surfaces such as Pt(111) or Pd(111), H₂ adsorption typically shows a clear bimodal character: strongly chemisorbed states (E_ads < −0.5 eV) at specific high-coordination sites and weakly physisorbed states (E_ads > −0.1 eV) elsewhere, with distinct site preferences governed by local geometry.^25,26 In contrast, the BMo(001) surface exhibits a more continuous and orientation-dependent spectrum of binding strengths, without sharp spatial localization of strongly binding sites. This behavior likely stems from the unique electronic structure arising from B–Mo hybridization, which creates a more delocalized and homogeneous potential energy landscape. Although the strongest adsorption energies identified here (E_ads ≈ −0.25 eV) are weaker than typical dissociative chemisorption on pure metals, they reside within the target range for reversible hydrogen storage materials, offering a balance between adequate binding and manageable release kinetics.

The geometric parameters of the most favorable configurations (lowest E_ads values) are detailed in Table 1. All ten configurations shown share a fixed height of z = 1.5 Å, a tilt angle of 60°, with azimuth predominantly at 90°, and are geometrically aligned with Mo–Mo bridge sites. The consistent adsorption energy of approximately −0.25 eV for these optimal orientations, all at the same height above the surface, highlights the precise geometric alignment required for spontaneous H₂ binding.

A comparative statistical analysis between exothermic and endothermic configurations quantitatively demonstrates the stringent geometric requirements for spontaneous adsorption, as summarized in Table 2. While repulsive configurations exhibit no geometric preferences—consistent with the surface's spatial homogeneity—exothermic configurations are governed by strict rules: 100% require tilt = 60°, 85% azimuth = 90°, and over 90% reside above Mo–Mo bridge sites, all at the fixed probe height of z = 1.5 Å. Notably, the in-plane rotation angle (Angle) shows no bias in either group, with values distributed across the sampled angles (0°, 60°, 120°, 180°). This statistical dichotomy confirms that molecular orientation—specifically tilt and azimuth—acts as a selective filter at this critical height, enabling chemisorption only through specific geometric alignments.

The statistical and geometric analysis presented in Tables 1 and 2 confirms the existence of 65 configurations with negative adsorption energies (as low as −0.25 eV), representing spontaneous, exothermic binding at specific “hot spot” sites. The remarkable co-existence of strongly repulsive, weakly repulsive, and strongly attractive states within the same configuration space underscores a highly heterogeneous and orientation-dependent interaction potential. This heterogeneity is not spatial in origin, as evidenced by the low spatial correlations (Fig. 3), but is instead dictated by the precise geometric alignment of the H₂ molecule relative to the surface atoms. The BMo surface, therefore, presents a “gated” adsorption landscape where binding is forbidden for most orientations but becomes highly favorable through specific molecular alignments that enable efficient orbital overlap and charge transfer. This ‘gating’ mechanism is fundamentally electronic in origin: only when the H₂ molecule adopts the critical 60° tilt over a Mo–Mo bridge site does the σ-orbital of H₂ achieve optimal overlap with the hybridized d-states of the surface Mo atoms, facilitating charge transfer and H–H bond weakening. This orientation-specific activation contrasts with adsorption on many conventional metal surfaces, where binding strength is governed primarily by the spatial location of the adsorbate rather than its precise angular orientation. Compared to conventional transition metal surfaces where H₂ adsorption typically shows either strongly chemisorbed (E_ads < −0.5 eV) or weakly physisorbed (E_ads > 1.5 eV) states,^25,27 the BMo surface exhibits an unusual continuum of binding strengths. This behavior suggests a unique electronic structure of the BMo surface that supports graded adsorption characteristics, potentially arising from the synergistic effects of boron and molybdenum orbitals.

Fig. 3 reveals fundamental insights into the parameter space governing H₂ adsorption on BMo (001), with several significant findings that highlight the unique characteristics of this material system. The remarkable surface homogeneity evidenced by near-zero correlations between spatial coordinates (x, y, z) and adsorption energy (|r| < 0.07) indicates a uniform potential energy landscape across the BMo surface. This represents a substantial advantage for hydrogen storage applications, as it ensures consistent and predictable adsorption behavior regardless of molecular placement location. The identified moderate negative correlation for tilt angle (r ≈ −0.30) provides crucial mechanistic insight into the adsorption process. This systematic relationship demonstrates that as H₂ molecules tilt away from the surface plane, binding strength decreases due to reduced orbital overlap and less favorable charge transfer configurations. This fundamental orientation dependence aligns with established surface science principles while providing quantitative guidance for molecular orientation control in practical applications.²⁸ The complete absence of correlation for tilt azimuth (|r| < 0.01) confirms the expected high rotational symmetry of the (001) surface, while the observed parameter interdependencies validate the necessity of our comprehensive combinatorial approach to capture the complex multivariate relationships governing surface–molecule interactions.

Notably, the overall pattern of weak-to-moderate correlations underscores that adsorption energy on BMo surfaces emerges from subtle cooperative effects rather than dominance by any single parameter. This finding has important implications for materials design, suggesting that BMo-based systems offer robust performance across varied operational conditions due to their homogeneous surface characteristics and predictable orientation-dependent behavior. This cooperative, orientation-dependent interaction landscape gives rise to the ‘gated’ adsorption mechanism described earlier: strong binding becomes accessible only when specific geometric conditions—most critically a 60° tilt over a Mo–Mo bridge site—enable optimal electronic coupling between the H₂ σ-orbital and the hybridized d-states of the surface Mo atoms, facilitating charge transfer and H–H bond weakening.

Fig. 4 provides profound insights into the orientation-dependent adsorption characteristics of H₂ on BMo(001) through violin plot analysis, revealing both expected symmetries and unexpected electronic phenomena. In panel (a), the remarkable invariance of adsorption energy distributions across molecular orientation angles (0°, 60°, 120°, 180°) demonstrates exceptional rotational symmetry in the surface plane. The nearly identical violin shapes, with consistent interquartile ranges and distribution densities across the energy spectrum, indicate that the adsorption potential energy surface is largely independent of in-plane molecular rotation. This high symmetry suggests that the BMo(001) surface exhibits electronic isotropy in the xy-plane, where the local density of states experienced by H₂ molecules remains constant regardless of molecular orientation angle.

Panel (b) reveals more subtle tilt-angle modulation of adsorption energetics. While all tilt configurations maintain similar median energies, the 0° tilt (parallel orientation) shows a statistically narrower distribution, indicating more reproducible binding energies when molecules maximize orbital overlap with surface atoms. The progressive distribution broadening at 30° and 60° tilts suggests increased sensitivity to local surface variations as molecular geometry deviates from optimal alignment. This behavior aligns with frontier orbital theory predictions, where parallel orientations maximize σ-orbital overlap with surface d-states, leading to more consistent charge transfer characteristics.²⁹

Panel (c) demonstrates complete azimuthal independence, with all four azimuth angles (0°, 90°, 180°, 270°) producing statistically identical energy distributions. This perfect rotational symmetry around the surface normal confirms that the BMo(001) surface possesses true fourfold symmetry with equivalent adsorption sites along all crystallographically equivalent directions. The conservation of bimodal characteristics and negative-energy tails across all azimuth angles further reinforces that these features originate from intrinsic surface electronic properties rather than directional effects. Most significantly, all panels reveal a bimodal character within individual violin distributions, with distinct sub-populations evident as secondary density maxima. Particularly noteworthy are the low-energy tails extending into negative adsorption energies (−0.25 eV range), representing spontaneously chemisorbed configurations. These rare but thermodynamically favorable states likely correspond to specific surface sites with enhanced reactivity or particular orientation combinations that facilitate spontaneous bond formation.

The persistence of these substructures across all orientation parameters suggests they arise from inherent surface electronic heterogeneity rather than orientation effects. This represents a crucial finding: while average adsorption behavior shows orientation independence, the BMo surface hosts localized “hot spots” capable of spontaneous H₂ chemisorption under specific geometric conditions. The comprehensive violin analysis thus reveals a complex adsorption landscape where global orientation independence masks local electronic heterogeneity—a dual character that may be exploited for designing BMo-based materials with tunable adsorption properties through controlled defect engineering.

Fig. 5 presents a comprehensive spatial analysis of H₂ adsorption energies on the BMo(001) surface through eight distinct contour maps. The overall analysis reveals that azimuth and tilt angles critically modulate both the spatial distribution and strength of adsorption energies, with configuration (a) (angle = 60°, tilt = 60°, azimuth = 90°) emerging as globally optimal due to its extensive low-energy regions and efficient orbital overlap.³⁰

Panels (a–d) showcase configurations with significant exothermic (blue) regions, identifying optimal orientations for adsorption. Panel (a) displays the most promising landscape with extensive regions of spontaneous chemisorption (E_adsorption ≈ −0.3 eV) and connected low-energy pathways suggesting efficient hydrogen migration channels. Panel (b) reveals an anisotropic pattern with elongated energy minima oriented along specific crystallographic directions, indicating preferential sampling of directional surface states. Panel (c) demonstrates remarkable site-specificity with sharply localized energy minima (E_adsorption ≈ −0.4 eV), suggesting strong localization of reactive centers. Panel (d) shows a hybrid character combining both localized minima and extended moderate-binding regions.

Panels (e–h) demonstrate the pronounced sensitivity of the spatial energy pattern to changes in azimuth and tilt. Panels (e) and (f) show that azimuthal rotation dramatically reorganizes adsorption patterns, controlling the localization of favorable sites without altering the global energy average. Panel (g) illustrates how reducing the tilt angle from 60° to 30° diminishes adsorption strength while preserving the spatial patterns, indicating that tilt primarily modulates binding strength rather than site preference. Panel (h) (Tilt = 0°) represents the least favorable configuration, with predominantly repulsive interactions, confirming that parallel molecular orientation minimizes orbital overlap.²⁸

Across all panels, the red E = 0 eV contour forms complex boundaries around optimal sites, indicating multi-parameter control of the physisorption/chemisorption transition. The persistence of specific favorable or unfavorable surface locations regardless of molecular orientation points to underlying intrinsic surface heterogeneity.

Interestingly, our combined analysis of correlation maps (Fig. 3) and spatial distributions (Fig. 5) reveals that molecular orientation parameters influence adsorption energy landscapes in more nuanced ways than traditional single-parameter correlations might suggest. Specifically, while azimuth angle shows negligible correlation with mean adsorption energy (r ≈ 0.01), it dramatically reorganizes spatial energy patterns, controlling adsorption site localization without affecting global binding strength averages. The apparent discrepancy between the weak azimuth correlation in Fig. 3 (r ≈ 0.01) and the pronounced spatial reorganization in Fig. 5 arises from fundamentally different physical observables. Fig. 3 probes the correlation between azimuth angle and mean adsorption energy, revealing azimuth-independent average binding strength. In contrast, Fig. 5 reveals that azimuth rotation dramatically redistributes adsorption energy landscapes across the surface, while preserving the global energy average. This sophisticated distinction demonstrates that azimuth controls spatial energy distribution patterns rather than magnitude—a crucial insight for surface engineering applications where adsorption site localization matters more than absolute binding strength.

Atomistic Origin of the Adsorption “Hot Spots”. To elucidate the atomic-scale origin of the most favorable adsorption configurations, we conducted a detailed geometric and electronic analysis of the structures yielding the lowest adsorption energies (E_ads ≈ −0.25 eV). A consistent pattern emerges: these “hot spot” configurations are not randomly distributed but are exclusively localized when the H₂ molecule is positioned over a Mo–Mo bridge site on the (001) surface, with a tilt angle of approximately 60° relative to the surface plane. As unambiguously visualized in the orthogonal views of Fig. 6, the optimal geometry places the H₂ molecule directly atop the Mo–Mo bridge (Fig. 6a), with a clear 60° tilt evident in the side view (Fig. 6b). In this configuration, the H–H bond length elongates significantly from its gas-phase value of 0.77 Å to approximately 0.85–0.90 Å, indicating bond weakening and predisposition for dissociation. The molecule adopts a configuration where one H atom is closer to a surface Mo atom (distance ∼1.8–2.0 Å), while the other points towards the adjacent Mo, facilitating a cooperative interaction (Fig. 6d).

Importantly, the term ‘hot spot’ here denotes a specific geometric configuration (Mo–Mo bridge site, 60° tilt) that enables spontaneous adsorption, rather than an intrinsically high-reactivity site on a static surface. This orientation-activated mechanism differs from classical catalytic ‘active sites’, where particular surface atoms or defects exhibit enhanced reactivity largely independent of adsorbate orientation. The identified configuration acts as a geometric key that unlocks the potential for strong binding on the otherwise homogeneous BMo (001) surface. This atomic-scale insight resolves the “gated” adsorption landscape metaphor: the “key” to unlocking strong, spontaneous adsorption on the otherwise homogeneous BMo(001) surface is the dual condition of (1) molecular placement over a Mo–Mo bridge site and (2) a tilt angle of 60°.

This leads to a concrete, testable design principle: to enhance the volumetric density of high-affinity H₂ binding sites in BMo-based materials, one should (i) maximize the density of exposed Mo–Mo bridge sites—achievable through nanostructuring (e.g., creating stepped surfaces or nanoparticles) or controlled defect engineering—and (ii) promote molecular orientations with a tilt angle of 60°, possibly via pre-adsorbed steering molecules or electric field alignment.

The electronic fingerprints of this activated state and the proposed ‘gating’ mechanism are quantified in Fig. 7. The Mulliken charge analysis shows a net electron transfer from the H₂ molecule to the surface Mo₄ atom (H₁: −0.128e, H₂: −0.115e, Mo₄: +0.660e), signifying the onset of a polar covalent interaction. Concurrently, the bond-population analysis reveals the dual signature of H₂ activation: the H–H bond population is reduced to 1.07 a.u., a substantial weakening from the gas-phase value (≈ 2.0 a.u.), while the newly formed Mo–H interactions exhibit negative populations (−0.33 and −0.21 a.u.), characteristic of incipient bonds in a transition-state-like configuration. This electronic picture demonstrates that the specific 60° tilt orientation does not merely place the molecule closer to the surface; it geometrically aligns the H–H σ-orbital to enable favorable overlap with metal d-states. This alignment facilitates charge transfer and weakens the H–H bond, acting as an orientation-activated ‘gate’ that lowers the kinetic barrier for dissociation. Thus, the identified ‘hot spot’ is not a geometric accident but an orientation-activated electronic configuration that serves as a precursor to dissociative chemisorption.

This study provides a fundamental, atomic-scale blueprint for H₂ interaction with BMo surfaces. While a comprehensive assessment of practical storage performance requires investigation of higher coverages, temperature/pressure effects, and direct comparison with other materials—which are crucial next steps—our work establishes the essential structure–property relationship: optimal adsorption occurs via a specific orientation-activated mechanism at Mo–Mo bridge sites. This mechanistic insight provides a concrete design principle: future efforts to enhance BMo's hydrogen storage capacity should focus on engineering surfaces that maximize the density or accessibility of these specific bridge sites, potentially through doping, defect creation, or nanostructuring. The identified −0.25 eV adsorption energy for the best configurations positions this material within the promising range for reversible storage, warranting the more applied investigations suggested by the reviewer. These findings suggest that engineering BMo-based materials for hydrogen storage should focus not only on increasing surface area but also on creating surfaces that promote the favorable 60° tilt orientation. This could be achieved through surface patterning, stepped structures, or external field alignment.

4. Conclusion

This comprehensive high-throughput DFT study unravels the complex adsorption behavior of H₂ on BMo(001) surfaces through combinatorial screening of 1728 configurations, yielding fundamental insights for hydrogen storage material design. The remarkable surface homogeneity evidenced by negligible spatial correlations (|r| < 0.07) ensures consistent adsorption regardless of molecular placement, representing a significant advantage for practical applications.

Our analysis reveals sophisticated orientation-space relationships where tilt angle moderately influences adsorption energy (r ≈ −0.30) while azimuth angle controls spatial distribution patterns without affecting global averages. This resolves the apparent contradiction between weak azimuth correlations and dramatic spatial reorganization, with configuration (angle = 60°, tilt = 60°, azimuth = 90°) identified as globally optimal for H₂ adsorption.

The complex energy distributions across all orientation parameters reveal a “gated” adsorption landscape: while the majority of configurations are energetically unfavorable, specific molecular orientations create “hot spots” capable of spontaneous chemisorption with energies as low as −0.25 eV. Our atomic-scale analysis reveals that these sites correspond specifically to configurations where the H₂ molecule is positioned over a Mo–Mo bridge with a tilt angle of 60°, facilitating optimal orbital overlap and charge transfer. While the precise atomic-scale origin of these sites—such as specific surface defects, local Mo/B atom arrangements, or enhanced electronic density of states—requires further investigation, including more advanced electronic structure analyses (e.g., projected density of states or charge density difference) to elucidate the underlying orbital interactions, this mechanistic understanding provides a concrete design principle: enhancing the density of accessible Mo–Mo bridge sites (e.g., via nanostructuring or controlled defect creation) and promoting favorable molecular orientations should directly increase the hydrogen storage capacity of BMo-based materials.

These findings establish strategic design principles for BMo-based hydrogen storage: homogeneous surfaces enable predictable performance, specific molecular orientations optimize adsorption strength, and azimuth control allows spatial patterning without altering overall capacity. The systematic combinatorial screening approach presented here offers a robust framework for mapping orientation-dependent adsorption landscapes, which can be extended to other gas–surface systems.

Looking forward, these computational predictions invite experimental validation and exploration of BMo surface functionalization strategies to enhance high-affinity site density. The fundamental insights from this study lay groundwork for rational design of advanced hydrogen storage materials based on BMo and related transition metal compounds, positioning this material family as promising candidates for sustainable energy applications.

While this study maps the orientation-dependent adsorption landscape at low coverage, several logical extensions follow. The immediate next step is full geometry optimization of the identified “hot-spot” configurations to confirm the predicted energy minima and characterize the fully relaxed chemisorbed state. Beyond this, future investigations should address higher H₂ coverages, temperature and pressure effects, and adsorption/desorption dynamics. Methodologically, future high-throughput screenings could enhance computational efficiency by employing symmetry-adapted grids that avoid redundant sampling of equivalent surface points. Finally, experimental validation of these predictions and comparison with other boride materials are essential to establish BMo's practical potential for hydrogen storage.

Author contributions

Cengiz Soykan: supervision, writing – original draft, methodoloy, investigation, writing – python code, dft calculation, formal analysis, validation, writing – review & editing, conceptualization.

Conflicts of interest

There are no conflicts to declare.

Data availability

Reasonable requests for the data supporting this study, including the optimized structures and the adsorption energy dataset for all 1728 configurations, will be fulfilled by the corresponding author. The Python scripts developed for configuration generation and data analysis are also available upon reasonable request.

Supplementary information (SI) is available. See DOI: https://doi.org/10.1039/d5cp04390e.

Acknowledgements

The open access publication fee was covered by TUBITAK ULAKBIM through the EKUAL agreement.

References

1. J. D. Holladay, J. Hu, D. L. King and Y. Wang, An overview of hydrogen production technologies, Catal. Today, 2009, 139(4), 244–260,  DOI:10.1016/j.cattod.2008.08.039.
2. Agency, I. E. Global Hydrogen Review 2021. 2021. https://www.iea.org/reports/global-hydrogen-review-2021 (accessed).
3. M. S. Dresselhaus and I. L. Thomas, Alternative energy technologies, Nature, 2001, 414(6861), 332–337,  DOI:10.1038/35104599.
4. A. Midilli, M. Ay, I. Dincer and M. A. Rosen, On hydrogen and hydrogen energy strategies: I: current status and needs, Renewable Sustainable Energy Rev., 2005, 9(3), 255–271,  DOI:10.1016/j.rser.2004.05.003.
5. L. Schlapbach and A. Züttel, Hydrogen-storage materials for mobile applications, Nature, 2001, 414(6861), 353–358,  DOI:10.1038/35104634.
6. S. S. Makridis, Hydrogen storage and compression, Methane and Hydrogen for Energy Storage, 2024, pp. 1–28 DOI:10.1049/PBPO101E_ch1.
7. A. R. Oganov, J. Chen, C. Gatti, Y. Ma, Y. Ma, C. W. Glass, Z. Liu, T. Yu, O. O. Kurakevych and V. L. Solozhenko, Ionic high-pressure form of elemental boron, Nature, 2009, 457(7231), 863–867,  DOI:10.1038/nature07736.
8. C. S. Rout, P. V. Shinde, A. Patra and S. M. Jeong, Recent Developments and Future Perspectives of Molybdenum Borides and MBenes, Adv. Sci., 2024, 11(21), 2308178,  DOI:10.1002/advs.202308178.
9. L. Zhang, P. Liang, H.-B. Shu, X.-L. Man, F. Li, J. Huang, Q.-M. Dong and D.-L. Chao, Borophene as Efficient Sulfur Hosts for Lithium–Sulfur Batteries: Suppressing Shuttle Effect and Improving Conductivity, J. Phys. Chem. C, 2017, 121(29), 15549–15555,  DOI:10.1021/acs.jpcc.7b03741.
10. T. Oku, Hydrogen Storage in Boron Nitride and Carbon Nanomaterials, Energies, 2015, 8(1), 319–337.
11. M. Zhang, H. Wang, H. Wang, T. Cui and Y. Ma, Structural Modifications and Mechanical Properties of Molybdenum Borides from First Principles, J. Phys. Chem. C, 2010, 114(14), 6722–6725,  DOI:10.1021/jp100225c.
12. S. J. Clark, M. D. Segall, C. J. Pickard, P. J. Hasnip, M. I. J. Probert, K. Refson and M. C. Payne, First principles methods using CASTEP, Z. Kristallogr. - Cryst. Mater., 2005, 220(5–6), 567–570,  DOI:10.1524/zkri.220.5.567.65075.
13. 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(1), 15–50,  DOI:10.1016/0927-0256(96)00008-0.
14. 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(16), 11169–11186,  DOI:10.1103/PhysRevB.54.11169 Scopus.
15. G. Kresse and J. Hafner, Norm-conserving and ultrasoft pseudopotentials for first-row and transition elements, J. Phys.: Condens. Matter, 1994, 6(40), 8245,  DOI:10.1088/0953-8984/6/40/015.
16. J. P. Perdew, K. Burke and M. Ernzerhof, Generalized gradient approximation made simple, Phys. Rev. Lett., 1996, 77(18), 3865.
17. S. Grimme, Semiempirical GGA-type density functional constructed with a long-range dispersion correction, J. Comput. Chem., 2006, 27(15), 1787–1799,  DOI:10.1002/jcc.20495.
18. E. R. McNellis, J. Meyer and K. Reuter, Azobenzene at coinage metal surfaces: Role of dispersive van der Waals interactions, Phys. Rev. B:Condens. Matter Mater. Phys., 2009, 80(20), 205414,  DOI:10.1103/PhysRevB.80.205414.
19. H. J. Monkhorst and J. D. Pack, Special points for Brillouin-zone integrations, Phys. Rev. B, 1976, 13(12), 5188–5192,  DOI:10.1103/PhysRevB.13.5188.
20. T. H. Fischer and J. Almlof, General methods for geometry and wave function optimization, J. Phys. Chem., 1992, 96(24), 9768–9774,  DOI:10.1021/j100203a036.
21. J. K. Nørskov, J. Rossmeisl, A. Logadottir, L. Lindqvist, J. R. Kitchin, T. Bligaard and H. Jónsson, Origin of the Overpotential for Oxygen Reduction at a Fuel-Cell Cathode, J. Phys. Chem. B, 2004, 108(46), 17886–17892,  DOI:10.1021/jp047349j.
22. G. Wu, J. Zhang, Y. Wu, Q. Li, K. Chou and X. Bao, First-Principle Calculations of the Adsorption, Dissociation and Diffusion of Hydrogen on the Mg(0001) Surface, Acta Phys.-Chim. Sin., 2008, 24(1), 55–60,  DOI:10.1016/S1872-1508(08)60006-6.
23. A. Züttel, Materials for hydrogen storage, Mater. Today, 2003, 6(9), 24–33,  DOI:10.1016/S1369-7021(03)00922-2.
24. S. K. Bhatia and A. L. Myers, Optimum Conditions for Adsorptive Storage, Langmuir, 2006, 22(4), 1688–1700,  DOI:10.1021/la0523816.
25. E. Skúlason, G. S. Karlberg, J. Rossmeisl, T. Bligaard, J. Greeley, H. Jónsson and J. K. Nørskov, Density functional theory calculations for the hydrogen evolution reaction in an electrochemical double layer on the Pt(111) electrode, Phys. Chem. Chem. Phys., 2007, 9(25), 3241–3250,  10.1039/B700099E.
26. J. Greeley and M. Mavrikakis, Alloy catalysts designed from first principles, Nat. Mater., 2004, 3(11), 810–815,  DOI:10.1038/nmat1223.
27. A. Michaelides, Z. P. Liu, C. J. Zhang, A. Alavi, D. A. King and P. Hu, Identification of General Linear Relationships between Activation Energies and Enthalpy Changes for Dissociation Reactions at Surfaces, J. Am. Chem. Soc., 2003, 125(13), 3704–3705,  DOI:10.1021/ja027366r.
28. G. S. Karlberg, J. Rossmeisl and J. K. Nørskov, Estimations of electric field effects on the oxygen reduction reaction based on the density functional theory, Phys. Chem. Chem. Phys., 2007, 9(37), 5158–5161,  10.1039/B705938H.
29. J. K. Nørskov, T. Bligaard, J. Rossmeisl and C. H. Christensen, Towards the computational design of solid catalysts, Nat. Chem., 2009, 1(1), 37–46,  DOI:10.1038/nchem.121.
30. J. K. Nørskov, F. Abild-Pedersen, F. Studt and T. Bligaard, Density functional theory in surface chemistry and catalysis, Proc. Natl. Acad. Sci. U. S. A., 2011, 108(3), 937–943,  DOI:10.1073/pnas.1006652108.

---