From Aromatic Motifs to Cluster-Assembled Materials: Silicon–Lithium Nanoclusters for Hydrogen Storage Applications

1. Introduction

Hydrogen is widely regarded as a promising energy carrier due to its high gravimetric energy density, environmental compatibility, and potential role in decarbonizing multiple sectors of the global energy system. However, its practical implementation is limited by the lack of storage technologies that simultaneously ensure safety, reversibility, and efficiency under ambient or near-ambient conditions. Conventional storage methods such as high-pressure compression and cryogenic liquefaction suffer from energy inefficiency, safety risks, and low volumetric densities. In contrast, material-based hydrogen storage—where hydrogen is stored via physisorption, chemisorption, or a combination of both—has emerged as a compelling alternative. For such systems to be practical, adsorption energies must typically fall within the range of −0.1 to −0.8 eV per H₂ molecule, striking a balance between sufficient binding strength and reversible release under operational conditions[1,2,3]. Theoretical studies have shown that nanostructured materials—including metal-organic frameworks, functionalized carbon materials, and transition-metal-decorated systems—can meet these criteria by fine-tuning their surface electronic structure, pore architecture, and active site polarity[4,5,6]. Within this landscape, atomically precise clusters offer several intrinsic advantages, including high surface-to-volume ratios, tunable reactivity, and well-defined sorption sites, positioning them as promising candidates for high-performance hydrogen storage applications.

Various silicon–lithium (Si–Li) clusters have been proposed for hydrogen storage, supported by computational predictions of high gravimetric capacities and adsorption energies compatible with reversible, non-dissociative adsorption[7,8,9,10,11]. Particular attention has focused on aromatic clusters combining enhanced stability with favorable interaction profiles. Si₅Li₆ (C_2v), studied by Jena et al. in 2006, was predicted to adsorb up to 14 H₂ molecules, though steric constraints limit effective uptake to around 10, with adsorption energies of 0.11–0.16 eV per H₂ range[7]. Shortly thereafter, in 2012, Pan, Merino, and Chattaraj reported Si₅Li₇⁺ (D_5h) and Si₄Li₄ (T_d) as viable hydrogen hosts, capable of storing 10–12 H₂ molecules with gravimetric capacities of 15.25 wt% and 14.7 wt%, and adsorption energies in the 0.10–0.20 eV per H₂ range[9]. Guo and Wang, in 2020, investigated SiLi₄⁺, composed of a central Si atom tetrahedrally coordinated by Li atoms, which adsorbs 12 H₂ molecules at 30.2 wt% and ~0.12 eV per H₂, with dynamic stability confirmed at 300 K via molecular dynamics simulations[10]. All these clusters—Si₄Li₄, Si₅Li₆, Si₅Li₇⁺, and SiLi₄⁺—were confirmed as global minima (GM) through systematic potential energy surface (PES) exploration, reinforcing the reliability of their predicted properties[10,12,13,14]. In contrast, larger clusters such as Si₆Li₆ (D_2h) and Si₁₀Li₁₀ (C_s), proposed by Jaiswal and Sahu in 2022[8], were constructed without global optimization. While they were predicted to adsorb 18 and 40 H₂ molecules, respectively, with gravimetric capacities of 14.7 wt% and 18.7%, and adsorption energies ranging from 0.059 to 0.141 eV per H₂, their thermodynamic relevance remains uncertain. Larger silicon-based assemblies have also been proposed, including Li-functionalized Si₂₀H₂₀ frameworks with lithium-containing organic groups (e.g., CN₂HLi, CONHLi), reported in 2015 to adsorb up to 60 H₂ molecules (12.5 wt%) [15], and Li₁₂Si₆₀H₆₀, a 2009 silicon analog of a decorated fullerene, predicted to bind 30 H₂ molecules (7.46 wt%)[11] . However, these extended systems have not undergone global optimization, leaving their stability and practical viability unverified.

