Mechanistic Study of the Electrocatalytic Carbon Dioxide Reduction Reaction Over Boron/Nitrogen Co-Doped Graphene-Supported Single-Atom Catalysts

1. Introduction

Excessive consumption of fossil fuels has led to a rapid increase in atmospheric CO₂ concentration, thereby intensifying global climate change. In this context, the electrocatalytic carbon dioxide reduction reaction (CO₂RR) has attracted extensive attention as a sustainable technology for converting CO₂ into value-added chemicals[1]. However, the high thermodynamic stability of the C=O bonds in CO₂ imposes substantial energy barriers on its chemical conversion, leading to sluggish CO₂ reduction kinetics. In addition, competition among multiple possible products, together with interference from the hydrogen evolution reaction (HER), gives rise to major challenges, including poor product selectivity and low Faradaic efficiency (FE) [1,2]. In recent years, metal-based materials, including pure metals, alloys, and metal oxides, have been widely investigated in CO₂RR catalytic systems and have served as key catalysts for improving reaction performance. Although notable progress has been achieved, conventional metal-based catalysts still suffer from relatively high overpotentials, which not only lower energy utilization efficiency but also promote particle aggregation during the reaction, thereby compromising catalyst stability [45]. Therefore, the development of new high-efficiency CO₂RR catalytic materials remains of great significance.

Single-atom catalysts (SACs) have shown remarkable advantages in electrocatalysis, particularly in CO₂RR, owing to their exceptionally high atomic utilization efficiency and precisely tunable active-site configurations. The catalytic activity and product selectivity of SACs are governed primarily by the type of metal center and its coordination environment. Accordingly, catalytic performance can be effectively regulated through rational design of the coordination structure. In recent years, N-doped graphene and its derivatives have emerged as ideal supports for SACs because they provide abundant coordination sites for firmly anchoring isolated metal atoms, while also offering high specific surface area, excellent electrical conductivity, and outstanding chemical stability [4,10].

Significant progress has recently been made in CO₂RR studies based on such systems. For example, the five-coordinated Fe-N₅/C catalyst prepared by Wang et al. [7] exhibited outstanding performance in regulating syngas composition. Through theoretical calculations, Yang et al. [48] revealed the critical role of spin polarization in determining the selectivity of M-N₄@Gr structures and found that Ni-N₄@Gr shows high selectivity toward HCOOH. Bai et al. [49] proposed a strategy that simultaneously regulates the central metal, axial coordination, and nonmetal coordination environment, successfully breaking the scaling relationships of carbon-based intermediates in CO₂RR. The screened systems, such as Ti-N₄-B and Ru-C₂NO-C, reduced the required applied voltage by 45%–71% relative to the current state-of-the-art level. Nevertheless, despite the great potential demonstrated by N-doped SACs, most existing studies remain limited to single planar coordination structures, such as M-N₃, M-N₄, and M-N₅. This limitation often results in insufficient control over the adsorption strength of key intermediates, such as COOH and OCHO, making it difficult to balance CO₂ activation with the synergistic optimization of multistep reaction intermediates. As a result, further improvement in catalytic performance remains restricted. Therefore, developing synergistic regulation strategies beyond the traditional single-coordination mode is essential for the design of a new generation of highly efficient CO₂RR catalysts.

To address these limitations, heteroatom doping has provided a new route for optimizing CO₂RR catalyst performance [50]. Previous studies have shown that the introduction of heteroatoms such as nitrogen (N) and boron (B) can regulate the adsorption behavior of key intermediates by modifying the local electronic structure and surface chemical environment of the catalyst, thereby affecting CO₂RR activity and selectivity. Specifically, N doping can alter the coordination environment and electronic distribution of single-atom active centers, thereby tuning product selectivity [48]. In contrast, B doping can optimize the adsorption behavior of CO₂ and its intermediates by changing the charge distribution and electronic structure around the metal sites, thus enhancing catalytic activity and stability [51]. In related work, Yan et al. [53] prepared an asymmetric Ni-N3-B catalyst supported on B, N-codoped carbon (Ni-SACs@BNC) by calcination, using Ni-SACs@NC with a symmetric Ni-N4 structure as the reference. The results showed that, under visible-light irradiation, the catalyst achieved a CO production rate of 37.1 μmol/h with a selectivity of 93.2% and operated stably for 24 h. Experimental characterization and DFT calculations further demonstrated that B doping broke the charge symmetry of the Ni center, enhanced electron-capturing capability, and shifted the d-band center upward, thereby strengthening CO₂ adsorption and lowering the energy barrier for *COOH formation. Liu et al. [52] constructed an Fe single-atom catalyst supported on B, N-codoped carbon and applied it to the electrocatalytic reduction of CO₂. Their study showed that the introduction of B regulated the coordination environment and electronic structure of the Fe active site, which enhanced CO₂RR activity while suppressing the competing HER. Overall, N doping, B doping, and B,N co-doped supports can effectively regulate the electronic structure of single-atom sites and the adsorption behavior of intermediates, further indicating that B, N-codoped graphene is a suitable support for SACs in CO₂RR.

