On the performance of pure and group 2B transition metal-doped metal oxide nanocages as single-atom catalysts for the hydrogen storage process: a DFT study

Abstract

For an efficient confrontation of the exhaustion of nonrenewable energy sources issue, the storage of hydrogen as an eco-friendly and renewable alternative energy source has received considerable attention. Herein, the performance of pure and group 2B transition metal-doped metal oxide nanocages (M₁₂O₁₂ and TM-M₁₁O₁₂; where M = Zn, Mg, and Be; TM = Zn, Cd, and Hg) as single-atom catalysts for the hydrogen dissociation reaction (HDR) was investigated using DFT calculations. Regarding step-I of the HDR, all the investigated catalysts exhibited remarkable potentiality to adsorb the H₂ molecule with negative BSSE-corrected adsorption energy values up to −5.22 kcal mol⁻¹. In step-II, further activation for the H₂ molecule over the surface of the M₁₂O₁₂ and TM-M₁₁O₁₂ catalysts occurred, and hence the transition state (TS) structure was obtained. Upon the energetic results, the Zn₁₂O₁₂-based catalysts exhibited higher performance toward the HDR compared to the Mg₁₂O₁₂- and Be₁₂O₁₂-based candidates. Furthermore, the Cd-Zn₁₁O₁₂ catalyst demonstrated the most promising catalytic activity with an activation energy of 9.58 kcal mol⁻¹ for the H₂⋯Cd-Zn₁₁O₁₂ complex. In step-III, one of two activated H atoms (H1) shifted to the Zn atom, whereas the other hydrogen atom (H2) migrated to the O atom. Analysis of natural bond orbitals and electron density difference outlined the charge transfer from M/TM atoms to their interacting hydrogen atom (H1) and from the O atom to the corresponding hydrogen atom (H2). Quantum theory of atoms in molecules outcomes demonstrated the partial covalent nature of the interactions within the TS structures, pinpointing the optimum catalytic efficiency. The obtained results will provide a comprehensive picture of the behavior of metal oxide-based SACs for HDR catalysis, and hence their performance for the hydrogen storage process.

---

Introduction

Nowadays, catalysts are considered to be the fundamental component in numerous commercial and industrial energy conversion methods.¹ The technology of catalysts is essential in managing the annual production of diverse goods worth one trillion dollars in the chemical, petroleum, food, and power industries.² Notably, noble metals (i.e., Pd, Pt, Rh, Ru, and Au) exhibit remarkable catalytic activity for energy production.^3–7 In this respect, there are a few obstacles that are linked with noble metal-based catalysts, such as their high cost, elevated operation temperature, and restricted availability.^8–10 In turn, the noble metal-based catalysts are deemed economically unfeasible.^11,12 Replacement of noble metals with low-cost ones is extremely intriguing. In this regard, the first-row transition metals have gained considerable attention as a result of their considerable availability and affordability.^13,14

In the realm of catalysis, single-atom catalysis has emerged as an innovative method for improving performance and minimizing the usage of metals.¹⁵ The single-atom catalysts (SACs) are materials with a catalytic active metal site anchored at the atomic scale on precise supports with a noteworthy surface area.¹⁶ In addition to their outstanding catalytic activity, SACs are characterized by their cost-effective nature and optimal utilization efficiency, which in turn reinforces their significant contributions to multilateral applications.^16–18 In SACs, two-dimensional layered materials, metal–organic frameworks (MOFs), or metal oxides are the predominant supports due to their propensity to enhance SACs’ stability.^19,20 The significant physical and chemical stability, along with a high surface area, of metal oxides enabled them to be frequently utilized as supporters.²¹ Recently, studying hydrogen storage capacity was limited to adsorption energy and bulk module concepts.²² Nowadays, the hydrogen storage vein has been extended to include the catalyzed hydrogen dissociation reaction (HDR) as an elementary process.²³ Notably, the preceding reports pinpointed that the HDR entails three consecutive steps,^24–28 dubbed hydrogen adsorption, hydrogen cleavage, and hydrogen migration, in turn forming intermediate (I), transition state (TS), and product (P) structures, respectively. Basically, the activation energy (E_activation) was considered as the most crucial factor in determining the reactivity of the catalysts.^29,30 Although a recent study highlighted the efficacy of metal oxide-based SACs for hydrogen storage applications from the formation/adsorption energy perspective and bulk modules,²² their HDR mechanism is still obscure. Additionally, the impact of the doping process using group 2B transition metals on metal oxide-based SACs’ performance through the HDR mechanism has not been investigated yet.

The current study was designed to uncover the efficiency of pure and group 2B transition metal-doped metal oxide nanocages as SACs for the HDR catalysis. The HDR mechanism is demonstrated in Fig. 1. First, geometrical optimization was performed for the SACs under study, followed by the electrostatic potential (EP) analysis to uncover the regions with electrophilic and nucleophilic natures over their surfaces. Frontier molecular orbital (FMO) theory was applied to investigate the changes in the electronic features after the doping process. Furthermore, natural bond orbitals (NBO) charges, electron density difference (EDD), and quantum theory of atoms in molecules (QTAIM) computations were accomplished for the TS structures of the H₂⋯M₁₂O₁₂ and ⋯TM-M₁₁O₁₂ complexes. The presented study would be a trustworthy reference for the forthcoming works concerned with designing cost-effective and high-performance single-atom catalysts for HDR catalysis.

Fig. 1 (a) Schematic representations for the M₁₂O₁₂ and TM-M₁₁O₁₂ catalysts, (b) relevant structures for the HDR mechanism, including reactants, intermediate, transition state, and product structures.

Computational methods

The performance of pure and group 2B transition metal-doped metal oxide nanocages (M₁₂O₁₂ and TM-M₁₁O₁₂; where M = Zn, Mg, and Be; TM = Zn, Cd, and Hg) as SACs for the HDR was investigated using DFT computations. All the considered calculations were performed with the help of Gaussian 09 software³¹ using the M06-2X method due to its well-documented accuracy for thermochemistry and kinetics.³² Moreover, Los Alamos National Laboratory 2-Double-zeta (LANL2DZ) basis set was utilized for the Zn, Cd, and Hg atoms³³ and the 6-311+G** one for the remaining atoms.³⁴ Based on the literature, LANL2DZ is a plausible basis set for describing the transition metals due to its ability to thoroughly describe their electronic nature and core potential,^35–38 while 6-311+G** offers a promising balance between computational efficiency and accuracy.^39,40

Geometrical optimization computations were first carried out at the singlet state for the M₁₂O₁₂ catalysts and at the lowest four possible spin states for the TM-M₁₁O₁₂ ones. For TM-M₁₁O₁₂ catalysts, relative energy (E_relative) computations were accomplished to pinpoint the most favorable spin state using eqn (1). All the forthcoming calculations for the isolated TM-M₁₁O₁₂ catalysts were conducted at the most favorable spin state. Spin contamination was also computed for the studied catalysts by computing the expected value of the <S²> operator through eqn (2).

E_relative = E_TM-M11O12 − E_{TM-M11O12 at the most favorable spin state} (1)