Cluster-assembled materials (CAMs)[16], constructed from discrete, atomically defined units that retain their structural and electronic identity upon aggregation, provide a robust framework for the modular design of functional nanomaterials[17]. In the realm of Si–Li systems, our computational studies have shown that Li₄Si₄ (T_d) and Li₆Si₅ (C_2v)—both global minima—are promising building blocks stabilized by distinct forms of aromaticity. Li₄Si₄ exhibits spherical σ-aromaticity, while Li₆Si₅ features both σ- and π-aromatic delocalization. These clusters assemble into stable oligomers such as Li₈Si₈, Li₁₀Si₉, and Li₁₂Si₁₀, which preserve their local bonding environments and remain dynamically stable even at elevated temperatures[12]. Notably, these systems' Si₄^4– and Si₅^6– motifs are also present in experimental Zintl phases such as Li₁₂Si₇, Li₈MgSi₆, and Li₂₁Si₅, supporting their chemical viability[18,19,20,21]. Among these assemblies, (Li₄Si₄)_n oligomers are particularly attractive for hydrogen storage due to their high density of surface-accessible Li⁺ centers, which offer multiple binding sites for physisorption.

This study examines a series of structurally diverse Si–Li clusters selected for their potential to enable hydrogen storage via polar Li⁺ adsorption centers, favorable charge distribution, and electronically stable, modular architectures. Although distinct in topology, all analyzed systems share a common underlying motivation: they are either based on aromatic Si–Li motifs or serve as computational models of CAMs with accessible surfaces for physisorption. Firstly, we revisit the Si₆Li₆ and Si₁₀Li₁₀ clusters, previously proposed as high-capacity sorbents due to their resemblance to benzene and naphthalene, respectively[8]. While these structures suggest delocalized bonding frameworks, they were modeled without global optimization, and their thermodynamic viability remains unresolved. We, therefore, carry out a comprehensive exploration of their potential energy surfaces (PES) to locate the true global minima and reassess their hydrogen adsorption properties. In addition, we analyze the Li₁₂Si₅ (D₅h) cluster[13], recently reported by our group as a global minimum sandwich-type system, built from a Si₅¹⁰⁻ ring flanked by two Li₆⁵⁺ units[13]. This compact, highly polarizable structure features σ-aromatic delocalization and a high density of Li⁺ sites favorable for hydrogen binding. Finally, we investigate (Li₄Si₄)_n oligomers (n = 1–3), which model CAMs based on the Li₄Si₄ (T_d) global minimum and exhibit a large number of accessible Li⁺ adsorption sites. Together, these systems span a range of structural motifs and degrees of modularity, allowing us to examine the interplay between PES stability, electronic structure, and hydrogen uptake capacity in Si–Li nanocluster design.

2. Computational Details

The potential energy surfaces (PES) of the Li₆Si₆, Li₁₀Si₁₀, and (Li₄Si₄)n (n = 1–3) clusters were explored using the AUTOMATON program[22], which combines a probabilistic automata-based framework with genetic algorithms, for the (Li₄Si₄)n oligomers, PES exploration was also carried out using the guided Kick-MEP method[23]; full methodological details are provided in the Supporting Information. Initial structure screening for all systems was performed for singlet and triplet multiplicities at the PBE0[24]/SDDALL[25] level. Low-energy isomers within 20.0 kcal·mol⁻¹ of the putative global minimum were re-optimized at the PBE0-D3[26]/def2-TZVP[27] level. Harmonic vibrational frequency calculations confirmed that all reported minima are true stationary points. Final relative energies were obtained from single-point refinements at the DLPNO-CCSD(T)[28]/CBS[29] //PBE0-D3/def2-TZVP level, carried out using ORCA 5.0.3; only these values are considered in the energetic discussion. Geometry optimizations and frequency calculations were performed using Gaussian 16[30]. All computational details and representative low-energy structures are included in the Supporting Information (Sections S1–S3).

Hydrogen adsorption energetics were evaluated by computing adsorption energies (E_ads) for all nH₂–cluster complexes using the expression:

$${\mathrm{E}}_{\mathrm{ads}}={\displaystyle \frac{\left[\mathrm{E}\left(\mathrm{complex}\right)-\mathrm{nE}({\mathrm{H}}_2\right)-\mathrm{E}\left(\mathrm{cluster}\right)]}{{\mathrm{nH}}_2}}$$

(1)

