Aromatic M₁₂B₆₀ (M = Y, Lu) metallo-borospherenes for reversible hydrogen storage

Abstract

Metallo-borospherenes, a rapidly evolving class of boron-based nanostructures, exhibit intriguing diversity in both geometric architectures and electronic properties. Recently, the experimental identification of the D_3h-symmetric Ln₃B₁₈⁻ (Ln = La, Tb) clusters marked a significant milestone in this emerging field. In this study, a new family of metallo-borospherenes M₁₂B₆₀ (M = Y, Lu) was predicted using first-principles calculations. Ab initio molecular dynamics simulations and vibrational frequency analyses confirm their thermodynamic and kinetic stability. Detailed bonding analysis reveals the presence of 108 delocalized multi-center two-electron bonds within the cage, accompanied by pronounced aromatic character, which collectively account for their exceptional stability. Notably, the Y₁₂B₆₀ cage exhibits potential as a promising hydrogen storage material. Its curved surface provides multiple favorable adsorption sites, enabling H₂ molecules to form layered distributions. The Y₁₂B₆₀ cage can adsorb up to 89 H₂ molecules with an average adsorption energy of −0.201 eV per H₂, corresponding to a gravimetric density of 9.44 wt%, meeting practical adsorption requirements. Ab initio molecular dynamics simulations further validate the structural robustness and adsorption capacity of Y₁₂B₆₀ under near-ambient conditions, underscoring its promise for reversible hydrogen storage applications.

---

Introduction

The electron deficiency of boron facilitated the formation of multi-center two-electron bonds (mc-2e), resulting in a wide variety of structural motifs and electronic properties in boron-based nanostructures. Over the past two decades, a combination of photoelectron spectroscopy and first-principles calculations has revealed a rich structural diversity among boron clusters, encompassing planar, quasi-planar, cage-like, bilayer, and core–shell configurations.^1–4

Metal atom doping has emerged as a powerful strategy to compensate for boron's electron deficiency, significantly altering geometries and introducing novel chemical bonding patterns and properties.^5–7 In 2013, inspired by the perfect molecular wheels B₈⁻ and B₉⁻,^8,9 a design principle for metal-centered boron rings (M©B_n) was proposed, where the central boron atom was replaced with a transition metal. This led to the prediction of several dually aromatic species, including D_8h Co©B₈⁻, D_9h Ru©B₉⁻ and D_10h Ta©B₁₀⁻.¹⁰ More recently, this concept has been extended to M^(m)©B_n^k− clusters, culminating in the discovery of the D_11h symmetry YB₁₁²⁻ structure.¹¹ Representative examples of single-metal doping included dually anti-aromatic compound BiB₆⁻ with planar configuration, dually aromatic compound BiB₇⁻ and BiB₈⁻,¹² half-sandwich configuration LaB_n^−/0 (n = 14–17)¹³ and twenty-coordinate boron molecular drums An@B₂₀ (An = U, Np, and Pu).¹⁴

By doping two metal atoms, the icosahedral boron-based clusters B₁₀M₂ (M = Rh, Ir)¹⁵ were first predicted. Additionally, it has been observed that B_n⁻ (n = 10, 11) doped by 2 La atoms could form the anti-sandwich complex La₂B_n⁻ (n = 10, 11), offering a novel pathway for the structural evolution of larger lanthanide boron clusters.¹⁶ A series of bimetal-doped inverse sandwich boron clusters, including Ta₂B₆^−/0, Pr₂B₈⁻, and La₂B_n⁻ (n = 7–9), as well as B₁₀La₂⁻ and B₁₁La₂⁻ extended inverse-sandwich complexes, have also been experimentally observed.^16–18

As the number of doped metal and boron atoms increased, the geometric configuration gradually evolved to metallo-borospherenes. The experimental discovery of trimetal-doped metallo-borospherenes D_3h La₃B₁₈⁻ marked a breakthrough advancement in the field of metallo-borospherenes.¹⁹ Subsequent studies revealed its potential as a hydrogen storage material, capable of adsorbing 18 H₂ molecules with a gravimetric density of 5.6 wt%.²⁰ Inspired by the D_3h La₃B₁₈⁻ structure, a core–shell metallo-borospherene, D_3h La₃&[B₂@B₁₈]⁻, sharing the same shell configuration was proposed.²¹ Further exploration with increased metal doping led to the proposal of tetrahedral cage T_d La₄B₂₄ (exhibiting spherical aromaticity) and core–shell structures T_d La₄B₂₉^0/+/−, and La@[La₅&B₃₀]^0/−/2−, enriching the family of metallo-borospherenes.^22,23 Most recently, spherical trihedral metallo-borospherenes B₂₀TM_n (TM = Sc, Y; n = 3, 4) were reported. Their B₂₀ framework consisted of three B₂ units bridging two equivalent B₇ triangles, providing a B₇ triangular motif for the B-framework.²⁴