Inspired by these studies, this work uses transition-metal single atoms supported on B, N -codoped graphene (TM@BNG) as a model system and systematically screens the CO₂RR catalytic performance of different transition-metal active site. Based on density functional theory (DFT) calculations, the free-energy evolution and reaction-path characteristics of a series of systems during CO₂ reduction are analyzed in detail. Active sites including Pd, Ag, Sc, Cu, Co, Cd, and Y are identified and confirmed to exhibit excellent activity and selectivity toward the formation of HCOOH or CO. Furthermore, the key descriptor φ is extracted from intrinsic material parameters using the GBR algorithm and the SISSO method, and a volcano-type relationship between φ and the limiting potential (UL) is revealed. These results clarify the regulatory role of the electronic structure of metal active sites in catalytic performance and provide a theoretical reference for the design of highly efficient CO₂RR catalysts.

2. Computational Details

All spin-polarized DFT calculations in this study were performed using the Vienna Ab initio Simulation Package (VASP) [1,11]. The exchange-correlation energy was described by the Perdew-Burke-Ernzerhof (PBE) functional within the generalized gradient approximation (GGA) framework [12], and the plane-wave cutoff energy was set to 450 eV. During structural optimization, a 3×3×1 Monkhorst-Pack k-point mesh was adopted, whereas a denser 9×9×1 k-point mesh was used for electronic structure calculations to ensure computational accuracy [13]. The convergence thresholds for energy and force were set to 1×10^-5 eV and 0.05 eV/Å, respectively. Van der Waals interactions were corrected using the DFT-D3 method [14], and solvent effects were described by an implicit solvation model implemented in the VASPsol program [15,16,17], with the dielectric constant set to 78.4. Bader charge analysis was carried out using the tools developed by the Henkelman group. To eliminate interference from interlayer interactions, a vacuum layer of 20 Å was introduced along the z direction. To evaluate the thermal stability of the catalysts, ab initio molecular dynamics (AIMD) simulations were performed at 400 K using the Nosé-Hoover thermostat, with a time step of 2.0 fs and a total simulation time of 10 ps [18].

The Gibbs free energy change (ΔG) of each elementary step in CO₂RR was calculated according to the following equation[22]:

ΔG = ΔE + ΔZPE – TΔS + ΔGU + ΔGpH

where ΔE represents the energy difference between the products and reactants. ZPE and ΔS denote the zero-point energy correction and entropy change at room temperature (T = 298.15 K), respectively. ΔGU = -neU, where n is the number of transferred electrons and U is the electrode potential. ΔGpH = kBT × ln10 × pH, where kB is the Boltzmann constant. In this work, the pH was set to 0 to simulate an acidic environment.

The limiting potential of CO₂RR was calculated using the following equation:

UL = -ΔGPDS/e

where ΔGPDS is the Gibbs free-energy change of the potential-determining step (PDS) in CO₂RR.

3. Results and Discussion

3.1. Structure and Stability of TM@BNG

The geometric configurations of the SAC models constructed in this study are shown in Figure 1. Because Sc, Ti, V, Cr, Mn, Fe, Co, Ni, Cu, Zn, Y, Zr, Nb, Mo, Ru, Rh, Pd, Ag, Cd, Hf, Ta, W, Re, Os, Ir, Pt, and Au have been widely investigated in the design of CO₂RR catalysts, these 27 metal atoms were selected as candidate active sites. Structural optimization showed that all systems retained a planar configuration after relaxation, and the optimized structures are presented in Fig.S1. To evaluate catalyst stability, both thermodynamic stability and electrochemical stability were considered in this study, and the formation energy (Ef) [20] and dissolution potential (Udiss) [20] were accordingly calculated.