<S²> = S(S + 1) (2)

where E_TM-M11O12 refers to the energy of TM-M₁₁O₁₂ catalysts at a definite spin state, whereas the E_{TM-M11O12 at the most favorable spin state} represents the energy of TM-M₁₁O₁₂ catalysts at the most favorable spin state. All the forthcoming calculations for the isolated TM-M₁₁O₁₂ catalysts were conducted at the most favorable spin state. S is the half number of unpaired electrons.

The formation energy (E_form) of M₁₂O₁₂/TM-M₁₁O₁₂ catalysts was calculated with respect to the isolated atoms as illustrated in eqn (3). Moreover, relative formation energy (E^relative_form) for the TM-M₁₁O₁₂ catalysts with reference to the corresponding pure nanocages was evaluated through eqn (4).^41,42

E_form = (E_{M12O12/TM-M11O12} − aE_M − bE_TM − 12E_O)/24 (3)

E^relative_form = E_TM-M11O12 − E_M12O12 ± n_iµ_i (4)

where E_{M12O12/TM-M11O12} refers to the energy of M₁₂O₁₂/TM-M₁₁O₁₂ catalyst. E_M, E_TM, and E_O represent the energies of a single M, TM, and O atom, respectively. a and b stand for the number of M and TM atoms within the investigated SACs, respectively. The n_i and µ_i represent the number and chemical potential of atoms, respectively. For one-to-one metal substitution in the M₁₂O₁₂ nanocages, eqn (4) could be simplified to be as follows:where E_TM12O12 refers to energies of the optimized TM₁₂O₁₂ nanocages.

Electrostatic potential (EP) analysis was carried out, aiming to pinpoint the electrophilic or nucleophilic surfaces over the SACs under investigation.⁴³ In this context, the molecular electrostatic potential (MEP) maps were built for the M₁₂O₁₂ and TM-M₁₁O₁₂ catalysts at the most favorable spin state with the implementation of a 0.002 au electron density envelope. Additionally, FMOs were conducted to study the TM-doping impact on the electronic features of the M₁₂O₁₂ catalysts at the most favorable spin state. In this vein, the distributions of the highest occupied/lowest unoccupied molecular orbitals (HOMO/LUMO) were created for the M₁₂O₁₂ and TM-M₁₁O₁₂ catalysts. Additionally, HOMO/LUMO energies (i.e., E_HOMO/E_LUMO) were calculated. Accordingly, the Fermi level (E_FL) and energy gap (E_gap) were assessed via eqn (6) and (7), respectively. The TM-doping impact on the reactivity of SACs was studied by evaluating the % ΔE_gap using eqn (8).

E_gap = E_LUMO − E_HOMO (7)

Moreover, the ionization potential (IP) and electron affinity (EA) were calculated from eqn (9) and (10), respectively.

IP ≈ −E_HOMO (9)

EA ≈ −E_LUMO (10)

In the power of Koopman’s theorem,⁴⁴ the chemical reactivity descriptors of the employed SACs, including global hardness (η), work function (Φ), electrophilicity index (ω), global softness (S), and chemical potential (µ) were computed through eqn (11)–(15). In the work function calculation, the vacuum-level electrostatic potential (V_el(+∞)) was suggested to be nearly 0.

Φ = V_el(+∞) − E_FL (15)

For the HDR mechanism, three consecutive steps were defined along the reaction pathway. In step-I, the propensity of the M₁₂O₁₂ and TM-M₁₁O₁₂ catalysts to adsorb the H₂ molecule was estimated by evaluating the adsorption energy (E_ads) as follows:

E_ads = E_intermediate − E_reactants (16)

where E_intermediate identifies the energy of the intermediate structure of H₂⋯M₁₂O₁₂/TM-M₁₁O₁₂ complexes. E_reactants represents the sum of energies relevant to the isolated M₁₂O₁₂/TM-M₁₁O₁₂ catalysts and H₂ molecule. Further, the corrected adsorption (E^corrected_ads) and interaction (E^corrected_int) energies of the H₂⋯M₁₂O₁₂/TM-M₁₁O₁₂ complexes were calculated via eqn (17) and (18), respectively. The annihilation of the basis set superposition error (BSSE) from the E^corrected_ads and E^corrected_int energies was carried out using the counterpoise-corrected method.⁴⁵ Furthermore, the solvent effect upon the adsorption process was investigated by assessing the solvation energy (ΔE_solv) values using eqn (19).

E^corrected_ads = E_{H2⋯M12O12/TM-M11O12} − (E_H2 + E_{M12O12/TM-M11O12}) + E_BSSE (17)

E^corrected_int = E_{H2⋯M12O12/TM-M11O12} − (E_{H2 in complex} + E_{M12O12/TM-M11O12 in complex}) + E_BSSE (18)

ΔE_solv = E_solvent − E_gas (19)

where E_{H2⋯M12O12/TM-M11O12} is the energy of the considered complexes. E_{M12O12/TM-M11O12} and E_H2 are energies of the isolated M₁₂O₁₂/TM-M₁₁O₁₂ catalysts and H₂ molecule, respectively. The E_{M12O12/TM-M11O12 in complex} and E_{H2 in complex} are the energies of the M₁₂O₁₂/TM-M₁₁O₁₂ catalysts and H₂ molecule with respect to their coordinates in the optimized complexes. E_BSSE is the energy of basis set superposition error. E_solvent and E_gas represent energies of the studied complexes in the gas and solvent phases, respectively.

In step-II, the TS structures were located at the lowest four possible spin states. Frequency calculations were performed to corroborate the TS with one imaginary frequency along with the absence of imaginary frequency for the rest geometries (i.e., true minima nature). The intrinsic reaction coordinate (IRC) calculations were performed for the TS structures to verify that all the procured structures were adequately connected to their respective minima along the imaginary vibration mode. In addition, the activation energy (E_activation) was computed as the energy difference between the TS structure (E_TS) and the intermediate (E_intermediate) one via eqn (20).

E_activation = E_TS − E_intermediate (20)

In step-III, the hydrogen migration step was studied. In this vein, the reaction energy (E_reaction) and released energy (E_released) along the reaction pathway were computed via eqn (21) and (22), respectively.

E_reaction = E_product − E_intermediate (21)

E_released = E_product − E_TS (22)

where the E_TS and E_product refer to the energies of the TS and product, respectively.

In this context, the literature states that the charge transfer was considered the most crucial factor beyond the occurrence of the investigated HDR catalysis.^24–27 Therefore, analyses of NBO and EDD were performed for the TS structures of H₂⋯M₁₂O₁₂/TM-M₁₁O₁₂ complexes at the most favorable spin state to provide a comprehensive explanation of the activation and dissociation processes pertinent to the H₂ molecule over the examined catalysts. In order to examine the interactions within the TS structures of H₂⋯ M₁₂O₁₂ and ⋯TM-M₁₁O₁₂ complexes at the most favorable spin state, QTAIM analysis was accomplished by extracting the bond critical points (BCPs) and bond paths (BPs). Moreover, the topological parameters at the desired BCP were evaluated. EP and QTAIM analyses were conducted using Multiwfn 3.7 software.⁴⁶ Using Visual Molecular Dynamics program, QTAIM schemes were created.⁴⁷