Metal-doped boron nanomaterials have been extensively studied for hydrogen storage. Commonly investigated metals included alkali (Li, Na, and K),^25–28 alkaline earth (Be, Mg, and Ca),^29–32 transition metals (Sc, Ti, Y, Fe, Co, and Ni),^33–38 and rare earth elements (La, Ce, Nd, Gd, and Er).^20,39–41 Notably, Ti-decorated B₃₈ cages can adsorb up to 6 H₂ molecules per Ti atom, achieving a gravimetric capacity of 7.44 wt%.⁴² Recently, Y₃C₃B₉H₁₂ was found to bind 14 H₂ molecules, yielding a storage capacity of 6.4 wt%.³⁷ Theoretical predictions indicated that C₄₈B₁₂ decorated with 12 Y atoms is an excellent hydrogen storage material. Each Y atom adsorbed 6 H₂ molecules primarily through Kubas interactions, with an average adsorption energy of −0.64 eV per H₂ and a corresponding gravimetric density of 7.51 wt%.⁴³ The hydrogen storage behavior of multi-yttrium-doped B₃₈ systems has been explored using first-principles calculations. These studies revealed that Y atoms preferentially occupied the hexagonal cavities of B₃₈ rather than forming clusters, effectively preventing Y atoms from diffusion. These results confirmed that Y₄@B₃₈ is a promising hydrogen storage material, with each Y atom adsorbing up to 6 H₂ molecules, corresponding to a gravimetric density of 4.96 wt%.³⁸

While prior work on metallo-borospherenes has largely concentrated on endohedral and exohedral motifs, investigations of hollow, cage-like structures remained relatively scarce. Notably, the experimentally characterized La₃B₁₈⁻ cage-like motif—conceived from three shared La©B₁₀ units—provided an important proof of concept for constructing larger cage architectures.¹⁹ In addition, previous theoretical studies of yttrium-decorated boron nanostructures (e.g., Y©B₁₁²⁻ molecular wheel,¹¹ B₂₀Y₃/Y₄ metallo-borospherenes,²⁴ and Y₄@B₃₈ hydrogen-storage candidates³⁸) demonstrated that Y can stabilize a variety of boron frameworks and exhibit favorable properties for hydrogen adsorption. These findings motivated us to explore whether a Y©B₁₀ building block can be extended to form a highly symmetric D_2h Y₁₂B₆₀ cage assembled from twelve shared Y©B₁₀ units; subsequent calculations indicated that this structure was energetically and kinetically viable and possessed promising hydrogen storage characteristics.

To probe the role of the metal center, we selected Lu as a comparative element. Y and Lu shared similar valence electron counts (both commonly trivalent), but Lu had a filled 4f shell and experienced a pronounced lanthanide contraction, resulting in a smaller atomic radius. This combination—similar valence electronic structure but different ionic/atomic sizes—provided a controlled way to disentangle size (steric) effects from purely electronic effects on the cage geometry and stability. Consequently, comparing Y₁₂B₆₀ and Lu₁₂B₆₀ yields direct insight into how metal-center size influences structural parameters and hydrogen adsorption behavior.

Methods

Geometry optimizations for the metallo-borospherenes M₁₂B₆₀ (M = Y, Lu) were conducted at the PBE0 theory level using the Gaussian 16 package. In these calculations, the 6-311+G* basis set was applied for B, while the Stuttgart/Dresden (SDD) effective core potential was used for Y, and the ECP60MWB pseudopotential for Lu.^44,45 Vibrational frequency analyses were carried out to confirm that the optimized structures correspond to true local minima on the potential energy surface. To further verify the robustness of the results, complementary calculations were performed using the TPSSh hybrid meta-GGA functional,⁴⁶ thereby minimizing potential bias arising from the limitations of any single functional.

The hydrogen storage properties were evaluated employing the WB97XD functional⁴⁷ in the Gaussian 16 package,⁴⁸ which included empirical dispersion corrections and provided a more reliable description of weak interactions.

Thermodynamic stability was further assessed by performing ab initio molecular dynamics (AIMD) simulations within the canonical (NVT) ensemble at 300 K, as implemented in the DMol³ package.⁴⁹ The simulations were run for a total of 5.0 ps using a time step of 1.0 fs.

Additional property calculations, including Raman spectra, photoelectron spectroscopy (PES), adaptive natural density partitioning (AdNDP) calculations, and interaction region indicator (IRI) were conducted using the Gaussian 16 package.⁴⁸ Data processing and visualization were conducted using Multiwfn⁵⁰ and VMD,⁵¹ respectively.

Results and discussion

Configurations

The predicted M₁₂B₆₀ (M = Y, Lu) structures comprised a B₆₀ cage coordinated with twelve metal atoms. The B₆₀ framework adopted a truncated dodecahedral geometry, consisting of twelve interconnected B₁₀ rings sharing B–B bonds between adjacent units. Twelve metal atoms occupied positions near the centroids of these B₁₀ rings, forming the geodesic truncated dodecahedron. Energy minimization revealed that four boron atoms at opposing poles exhibited slight inward deformation, resulting in a distinctive drum-shaped morphology (Fig. 1).

Fig. 1 Optimized structures (a) Y₁₂B₆₀ and (b) Lu₁₂B₆₀. Boron atoms were shown in pink; yttrium atoms, in cyan; and lutetium atoms, in green.