Ef = (Ecatalyst − ETM-bulk − EBNG)

Udiss = U0diss(metal,bulk) − Ef/ne

where Ecatalyst is the total energy of the catalyst system, ETM-bulk is the energy of a single metal atom in the bulk metal, EBNG is the energy of the substrate after removal of the metal atom, U0diss(metal, bulk) is the standard dissolution potential of the bulk metal, and n is the number of electrons transferred in the dissolution reaction [21].

As shown in Figure 2a all single-atom catalysts except Au-SAC exhibit negative formation energies, ranging from -8.42 to -0.13 eV, indicating favorable thermodynamic stability. Meanwhile, the dissolution potentials (Udiss) of these catalysts range from 0.21 V for Hf@BNG to 1.57 V for Rh@BNG, suggesting high structural stability under electrochemical conditions. Therefore, 26 SACs that are stable under both thermodynamic and electrochemical conditions were selected for subsequent evaluation of catalytic performance. Taking Cu@BNG as a representative example, the total density of states analysis (Figure 2band) shows a pronounced distribution of electronic states near the Fermi level, mainly contributed by the Cu 3d orbitals. This feature is favorable for electron transport and for the adsorption and activation of CO₂RR intermediates. Among these 26 SACs, the bond lengths between the transition-metal atoms and N atoms range from 1.78 Å to 2.25 Å, indicating that the metal atoms can be firmly anchored on the BNG substrate. Furthermore, electron localization function analysis (Figure 2cand), again using Cu@BNG as an example, shows that the ELF values in the Cu-N bond region are concentrated in the range of 0.70-0.80, displaying clear covalent-bonding characteristics [23,24]. This result further confirms the structural stability of the TM@BNG systems from the perspective of electronic structure.

3.2. Reaction Mechanism of CO2RR on TM@BNG

As shown in Figure 3, CO₂RR on the surface of SACs proceeds through a series of proton-coupled electron transfer (PCET) steps. Through different reaction pathways, this process can generate a variety of C1 products. The 2e- reduction pathway mainly yields CO and HCOOH, whereas further reduction can produce HCHO, CH₃OH, and CH4, corresponding to 4e-, 6e-, and 8e- processes, respectively. At the initial stage of CO₂RR, the first PCET step is particularly important because it determines the subsequent reaction route. Protonation at different atomic sites leads to distinct intermediates and thus different product selectivities. Specifically, preferential protonation at an O atom forms the *COOH intermediate, whereas protonation at the C atom generates *HCOO, thereby directing the subsequent product distribution. According to the formation of the final C1 products, the 26 SACs investigated in this work can be classified into two categories: one tends to yield only a single product, whereas the other can generate multiple reduction products. The corresponding results are summarized in Figure 4 and Figure 5.

Taking Pd@BNG as an example, the reaction mechanism responsible for its tendency to produce a single product during CO₂RR was systematically analyzed. As shown in Figure 4a, after adsorption of the CO₂ molecule, the *CO₂ intermediate is first formed, followed by an initial hydrogenation step that can proceed through two competing pathways. The free-energy change for *HCOO formation is -0.47 eV, whereas only 0.06 eV is required for *COOH formation, indicating that both intermediates are thermodynamically accessible. In the subsequent steps, dehydration of *COOH to form *CO is spontaneous, whereas *HCOO must undergo an additional hydrogenation step to produce HCOOH, with a corresponding ΔG of 0.57 eV. Alternatively, *HCOO may be further converted into the higher-energy *OCH₂O intermediate (1.21 eV), which is thermodynamically unfavorable. Notably, the desorption energy of *CO on the Pd@BNG surface is low, indicating facile release from the active site. Therefore, the conversion of *CO₂ to *COOH is identified as the potential rate-controlling step for CO formation, corresponding to a limiting potential of -0.06 V. In contrast, for the HCOOH pathway, the potential-determining step is the conversion of *HCOO to HCOOH, which requires a limiting potential of -0.57 V and is thus unfavorable. As hydrogenation proceeds, the system preferentially follows the *CO → *OCH and *OCH₂O → *OCH₂OH pathways. Among these steps, the formation of *OCH requires overcoming a free-energy barrier of 0.18 eV, whereas the formation of *OCH₂O is thermodynamically spontaneous. Subsequently, *OCH₂OH undergoes dehydration to form *OCH₂ (-0.20 eV). In parallel, *OCH can be further hydrogenated to generate *OCH₂ and *CHOH, requiring free-energy changes o 0.21 eV and 0.70 eV, respectively. Here, *OCH₂ corresponds to the adsorbed state of formaldehyde, and its desorption requires only 0.20 eV. During further hydrogenation, *OCH₂ can spontaneously transform into *OCH₃ and *CH₂OH, whereas *CHOH preferentially converts into *CH₂OH (-0.69 eV) rather than following the higher-energy *CH pathway (1.06 eV). For the 6e- product CH₃OH, formation is thermodynamically favorable through either hydrogenation of *OCH₃ (-0.35 eV) or conversion of *CH₂OH (-0.47 eV), and the CH₃OH molecule can also desorb readily from the surface (-0.21 eV). On this basis, *CH₃OH can be further converted into the 8e- product CH₄ through dehydration and hydrogenation, and this process is spontaneous. It should be noted that the formation of HCHO, CH₃OH, and CH₄ is limited by the same rate-controlling step, namely the conversion of *OCH to *OCH₂, with a limiting potential of -0.21 V. Overall, comparison of the limiting potentials for the different product pathways indicates that Pd@BNG preferentially produces CO at low applied potentials (UL = -0.06 V), whereas the formation of other competing products is suppressed. Therefore, this catalyst exhibits pronounced selectivity toward CO.