Results and discussion

Spin state study

Geometrical optimization computations were accomplished for all TM-M₁₁O₁₂ nanocages at the lowest four possible spin states to investigate the spin state effect on the stability of the studied nanocages. To determine the most preferable structure, the E_relative and the expected values of <S²> operator are listed in Table S1. Frequency calculations were also conducted to verify the true minima of the obtained geometries.

From Table S1, favorable structures were procured for all TM-M₁₁O₁₂ catalysts at the singlet state that could be illustrated as an upshot of filling up all d orbitals. Accordingly, the expected values of the <S²> operator were found to be zero for all the studied catalysts, highlighting the absence of spin contamination. Frequency calculations demonstrated that there was no imaginary frequency for all TM-M₁₁O₁₂ catalysts, confirming the true minima character.

Geometric structure and EP analyses

To pictorially illustrate the nucleophilic and electrophilic regions on the surfaces of the studied catalysts, EP analysis was employed.^48,49 The optimized structures and MEP maps of M₁₂O₁₂ and TM-M₁₁O₁₂ catalysts at the most favorable spin state are pictured in Fig. 2. The geometrical features of the optimized catalysts are tabulated in Table 1.

Fig. 2 Optimized structures and MEP maps of the M₁₂O₁₂ and TM-M₁₁O₁₂ catalysts at the most favorable spin state.

Table 1 Geometrical features of the studied M₁₂O₁₂ and TM-M₁₁O₁₂ catalysts at the most favorable spin state

Catalyst	Distance (Å)		Eform (kcal mol^−1)	E^relativeform (kcal mol^−1)
	d1	d2
Zn12O12	1.91	1.98	−102.28	—
Cd-Zn11O12	2.10	2.18	−101.20	−0.29
Hg-Zn11O12	2.17	2.28	−99.36	0.04
Mg12O12	1.88	1.94	−125.10	—
Zn-Mg11O12	1.91	1.98	−123.28	−1.87
Cd-Mg11O12	2.09	2.17	−122.12	−0.43
Hg-Mg11O12	2.14	2.27	−120.26	0.48
Be12O12	1.52	1.58	−159.86	—
Zn-Be11O12	1.91	1.98	−155.12	−1.33
Cd-Be11O12	2.09	2.17	−154.05	−2.00
Hg-Be11O12	2.17	2.27	−152.13	0.37

---

Inspecting the optimized structures of the studied catalysts unveiled that the M₁₂O₁₂ and TM-M₁₁O₁₂ catalysts consisted of six four- (tetragons) and eight six-membered (hexagons) rings. Moreover, two distances of the M/TM–O bond were disclosed; the first was denoted within the four-membered rings (d₁), and the other within the six-membered rings (d₂). Furthermore, the MEP maps illustrated the localization of electrophilic and nucleophilic areas on the M/TM and O atoms, respectively.

From Table 1, values of d₁ and d₂ relevant to the M₁₂O₁₂ catalysts were in line with the previous reports.^50–52 Upon the doping process, d₁ and d₂ were greatly lengthened; for instance, d₁/d₂ were 1.91/1.98 and 2.10/2.18 Å for the Zn₁₂O₁₂ and Cd-Zn₁₁O₁₂ catalysts, respectively. Clearly, remarkable negative E_form values were perceived for all studied M₁₂O₁₂ and TM-M₁₁O₁₂ catalysts, in the range from −99.36 to −159.86 kcal mol⁻¹. Moreover, the lack of any structural deformation following the optimization within the TM-M₁₁O₁₂ catalysts highlighted the facile character of the doping process on the M₁₂O₁₂ catalysts. Moreover, E^relative_form for the TM-M₁₁O₁₂ catalysts were predominantly negative, with only small positive values in a few cases, indicating the facile nature of the doping process and the energetic stability of the resulting structures.⁴¹ Overall, the outcomes of E_form and E^relative_form, along with the absence of imaginary frequency of all M₁₂O₁₂ and TM-M₁₁O₁₂ catalysts, demonstrated that the investigated catalysts were stable nanocages rather than aggregates. Such findings were in line with the molecular dynamics-based results relevant to Zn₁₂O₁₂, Mg₁₂O₁₂, and Be₁₂O₁₂ nanocages.^22,53

Electronic parameters

In order to investigate the electronic features of chemical systems, FMOs theory is employed.^54,55 Using the FMOs, the changes in the electronic characteristics of the studied catalysts after the doping process were studied. Fig. 3 demonstrates the HOMO/LUMO electron densities for the M₁₂O₁₂ and TM-M₁₁O₁₂ catalysts. Table 2 compiles the electronic parameters, including E_HOMO, E_LUMO, E_FL, and E_gap values.

Fig. 3 HOMO and LUMO distribution patterns of the M₁₂O₁₂ and TM-M₁₁O₁₂ catalysts at the most favorable spin state.

Table 2 Electronic features (in eV) of the M₁₂O₁₂ and TM-M₁₁O₁₂ catalysts at the most favorable spin state

Catalyst	EHOMO	ELUMO	EFL	Egap	% ΔEgap
Zn12O12	−8.815	−1.952	−5.384	6.864	—
Cd-Zn11O12	−8.652	−2.023	−5.338	6.629	−3
Hg-Zn11O12	−8.504	−2.180	−5.342	6.325	−5
Mg12O12	−8.485	−1.093	−4.789	7.392	—
Zn-Mg11O12	−8.461	−1.152	−4.807	7.309	−1
Cd-Mg11O12	−8.266	−1.274	−4.770	6.992	−4
Hg-Mg11O12	−8.077	−1.448	−4.763	6.629	−5
Be12O12	−10.611	−0.671	−5.641	9.940	—
Zn-Be11O12	−10.062	−1.676	−5.869	8.385	−16
Cd-Be11O12	−9.714	−2.508	−6.111	7.207	−14
Hg-Be11O12	−9.556	−3.035	−6.295	6.521	−10

---

Distributions of the HOMO and LUMO patterns illustrated in Fig. 3 highlighted the localization of nucleophilic and electrophilic sites over the O and M/TM atoms, respectively. Clearly, the redistribution of molecular orbitals of the TM-M₁₁O₁₂ nanocages was indicated in comparison to the M₁₂O₁₂ ones, highlighting the significant impact of the doping process on altering the electronic properties of the M₁₂O₁₂ catalysts.