Following geometry optimization, average bond lengths for the Y₁₂B₆₀ structure were determined. The average Y–B bond length was 2.68 Å. The average B–B bond length was 1.63 Å, intermediate between typical B=B double bonds (about 1.56 Å) and B–B single bonds (about 1.70 Å). For Lu₁₂B₆₀, the average Lu–B bond length was 2.63 Å, while the average B–B bond length was 1.63 Å. Additionally, the average binding energy per atom of Lu₁₂B₆₀ was −5.48 eV per atom, slightly lower than that of Y₁₂B₆₀ (−5.58 eV per atom), indicating lower stability in Lu₁₂B₆₀. This originated from Lu's larger atomic radius, which led to a less optimal match with the B₁₀ rings in the B₆₀ cage. Consequently, Lu₁₂B₆₀ adopted a lower C_s symmetry compared to the D_2h-symmetric Y₁₂B₆₀. Given its superior stability and symmetry, the D_2h Y₁₂B₆₀ structure was the focus of subsequent analysis.

Stability

The thermal stability of the Y₁₂B₆₀ and Lu₁₂B₆₀ structures was evaluated through ab initio molecular dynamics simulations at 800 K, 1000 K, and 1200 K. The simulations employed the NVT ensemble with an NH Chain thermostat.⁵² The total simulation time was set to 5.0 ps, with a time step of 1.0 fs. The results displayed in Fig. 2 showed that at 800 K and 1000 K, the Y₁₂B₆₀ structure underwent slight distortion but maintained a well-preserved cage-like framework. The analysis showed that the total energy and temperature stabilized during the simulation periods without systematic drift (see Fig. S3(a, b, d and e)), supporting the numerical stability of the AIMD trajectories. The RDF computed at 1000 K for Y₁₂B₆₀ displayed clear Y–B and B–B peaks at approximately 2.65 Å and 1.65 Å, respectively, which were in close agreement with the corresponding distances in the optimized Y₁₂B₆₀ geometry (2.68 Å for Y–B and 1.63 Å for B–B), with deviations consistent with expected thermal broadening (Fig. S2). Similar behaviour was observed for Lu₁₂B₆₀ at 800 K. The average root-mean-square deviation (RMSD) of the Y₁₂B₆₀ structure was 0.28 Å at 1000 K, indicating that all atoms were fluctuating around their equilibrium positions, as illustrated in Fig. S3(c). However, at 1200 K, the Y₁₂B₆₀ structure collapsed, accompanied by the destruction of the Y©B₁₀ structural unit on the structure surface. These findings demonstrated that the Y₁₂B₆₀ structure remained stable up to approximately 1000 K, exhibiting good thermal stability. Similarly, the Lu₁₂B₆₀ structure demonstrated good thermal stability, maintaining its structural integrity at about 800 K.

Fig. 2 Representative snapshots from AIMD simulations of the Y₁₂B₆₀ structure at 800 K, 1000 K and 1200 K under the NVT ensemble. The total simulation time was 5.0 ps with a time step of 1.0 fs.

Subsequently, to further assess their relative stability, a set of low-energy isomers extracted from high-temperature AIMD trajectories were fully reoptimized. The energy differences (ΔE) between these isomers and the primary structures are summarized in Fig. S4 and S5 of the SI. For Y₁₂B₆₀, isomers 1 and 2 exhibited slightly lower total energies than the primary D_2h structure, accompanied by reduced symmetry and partial disruption of surface Y©B₁₀ units, which evolved into B₆ nanoribbon-like motifs. In contrast, isomers 3–6, which retained intact Y©B₁₀ units, were slightly higher in energy. A similar energetic ordering was observed for Lu₁₂B₆₀. Although the structures discussed in the manuscript were not the absolute global minima, they represented low-lying metastable isomers characterized by intact Y©B₁₀/Lu©B₁₀ motifs, high symmetry, and pronounced kinetic stability—features that were crucial for potential experimental realization and functional applications.

The Raman spectra and infrared (IR) of the Y₁₂B₆₀ and Lu₁₂B₆₀ structures were simulated at the PBE0 level. While the frequency calculation provided all normal vibrational modes, only those that are IR- or Raman-active appeared as peaks in the corresponding spectra. The simulated spectra therefore served not only as a theoretical reference for future experimental characterization but also as a complementary verification of the stability of the predicted structures. Raman spectra (Fig. S6) revealed the frequencies of Y₁₂B₆₀ and Lu₁₂B₆₀ structures spanning 39.7 to 1405.1 cm⁻¹ and 39.6 to 1427.1 cm⁻¹, respectively, with no imaginary frequencies presented, confirming the good kinetic stability. For the Y₁₂B₆₀ structure, major vibration modes corresponding to specific spectral features are detailed in Fig. S7. The dominant peak at 61.4 cm⁻¹ corresponded to the global breathing mode of the entire cage structure, while prominent bands observed at 279.6, 603.4, 937.9, and 1346.4 cm⁻¹ originated from localized stretching vibrations of B₃ triangles. For infrared (IR) spectra (Fig. S6), the Y₁₂B₆₀ structure exhibited several major asymmetric infrared activity peaks. Among them, the first infrared peak at 104.0 cm⁻¹ and the highest peak at 845.6 cm⁻¹ mainly corresponded to localized stretching vibrations of B₃ triangles.