As shown in Figure 4b, the initial protonation of CO₂ on the Ag@BNG surface,can generate two intermediates, namely *COOH and *HCOO. Among them, the formation of *HCOO requires only 0.09 eV, whereas *COOH formation requires 0.68 eV, indicating that *HCOO is thermodynamically more favorable at the initial stage of the reaction. Continued hydrogenation along the *HCOO pathway yields *HCOOH with a free-energy change of 0.11 eV. Therefore, the conversion of *HCOO to *HCOOH is identified as the potential-determining step (PDS) for HCOOH formation. By contrast, CO formation proceeds through the CO₂ → *COOH → *CO pathway, in which the CO₂ → *COOH step has a ΔG of 0.68 eV and serves as the PDS of this route. The further conversion of *HCOO to *OCH₂O must overcome a free-energy barrier as high as 1.73 eV, making this pathway thermodynamically unfavorable. In the subsequent PCET steps, hydrogenation of *HCOOH to *OCH requires a free-energy input of 0.34 eV, and the resulting *OCH can then be spontaneously hydrogenated to form *OCH₂, *CH₂OH, and *OCH₃ intermediates. Accordingly, the formation of both HCHO and CH3OH is limited by the *HCOOH → *OCH step, with a ΔG of 0.34 eV for each pathway. For the 8e- product CH₄, the optimal reaction pathway is CO₂ → *HCOO → *HCOOH → *OCH → *OCH₂ → *CH₂OH → *CH₂ → *CH₃ → CH₄, in which the conversion of *CH₂OH to *CH₂ is the PDS, corresponding to a ΔG of 0.58 eV. Notably, product desorption on the Ag@BNG surface involves relatively low free-energy changes (ΔG < 0.70 eV) and therefore does not impose additional limitations on the reaction. Overall, the energetic features of the different pathways indicate that HCOOH formation is favored under lower free-energy conditions, making HCOOH the major reduction product on Ag@BNG.

As shown in Figure 4c after CO₂ is adsorbed on the Os@BNG surface to form the *CO₂ intermediate, it can spontaneously convert into *COOH and *HCOO, with ΔG values of -0.13 eV and -0.36 eV, respectively. Along the *COOH pathway, the subsequent formation of *CO is also thermodynamically favorable, giving the sequence CO₂ → *CO₂ → *COOH → *CO. However, *CO desorption requires overcoming a high energy barrier of 2.02 eV, which makes this pathway unfavorable under practical reaction conditions. By contrast, further hydrogenation of *HCOO to *HCOOH requires a free-energy input of 0.91 eV, indicating that this step is energetically demanding within that branch. Considering the unfavorable kinetic or desorption characteristics of these two pathways, the reaction is more likely to proceed through the *OCH₂O intermediate, which involves a lower energy barrier. The hydrogenation of *OCH₂O to *OCH₂OH is exothermic, with a ΔG of -0.26 eV, and the subsequent conversion of *OCH₂OH to *OCH₂ is also thermodynamically favorable, with a ΔG of -0.24 eV. In addition, the reduction of *OCH₂ to *CH₂OH proceeds spontaneously, with a ΔG of -0.30 eV. Further hydrogenation of *CH₂OH to *CH₃OH requires an energy input of 0.48 eV, after which *CH₃OH can desorb readily to yield the final stable product. Throughout the overall CO₂RR process, the formation of other possible products is suppressed because the corresponding pathways involve either high limiting potentials for key intermediates or large desorption energies. Overall, the free-energy profiles indicate that Os@BNG preferentially produces CH₃OH at relatively low applied potentials. The potential-determining step (PDS) is the conversion of *HCOO to *OCH₂O, corresponding to a limiting potential UL of -0.60 V.