As compiled in Table 2, the electronic parameters of the M₁₂O₁₂ nanocage were generally changed upon the doping process. For example, the E_HOMO/E_LUMO/E_FL/E_gap values were −10.611/−0.671/−5.641/9.940 and −9.556/−3.035/−6.295/6.521 eV for the Be₁₂O₁₂ and Hg-Be₁₁O₁₂ catalysts, respectively. Clearly, values of E_gap were observed to increase in the order Zn₁₂O₁₂ < Mg₁₂O₁₂ < Be₁₂O₁₂, demonstrating the high reactivity of the Zn₁₂O₁₂ catalyst compared to the Mg₁₂O₁₂ and Be₁₂O₁₂ candidates. Further, the value of E_gap demonstrated the great potential of Zn₁₂O₁₂-based catalysts toward the charge transfer process, which in turn supported the remarkable catalytic activity of the Zn₁₂O₁₂-based catalysts compared to the Mg₁₂O₁₂- and Be₁₂O₁₂-based ones. Moreover, negative % ΔE_gap values were disclosed upon the doping process using group 2B transition metals. Accordingly, such an annotation pinpointed the effective effect of the doping process on upgrading the reactivity of the considered nanocages.

Global indices of reactivity

In the vein of electronic parameters, global indices of reactivity were calculated for the M₁₂O₁₂ catalysts before and after the doping process. Accordingly, various global indices of reactivity, including IP, EA, µ, η, S, ω, and Φ, were assessed and are listed in Table S2.

As could be observed in Table S2, the Zn₁₂O₁₂ catalyst exhibited a higher reactivity character compared to the Mg₁₂O₁₂ and Be₁₂O₁₂ ones via observing a more positive S and less positive η values. This phenomenon might support the higher preferability of the Zn₁₂O₁₂ catalyst toward the charge transfer process and hence higher catalytic activity. Numerically, S values were 0.291, 0.271, and 0.201 eV⁻¹ for the Zn₁₂O₁₂, Mg₁₂O₁₂, and Be₁₂O₁₂ catalysts, respectively. After the doping process, less positive IP and η, along with more positive EA, S, and ω values, were observed. Evidently, values of IP/EA/µ/η/ω/Φ were 10.611/0.671/−5.641/4.970/3.201/5.641 eV and shifted to 9.556/3.035/−6.295/3.260/6.078/6.295 eV for the Be₁₂O₁₂ and Hg-Be₁₁O₁₂ catalysts, respectively. Accordingly, this finding exhibited the ample effects of the doping process on enhancing the electronic features of the M₁₂O₁₂ catalyst.

Adsorption study

In this study, the HDR was processed through three steps. First, the adsorbing ability of M₁₂O₁₂ and TM-M₁₁O₁₂ catalysts toward the H₂ molecule was investigated. The optimized structures of H₂⋯M₁₂O₁₂ and ⋯TM-M₁₁O₁₂ complexes at the most favorable spin state are represented in Fig. 4 and S1, and the corresponding complexation parameters are compiled in Table 3.

Fig. 4 Optimized structures of the H₂⋯M₁₂O₁₂ complexes.

Table 3 Complexation parameters of the H₂⋯M₁₂O₁₂ and ⋯TM-M₁₁O₁₂ complexes at the most favorable spin state throughout the three hydrogen dissociation reaction steps. Distances and energies are given in Å and kcal mol⁻¹, respectively

Step-I: hydrogen adsorption
Complex	Distance			Eads	E^correctedads	E^correctedint	ΔEsolv
	dM/TM–H	dO–H	dH–H^a
a The dH–H was 0.74 Å for the isolated H2 molecule.
H2⋯Zn12O12	2.23	2.88	0.75	−3.65	−3.18	−3.02	−143.38
H2⋯Cd-Zn11O12	2.47	2.80	0.75	−3.27	−3.09	−3.02	−149.80
H2⋯Hg-Zn11O12	2.86	2.50	0.75	−2.27	−2.00	−1.96	−160.84
H2⋯Mg12O12	2.28	2.77	0.75	−3.58	−3.21	−3.21	−111.31
H2⋯Zn-Mg11O12	2.28	2.93	0.75	−3.05	−2.77	−2.77	−113.24
H2⋯Cd-Mg11O12	2.45	3.44	0.75	−2.82	−2.67	−2.67	−119.22
H2⋯Hg-Mg11O12	3.06	2.50	0.75	−1.93	−1.68	−1.75	−131.50
H2⋯Be12O12	2.03	2.58	0.75	−3.04	−2.80	−3.17	−39.85
H2⋯Zn-Be11O12	2.14	3.15	0.75	−5.59	−5.22	−5.40	−49.43
H2⋯Cd-Be11O12	2.31	3.63	0.75	−5.13	−5.00	−5.14	−62.41
H2⋯Hg-Be11O12	2.42	3.91	0.75	−3.97	−3.78	−4.03	−74.65

---

Step-II: hydrogen cleavage
Complex	Distance			Eactivation
	dM/TM–H	dO–H	dH–H^a
H2⋯Zn12O12	1.91	1.32	0.95	13.67
H2⋯Cd-Zn11O12	2.10	1.37	0.93	9.58
H2⋯Hg-Zn11O12	2.18	1.43	0.89	9.84
H2⋯Mg12O12	1.95	1.24	1.03	14.82
H2⋯Zn-Mg11O12	1.92	1.32	0.97	14.80
H2⋯Cd-Mg11O12	2.11	1.37	0.93	10.76
H2⋯Hg-Mg11O12	2.19	1.43	0.90	11.29
H2⋯Be12O12	1.47	1.15	1.12	32.90
H2⋯Zn-Be11O12	1.83	1.30	0.99	18.76
H2⋯Cd-Be11O12	2.01	1.34	0.96	14.94
H2⋯Hg-Be11O12	2.08	1.43	0.91	13.25

---

Step-III: hydrogen migration
Complex	Distance		Ereaction	Ereleased
	dM/TM–H	dO–H
H2⋯Zn12O12	1.66	0.96	−3.81	−17.48
H2⋯Cd-Zn11O12	1.79	0.96	−12.39	−21.97
H2⋯Hg-Zn11O12	1.72	0.96	−26.48	−36.32
H2⋯Mg12O12	1.78	0.98	8.41	−6.41
H2⋯Zn-Mg11O12	1.66	0.96	−3.12	−17.92
H2⋯Cd-Mg11O12	1.79	0.96	−11.71	−22.47
H2⋯Hg-Mg11O12	1.73	0.96	−25.92	−37.21
H2⋯Be12O12	1.37	0.98	24.27	−8.63
H2⋯Zn-Be11O12	1.63	0.96	1.18	−17.58
H2⋯Cd-Be11O12	1.76	0.96	−6.68	−21.62
H2⋯Hg-Be11O12	1.70	0.96	−24.15	−37.40

---

Obviously, in step-I, the adsorption process of the H₂ molecule over all M₁₂O₁₂ and TM-M₁₁O₁₂ catalysts was distinguished by the lack of any structural distortion. In addition, the d_M/TM–H and d_O–H intermolecular distances were observed to be aligned within the range of 2.03–3.06 and 2.50–3.91 Å, respectively (Fig. 4 and S1). Moreover, it was seen that the length of the H–H bond had a relatively small extension to 0.75 Å, compared to its counterpart for the isolated H₂ molecule (d_H–H = 0.74 Å). Such variation was attributed to the emerging interactions between the M/TM and O atoms with the H₂ molecule. This process entirely resulted in a weakening of the H–H bond, which led to a slight activation for the H₂ molecule over the utilized M₁₂O₁₂ and TM-M₁₁O₁₂ catalysts.