Photoelectron spectroscopy (PES) is a powerful analytical technique that probes the elemental composition, chemical state, and electronic structure of matter by measuring the kinetic energies of electrons ejected upon irradiation with photons, with particular sensitivity to surface properties. The simulated PES offers valuable predictive insights for subsequent experimental studies. The simulated photoelectron spectrum of the Y₁₂B₆₀⁻ structure (Fig. S6) displayed distinct features that were expected to be readily detectable in experimental measurements. The first vertical detachment energy (VDE) was calculated to be 1.56 eV, which corresponded to the energy difference between the anionic ground state and the neutral species at the anion geometry. The first adiabatic detachment energy (ADE) was found to be 1.44 eV, representing the energy difference between the lowest potential energy surfaces of the anion and the neutral molecule. The small difference between the VDE and ADE (0.12 eV) suggested that the geometric structures of the neutral molecule and the anion in their ground states were likely similar. Therefore, it was inferred that the first ionization band may exhibit a relatively sharp profile, possibly with vibrational fine structure. Furthermore, the noticeable intervals (>0.2 eV) among the computed VDE values, including 1.82, 2.10, 2.70, 3.02, 3.32, 3.67, 4.11, 4.87, 5.14, 5.42, 5.65 and 6.06 eV, indicated that the corresponding electronic states were expected to be resolved as distinct bands in the experimental spectrum. These predicted spectral features may provide useful theoretical guidance for future characterization of the Y₁₂B₆₀⁻ structure using photoelectron spectroscopy in experimental studies.

Recent studies have demonstrated that the spherical trihedral metallo-borospherene La₃B₁₈⁻ has been successfully realized through a combination of theoretical prediction and experimental photoelectron spectroscopy, where three La©B₁₀ structural units formed the stable cage.¹⁹ In an analogous manner, the Y₁₂B₆₀ (Lu₁₂B₆₀) cluster can be regarded as being constructed from twelve such Y©B₁₀ (Lu©B₁₀) subunits, thereby exhibiting both good thermodynamic stability and strong structural similarity to the experimentally confirmed La₃B₁₈⁻.¹⁹

Based on structural analogies to previously reported metallo-borospherenes, it was conceivable that M₁₂B₆₀ (M = Y, Lu) could potentially be generated experimentally through laser ablation of a Y/B (Lu/B) mixed target, followed by supersonic expansion in a helium carrier gas—the same method employed in the synthesis of La₃B₁₈⁻.¹⁹

Electronic properties

The partial density of states (PDOS) and molecular frontier orbitals of the Y₁₂B₆₀ structure provided detailed insights into the contributions of individual atoms or atomic orbitals to molecular orbitals. As shown in Fig. 3, the HOMO–LUMO gap of 1.18 eV in the Y₁₂B₆₀ structure indicated good chemical inertness. These molecular orbitals primarily originated from hybridization of B 2s, 2p orbitals and Y 4d, 5s orbitals. The HOMO level was primarily contributed by B 2p orbitals, followed by Y 4d and B 2s orbitals. The LUMO showed dominant contributions from B 2p orbitals, alongside participation from Y 4d, 5s orbitals. These characteristic features were clearly visualized in the molecular frontier orbitals in Fig. S8. Notably, the Y 4d orbital contribution increased above the LUMO level, suggesting enhanced participation of Y orbitals in higher-energy unoccupied states.

Fig. 3 Total density of states (TDOS) and projected density of states (PDOS) of the Y₁₂B₆₀ structure. The HOMO and LUMO energy levels were indicated by pink and orange dashed lines, respectively.

The deformation electron density was analyzed to clarify the electronic properties of the Y₁₂B₆₀ structure (Fig. 4). Regions of charge accumulation and depletion were colored in blue and yellow, respectively. The results revealed significant charge depletion around the metal atoms and B atoms. Charge accumulation occurred between the Y and B atoms, as well as among the B atoms, indicating the presence of multi-center two-electron bonds. This finding was consistent with the subsequent AdNDP analysis. The Hirshfeld charge analysis further confirmed an average transfer of 0.69e from the Y atom to the B atom.

Fig. 4 Deformation electron density for the Y₁₂B₆₀ structure, with an isosurface of 0.025 e Å⁻³. The blue and yellow represented the accumulation and depletion of electrons, respectively.

Adaptive natural density partitioning (AdNDP) was a method of chemical bonding analysis that described the electronic structure of a molecule by decomposing its electron density (typically derived from computations) into a set of chemically meaningful units, such as lone pairs, two-center two-electron bonds, and multi-center bonds. Therefore, to further investigate the bonding properties of the M₁₂B₆₀ (M = Y, Lu) structures, AdNDP bonding analysis was conducted. A total of 108 multi-center two-electron bonds were identified within the Y₁₂B₆₀ framework. As shown in Fig. 5, Y₁₂B₆₀ possessed thirty localized 2c–2e σ bonds on the B–B bonds shared by adjacent B₁₀ rings. Twenty 3c–2e σ bonds were located on the B₃ triangular units, and another forty 3c–2e σ bonds on the two B atoms and one Y atom, where the B atoms originated from the two B atoms in the B₃ triangles. Fourteen 8c–2e σ bonds are located between the six B atoms and the two adjacent Y atoms. At the drum-face positions of the Y₁₂B₆₀ structure, there were four 9c–2e σ bonds between the six B atoms and the three Y atoms. The existence of these covalent bonds enabled the system to maintain higher stability. Previous studies have indicated that the Lu 4f orbitals were fully occupied and their energy levels were significantly lower than those of the 5d orbitals, rendering them largely non-participatory in bonding interactions.⁵³ Therefore, the Lu₁₂B₆₀ structure has similar bonding characteristics, as shown in Fig. S12.