As shown in Figure 4d, after adsorption on the Re@BNG surface, CO₂ can form two stable initial intermediates, *COOH and *HCOO, with free-energy changes of -0.85 eV and -1.39 eV, respectively, indicating that both protonated configurations are thermodynamically favorable. Along the *COOH pathway, the subsequent formation of *CO is spontaneous. However, CO desorption from the catalyst surface requires overcoming a relatively high energy barrier of 1.85 eV, which limits the feasibility of this route under practical reaction conditions. Although the *HCOO configuration is more stable, its further hydrogenation to *HCOOH still requires a considerable energy input. Therefore, the reaction is more likely to proceed through the lower-energy *OCH₂O intermediate and then enter the subsequent reduction steps. In the following PCET steps, the hydrogenation of *OCH₂O to *OCH₂OH is exothermic, with a ΔG of -0.41 eV, whereas the subsequent conversion of *OCH₂OH to *OCH₂ requires a free-energy input of 0.31 eV. As the reaction proceeds, *OCH₂ can be further reduced to *CH₂OH with the same energy requirement of 0.31 eV. Thereafter, *CH₂OH undergoes dehydration to form *CH₂, which is thermodynamically favorable, with a ΔG of -0.49 eV. The resulting *CH₂ is then hydrogenated to *CH₃, and the final conversion of *CH₃ to CH₄ requires overcoming an energy input of 0.67 eV. Because other possible pathways involve either relatively high limiting potentials for key intermediates or hindered desorption of certain species, the formation of alternative reduction products on the Re@BNG surface is suppressed. Overall, the energetic features of the reaction pathways indicate that Re@BNG preferentially produces CH₄ at relatively low applied potentials. The final PCET step (*CH₃ → CH₄) is the PDS, corresponding to a limiting potential UL of -0.67 V.

As shown in Figure 5a, after CO₂ is adsorbed on the Pt@BNG surface, it can be protonated to form two intermediates, *COOH and *HCOO. Among them, *COOH can spontaneously convert into *CO, whereas the conversion of *HCOO to *HCOOH requires overcoming an energy barrier of 0.88 eV. Therefore, this branch is suppressed at the initial stage of the reaction. The resulting *CO can then be further hydrogenated to *OCH with a free-energy input o 0.36 eV, followed by subsequent conversion to *OCH₂. It should be noted that the desorption energies of CO and HCHO on the Pt@BNG surface are as high as 1.53 eV and 1.01 eV, respectively, indicating that these two products are difficult to release from the active sites and are thus kinetically disfavored in practice. In contrast, *OCH₂ can undergo successive hydrogenation to form CH₃OH and CH₄. The formation of these two C1 products is governed by the same PDS, namely the conversion of *CO to *OCH, with a corresponding ΔGPDS of 0.36 eV. Therefore, Pt@BNG mainly tends to produce CH₃OH and CH₄ with a limiting potential UL of -0.36 V. Similar reaction characteristics are observed for the W@BNG and Ta@BNG systems. In these systems, the formation of CH₃OH and CH₄ is limited by the *HCOO → *HCOOH and *OCH₂ → *CH₂OH steps, with corresponding limiting potentials of -0.55 V and -0.84 V, respectively. By contrast, the reaction pathway on Ru@BNG is slightly different, with the conversion of *CO to *OCH serving as the PDS (UL = -0.39 V). Under these conditions, the system can simultaneously generate multiple reduction products, including HCHO, CH₃OH, and CH₄.