Here, E(complex) is the total energy of the hydrogen–cluster complex, nE(H₂) is the energy of n isolated H₂ molecules, and E(cluster) is the energy of the bare cluster. To improve accuracy, adsorption energies were corrected for basis set superposition error (BSSE)[31] arising from artificial stabilization due to basis function overlap in weakly bound systems. Boys and Bernardi's counterpoise (CP) correction method was employed to account for BSSE[32]. The CP-corrected interaction energy (E_CP) is given by:

$${\mathrm{E}}_{\mathrm{CP}}={\mathrm{E}}_{\mathrm{Int}}-{\mathrm{E}}_{\mathrm{BSSE}}$$

(2)

Here, Eint is the uncorrected interaction energy, and EBSSE is the basis set superposition error correction. All cluster geometries and their hydrogenated analogues were optimized using Gaussian 16 at the M06[33]/6-311+G(d,p)[34] level, chosen for its proven accuracy in describing dispersion-driven interactions relevant to hydrogen storage[8]. Practical storage capacity was estimated through gravimetric density calculations at approximate saturation using the expression:

$$\mathrm{wt}\%={\displaystyle \frac{\mathrm{M}({\mathrm{nH}}_2)}{\mathrm{M}({\mathrm{nH}}_2)+\mathrm{M}\left(\mathrm{cluster}\right)}}\ast 100$$

(3)

Here, M(nH₂) is the total mass of the adsorbed hydrogen molecules, and M(cluster) is the molecular mass of the bare cluster.

Non-covalent interactions were analyzed using the Independent Gradient Model based on Hirshfeld partitioning (IGMH)[35]. This method improves upon the original IGM method by using Hirshfeld-derived atomic densities, offering a more physically grounded and higher-resolution depiction of weak interactions. Compared to conventional NCI plots, IGMH provides more detailed insights into dispersion-driven adsorption. All analyses were performed at the M06/6-311+G(d,p) level using Multiwfn[36], with visualizations rendered in VMD[37].

To evaluate the structural resilience and reversibility of hydrogen adsorption under thermal conditions, Born–Oppenheimer molecular dynamics (BOMD)[38] simulations were conducted at 300 K and 400 K, temperatures relevant to practical hydrogen storage. The ADMP method[39], a Lagrangian formulation using Gaussian basis sets, was employed as a computationally efficient alternative to conventional BOMD, offering comparable accuracy for trajectory prediction. Simulations were run for 10 ps to assess structural stability and the dynamic retention of adsorbed hydrogen over time.

3. Results and Discussions

3.1. Confirming Li-Si Lowest Energy Structures

Figure 1 shows the global minima and the previously proposed high-symmetry structures for the Li₆Si₆ and Li₁₀Si₁₀ clusters[8], as identified through comprehensive PES exploration. For Li₆Si₆ (panel a, left), the global minimum exhibits Cₛ symmetry (¹A′ state) and is composed of a Si₄ tetrahedral unit (T_d) and a Si₂ dimer, surrounded by asymmetrically distributed Li atoms. This compact three-dimensional configuration is significantly more stable—by 20.4 kcal·mol⁻¹—than the previously proposed D_2h-symmetric structure (panel a, right), which features a planar Si₆ ring and had been modeled as a benzene analog in earlier hydrogen storage studies. For Li₁₀Si₁₀ (panel b), the global minimum adopts C₁ symmetry (¹A′ state) and consists of two Si₄ tetrahedra and one Si–Si dimer, coordinated by a spatially dispersed set of lithium atoms. This low-symmetry arrangement lies 52.0 kcal·mol⁻¹ below the Cₛ-symmetric isomer (panel b, right), previously proposed as a naphthalene-like candidate for molecular hydrogen adsorption. These results indicate that the earlier high-symmetry, π-aromatic-like structures do not correspond to thermodynamically favored forms. Instead, the emergence of Si₄-based building units in both global minima highlight the stabilizing role of such substructures in lithium–silicon chemistry and supports their relevance for the design of hydrogen storage clusters. Additional low-energy isomers within 20 kcal·mol⁻¹ of the global minima are reported in Figure S1-S2 (Supporting Information), further illustrating the structural diversity of these systems.