Fig. 5 AdNDP analysis of the Y₁₂B₆₀ structure. The occupation numbers were indicated for each bonding pattern, and the isosurface was set to 0.05 e Å⁻³.

The nucleus-independent chemical shifts (NICSs) provided a reliable assessment of molecular aromaticity. The NICS(0) value of the Y₁₂B₆₀ structure, calculated at the PBE0 level, was −18.19 ppm, confirming the Y₁₂B₆₀ structure's aromatic character. The calculated ICSS_ZZ of the Y₁₂B₆₀ structure reflected the isosurface of the ZZ component of the NICS, providing a better understanding of aromaticity. The isochemical shielding surface (ICSS) method can be considered as a three-dimensional extension of the nucleus-independent chemical shift (NICS) concept. By calculating the magnetic shielding values at grid points in space and visualizing the corresponding isosurfaces, this approach offered a more intuitive depiction of the spatial distribution of magnetic response. In particular, ICSS_ZZ denoted the isosurface associated with the ZZ component of the NICS. A principal advantage of the ZZ component lies in its enhanced capacity to differentiate between long-range magnetic effects contributed by delocalized π-electrons and local magnetic anisotropy arising from the σ-framework. This distinction facilitated the identification of magnetically shielded regions that were indicative of aromatic character. As shown in Fig. 6, distinct green regions represented the chemical shielding area, while the yellow regions indicated the chemical de-shielding area. The analysis revealed a chemical shielding region along the radial z direction and a chemical de-shielding area in the horizontal direction surrounding the structure. These results collectively confirmed the aromaticity of the Y₁₂B₆₀ structure, which represents a key factor underpinning the good stability.

Fig. 6 Side view (a) and top view (b) of the ICSS_ZZ for the Y₁₂B₆₀ structure. Green and yellow represented chemical shielding and de-shielding areas, respectively.

Hydrogen storage properties

To evaluate the hydrogen storage capabilities of the Y₁₂B₆₀ structure, four distinct adsorption sites were selected: above the Y atom (A), above the Y–B bond (B₁), above the triangular vacancy formed by Y–B–B (B₂), and above the B₃ triangular vacancy (T), as shown in Fig. S14. Calculations yielded adsorption energies of −0.243, −0.298, −0.262 and −0.192 eV for sites A, B₁, B₂ and T, respectively, with corresponding desorption temperatures of 310 K, 380 K, 334 K, and 245 K based on Van't Hoff equation, indicating potential suitability for hydrogen storage applications.

To determine the hydrogen storage capacity of the Y₁₂B₆₀ structure, 150 H₂ molecules were placed around it and structural optimization was conducted. The results revealed that the H₂ molecules formed distinct adsorption layers (Fig. 7 and 8). The first adsorption layer, located 6.60–7.00 Å from the structure center, contained about 38 H₂ molecules with an average adsorption energy of −0.253 eV per H₂. The second adsorption layer was positioned 7.30–7.81 Å from the Y₁₂B₆₀ structure center (corresponding to an interlayer distance of 0.70–0.81 Å relative to the first layer), which was significantly shorter than the 2.7 Å equilibrium distance between free H₂ molecules, thereby confirming the adsorption of H₂ molecules by the Y₁₂B₆₀ structure. Two layers accommodated about 89 H₂ molecules with an adsorption energy of −0.201 eV per H₂, achieving a gravimetric hydrogen storage density of 9.44 wt% and a corresponding desorption temperature of 256 K based on the Van’t Hoff equation. Furthermore, the structural schematic diagram of the Y₁₂B₆₀-38 H₂ (the first layer) and the Y₁₂B₆₀-89 H₂ systems has been provided to more intuitively display the molecular arrangements within each layer (Fig. S15).

Fig. 7 The distributions of the number density of B, Y, and H atoms, and the total number of H₂ molecules, as a function of distance from the Y₁₂B₆₀ structure center.

Fig. 8 The distributions of average binding energy of per H₂ molecule and total number of H₂ molecules, as a function of distance from the structure center. Here, the average binding energy of per H₂ was calculated as E_b = {[E(Y₁₂B₆₀)-(H₂)_n] − E(Y₁₂B₆₀) − E(H₂) × n}/n.

Quantum theory of atoms in molecules (QTAIM) provided a robust framework for analysing chemical bonding. To investigate hydrogen adsorption between H₂ molecules and the Y₁₂B₆₀ structure, we performed QTAIM topological analysis at four adsorption sites listed in Table S3. For non-covalent interactions, the Laplacian of electron density ∇²ρ(r) exhibited positive values. Higher electron density values ρ(r) correlated with more negative potential energy density V(r), signifying stronger interactions. Obviously, positive ∇²ρ(r) values were characteristic of non-covalent interactions between H₂ and the Y₁₂B₆₀ framework. The adsorption site above the Y–B bond showed the strongest binding, consistent with the calculated adsorption energies.