As shown in Figure 5b, after CO₂ is adsorbed on the Rh@BNG surface to form the *CO₂ intermediate, the optimal reaction pathway for CO formation is *CO₂ → *COOH → *CO → CO, in which the conversion of *CO₂ to *COOH is the PDS, corresponding to a limiting potential of -0.17 V. By comparison, the formation of HCOOH, HCHO, CH₃OH, and CH₄ is all limited by the conversion of *HCOO to *HCOOH, with a corresponding limiting potential of -0.31 V. This result indicates that multiple C1 products on Rh@BNG share the same limiting step and may therefore be activated simultaneously during the reaction, leading to reduced product selectivity. In terms of desorption behavior, all C1 products except CO can desorb relatively easily from the catalyst surface, whereas CO exhibits a comparatively high desorption free energy (ΔGCO = 0.98 eV). Therefore, Rh@BNG tends to produce multiple products simultaneously, including COOH, HCHO, CH₃OH, and CH₄, rather than showing selectivity toward a single product. Similar reaction characteristics are also found for the Fe@BNG, Ni@BNG, and Ir@BNG catalysts. In all three systems, the PDS is *HCOO → *HCOOH, with corresponding limiting potentials UL of -0.43 V, -0.35 V, and -0.56 V, respectively, indicating poor product selectivity. For the V@BNG catalyst, the reaction is limited by the conversion of *CO₂ to *COOH (UL = -0.07 V), and at this potential the system can simultaneously generate multiple C1 products, such as CO, HCHO, CH₃OH, and CH₄. Taken together, the pathway analysis indicates that these catalysts, because of their broad product distributions, do not satisfy the screening criteria for highly selective catalysts.

Based on the ability of the SACs in Figure 6 to produce C1 products, the reaction behaviors of the 26 SACs were systematically classified into five categories in this work. Pd@BNG, Ti@BNG, Mn@BNG, Zr@BNG, Mo@BNG, and Hf@BNG mainly produce CO. Ag@BNG, Sc@BNG, Cr@BNG, Co@BNG, Cu@BNG, Zn@BNG, Y@BNG, Nb@BNG, and Cd@BNG show a stronger tendency to produce HCOOH. Os@BNG mainly produces CH₃OH, whereas Re@BNG mainly produces CH₄. In contrast, while Rh@BNG, Pt@BNG, V@BNG, Fe@BNG, Ni@BNG, Ru@BNG, Ta@BNG, W@BNG, and Ir@BNG can generate multiple C1 products simultaneously under similar reaction conditions and therefore exhibit multi-product characteristics. Previous studies have identified the Cu(211) surface as a high-performance metal catalyst for CO₂RR. Therefore, the limiting potentials of Cu(211) for different products were used in this work as the screening benchmarks, namely UL(CO) = -0.46 V, UL(HCOOH) = -0.32 V, and UL(CH4) = -0.61 V, while candidate catalysts were also required to exhibit selectivity toward a single product [22]. According to these criteria, 12 representative catalysts were screened from the 26 SACs, including Pd@BNG, Ti@BNG, Mn@BNG, Zr@BNG, Mo@BNG, Hf@BNG, Ag@BNG, Sc@BNG, Co@BNG, Cu@BNG, Y@BNG, and Cd@BNG, with corresponding UL values of -0.06 V, 0.00 V, -0.10 V, 0.00 V, -0.45 V, 0.00 V, -0.11 V, -0.22 V, -0.24 V, -0.22 V, -0.32 V, and -0.31 V, respectively. All of these catalysts outperform the Cu(211) benchmark in the production of a single C1 product. The remaining systems were excluded because of either excessively high limiting potentials or the simultaneous formation of multiple products at the same potential. The dominant product type, corresponding rate-controlling step, and limiting potential are summarized in Table 1.

3.3. Selectivity of CO2RR on TM@BNG

Because CO₂RR and HER have similar thermodynamic reaction potentials, and HER usually occurs concurrently with CO₂RR in aqueous electrolytes, the competing HER can lower the Faradaic efficiency and selectivity toward the target products [54]. Therefore, high-performance CO2RR catalysts should not only exhibit high activity but also effectively suppress HER to enhance selectivity toward the desired products [55]. To evaluate hydrogen evolution activity, this work calculated the Gibbs free-energy changes of HER on Pd@BNG, Ti@BNG, Mn@BNG, Zr@BNG, Mo@BNG, Hf@BNG, Ag@BNG, Sc@BNG, Co@BNG, Cu@BNG, Y@BNG, and Cd@BNG, and the corresponding results are shown in Figure 7a. The calculated limiting potentials for HER are -0.30 V, -0.54 V, -0.36 V, -1.26 V, 0.12 V, -0.71 V, -0.73 V, -2.27 V, -0.72 V, -0.91 V, -2.32 V, and -0.56 V, respectively. A further comparison between the limiting potentials of HER and CO₂RR (Figure 7b) shows that, except for Mo@BNG, the HER limiting potentials of the other 11 catalysts are all higher than their corresponding CO₂RR limiting potentials. This result indicates that Pd@BNG, Ti@BNG, Mn@BNG, Zr@BNG, Hf@BNG, Ag@BNG, Sc@BNG, Co@BNG, Cu@BNG, Y@BNG, and Cd@BNG can effectively suppress the competing HER during the reaction process and therefore exhibit high selectivity toward CO₂RR.