From an energetic perspective, the calculated E_ads and E^corrected_ads values ranged from −1.93 to −5.59 and from −1.68 to −5.22 kcal mol⁻¹, respectively, pinpointing the potentiality of M₁₂O₁₂ and TM-M₁₁O₁₂ catalysts to adsorb the H₂ molecule (Table 3). Furthermore, negative E^corrected_int for all complexes with values up to −5.40 kcal mol⁻¹ were also observed. Notably, for the H₂⋯M₁₂O₁₂ complexes, the adsorbing ability was denoted to increase in the following order: H₂⋯Zn₁₂O₁₂ ≈ ⋯Mg₁₂O₁₂ > ⋯Be₁₂O₁₂ complexes. Regarding the H₂⋯TM-M₁₁O₁₂ complexes, negative values of E_ads/E^corrected_ads were detected to decrease on going from TM = Zn to Cd and Hg atoms. For example, E_ads/E^corrected_ads values were −3.65/−3.18, −3.27/−3.09, and −2.27/−2.00 kcal mol⁻¹ for H₂⋯Zn-, Cd-, and Hg-Zn₁₁O₁₂ complexes, respectively. Furthermore, negative ΔE_solv values were observed for all H₂⋯M₁₂O₁₂ and ⋯TM-M₁₁O₁₂ complexes, announcing the crucial role of the solvent effect in enhancing the investigated interactions.

Hydrogen dissociation reaction (HDR)

In order to evaluate the catalytic performance of the M₁₂O₁₂ and TM-M₁₁O₁₂ nanocages as SACs toward the hydrogen storage process, their reaction pathways were generated with the help of IRC calculations. Initially, the TS structures for H₂⋯M₁₁O₁₂ and ⋯TM-M₁₁O₁₂ complexes at the lowest four possible spin states were located. In the IRC vein, all the localized TS structures were authenticated to connect their corresponding minima along the imaginary vibration mode adequately. The calculated E_relative values for H₂⋯TM-M₁₁O₁₂ complexes at the lowest four possible spin states are listed in Table S3. The reaction pathways of the H₂⋯Cd-Zn₁₁O₁₂, ⋯Cd-Mg₁₁O₁₂, and ⋯Hg-Be₁₁O₁₂ complexes, as the most proficient ones, are displayed in Fig. 5. The reaction pathways of the remaining complexes are delineated in Fig. S2. Table 3 compiles the complexation parameters of all the studied complexes.

Fig. 5 Reaction pathway of the catalyzed HDR for the (a) H₂⋯Cd-Zn₁₁O₁₂, (b) H₂⋯Cd-Mg₁₁O₁₂, and (c) H₂⋯Hg-Be₁₁O₁₂ complexes at the most favorable spin state. R, I, TS, and P symbols refer to reactants, intermediate, transition state, and product structures, respectively.

From Table S3, the most favorable spin state for the H₂⋯M-M₁₁O₁₂ complexes was selected as the one with the lowest activation energy barrier (E_activation). Accordingly, the singlet spin state was the most proficient one for all complexes. This observation coincides with its analog for the isolated TM-Zn₁₂O₁₂ catalysts (Table S1).

In step-II, the TS structures of all H₂⋯M₁₂O₁₂ and ⋯TM-M₁₁O₁₂ complexes pinpointed the occurrence of heterolytic hydrogen cleavage from the molecular form to the atomic one (i.e., ) (Fig. 5 and S2). As a consequence, a remarkable increase in the d_H–H value was noticed to range from 0.89 to 1.12 Å, announcing that the H₂ was entirely activated over the M₁₂O₁₂ and TM-M₁₁O₁₂ catalysts. Further, the d_M/TM–H and d_O–H ones declined till they ranged in the 1.47–2.19 and 1.15–1.43 Å distances. Generally, the E_activation was aligned within the 9.58–32.90 kcal mol⁻¹ energetic scope. Comparatively, lower E_activation values were noticed for the H₂⋯TM-complexes compared to the H₂⋯M₁₂O₁₂ ones, announcing the significant role of the doping process in reducing the E_activation values. Moreover, E_activation values were detected to diminish in the order H₂⋯Be₁₂O₁₂/⋯TM-Be₁₁O₁₂ > H₂⋯Mg₁₂O₁₂/⋯TM-Mg₁₁O₁₂ > H₂⋯Zn₁₂O₁₂/⋯TM-Zn₁₁O₁₂, announcing the remarkable catalytic behavior of Zn₁₂O₁₂ and TM-Zn₁₁O₁₂ catalysts compared to their analogs. Among all the investigated complexes, the minimum E_activation value was observed in the case of the H₂⋯Cd-Zn₁₁O₁₂ complex, highlighting that the Cd-Zn₁₁O₁₂ catalyst produced the best catalytic performance. Notably, the significant E_activation value relevant to the H₂⋯Be₁₂O₁₂ complex (i.e., poor catalytic performance) could be interpreted as owing to the low chemical reactivity of the Be₁₂O₁₂ catalyst with a high E_gab value.

In step-III, one of the dissociated 2H* atoms migrated to the M/TM atoms and the other to the O atom. For the product structures, values of d_M/TM–H and d_O–H were arranged from 1.37 to 1.79 and from 0.96 to 0.98 Å, respectively. This hydrogen migration step was accompanied by releasing a significant amount of energy in the scope from −6.41 to −37.40 kcal mol⁻¹. Overall, the obtained findings showed the considerable potentiality of M₁₂O₁₂ and TM-M₁₁O₁₂ catalysts toward HDR catalysis.

As an exemplary model, the complexation parameters of the H₂⋯Cd-Zn₁₁O₁₂ complex were demonstrated due to the ample catalytic performance of the Cd-Zn₁₁O₁₂ catalyst. As shown in Fig. 5, the H–H bond distance was prolonged from 0.74 to 0.75 Å after the hydrogen adsorption step. Subsequently, TS structures manifested the existence of the H₂ molecule heterolytic cleavage to the atomic form (i.e., ), causing an increase in the H–H bond from 0.75 to 0.93 Å via an E_activation of 9.58 kcal mol⁻¹. In addition, the d_TM–H and d_O–H decreased from 2.47 to 2.10 Å and from 2.80 to 1.37 Å, respectively. These findings delineated the full activation of the H₂ molecule over the Cd-Zn₁₁O₁₂ catalyst. Finally, the two activated hydrogen atoms independently migrated to the Cd and O atoms with d_TM–H and d_O–H of 1.79 and 0.96 Å, respectively. Upon the product formation, an energy of −21.97 kcal mol⁻¹ was released. It is worth noting that the product structure was more favorable than the intermediate one (E_reaction = −12.39 kcal mol⁻¹).

Overall, for all H₂⋯M₁₂O₁₂ and ⋯TM-M₁₁O₁₂ complexes, the calculated E_bind and E_activation values were generally observed to be close to the accessible ranges for reversible hydrogen storage, as previously announced.^40,56–58 This observation demonstrates the reversibility of the hydrogen storage process and hence its experimental feasibility.