The interaction region indicator (IRI) was a real-space function derived from the electron density and its derivatives. It was employed to visualize and qualitatively analyse diverse chemical interactions in molecular systems, encompassing covalent bonds, hydrogen bonds, van der Waals forces, and steric repulsions. Accordingly, to further investigate the hydrogen adsorption mechanism between H₂ molecules and the Y₁₂B₆₀ structure, we conducted interaction region indicator (IRI) and partial density of states (PDOS) analysis for the Y₁₂B₆₀-89 H₂ system (Fig. 9 and Fig. S16, S17). Weak interactions in Fig. 9 regions were coloured in green. Between the Y₁₂B₆₀ structure and the adsorbed 89 H₂ molecules were surrounded by green regions, confirming that hydrogen adsorption was primarily mediated by weak interactions. This finding was supported by PDOS analysis (Fig. S16), which revealed weak hybridization between the H 1s orbitals and Y 5d orbitals within the −14.5 to −11 eV energy range. Furthermore, no significant hybridization peaks (Fig. S17) were observed between the Y 5d orbitals and H 1s orbitals near the Fermi level, providing additional evidence that the adsorption between Y₁₂B₆₀ and H₂ molecules occurred via van der Waals interactions. Most importantly, the H–H bond lengths for all adsorbed H₂ molecules in both layers remained in the range of 0.75–0.76 Å. Although slightly elongated compared to that of the free H₂ molecule (0.74 Å), no significant bond cleavage occurred, providing key evidence for a mechanism dominated by physical adsorption.

Fig. 9 The geometric structure (a) and IRI analysis (b) of the Y₁₂B₆₀-89 H₂ complex. The weak interaction regions were coloured green.

To elucidate the desorption mechanism of hydrogen molecules, we performed ab initio molecular dynamics simulations on the Y₁₂B₆₀-89 H₂ system at different temperatures (250 K, 300 K, 350 K, 400 K, 500 K, 600 K), with corresponding statistical results presented in Fig. 10. The simulation duration was set to 5.0 ps with a time step of 0.5 fs to accurately capture the high frequency vibrations of hydrogen molecules. Besides, the analysis of the energy and temperature evolution confirmed that the system reached a well-equilibrated state (Fig. S18), indicating the reliability of the computational results.

Fig. 10 The number of adsorbed H₂ molecules as a function of AIMD simulation temperature for the Y₁₂B₆₀-89 H₂ complex.

The simulations revealed a clearly temperature-dependent adsorption/desorption behavior. At 250 K, the system exhibited a low desorption rate, retaining ∼84 H₂ molecules (8.97 wt%). Increasing the temperature to 300 K reduced the adsorption to ∼66 H₂ molecules (7.19 wt%), while at 350 K the number decreased further to ∼53 H₂ molecules (5.85 wt%). At 600 K, nearly complete desorption occurred, with only ∼8 H₂ molecules remaining adsorbed. These results indicated that the desorption process was effectively activated with increasing temperature.

Importantly, according to the U.S. Department of Energy (DOE) target, a practical hydrogen storage material should achieve a gravimetric capacity of ≥5.5 wt% near ambient temperature (∼300 K). The Y₁₂B₆₀ cluster maintains a capacity above this threshold within the 250–350 K range, demonstrating its potential for practical hydrogen storage applications. The temperature-activated desorption further supported the material's applicability.

Conclusions

In summary, metallo-borospherenes M₁₂B₆₀ (M = Y, Lu) have been designed and investigated through comprehensive first-principles calculations. Vibrational frequency analysis confirmed their kinetic stability. Molecular dynamics simulations within the NVT ensemble demonstrated that both Y₁₂B₆₀ and Lu₁₂B₆₀ maintained structural integrity at 1000 K and 800 K, respectively. The high stability of these cages was attributed to the presence of 108 delocalized multi-center two-electron covalent bonds, together with pronounced aromaticity. Furthermore, the research revealed multiple potential hydrogen adsorption sites on the Y₁₂B₆₀ surface, enabling the adsorption of about 89 H₂ molecules. These adsorbed H₂ molecules exhibited a layered distribution around the structure, with a corresponding adsorption energy of −0.201 eV per H₂ and a gravimetric hydrogen storage density of 9.44 wt%, surpassing that of La₃B₁₈-18 H₂ (5.6 wt%)²⁰ and Y₄@B₃₈-24 H₂ (4.96 wt%).³⁸Ab initio molecular dynamics simulations demonstrated that the Y₁₂B₆₀ structure may be a feasible hydrogen storage material under near-ambient conditions.

Author contributions

Conceptualization: Hui-Yan Zhao and Ying Liu; writing – original draft: Yi-Sha Chen; writing – review & editing: Yi-Sha Chen, Jing-Jing Guo, Peng-Bo Liu, Hui-Yan Zhao, and Jing Wang; validation: Hui-Yan Zhao and Ying Liu; investigation: Hui-Yan Zhao and Ying Liu; methodology: Yi-Sha Chen, Hui-Yan Zhao and Ying Liu; funding acquisition: Yi-Sha Chen, Jing-Jing Guo, Peng-Bo Liu and Ying Liu.

Conflicts of interest

The authors declare no conflicts of interest.

Data availability