The chemical state of the electrode surface plays a crucial role in electrochemical reaction processes. Under specific potential and pH conditions, adsorbed species in solution, such as *O and *OH, may accumulate on the electrode surface, thereby changing the availability of active sites and consequently affecting the reaction rate and product selectivity [56]. To clarify the surface coverage behavior under CO₂RR conditions, surface Pourbaix diagrams were constructed for Pd@BNG, Ti@BNG, Mn@BNG, Zr@BNG, Hf@BNG, Ag@BNG, Sc@BNG, Co@BNG, Cu@BNG, Y@BNG, and Cd@BNG, and the corresponding results are presented in Figure 8. When the standard hydrogen electrode potential (USHE) exceeds the onset potential (Upristine, namely the critical potential for maintaining the pristine surface at pH = 0), the surfaces of these catalysts tend to be covered by *OH. However, as the potential increases further, all systems can revert to the bare-surface state, thus ensuring effective exposure of the single-atom active sites. Further analysis shows that the minimum potentials required to maintain an uncovered surface for Pd@BNG, Ag@BNG, Sc@BNG, Co@BNG, Cu@BNG, Y@BNG, and Cd@BNG are all lower than their corresponding limiting potentials for CO₂RR, indicating that surface oxidation of these catalysts can be effectively suppressed under actual reaction conditions. To further evaluate catalyst stability, ab initio molecular dynamics simulations were performed for these seven catalysts at 400 K for 10 ps. As shown in Figure 9, throughout the entire simulation, the temperature and total energy of each system fluctuate only within a narrow range, while the structures remain intact without obvious distortion. These results indicate that the catalysts possess good thermodynamic structural stability and high experimental feasibility.

As shown in Figure 10, the projected density of states (PDOS) of Pd@BNG, Ag@BNG, Sc@BNG, and Co@BNG after *HCOO adsorption is presented, with the dashed line denoting the Fermi level. A comparison of these systems reveals clear differences in the d-state intensity near the Fermi level, which usually governs the electron-coupling efficiency between the adsorbed species and the active site. In Pd@BNG, the Pd 4dxz orbital exhibits a sharp peak near 0 eV, while the C-2p orbital shows pronounced overlap within a similar energy range, indicating strong orbital hybridization between the Pd site and *HCOO. The high density of states near the Fermi level is also favorable for interfacial charge redistribution, thereby stabilizing key intermediates and facilitating subsequent charge transfer. Similarly, Co@BNG shows a significant contribution from Co-3d states near the Fermi level, and the corresponding peak profile varies consistently with that of C-2p over a neighboring energy range, indicating that the Co site can form stable electronic interactions with *HCOO and thus promote intermediate activation. By contrast, in Sc@BNG, the Sc-3d peak is located close to the Fermi level but is more broadly distributed, and the coupling signal extends over a wider energy range. This feature suggests that the regulation of *HCOO by the Sc site is characterized by more continuous and moderate electronic participation. Ag@BNG is more representative of a moderately adsorbed state. Its 4d orbitals are mainly distributed at lower energy levels, and the overlap with C-2p near the Fermi level is less concentrated, although a clear response to the adsorbate orbitals is still observed. This result indicates relatively mild participation of near-Fermi electrons, which is favorable for intermediate desorption after conversion and reduces the risk of long-term occupation of active sites. Overall, Pd@BNG is characterized by strong orbital coupling and high electronic accessibility, Ag@BNG by moderate adsorption strength and favorable desorption behavior, Sc@BNG by more continuous electronic participation, and Co@BNG by more pronounced adsorption and activation capability. These four catalysts highlight different aspects of electronic-structure characteristics and adsorption-regulation behavior, and all can therefore be regarded as reasonable candidates for effective active sites. In addition, the PDOS results of Cd@BNG, Cu@BNG, and Y@BNG are provided in the Appendix for horizontal comparison and supplementary analysis.