In parallel with our identification of new global minima, we re-evaluated the potential energy surfaces of two previously proposed systems—Li₁₂Si₅ (D₅h, ¹A₁) and the (Li₄Si₄)ₙ oligomers (n = 2–3)—reported initially as global minima and selected here for their relevance to hydrogen storage. Our independent PES analysis confirms that these structures correspond to ground-state configurations. As shown in Figure 2, the Li₄Si₄ (T_d) monomer maintains its structural integrity upon oligomerization, with Li₈Si₈ and Li₁₂Si₁₂ preserving the local Si₄-based connectivity and Li⁺ coordination. The Li₁₂Si₅ cluster was likewise verified to occupy the lowest point on the potential energy surface. These results substantiate the thermodynamic stability of the selected clusters and support their role as structurally persistent, modular building units for hydrogen-rich Si–Li assemblies.

3.2. Structural and Electronic Features of Bare and Hydrogen-Adsorbed Si–Li Clusters

Table 1 depicts the computed interatomic distances for the Si–Li clusters examined in this study and their corresponding hydrogen-adsorbed complexes. Si–Si bond lengths remain largely invariant upon H₂ adsorption, spanning 2.12–2.57 Å, reflecting the structural integrity of the silicon frameworks. Si–Li distances in the bare clusters range from approximately 2.40 to 2.76 Å and undergo modest elongation upon hydrogenation, particularly in larger assemblies such as Li₈Si₈ and Li₁₂Si₅, where maximum distances approach ~2.89 Å. Li–H separations offer insight into adsorption, with the shortest distances (2.08–2.30 Å) corresponding to favorable electrostatic interactions with exposed Li⁺ centers. As hydrogen coverage increases, longer Li–H contacts—up to ~3.9 Å—are observed, indicating weaker physisorption at more peripheral or sterically hindered sites. This trend appears in both the global minimum (GM) structures identified through PES exploration and in the previously proposed high-symmetry local minima, Li₆Si₆* and Li₁₀Si₁₀*. While both structures support molecular hydrogen adsorption, the GM isomers generally present more spatially accessible and topologically diverse Li⁺ coordination, which can enhance adsorption site availability. In all cases, H–H bond distances remain close to 0.75 Å, consistent with nondissociative molecular adsorption. These structural characteristics underscore the relevance of using thermodynamically validated GM structures for accurately predicting hydrogen storage performance in Si–Li clusters.

Table 2 compiles the HOMO–LUMO energy gaps (ΔE_H–L) of the studied lithium–silicon clusters, both in their bare and hydrogen-adsorbed configurations. These values offer insight into electronic stability and chemical hardness. Among the bare clusters, Li₄Si₄ exhibits the highest gap (3.1 eV), which increases to 3.4 eV upon adsorption of 8 H₂ molecules and remains relatively high (3.2 eV) for 12 H₂, confirming its closed-shell character and resilience to electronic perturbation. The global minimum structures of Li₆Si₆ and Li₁₀Si₁₀ show gaps of 2.8 and 2.6 eV, respectively, which are also maintained or even slightly enhanced upon hydrogenation, reaching 2.9 eV for both 12H₂@Li₆Si₆ and 30H₂@Li₁₀Si₁₀. In contrast, the previously proposed high-symmetry structures—Li₆Si₆* (D_2h) and Li₁₀Si₁₀* (C_s)—exhibit narrower gaps of 2.2 eV and 1.8 eV, respectively, with negligible changes upon hydrogen loading. These values match closely with those reported by Jaiswal et al., who found ΔE_H–L values of 2.33 eV for Li₆Si₆ and 1.81 eV for Li₁₀Si₁₀ using the B3LYP/6-31G(d,p) level of theory. Although method-dependent differences are expected, the trend is consistent: higher energy isomers on PES tend to display smaller HOMO–LUMO gaps and reduced electronic stability. Additional insights are observed in other systems. Li₈Si₈ displays an increase from 2.7 to 3.1 eV as H₂ adsorption progresses up to 16 molecules, followed by a decline to 2.1 eV at 24H₂, suggesting a saturation threshold for electronic stabilization. Li₁₂Si₁₂ maintains gaps above 2.6 eV, while Li₁₂Si₅, the most compact and polarizable structure, shows the lowest gaps (1.3–1.7 eV), which increase modestly with hydrogen coverage.

3.3. Hydrogen Adsorption Energetics, Charge Redistribution, and Storage Capacity