The data supporting this article have been included as part of the supplementary information (SI). Supplementary information: the results of AIMD simulations of the Lu₁₂B₆₀ structure at different temperatures; radial distribution functions (RDF), energy, temperature, and the root-mean-square deviation (RMSD) varied with the simulation time for M₁₂B₆₀ (M = Y, Lu) structures; the initial and optimized M₁₂B₆₀ (M = Y, Lu) isomers; PES, Raman, IR, vibration modes and molecular frontier orbitals of the M₁₂B₆₀ (M = Y, Lu) structures; TDOS, PDOS, AdNDP and ICSS_ZZ for Lu₁₂B₆₀ structure; several possible locations of absorbed H₂ molecules on the Y₁₂B₆₀ structure; schematic diagram of Y₁₂B₆₀-38 H₂ system and Y₁₂B₆₀-89 H₂ system; TDOS and PDOS of the Y₁₂B₆₀ structure and Y₁₂B₆₀-89 H₂ system; energy and temperature varied with the simulation time for Y₁₂B₆₀-89 H₂ system; the structural and electronic properties of M₁₂B₆₀ (M = Y, Lu) structures; the results of QTAIM topological analysis, the distance of H₂ to the nearest atoms (d_H-M), adsorption energies (E_ad) and the desorption temperature (T_d) of H₂ molecule on the four specific adsorption locations.  See DOI: https://doi.org/10.1039/d5cp03352g.

Acknowledgements

This work is supported by the National Natural Science Foundation of China (Grant No. 12174084 and 12404338), the Doctoral Research Start-up Foundation of Hebei Normal University (Grant No. L2023B08 and L2025B12), and the Innovation Funding Project for Doctoral Students of Hebei Province (Grant No. CXZZBS2025144).

References

