## 2 Methodology  
This research utilizes the CASTEP module within Materials Studio, leveraging spin-polarized Density Functional Theory (DFT) and the Ultrasoft Pseudopotentials (USPPs) method to accurately simulate the adsorption behavior of gaseous hydrogen on iron surfaces [31; 32]. DFT, serving as a core computational tool, streamlines the complexity of electron systems, thus reducing computational difficulty and establishing it as an ideal choice for studying molecular structures and properties [33]. The USPPs method achieves an accurate simulation of complex electronic structures and enhances computational efficiency by describing the interaction between nuclei and electrons with a softer pseudopotential [34; 35]. Meanwhile, the BFGS algorithm exhibits its efficiency and accuracy in optimizing the BCC lattice structure of iron, effectively accelerating the prediction process of lattice parameters and physical properties by approximating the Hessian matrix, thus becoming a preferred tool for analyzing complex lattice structures [36]. The simple 1s1 electron configuration of hydrogen atoms enables them to form chemical bonds with other atoms through electron sharing. For Fe, its electron configuration is Fe:1s2s22p63s23p63d64s2, with the valence electrons primarily located in the 3d and 4s orbitals. In studying the electronic interaction between hydrogen and iron, the USPPs method effectively manages this complex interaction [37]. This study selected the stable body-centered cubic structure of \(\alpha\)-Fe as the subject material, specifically focusing on the atomic arrangement characteristics of the Fe(110) crystal plane. To simulate the dissociative adsorption process of hydrogen molecules (H\({}_{2}\)) on the Fe(110) surface, a periodic cell model was employed. Computational parameters were carefully chosen, including a plane wave cutoff energy of 350 eV and a k-point grid size of 12 \(\times\) 12 \(\times\) 12. Brillouin zone integration was executed utilizing the Monkhorst-Pack scheme, with the aim of ensuring the high accuracy of the computational findings. The convergence criterion for self-consistent field (SCF) calculations was set at a rigorous threshold of 2.0 \(\times\) 10\({}^{-6}\) eV/atom, underscoring the precision of the calculation efforts. Tosimulate the real surface environment more accurately, a 15 A vacuum layer was introduced above the Fe(110) surface, mirroring conditions closer to those in experimental settings. The geometric optimization results indicated that the optimized lattice constant a of \(\alpha\)-Fe is 2.844 A, representing only a 0.77% deviation from the experimental value of 2.866 A, thereby verifying the high accuracy of the model setup [38]. Furthermore, a model comprising 7 layers of Fe(110) crystal planes was constructed, with a 15 A vacuum layer set to accurately simulate the surface environment. In this model, the positions of the bottom four layers of Fe atoms were fixed, whereas the top three layers of Fe atoms and the atoms adsorbed on them were permitted to relax, to more closely simulate real-world conditions. All simulation parameters underwent convergence tests to guarantee the reliability of this study's results. The computational workflow consists of five sequential stages. First, an initial geometry is constructed by placing a physisorbed H\({}_{2}\) molecule above the Fe(110) surface, followed by a pre-relaxation step at the Density Functional Theory (DFT) level. Second, static adsorption energies at 0 K are calculated for all high-symmetry adsorption sites--including top, bridge, and hollow positions--to identify the most thermodynamically favorable configurations. Third, the lowest-energy structures are subjected to full ionic relaxation to obtain fully optimized adsorbate-surface geometries. Fourth, finite-temperature and finite-pressure corrections are incorporated by combining the 0 K adsorption enthalpy with the temperature- and pressure-dependent chemical potential of H\({}_{2}\), yielding operational Gibbs free energies.