3.4. Machine Learning

When constructing a machine-learning model to correlate catalytic performance with the intrinsic properties of materials, it is first necessary to determine whether redundancy exists among the input variables, because linear dependence or strong correlation between features may lead to multicollinearity and thus compromise model stability [57,58]. On this basis, pairwise correlation analysis was performed for all candidate parameters, and the correlation coefficient R was used to quantify the linear dependence between variables, with the results visualized in a correlation matrix [37]. As shown in Figure 11, the overall correlation distribution is relatively scattered, and the relationships among the features are clearly distinguishable. From a physical perspective, the nine input variables can be classified into three categories: the d-band center (εd), p-band center (εp), and Fermi level (Efermi), which describe band-position information; the first ionization energy (Ie), electronegativity (χ), and electron affinity (Ea), which reflect electron-donating and electron-accepting ability; and the magnetic moment (Mag), covalent radius (r), and atomic charge (Q), which represent local properties. The statistical results show that the absolute values of the correlation coefficients for all feature combinations remain outside the high-correlation range, indicating the absence of significant coupling among the variables. Such a data structure is beneficial for reducing the influence of multicollinearity on the model. Therefore, the nine parameters were collectively used as model inputs for subsequent performance prediction and trend analysis.

This chapter also employs the GBR method to train and predict the sample data [38,39,40,41,42]. As shown in Figure 12a, the model exhibits good fitting performance on the training set, with predicted values generally consistent with the corresponding actual values and only slight dispersion. The coefficient of determination (R₂) for the training set reaches 0.97, while the root mean square error (RMSE) is 0.04, indicating low prediction error and confirming that the model provides sufficient accuracy and stability for the present study. On this basis, the contributions of the first nine input variables were further evaluated (F Figure 12). The results show that different features contribute to UL to different extents, among which Q, Ie, and Ea are particularly important, accounting for approximately 22%, 20%, and 18%, respectively. These three features therefore play dominant roles in the model decision process. Combined with the feature-ranking results, it is evident that Q, Ie, and Ea are closely associated with UL and may jointly affect the reaction-limiting step at the catalytic site through synergistic effects, thus providing a basis for subsequent mechanistic analysis.

To further reveal the relationship between the intrinsic parameters of single-atom catalysts and UL in mathematical terms, the SISSO method was introduced to identify an appropriate descriptor [43]. Through sure independence screening and sparsifying operator construction, this method selected, from candidate expressions associated with Ie and Ea, the descriptor that best characterizes UL. The resulting descriptor is denoted as φ, with the expression φ=sin(abs((εd/Ea)-(Ie*x))). As shown in Figure 13, φ exhibits a volcano-type relationship with UL, indicating the presence of an optimal range associated with relatively high catalytic activity. Specifically, for the active sites in TM@BNG, when φ falls within the approximate range of 0.35-0.66, the predicted UL is generally higher than that of Cu(211), suggesting more favorable CO₂RR performance at these sites. In addition, this descriptor can be used to quantify the influence of the metal center and its coordination structure on catalytic activity, thereby providing guidance for the rational design of single-atom catalysts and the efficient screening of active sites in future studies.

4. Conclusions

In our study, the stability and catalytic performance of 27 TM@BNG single-atom systems toward CO₂RR were systematically investigated based on DFT calculations. Analyses of formation energies and dissolution potentials indicate that most candidate structures possess good thermodynamic and electrochemical stability under reaction conditions. On this basis, combined with analyses of reaction pathways and product selectivity, seven active sites with superior performance were further screened from the 27 candidate systems, namely Pd, Ag, Sc, Cu, Co, Cd, and Y. These sites exhibit high selectivity toward the formation of HCOOH or CO, with corresponding limiting potentials of -0.06 V, -0.11 V, -0.22 V, -0.22 V, -0.24 V, -0.31 V, and -0.32 V, respectively. Further machine-learning analysis shows that catalytic activity is closely correlated with intrinsic parameters such as atomic charge, ionization energy, and electron affinity. Grounded in these findings, a descriptor containing key parameters was constructed, and a volcano-type relationship between this descriptor and UL was established. These results not only clarify the distribution pattern of favorable active sites in the TM@BNG system, but also provide a theoretical basis for the screening and performance regulation of single-atom catalysts for CO₂RR.