1. A. P. Sergeeva, D. Y. Zubarev, H. J. Zhai, A. I. Boldyrev and L. S. Wang, J. Am. Chem. Soc., 2008, 130, 7244–7246.
2. W. L. Li, C. Romanescu, T. Jian and L. S. Wang, J. Am. Chem. Soc., 2012, 134, 13228–13231.
3. A. P. Sergeeva, Z. A. Piazza, C. Romanescu, W. L. Li, A. I. Boldyrev and L. S. Wang, J. Am. Chem. Soc., 2012, 134, 18065–18073.
4. Q. Chen, G. F. Wei, W. J. Tian, H. Bai, Z. P. Liu, H. J. Zhai and S. D. Li, Phys. Chem. Chem. Phys., 2014, 16, 18282–18287.
5. T. Jian, X. N. Chen, S. D. Li, A. I. Boldyrev, J. Li and L. S. Wang, Chem. Soc. Rev., 2019, 48, 3550–3591.
6. J. Barroso, S. Pan and G. Merino, Chem. Soc. Rev., 2022, 51, 1098–1123.
7. W. L. Li, X. Chen, T. Jian, T. T. Chen, J. Li and L. S. Wang, Nat. Rev. Chem., 2017, 1, 0071.
8. H. J. Zhai, A. N. Alexandrova, K. A. Birch, A. I. Boldyrev and L. S. Wang, Angew. Chem., Int. Ed., 2003, 42, 6004–6008.
9. L. L. Pan, J. Li and L. S. Wang, J. Chem. Phys., 2008, 129, 024302.
10. C. Romanescu, T. R. Galeev, W. L. Li, A. I. Boldyrev and L. S. Wang, Acc. Chem. Res., 2013, 46, 350–358.
11. X. R. Dong, J. X. Zhang, T. T. Chen, C. Q. Xu and J. Li, Inorg. Chem., 2024, 63, 6276–6284.
12. W. J. Chen, M. Kulichenko, H. W. Choi, J. Cavanagh, D. F. Yuan, A. I. Boldyrev and L. S. Wang, J. Phys. Chem. A, 2021, 125, 6751–6760.
13. X. Q. Lu, Z. H. Wei and S. D. Li, J. Cluster Sci., 2022, 33, 1467–1474.
14. J. Wang, N. X. Zhang, C. Z. Wang, Q. Y. Wu, J. H. Lan, Z. F. Chai, C. M. Nie and W. Q. Shi, Phys. Chem. Chem. Phys., 2021, 23, 26967–26973.
15. W. Y. Liang, J. Barroso, S. Jalife, M. Orozco-Ic, X. Zarate, X. Dong, Z. H. Cui and G. Merino, Chem. Commun., 2019, 55, 7490–7493.
16. Z. Y. Jiang, T. T. Chen, W. J. Chen, W. L. Li, J. Li and L. S. Wang, J. Phys. Chem. A, 2021, 125, 2622–2630.
17. W. L. Li, T. T. Chen, D. H. Xing, X. Chen, J. Li and L. S. Wang, Proc. Natl. Acad. Sci. U. S. A., 2018, 115, E6972–E6977.
18. T. T. Chen, W. L. Li, J. Li and L. S. Wang, Chem. Sci., 2019, 10, 2534–2542.
19. T. T. Chen, W. L. Li, W. J. Chen, X. H. Yu, X. R. Dong, J. Li and L. S. Wang, Nat. Commun., 2020, 11, 2766.
20. X. Y. Sun, P. F. Yin, Y. Zhang, C. Y. Zhang, X. Feng and G. Jiang, Int. J. Hydrogen Energy, 2023, 48, 7807–7813.
21. X. Y. Zhao, M. Yan, Z. H. Wei and S. D. Li, RSC Adv., 2020, 10, 34225–34230.
22. X. Q. Lu, C. Y. Gao, Z. H. Wei and S. D. Li, J. Mol. Model., 2021, 27, 130.
23. M. Z. Ao, X. Q. Lu, Y. W. Mu, W. Y. Zan and S. D. Li, Phys. Chem. Chem. Phys., 2022, 24, 3918–3923.
24. L. J. Yan, Inorg. Chem., 2022, 61, 10652–10660.
25. L. F. Z. Wang, X. F. Chen, H. Y. Du, Y. Q. Yuan, H. Qu and M. Zou, Appl. Surf. Sci., 2018, 427, 1030–1037.
26. H. Bai, B. Bai, L. Zhang, W. Huang, Y. W. Mu, H. J. Zhai and S. D. Li, Sci. Rep., 2016, 6, 35518.
27. L. Si and C. M. Tang, Int. J. Hydrogen Energy, 2017, 42, 16611–16619.
28. C. M. Tang, H. L. Liu, H. B. Yao and L. Fu, Int. J. Hydrogen Energy, 2020, 45, 21646–21654.
29. Y. J. Wang, L. Xu, L. H. Qiao, J. Ren, X. R. Hou and C. Q. Miao, Int. J. Hydrogen Energy, 2020, 45, 12932–12939.
30. S. Jeong, T. W. Heo, J. Oktawiec, R. P. Shi, S. Kang, J. L. White, A. Schneemann, E. W. Zaia, L. F. Wan, K. G. Ray, Y. S. Liu, V. Stavila, J. Guo, J. R. Long, B. C. Wood and J. J. Urban, ACS Nano, 2020, 14, 1745–1756.
31. Y. S. Liu, K. G. Ray, M. Jorgensen, T. M. Mattox, D. F. Cowgill, H. V. Eshelman, A. M. Sawvel, J. L. Snider, W. York, P. Wijeratne, A. L. Pham, H. Gunda, S. Li, T. W. Heo, S. Kang, T. R. Jensen, V. Stavila, B. C. Wood and L. E. Klebanoff, J. Phys. Chem. C, 2020, 124, 21761–21771.
32. E. Beheshti, A. Nojeh and P. Servati, Carbon, 2011, 49, 1561–1567.
33. P. P. Liu, H. Zhang, X. L. Cheng and Y. J. Tang, Int. J. Hydrogen Energy, 2017, 42, 15256–15261.
34. H. L. Dong, T. J. Hou, S. T. Lee and Y. Y. Li, Sci. Rep., 2015, 5, 09952.
35. Y. F. Zhang, X. Y. Han and X. L. Cheng, Chem. Phys. Lett., 2020, 739, 136961.
36. C. M. Tang and X. Zhang, Int. J. Hydrogen Energy, 2016, 41, 16992–16999.
37. C. Guo and C. Wang, Int. J. Hydrogen Energy, 2024, 60, 524–530.
38. M. D. Esrafili and S. Sadeghi, Int. J. Hydrogen Energy, 2022, 47, 11611–11621.
39. P. Jin, Q. H. Hou, C. C. Tang and Z. F. Chen, Theor. Chem. Acc., 2015, 134, 13.
40. Z. W. Zhang, W. T. Zheng and Q. Jiang, Int. J. Hydrogen Energy, 2012, 37, 5090–5099.
41. L. L. Zhang, B. Chen, S. Zhang and C. Lu, Int. J. Hydrogen Energy, 2025, 102, 1181–1187.
42. P. P. Liu, H. Zhang, X. L. Cheng and Y. J. Tang, Int. J. Hydrogen Energy, 2016, 41, 19123–19128.
43. T. W. Wang and Z. Y. Tian, Int. J. Hydrogen Energy, 2020, 45, 24895–24901.
44. M. Dolg, H. Stoll, A. Savin and H. Preuss, Theor. Chim. Acta, 1989, 75, 173–194.
45. M. Dolg, H. Stoll and H. Preuss, Theor. Chim. Acta, 1993, 85, 441–450.
46. V. Staroverov, G. Scuseria, J. Tao and J. Perdew, J. Chem. Phys., 2003, 119, 12129–12137.
47. J. D. Chai and M. Head Gordon, Phys. Chem. Chem. Phys., 2008, 10, 6615–6620.
48. M. J. Frisch, G. W. Trucks, H. B. Schlegel, G. E. Scuseria, M. A. Robb, J. R. Cheeseman, G. Scalmani, V. Barone, G. A. Petersson and H. Nakatsuji, Gaussian 16, Revision C.01, 2016.
49. B. Delley, J. Chem. Phys., 2000, 113, 7756–7764.
50. T. Lu and F. W. Chen, J. Comput. Chem., 2012, 33, 580–592.
51. W. Humphrey, A. Dalke and K. Schulten, J. Mol. Graphics, 1996, 14, 33–38.
52. G. J. Martyna, M. L. Klein and M. Tuckerman, J. Chem. Phys., 1992, 97, 2635–2643.
53. W. L. Li, C. Ertural, D. Bogdanovski, J. Li and R. Dronskowski, Inorg. Chem., 2018, 57, 12999–13008.

---