NBO and EDD analyses

NBO and EDD investigations are performed to provide further characterization for the activation and dissociation processes of the hydrogen molecule over the surfaces of the M₁₂O₁₂ and TM-M₁₁O₁₂ catalysts. Table 4 includes the calculated NBO charges of M/TM, O, H1, and H2 atoms for the TS structures of the H₂⋯M₁₂O₁₂ and ⋯TM-M₁₁O₁₂ complexes at the most favorable spin state. EDD isosurfaces for the investigated complexes are shown in Fig. 6.

Table 4 NBO charges (in e) of the TS structures relevant to the studied H₂⋯M₁₂O₁₂ and ⋯TM-M₁₁O₁₂ complexes at the most favorable spin state

Complex	M/TM	H1^a	O	H2^a
a The H1 and H2 represent the H atoms interacting with TM and O ones, respectively.
H2⋯Zn12O12	1.464	−0.287	−1.493	0.271
H2⋯Cd-Zn11O12	1.594	−0.280	−1.520	0.239
H2⋯Hg-Zn11O12	1.573	−0.243	−1.533	0.214
H2⋯Mg12O12	1.462	−0.364	−1.482	0.308
H2⋯Zn-Mg11O12	1.433	−0.296	−1.505	0.271
H2⋯Cd-Mg11O12	1.570	−0.288	−1.535	0.239
H2⋯Hg-Mg11O12	1.545	−0.253	−1.547	0.216
H2⋯Be12O12	0.937	−0.262	−1.180	0.433
H2⋯Zn-Be11O12	1.468	−0.287	−1.266	0.305
H2⋯Cd-Be11O12	1.590	−0.295	−1.282	0.275
H2⋯Hg-Be11O12	1.568	−0.241	−1.290	0.238

---

Fig. 6 EDD isosurfaces of the TS structures relevant to the studied H₂⋯M₁₂O₁₂ and ⋯TM-M₁₁O₁₂ complexes at the most favorable spin state.

Notably, the M/TM atoms had the maximum positive NBO charges within all H₂⋯M₁₂O₁₂ and ⋯TM-M₁₁O₁₂ complexes which revealed their electropositive nature, whereas the H1 atom had a negative NBO charge (Table 4). Besides, the O and H2 atoms had negative and positive NBO charges, respectively. These annotations affirmed the presence of charge transfer from the M/TM and H2 atoms to the H1 and O candidates, respectively. For instance, charges of the Cd, H1, O, and H2 atoms were 1.594, −0.280, −1.520, and 0.239e, respectively, in the case of the H₂⋯Cd-Zn₁₁O₁₂ complex. A great consistency between the NBO charge-based results and the E_activation values was noticed. The highest and lowest positive NBO charges of M/TM transition metals were observed in the case of the H₂⋯Cd-Zn₁₁O₁₂ (i.e., minimum E_activation) and H₂⋯Be₁₂O₁₂ (i.e., maximum E_activation) complexes, respectively. Numerically, NBO charges were 1.594 and 0.937e, along with E_activation of 9.58 and 32.90 kcal mol⁻¹ for the H₂⋯Cd-Zn₁₁O₁₂ and ⋯Be₁₂O₁₂ complexes, respectively. Furthermore, the amount of NBO charges relevant to the M/TM atoms within the investigated complexes was found to increase in the following sequence: H₂⋯Be₁₂O₁₂ < ⋯Mg₁₂O₁₂ < ⋯Zn₁₂O₁₂ complexes. This progressive increase in charge transfer facilitated H–H bond cleavage, thereby highlighting the significant catalytic activity of the Zn₁₂O₁₂ catalyst compared to the Mg₁₂O₁₂ and Be₁₂O₁₂ catalysts.

From Fig. 6, for all H₂⋯M₁₂O₁₂ and ⋯TM-M₁₁O₁₂ complexes, the EDD isosurfaces outlined the existence of sea green-coded patches on the surfaces of the TM and H2 atoms, suggesting their rich electron density character. The opposite finding was true in the case of O and H1 atoms that were adorned with gray-coded patches. Accordingly, the EDD findings highlighted the existence of M/TM → H1 and H2 → O charge transfer, coincident with the NBO claims. This phenomenon facilitated filling up the antibonding orbitals of the H₂ molecule and hence the hydrogen cleavage process over the surfaces of M₁₂O₁₂ and TM-M₁₁O₁₂ catalysts. Overall, these observations outlined the significant role of M₁₂O₁₂ and TM-M₁₁O₁₂ SACs in catalyzing the HDR.

QTAIM analysis

Analysis of QTAIM has been recognized as a reliable tool that provides a comprehensive understanding of the characteristics of intra- and inter-molecular interactions.^59,60 Portrayals of QTAIM for TS structures relevant to the H₂⋯M₁₂O₁₂ and ⋯TM-M₁₁O₁₂ complexes at the most favorable spin state are displayed in Fig. S3. QTAIM topological parameters are tabulated in Table S4.

The QTAIM diagrams of the TS structures shed light on the presence of two BPs within all the H₂⋯M₁₂O₁₂ and ⋯TM-M₁₁O₁₂ complexes (Fig. S3). The first BP was formed between the M/TM atom and its interacting H one, while the other BP was observed between the O atom and its interacting H counterpart. Such annotations announced the potency of the H₂ molecule to preferably interact with the M₁₂O₁₂ and TM-M₁₁O₁₂ catalysts.

From Table S4, the topological parameters relevant to the BCPs at the M/TM–H and O–H bonds within the H₂⋯M₁₂O₁₂ and ⋯TM-M₁₁O₁₂ complexes could be generally concluded as follows: low positive ρ_b, negative H_b, and −G_b/V_b < 1 values were noted. Numerically, values of ρ_b, H_b, and −G_b/V_b for the BCP at the Cd–H bond of the H₂⋯Cd-Zn₁₁O₁₂ complex were 0.0449, −0.0040, and 0.9087 au, whereas they were 0.1125, −0.0553, and 0.5858 au for the BCP at the O–H bond. The procured findings unveiled the partially covalent nature of the interactions between the H₂ molecule and the M₁₂O₁₂ and TM-M₁₁O₁₂ catalysts. From the literature, this character represented the most favorable circumstance for the catalytic interactions.^26,61

Conclusion