The hydrogen adsorption properties of the investigated lithium–silicon clusters were assessed through BSSE-corrected adsorption energies (E_ads), partial charges on Li centers, and gravimetric hydrogen capacities (wt%), as summarized in Table 3. In all systems, Li atoms initially exhibit partial charges between +0.75 and +0.89, consistent with their role as electropositive adsorption sites. Upon hydrogen adsorption, these charges decrease progressively—reaching as low as +0.30 in Li₁₂Si₅—reflecting charge redistribution driven by weak donor-acceptor interactions. The BSSE-corrected Eads values lie within the optimal range for reversible hydrogen storage (−0.11 to −0.16 eV per H₂), with slightly stronger binding at low coverage. The difference from uncorrected values (~0.01–0.02 eV) confirms the weakly bound nature of the interactions and the importance of applying BSSE corrections for reliable energetic estimates. All PES-validated GM structures—Li₄Si₄, Li₆Si₆, Li₈Si₈, and Li₁₀Si₁₀—reach 14.72 wt% through adsorbing 12, 18, 24, and 30 H₂ molecules, respectively. These values align with the upper limit of approximately three H₂ molecules per Li⁺, as Pan, Merino, and Chattaraj proposed[9]. Li₄Si₄, a symmetric and modular unit (T_d), and its oligomeric derivatives Li₈Si₈ and Li₁₂Si₁₂, maintain favorable adsorption behavior and extend capacity, with Li₁₂Si₁₂ reaching 14.72 wt% with 36 H₂. While the high-symmetry isomers of Li₆Si₆ (D_2h) and Li₁₀Si₁₀ (C_s) achieve comparable hydrogen uptake and adsorption energetics, they correspond to higher-energy local minima and show narrower Li charge distributions relative to the GMs. The Li₁₂Si₅ GM attains the highest capacity in the series (23.45 wt% with 34 H₂), enabled by its compact D_5h geometry, slightly stronger adsorption energies, and a broader Li charge range (+0.30 to +0.78), which reflects increased electrostatic heterogeneity across adsorption sites. For example, Figure 3 illustrates the stepwise hydrogen uptake over Li₈Si₈ and Li₁₂Si₅, showing progressive occupation of Li⁺ centers with minimal structural deformation. These results demonstrate that the Si–Li GMs studied here combine thermodynamic stability, favorable charge redistribution, and optimal adsorption energetics to enable efficient and reversible hydrogen storage.

3.3. Thermal Stability and Hydrogen Release Dynamics via BOMD Simulations

The reversibility and thermal resilience of hydrogen adsorption were assessed through Born–Oppenheimer molecular dynamics (BOMD) simulations on hydrogen-loaded Si–Li clusters: 12H₂@Li₄Si₄, 18H₂@Li₆Si₆, 24H₂@Li₈Si₈, 30H₂@Li₁₀Si₁₀, 34H₂@Li₁₂Si₅, and 36H₂@Li₁₂Si₁₂. Simulations were conducted for 10 ps at 300 K and 400 K to probe hydrogen release under operating conditions. As shown in Figure 4 and Figure 5, desorption behavior depends strongly on cluster size. At 300 K, smaller clusters such as Li₄Si₄ and Li₆Si₆ release the majority of their adsorbed hydrogen within the first few picoseconds, with only 1–2 H₂ molecules retained by the end of the simulation. In contrast, larger systems like Li₁₂Si₅ and Li₁₂Si₁₂ retain 6 and 12 H₂ molecules, respectively, under identical conditions. This size-dependent stability becomes even more pronounced at 400 K, where small clusters undergo complete or near-complete desorption, while Li₁₂Si₁₂ maintains nearly one-third of its original hydrogen load. These observations confirm that increasing cluster size and coordination density enhances hydrogen retention under elevated thermal conditions.