In the essence of searching for a sustainable alternative for nonrenewable energy sources, the current study was dedicated to studying the performance of M₁₂O₁₂ and TM-M₁₁O₁₂ catalysts (i.e., M = Zn, Mg, and Be; TM = Zn, Cd, and Hg) toward the hydrogen storage process. Notably, the most preferable structures for all the studied catalysts were denoted at the singlet spin state. For all the H₂⋯M₁₂O₁₂ and ⋯TM-M₁₁O₁₂ complexes, negative adsorption energies in the ambit from −1.68 to −5.22 kcal mol⁻¹ were noticed, outlining the adsorbing ability of all the studied M₁₂O₁₂ and TM-M₁₁O₁₂ catalysts towards the H₂ molecule. Remarkably, lower E_activation values for the Zn₁₂O₁₂-based catalysts were revealed rather than for their Mg₁₂O₁₂- and Be₁₂O₁₂-based counterparts, highlighting their elevated catalytic activity. Among the investigated catalysts, the Cd-Zn₁₁O₁₂ catalyst demonstrated the highest catalytic activity, where the H₂⋯Cd-Zn₁₁O₁₂ complex was characterized by the lowest activation energy barriers (E_activation) with a value of 9.58 kcal mol⁻¹. Results of the natural bond orbitals and electron density difference analysis unveiled the occurrence of M/TM → H1 and H2 → O charge transfer. Further, the quantum theory of atoms in molecules confirmed the partially covalent nature of the interactions between the interacting species at the TS, which was announced as the optimum catalytic conditions. These observations confirmed that the investigated M₁₂O₁₂ and TM-M₁₁O₁₂ nanocages have a high potential to act as single-atom catalysts (SACs) for H–H bond dissociation rather than serving merely as adsorption sites. The procured results will help experimentalists in designing highly effective SACs for hydrogen storage applications.

Author contributions

Mohammed N. I. Shehata: methodology, formal analysis, investigation, data curation, visualization, writing—original drafts. Lamiaa A. Mohamed: supervision, writing—review and editing. Hu Yang: methodology, writing—review and editing. Tamer Shoeib: supervision, conceptualization, software, writing—review and editing. Jabir H. Al-Fahemi: methodology, resources, writing—review and editing. Mahmoud A. A. Ibrahim: conceptualization, methodology, software, resources, project administration, supervision, writing—review and editing. All authors have read and agreed to the published version of the manuscript.

Conflicts of interest

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

Data availability

The data supporting this article have been included as part of the supplementary information (SI). Supplementary information: reaction pathway of the HDR catalysis for the H₂⋯M₁₂O₁₂ and ⋯TM-M₁₁O₁₂ complexes at the most favorable spin state. Optimized structures of the H₂⋯TM-M₁₁O₁₂ complexes at the most favorable spin state. QTAIM diagrams of the TS structures relevant to the H₂⋯M₁₂O₁₂ and ⋯TM-M₁₁O₁₂ complexes at the most favorable spin state. E_relative of TM-M₁₁O₁₂ nanocages as SACs. Global indices of reactivity for the M₁₂O₁₂ and TM-M₁₁O₁₂ catalysts at the most favorable spin state. E_relative of TM-M₁₁O₁₂ nanocages and the utilized H₂⋯TM-M₁₁O₁₂ complexes at the lowest four possible spin states. Topological parameters (in au) of the TS structures relevant to the H₂⋯M₁₂O₁₂ and ⋯TM-M₁₁O₁₂ complexes at the most favorable spin state. See DOI: https://doi.org/10.1039/d5ra06216k.

Acknowledgements

The authors extend their appreciation to Umm Al-Qura University, Saudi Arabia for funding this research work through grant number: 25UQU4200274GSSR07.

References

1. X. F. Yang, A. Wang, B. Qiao, J. Li, J. Liu and T. Zhang, Acc. Chem. Res., 2013, 46, 1740–1748.
2. J. Hagen, Industrial Catalysis: a Practical Approach, John Wiley & Sons, 2015.
3. B. Qiao, A. Wang, X. Yang, L. F. Allard, Z. Jiang, Y. Cui, J. Liu, J. Li and T. Zhang, Nat. Chem., 2011, 3, 634–641.
4. Y. L. Guo, R. Lang and B. T. Qiao, Catalysts, 2019, 9, 135.
5. M. Xiao, L. Gao, Y. Wang, X. Wang, J. Zhu, Z. Jin, C. Liu, H. Chen, G. Li, J. Ge, Q. He, Z. Wu, Z. Chen and W. Xing, J. Am. Chem. Soc., 2019, 141, 19800–19806.
6. T. K. Ghosh and N. N. Nair, ChemCatChem, 2013, 5, 1811–1821.
7. X. Liu, Y. Yang, M. M. Chu, T. Duan, C. G. Meng and Y. Han, Catal. Sci. Technol., 2016, 6, 1632–1641.
8. X. Li, X. Yang, Y. Huang, T. Zhang and B. Liu, Adv. Mater., 2019, 31, e1902031.
9. S. Zhang, X. Zhang, Y. Rui, R. H. Wang and X. J. Li, Green Energy Environ., 2021, 6, 458–478.
10. W. Zhao, Z. Chen, X. R. Yang, X. X. Qian, C. X. Liu, D. T. Zhou, T. Sun, M. Zhang, G. Y. Wei, P. D. Dissanayake and Y. S. Ok, Renewable Sustainable Energy Rev., 2020, 132, 110040.
11. N. C. Bigall, M. Reitzig, W. Naumann, P. Simon, K. H. van Pée and A. Eychmüller, Angew. Chem., Int. Ed., 2008, 47, 7876–7879.
12. G. S. Parkinson, Catal. Lett., 2019, 149, 1137–1146.
13. N. C. Cheng, L. Zhang, K. Doyle-Davis and X. L. Sun, Electrochem. Energy Rev., 2019, 2, 539–573.
14. J. X. Liang, X. F. Yang, A. Q. Wang, T. Zhang and J. Li, Catal. Sci. Technol., 2016, 6, 6886–6892.
15. A. B. Shah, S. Sarfaraz, M. Yar, N. S. Sheikh, H. H. Hammud and K. Ayub, Nanomaterials, 2022, 13, 29.
16. S. X. Liang, C. Hao and Y. T. Shi, ChemCatchem, 2015, 7, 2559–2567.
17. D. Q. Liu, A. Barbar, T. Najam, M. S. Javed, J. Shen, P. Tsiakaras and X. K. Cai, Appl. Catal., B, 2021, 297, 120389.
18. F. Chen, X. Z. Jiang, L. L. Zhang, R. Lang and B. T. Qiao, Chin. J. Catal., 2018, 39, 893–898.
19. M. B. Gawande, P. Fornasiero and R. Zboril, ACS Catal., 2020, 10, 2231–2259.
20. X. Wang, W. Chen, L. Zhang, T. Yao, W. Liu, Y. Lin, H. Ju, J. Dong, L. Zheng, W. Yan, X. Zheng, Z. Li, X. Wang, J. Yang, D. He, Y. Wang, Z. Deng, Y. Wu and Y. Li, J. Am. Chem. Soc., 2017, 139, 9419–9422.
21. N. J. O’Connor, A. S. M. Jonayat, M. J. Janik and T. P. Senftle, Nat. Catal., 2018, 1, 531–539.
22. M. Ghorbanzadeh Ahangari and A. Hamed Mashhadzadeh, Int. J. Hydrogen Energy, 2020, 45, 6745–6756.
23. T. Fujitani, I. Nakamura, T. Akita, M. Okumura and M. Haruta, Angew Chem. Int. Ed. Engl., 2009, 48, 9515–9518.
24. S. Sarfaraz, M. Yar, N. S. Sheikh, I. Bayach and K. Ayub, ACS Omega, 2023, 8, 14077–14088.
25. S. Sarfaraz, M. Yar, A. Hussain, A. Lakhani, A. Gulzar, M. Ans, U. Rashid, M. Hussain, S. Muhammad, I. Bayach, N. S. Sheikh and K. Ayub, ACS Omega, 2023, 8, 36493–36505.
26. S. Sarfaraz, M. Yar, R. Hussain and K. Ayub, Int. J. Hydrogen Energy, 2023, 48, 20633–20645.
27. A. Hussain, M. Yar, R. S. Alshareef, T. Mahmood, K. Ayub, R. Nawaz, M. Yasin and M. A. Gilani, Int. J. Hydrogen Energy, 2024, 58, 806–821.
28. I. Bayach, S. Sarfaraz, N. S. Sheikh, K. Alamer, N. Almutlaq and K. Ayub, Materials, 2023, 16, 2792.
29. W. Li, S. E. Madan, R. Reocreux and M. Stamatakis, ACS Catal., 2023, 13, 15851–15868.
30. Z. K. Han, D. Sarker, R. Ouyang, A. Mazheika, Y. Gao and S. V. Levchenko, Nat. Commun., 2021, 12, 1833.
31. G. W. T. J. Frisch, H. B. Schlegel, G. E. Scuseria, M. A. Robb, J. R. Cheeseman, G. Scalmani, V. Barone, B. Mennucci, G. A. Petersson, H. Nakatsuji, M. Caricato, X. Li, H. P. Hratchian, A. F. Izmaylov, J. Bloino, G. Zheng, J. L. Sonnenberg, M. Hada, M. Ehara, K. Toyota, R. Fukuda, J. Hasegawa, M. Ishida, T. Nakajima, Y. Honda, O. Kitao, H. Nakai, T. Vreven, J. A. Montgomery Jr, J. E. Peralta, F. Ogliaro, M. Bearpark, J. J. Heyd, E. Brothers, K. N. Kudin, V. N. Staroverov, R. Kobayashi, J. Normand, K. Raghavachari, A. Rendell, J. C. Burant, S. S. Iyengar, J. Tomasi, M. Cossi, N. Rega, J. M. Millam, M. Klene, J. E. Knox, J. B. Cross, V. Bakken, C. Adamo, J. Jaramillo, R. Gomperts, R. E. Stratmann, O. Yazyev, A. J. Austin, R. Cammi, C. Pomelli, J. W. Ochterski, R.L. Martin, K. Morokuma, V. G. Zakrzewski, G. A. Voth, P. Salvador, J. J. Dannenberg, S. Dapprich, A. D. Daniels, Ö. Farkas, J. B. Foresman, J. V. Ortiz, J. Cioslowski and D. J. Fox, Gaussian 09, Revision E.01, Gaussian, Inc., Wallingford CT, 2009.
32. Y. Zhao and D. G. Truhlar, Theor. Chem. Acc., 2008, 120, 215–241.
33. P. J. Hay and W. R. Wadt, J. Chem. Phys., 1985, 82, 270–283.
34. Y. Zhao and D. G. Truhlar, J. Chem. Theory Comput., 2008, 4, 1849–1868.
35. X. L. Ding, W. Xue, Y. P. Ma, Z. C. Wang and S. G. He, J. Chem. Phys., 2009, 130, 014303.
36. A. Irfan, F. K. Al-Zeidaneen, I. Ahmed, A. G. Al-Sehemi, M. A. Assiri, S. Ullah and G. Abbas, Bull. Mater. Sci., 2020, 43, 45.
37. K. W. Paul, J. D. Kubicki and D. L. Sparks, Eurasian J. Soil Sci., 2007, 58, 978–988.
38. T. H. Dunning and P. J. Hay, in Methods of Electronic Structure Theory, ed. H. F. Schaefer, Springer US, Boston, MA, 1977, ch. 1, pp. 1–27,  DOI:10.1007/978-1-4757-0887-5_1.
39. M. Bursch, J. M. Mewes, A. Hansen and S. Grimme, Angew Chem. Int. Ed. Engl., 2022, 61, e202205735.
40. C. M. Ramos-Castillo, J. U. Reveles, M. E. Cifuentes-Quintal, R. R. Zope and R. de Coss, J. Phys. Chem. C, 2016, 120, 5001–5009.
41. S. A. Aravindh, I. S. Roqan and H. Alawadhi, J. Cluster Sci., 2021, 32, 55–62.
42. S. A. Aravindh, U. Schwingenschlögl and I. S. Roqan, J. Appl. Phys., 2014, 116, 233906.
43. M. A. A. Ibrahim, J. Mol. Model., 2012, 18, 4625–4638.
44. M. Ernzerhof, J. Chem. Theory Comput., 2009, 5, 793–797.
45. S. F. Boys and F. Bernardi, Mol. Phys., 1970, 19, 553–566.
46. T. Lu and F. Chen, J. Comput. Chem., 2012, 33, 580–592.
47. W. Humphrey, A. Dalke and K. Schulten, J. Mol. Graph., 1996, 14, 33–38.
48. J. S. Murray and P. Politzer, Wiley Interdiscip. Rev. Comput. Mol. Sci., 2011, 1, 153–163.
49. P. K. Weiner, R. Langridge, J. M. Blaney, R. Schaefer and P. A. Kollman, Proc. Natl. Acad. Sci. U. S. A., 1982, 79, 3754–3758.
50. A. A. Peyghan, M. T. Baei and S. Hashemian, J. Cluster Sci., 2013, 24, 341–347.
51. H. Y. Ammar, K. M. Eid and H. M. Badran, Results Phys., 2022, 35, 105349.
52. Q. Ali, A. Shakoor, G. Rehman, M. Ur Rehman, M. Khan, R. Ahmad, I. Ahmad, A. F. AlAsmari and F. Alasmari, Sci. Rep., 2025, 15, 1020.
53. H. M. Badran, K. M. Eid, S. Baskoutas and H. Y. Ammar, Nanomaterials, 2022, 12, 1757.
54. K. Fukui, T. Yonezawa and H. Shingu, J. Chem. Phys., 2004, 20, 722–725.
55. K. Fukui, Acc. Chem. Res., 1971, 4, 57.
56. H. Cui, Y. Zhang, W. Tian, Y. Wang, T. Liu, Y. Chen, P. Shan and H. Yuan, RSC Adv., 2021, 11, 13912–13918.
57. C. I. Contescu, K. van Benthem, S. Li, C. S. Bonifacio, S. J. Pennycook, P. Jena and N. C. Gallego, Carbon, 2011, 49, 4050–4058.
58. A. Sharma and B. Chakraborty, Int. J. Hydrogen Energy, 2024, 51, 1306–1313.
59. J. Contreras-Garcia, E. R. Johnson, S. Keinan, R. Chaudret, J. P. Piquemal, D. N. Beratan and W. Yang, J. Chem. Theory Comput., 2011, 7, 625–632.
60. P. Panini, R. G. Gonnade and D. Chopra, New J. Chem., 2016, 40, 4981–5001.
61. D. R. Aireddy and K. L. Ding, ACS Catal., 2022, 12, 4707–4723.

---