A comparative analysis between PES-validated global minima (GM) and previously reported local minima (LM) structures for Li₆Si₆ and Li₁₀Si₁₀ (Figure 4) further underscores the role of structure optimization. While initial desorption rates are similar for GM and LM configurations, the GMs consistently retain more hydrogen at later times. For instance, at 400 K, 18H₂@Li₆Si₆(GM) retains ~2 H₂ molecules, whereas the LM counterpart undergoes full desorption. Similarly, 30H₂@Li₁₀Si₁₀(GM) retains ~4 H₂, while the LM form releases nearly all of its hydrogen content. These differences, though subtle in kinetics, reveal that GM structures offer more resilient binding sites capable of maintaining adsorbed hydrogen under thermal fluctuations. Altogether, the BOMD results support the conclusion that larger, PES-validated Si–Li clusters—particularly Li₁₀Si₁₀, Li₁₂Si₅, and Li₁₂Si₁₂—combine structural stability and dynamic retention, making them promising candidates for reversible hydrogen storage under realistic operating conditions.

3.4. Visualization of Non-Covalent Interactions via IGMH Analysis

To gain visual insight into the nature of hydrogen adsorption in lithium–silicon clusters, we applied the Independent Gradient Model based on Hirshfeld partitioning (IGMH) to a series of representative hydrogenated systems (Figure 6). As a qualitative method, IGMH enables high-resolution visualization of non-covalent interactions by partitioning electron density based on the actual molecular environment, offering superior clarity compared to traditional NCI plots—especially for dispersion-dominated systems. In all cases, the interactions between H₂ molecules and Li⁺ centers are characterized by green-colored isosurfaces located between the adsorbates and the cluster surface, confirming that physisorption is mediated primarily by weak van der Waals forces. No significant blue or red regions are observed, indicating the absence of strong electrostatic attractions or steric repulsion. The interaction fields are more extensive and homogeneously distributed in PES-validated structures such as 12H₂@Li₄Si₄, 24H₂@Li₈Si₈, and 34H₂@Li₁₂Si₅, reflecting their favorable electrostatic landscapes and high hydrogen uptake. In contrast, the local minima geometries of Li₆Si₆(D_2h) and Li₁₀Si₁₀(Cₛ) exhibit more fragmented interaction regions, correlating with reduced charge delocalization and lower retention in BOMD simulations. Overall, the IGMH results reinforce that dispersion-driven, non-dissociative physisorption governs hydrogen storage in these Si–Li clusters, complementing the energetic and dynamic analyses and validating their potential as reversible hydrogen carriers.

4. Conclusions

This study comprehensively investigates Si–Li clusters as hydrogen storage candidates, emphasizing the importance of accurate potential energy surface (PES) exploration. For the Li₆Si₆ and Li₁₀Si₁₀ systems, global optimization revealed that the true ground-state structures do not correspond to previously proposed benzene- and naphthalene-like motifs but rather consist of compact arrangements built from tetrahedral Si₄ (T_d) units and a Si₂ dimer, all stabilized by surrounding lithium atoms. Specifically, Li₆Si₆ features a single Si₄ unit and a Si₂ bridge, while Li₁₀Si₁₀ includes two Si₄ units linked via a Si–Si dimer. These findings highlight the necessity of PES validation in cluster-based materials design, ensuring that the predicted structures are electronically stable and experimentally feasible.

Guided by the recurrence of the Si₄–T_d motif, we evaluated oligomeric systems—dimers, trimers, and tetramers—of Li₄Si₄ as model cluster-assembled materials (CAMs). These oligomers preserve the local Si₄ environments and exhibit additive hydrogen storage behavior: each Li₄Si₄ unit binds 12 H₂ molecules (3 per Li⁺), with total storage capacities increasing proportionally with the number of monomers. Additionally, we explored the aromatic Li₁₂Si₅ sandwich cluster, which incorporates a planar Si₅ ring stabilized by two Li₆ caps. This system achieves the highest gravimetric capacity in the series (23.45 wt%), attributable to its high density of accessible Li⁺ sites. Collectively, our results establish clear structure–property relationships among Si–Li clusters and demonstrate the potential of PES-validated, modular architectures for designing thermally stable, high-capacity hydrogen storage materials.

Importantly, this work significantly advances the field of Li–Si clusters for hydrogen storage by delivering (i) a rigorous benchmark of global minimum structures for key systems previously proposed in the literature, (ii) new modular design principles based on Si₄–T_d building blocks, and (iii) the identification of high-performing sandwich-type and oligomeric clusters with verified thermodynamic and dynamic stability. These contributions provide a robust theoretical foundation for guiding future experimental development of silicon–lithium nanomaterials in next-generation hydrogen storage technologies.