Design of Graphdiyne and Holey Graphyne‐Based Single Atom Catalysts for CO2 Reduction With Interpretable Machine Learning

Using electrochemical CO2 reduction reaction (CO2RR) to synthesize value‐added hydrocarbons provides a useful solution for environmental issues and energy crisis. However, this process is impeded by the low activity and selectivity of electrocatalysts toward targeted products. Employing density functional theory computations, the graphdiyne and holey graphyne supported single‐atom catalysts (SACs, M/GDY and M/HGY) are demonstrated to be the promising candidates for the CO2RR. By taking full elemental diversity of metal sites across the periodic table, 25 catalysts are found to effectively activate CO2 and inhibit competitive hydrogen evolution, and 8 of them show higher activity for CH4 production than Cu(211). Remarkably, changing supports are found to greatly affect limiting potentials and reaction pathways, even leading to an “inert‐active” transition for some metal centers. The resulting SACs, including Mn/GDY, Co/HGY, Ru/GDY, and Os/GDY, can achieve high activity with low limiting potentials of ≈ −0.22 to −0.58 V. Machine learning enables to identify the critical role of the polarized charge and magnetic moment of metal atoms in affecting the activity. The built machine learning model also shows an interpretable capability to predict the activity of the other types of SACs, offering a great promise to quick screening of high‐performance SACs.


Introduction
The rapid growth of energy demand has caused the sharp consumption of fossil fuels and high emission of carbon dioxide (CO 2 ), resulting in the damage of the natural normal carbon cycle and the generation of a series of environmental problems. [1]o address the above issues and achieve the low-carbon economy and socially sustainable development, a lot of carbon capture and conversion strategies have been proposed.Among them, the electrochemical reduction of CO 2 into high value-added hydrocarbons has become one of the most promising solutions due to environmental friendliness, low cost, and high efficiency. [2]lthough electrochemical CO 2 reduction reaction (CO 2 RR) has a broad prospect, there are still some challenges that urgently need to be resolved.On the one hand, the CO 2 molecule with a linear structure has a strong molecular inertness.High energy input is required to activate/break the C═O bond.In addition, the electrochemical CO 2 RR is a complex process that involves multiple proton-electron pair transfer steps.Diverse products and very similar thermodynamic redox potentials of each reaction pathway make a poor catalytic selectively toward the target products, such as CO, HCOOH, CH 4 , and CH 3 OH. [3]On the other hand, hydrogen evolution reactions (HER) can easily occur as a competitive reaction in aqueous solutions, which decreases the Faradaic efficient of CO 2 RR and sometimes becomes the main barrier of reaction. [4]u-based materials showed a good performance for the CO 2 RR.Mild binding strength for the reaction intermediates enabled these catalysts to produce various products, while the separation and purification of these products were quite difficult. [5]herefore, the development of cost-effective electrocatalysts with superior activity and selectivity is of great significance for the efficient CO 2 RR.Recently, single-atom catalysts (SACs) with maximum atom efficiency and unique electronic properties have become one of the most attractive frontiers and are widely used as electrodes for the nitrogen reduction reaction (NRR), [6] oxygen reduction reaction (ORR), [7] HER, [8] and CO 2 RR. [9]Notably, the aggregation of single atoms on the catalytic support is a thorny problem for the preparation of SACs, which could reduce the catalytic performance of the catalyst.Appropriate supports are required to strongly bond with single atoms to prevent metal aggregation.Typical supports include chalcogenides, oxides, carbon materials, and so on. [10]Among them, 2D carbon-based materials are of great promise in terms of their good stability under the reaction conditions. [11]pecifically, Graphdiyne (GDY), a new type of 2D carbon material, has been successfully synthesized on the copper surface. [12]he unique structure and excellent electronic characteristics have made GDY gain increasing attention in recent years. [13]The sp hybridization of −C≡C− in GDY enables /* to vertically point toward the central metal, to potentially form a strong interaction with metal atoms and maximum prohibit their aggregation.Experimentally, single Cr, [14] Mn, [14] Fe, [15] Ni, [15a] Cu, [16] Mo, [14,17] Ru, [18] Pd, [19] W, [14] Re, [14] and Pt [20] atoms anchored on GDY (M/GDY) are successfully realized and used toward OER, HER, ORR, and NRR.Theoretically, Sc, Ti, and Ir supported by GDY have been reported to have excellent thermal stability and outstanding catalytic performance for CO oxidation. [21]Unfortunately, studies related to the M/GDY toward the CO 2 RR are still limiting, and most reports are limited to 3d transition metals. [22]The great potential of transition metal beyond 3d needs to be investigated.Besides GDY, a novel 2D carbon allotrope, namely holey graphyne (HGY), was also synthesized using Castro-Stephens coupling reaction. [23]HGY is alternately linked by the six-membered ring and eight-membered ring, which are composed of sp 2 and sp hybridized C atoms.The special -conjugated structure and nonlinear sp bonding make HGY also a potential support to form metal anchored SACs (M/HGY).However, HGY-based SACs have not been realized in the experiment and the CO 2 RR mechanism is unclear.To this end, comprehensive understanding the potential role of GDY and HGY in the design of SACs for the CO 2 RR is quite important to optimize the catalysts.
In this work, we explored the full elemental diversity of SACs by anchoring a series of transition metal atoms on the experimentally available GDY and HGY.Among 60 SAC models, 25 catalysts were demonstrated that could be served as efficient CO 2 RR catalysts due to effective CO 2 activation ability and high selectivity toward CO 2 RR.Eight catalysts could break the metal-based benchmark to achieve higher catalytic activity.Remarkably, the catalytic activity, product selectivity, as well as poisoning issue of metal sites could be well addressed through a rational selection of catalytic supports.Statistical analyses including machine learning enable to identify the key properties in affecting the overall performance and build the predictable framework to screen a large composition space of SACs.The improved activity of CO 2 RR was found to come from the tunable interaction between the reaction intermediates such as *OH and metal in SACs, wherein the central metal atom with smaller polarized charge is more likely to induce the moderate activation of CO 2 molecules to produce high CO 2 RR catalytic activity.These results contribute to a comprehensive understanding of the role of the metal elements and supports in CO 2 RR, which provides a guidance to enhance the CO 2 RR performance.

Computational Details
All calculations were performed by the spin-polarized DFT method that was implemented in the Vienna ab initio simulation package (VASP) with the projector-augmented-wave (PAW) potential. [24]The Perdew-Burke-Ernzerhof (PBE) functional within the generalized gradient approximation (GGA) [25] was used to model the exchange-correlation energy.An energy cutoff of 450 eV was adopted for the plane-wave basis-set, and the Brillouin zones were sampled by a Monkhorst-Pack -point grid of 3 × 3 × 1 for all geometry relaxation, while a denser mesh of 6 × 6 × 1 was chosen for the density of states (DOS) calculations.To fully relax the structure, thresholds of 10 −5 eV and 0.02 eV Å −1 were set for the energy and force convergence, respectively.The dispersion correction was carried out via the DFT+D3 method within Grimme's scheme to properly describe the van der Waals (vdW) interaction. [26]he catalytic supports were modeled by 2 × 2 GDY and HGY monolayer supercells, of which a vacuum layer of 20 Å was added to avoid interactions between periodic images.The thermal stability of metal atom anchored GDY and HGY was evaluated by ab initio molecular dynamics (AIMD) simulation under the NVT ensemble at 400 K for 10 ps.The binding energies (E b ) for the SACs were calculated as , where E SUB and E M/SUB are the total energies of the substrate and single metal atom anchored substrate.E M presents the total energy of an isolated metal atom.The cohesive energy (E c ) of metal bulk was defined as E c = E bulk ∕n − E M , where E bulk and E M are the total energies of metal bulk and the isolated metal atom in vacuum and n is the number of metal atom in the bulk.The relative Gibbs free energy change along the reaction pathways was computed on the basis of computational hydrogen electrode (CHE) model which was proposed by Nørskov et al. [27] More calculation details related to electrochemical reactions are shown in the Supporting Information (SI).

Structural Models and Stability of M/GDY and M/HGY
GDY and HGY are the two types of 2D materials that contain rich natural and uniform holes (see Figure 1a).The unique structural characteristic makes them promising substrates to capture metal atoms and form atomically dispersed single atom catalysts (i.e., M/GDY and M/HGY).In principle, various 3d, 4d, and 5d transition metal and the main group metal elements can act as the anchoring atoms in the M/GDY and M/HGY.Considering the toxic or radioactivity properties, 30 metal atoms, including 26 transition metal atoms (except Tc, Cd, and Hg) and four main group elements (Al, Ga, Sn, and Bi) were chosen and anchored on GDY and HGY monolayers (see Figure S1, Supporting Information), and a total of 60 SAC models were built.
Figure 1a illustrates the possible anchoring sites in the GDY and HGY.For GDY, there were three possible anchoring sites: 1) the acetylenic-ring center (namely, S1), 2) the acetylenic-ring corner (S2), and 3) the hexatomic-ring center (S3).For HGY, four different hole sites were taken into consideration, including: 1) the center of the large pore formed by the six acetylenic linkages (S1), 2) the corner of the large pore formed by the six acetylenic linkages (S2), 3) the hexatomic-ring center (S3), and 4) the center of eight-membered ring containing two acetylene (S4).After full relaxation, the most stable geometric configurations of M/GDY and M/HGY were determined, which are displayed in Figures S2  and S3, Supporting Information.
For GDY, most of the metal atoms can be anchored stably at the acetylenic-ring corner (S2), where Ti atom is bonded to five C atoms; Sc and Y are bonded to six C atoms; Au, Al, and Sn are bonded to two C atoms; and the rest are bonded to four C atoms.Particularly, Zr and Hf tend to be anchored at the side of the acetylenic ring with four C atoms, while Zn and Ga are in the center of the acetylenic ring (S1).Similar to GDY, a majority of the metal atoms can be anchored stably at the corner of the large pore of HGY (S2) by bonding two (M = Cu, Pd, Ag, Au, Al, Ga), three (M = V, Co, Ni), or four (M = Mn, Fe, Zr, Nb, Mo, Ru, Rh, Ta, W, Re, Os, Ir, Pt) C atoms; whereas for Ti, Cr, Y, Hf, Sn, and Bi, they are more easy to bond with the two acetylenic C. In addition, the favorable anchor site for Sc is S4 (Figure 1a) and Zn is located in the center of the large pore of HGY (S1).The calculated average M─C bond lengths in most M-GDY (≈1.90-2.48Å) and M-HGY (≈1.94-2.54Å) are less than 2.60 Å, and the corresponding shortest M─C bond lengths are ≈1.85-2.41Å and ≈1.82-2.54Å, respectively.However, for Zn/GDY, Ga/GDY, and Zn/HGY, the shortest distances between the metal atom and the adjacent C atom are 3.68, 2.96, and 4.15 Å. Partial structural parameters are listed in Table S1, Supporting Information.
Note that the effective capture of CO 2 on the electrode surface is the prerequisite for catalytic reactions.Among 60 SAC models, 32 catalysts, including 17 M/GDY (M = Ti, V, Cr, Mn, Fe, Y, Zr, Nb, Mo, Ru, Rh, Hf, Ta, W, Re, Os, Ir) and 15 M/HGY (M = Ti, V, Cr, Fe, Co, Nb, Mo, Ru, Rh, Hf, Ta, W, Re, Os, Ir) display strong capability of adsorbing CO 2 (usually accompanied by the deformation of the CO 2 molecule), and thus; are selected for the further stability analysis.
The metal atoms might agglomerate into clusters in some of the SACs due to the insufficient binding between the metal and substrate, which is harmful to the real application of catalysts.Therefore, we calculated the binding energy (E b ) of single metal atoms on GDY/HGY monolayer to evaluate their stability (see Figure 1b).For comparison, the cohesive energy (E c ) of the metal atoms in their bulk phase was also computed.According to the definition of E b , a more negative value means stronger binding between metal atoms and supports, which can effectively prevent the agglomeration of metal atoms and improve the stability of catalytic materials.Promisingly, the E b values of all single metal atoms anchored on GDY/HGY monolayer are negative than −2.00 eV (≈ −2.20 to −5.88 eV), indicating the strong binding between metal atoms and GDY/HGY support.7 M/GDY (M = Ti, V, Fe, Co, Zr, Hf and Os) show more negative E b values relative to M/HGY; whereas, for the other 11 metal atoms (Cr, Mn, Y, Nb, Mo, Ru, Rh, Ta, W, Re, and Ir), HGY displays a stronger binding than the GDY.Similar trend can be observed in the energy difference between E b and E c (see Figure 1b).
However, different from studied graphene supported SACs such as the M-N x -C system, [6a] the energy differences between E b and E c of M/GDY and M/HGY are positive and within 0.49 and 3.96 eV (Table S1, Supporting Information) as the negative value of E b − E c generally suggests the thermodynamic feasibility of formation of SACs from their bulk phase.The computed E b − E c values for the M/GDY and M/HGY indicate a difficulty of using the thermal decomposition of metal precursors to synthesize these two types of materials. [28]Other fabrication strategies, such as control of electrochemical potential window of metal atoms via electrochemical deposition or in situ coordination approach, might be carried out to obtain the high loading catalysts. [29]Indeed, Zou et al. have recently successfully synthesized five M/GDY (M = Cr, Mo, W, Mn, and Re) materials via the in situ coordination approach.Among them, the mass loading of W on GDY is ≈1.03 wt%, which is much higher than that through a heterogeneous reaction between GDY and [W(CO) 6 ] (0.15 wt%). [14]mong all the studied 32 SACs, the computed E b − E c values of M/GDY and M/HGY are lower than the experimentally available W/GDY (3.96 eV), which implies that these catalysts could also be synthesized under appropriate conditions.In particular, though Cr atoms exhibited the weakest bonding with the substrates (Cr/GDY [ E b = − 2.20 eV] and Cr/HGY [E b = − 2.33 eV]) compared with other metal atoms (see Figure 1b), AIMD simulations also demonstrate superior thermal stability (see Figure 1c).The Cr atoms in the Cr-GDY and Cr-HGY keep stable at their favorable anchoring sites, and the average Cr─C bond length shows a negligible fluctuation of 0.01 Å during 10 ps timesteps at a temperature of 400 K.

Adsorption and Activation of CO 2 on M/GDY and M/HGY
Following the above analyses, the detailed adsorption behavior of CO 2 on the 32 SAC surfaces was then explored, and the corresponding adsorption energies of most stable adsorption geometries are shown in Figure 2a.From Figure 2d; Figures S4 and  S5, Supporting Information, it can be seen that CO 2 can be adsorbed steadily on the surface of most catalysts by forming a closed ternary ring structure of M-O-C, of which the average C─O bond length and O─C─O angle in CO 2 are in the range of ≈1.21-1.30Å and ≈127.9-154.7 ○ , respectively.Differently, the CO 2 molecule can be adsorbed stably on Y/GDY through forming a Y─O─C─O four-membered ring, and the average C─O bond length and the O─C─O angle of CO 2 are 1.27 Å and 127.0 ○ , respectively.
Analysis of the charge density difference and Bader charge indicates that a large number of charges were transferred to the adsorbed CO 2 molecules by the M─C and M─O bonds (see Figure 2d; Figure S4 and S5, Supporting Information).Interestingly, the charge transferring to CO 2 molecules and the average C─O bond length or O─C─O angle of adsorbed CO 2 molecules all have good linear scaling relationships with R 2 of 0.86 and 0.88 (Figure 2b,c).With the increase of charge transfer, the average C─O bond length gradually elongates and O─C─O angle gradually decreases.This result indicates that the charge transfer between the SACs and CO 2 can effectively promote the activation of CO 2 , resulting in the elongated C─O bond length and the smaller O─C─O angle.
In order to further understand the activation mechanism of the CO 2 on the SAC surfaces, the partial density of states (PDOSs) and crystal orbital Hamilton populations (COHP) [30] were then computed, as shown in Figure 2e; Figures S6-S8, Supporting Information.A clear hybridization between the d-orbital of metal atom and the molecular orbitals of CO 2 was found.Specifically, COHP analysis indicates the  * u molecular orbital of CO 2 interacts with d-orbital to form the partial occupied states.Moreover, the integrated crystal orbital Hamilton population (ICOHP) is in the range of ≈ −0.44 to −2.07, which also confirms the strong interaction between CO 2 and SACs.
Taking Ti/GDY as an example, we further analyzed the dorbitals of Ti before and after CO 2 adsorption, as shown in Figure 3.In pristine Ti/GDY, Bader charge analysis indicates that Ti atom is polarized, holding 2.80 electrons after transferring 1.20 electrons to GDY support.The calculated PDOS illustrates that the five degenerate 3d orbitals are divided into different energy states, three of which, that is, d xy , d x 2 −y 2 , and d z 2 , have lower energy states occupied by 2.80 electrons, and the other two orbitals (d yz , d xz ) are almost empty and most of the density of states are located above the Fermi level (Figure 3a).Such a characteristic offers the metal atoms a great potential to interact with both the occupied and unoccupied orbital of CO 2 .
Indeed, after CO 2 adsorption, the adsorbed CO 2 molecule can be activated via "donation-backdonation" mechanism (Figure 2f).PDOS analysis indicates that the CO 2  g orbital donates electrons to Ti d yz empty orbital, forming the bonding states at −3.00 eV to enhance adsorption.On the other hand, the d yz and d z 2 orbitals back-donate electrons to  * u unoccupied orbitals of CO 2 , leading to the partial occupied orbital near the Fermi level.Such a donation-backdonation process can effectively activate CO 2 , which results in an elongation of the C─O bond length and a bent configuration of CO 2 molecule.Remarkably, Ti d yz orbital is found to participate in both charge donation and backdonation processes; and thus, acts as an "electron transfer station" to transfer electrons between metal atom and the CO 2 (see Figure 3c).As a result, the significant hybridization and charge transfer between molecular orbitals of CO 2 and metal d orbital can strengthen the interaction between CO 2 molecules and M/GDY and M/HGY and benefit to the subsequent hydrogenation steps.

Selectivity for CO 2 RR Versus HER
The electrocatalytic CO 2 RR is a process that constantly accompanies the proton-electron coupling.HER is considered as a competitive reaction due to the participation of proton-electron pair from the electrolyte solution, which seriously affects the Faradaic efficiency of CO 2 RR.To evaluate the catalytic selectivity of CO 2 RR, we first compared the adsorption energies of CO 2 and H atom on M/GDY and M/HGY, as given in Figure 2a; Table S2, Supporting Information.24 catalysts (14 M/GDY, M = Ti, V, Cr, Mn, Fe, Y, Zr, Nb, Mo, Ru, Hf, Ta, W, Os and 10 M/HGY, M = Ti, V, Cr, Fe, Co, Nb, Mo, Hf, Ta, and W) were found to preferentially adsorb CO 2 compared to H atom with the more negative adsorption energy for CO 2 , which makes the CO 2 RR more likely to occur.Re/GDY has the closest adsorption energy values for CO 2 (−1.10 eV) and H atom (−1.22 eV), which is also expected to be a CO 2 RR candidate.
Notably, with the increase of external electrode potential, the hydrogen adsorption process will be accelerated due to the participation of electrons, leading to the decreased coverage of *CO 2 .Efficiently converting CO 2 is also the key to suppress the hydrogen adsorption.Thus, we then evaluated the Gibbs free energy changes (ΔG) of the initial hydro-

Reaction Mechanism and Activity for CO 2 RR
Fundamentally, various carbon products such as C1 and C2 products can be obtained through the CO 2 RR.Due to the lack of suitable C─C coupling sites, the main products on the SACs are C1 products. [31]So, we focused only on the C1 pathways, and the possible elementary reaction steps for CO 2 RR to C1 products are summarized in Figure 5a.The limiting potential (U L ) was employed as an activity indicator to evaluate the activity of catalysts, which is obtained by the U L = − ΔG max /e, where ΔG max is the largest free energy change in a specific reaction pathway.
Figure 5b; Table S3, Supporting Information summarize the limiting potentials (U L ), potential determining steps (PDS), and the final products along the most favorable reaction pathway for CO 2 RR on the 15 M/GDY and 10 M/HGY, which are screened from catalytic selectivity analysis.Remarkably, all these catalysts show strong capability to produce CH 4 products, eight of which (V/GDY, Cr/GDY, Mn/GDY, Fe/GDY, Ru/GDY, Os/GDY, Fe/HGY, and Co/HGY) show excellent CO 2 RR activity compared to the metal-based benchmark Cu(211) ( U L = − 0.74 V), [32] with less negative limiting potential range from −0.13 to −0.72 V.
Interestingly, for the catalysts having same metal center but different substrates, the catalytic performance can be varying.Simply, the effects of supports on CO 2 RR activity can be roughly divided into three categories: 1) Changing the reaction energy barriers without affecting the reaction pathways, such as Ti, V, Nb, Mo, Hf, Ta, and W SACs. 2) Changing reaction pathway and catalytic activity, such as Cr and Fe SACs.3) Re-activation of the stuck CO 2 RR by effective adsorption and activation of CO 2 , such as Mn, Co, Y, Zr, Ru, Re, and Os SACs.In the following sections, the different role of substrates will be discussed separately.

Tuning Limiting Potential
Reaction mechanism analysis indicates that the theoretical limiting potential can be tuned by the substrates on the Ti, V, Nb, Mo, Hf, Ta, and W SACs. Taking V SACs as an example, the corresponding Gibbs free energy profiles of CO 2 RR on the V/GDY and V/HGY are shown in Figure 5c,d, respectively.For both materials, the H atom prefers to couple with the C atom of adsorbed CO 2 to form * OCHO reaction intermediate with the more negative Gibbs free energy changes (−1.05 and −1.18 eV for V/GDY and V/HGY) compared to * COOH (−0.21 and −0.12 eV).Next, H atom continues to attack the C site to get * H 2 COO intermediate with the change in free energy of 0.44 and −0.18 eV for V/GDY and V/HGY.In contrast, the generation of * HCOOH intermediate needs a higher energy input of 0.77 and 0.94 eV (see Figure 5c,d Similarly, the U L values for M/GDY and M/HGY (M = Ti, Nb, Mo, Hf, Ta, and W) can be obtained, as presented in Figure S9, Supporting Information.The potential difference of U L can reach to 0.34 V, demonstrating great impact of substrate on the catalytic activity.Among them, the U L value of V/GDY is lower than that of Cu(211). [32]These results indicate that the appropriate support helps to improve the CO 2 RR activity and even break the metalbased benchmark.

Tuning Reaction Pathway
CO 2 RR involves multiple reaction pathways and reaction species; effectively balancing adsorption strength of each reaction intermediate is the key to achieve high activity.Changing the reaction pathway is a potential method to achieve this goal for some catalysts, such as Cr and Fe SACs.Taking Fe SACs as an ex-ample, the detailed free energy profiles of CO 2 RR are shown in Figure 5e,f.For Fe/GDY, CO 2 readily protonates to * OCHO, and the formed * OCHO intermediate prefers to generate * HCOOH by the hydrogenation of the O site with small increase of free energy (0.03 eV).Compared to desorption as HCOOH product (0.3 eV for desorption energy), * HCOOH is easier to be further reduced to * CHO with the free energy change of 0. As for Fe/HGY, the most favorable reaction pathway is * OCHO → * HCOOH is the PDS with the U L value of 0.5 V.By comparison with Fe/GDY, the reaction pathway is changed from * OCHO → * HCOOH → * CHO to * OCHO → * HCOOH → * H 2 COOH, and the maximum free energy is increased by 0.37 eV.The role of GDY substrate is to help to weaken the adsorption of * OCHO, * HCOOH and * H 2 COOH intermediates so that the adsorption of * OCHO, * HCOOH and * CHO reaches an effective balance which minimizes the potential barrier of each step and improves the overall activity (see Figure 5e,f).Similarly, the limiting potential of Cr SACs changes to 0.08 V and the pathway is changed from * OCHO → * H 2 COO on Cr/HGY to * OCHO → * HCOOH on Cr/GDY (see Figure S10, Supporting Information).These results show that the substrate can help to selectively tune the adsorption of some reaction species, offering another freedom to optimize the reaction pathway and increase the activity.

Reactivating Stuck CO 2 RR
Effective activation of the CO 2 molecule is a prerequisite for the CO 2 RR.For SACs (M = Mn, Co, Y, Zr, Ru, Re, and Os), changing the support is found to significantly affect the adsorption and activation of CO 2 on metal centers.For example, the chemical adsorption of CO 2 on Co/GDY is almost prohibited, which will cause CO 2 RR to be stuck; while the CO 2 molecule can be effectively captured and activated on Co/HGY.As shown in Figure S11a, Supporting Information, the adsorbed CO 2 is preferentially reduced to * OCHO on the Co/HGY rather than to * COOH intermediate due to downhill free energy (−0.35 eV).Then, H atom attacks O atom in * OCHO to form * HCOOH with a small increase of free energy (0.09 eV).Further, the larger desorption energy (0.62 eV) makes the desorption of * HCOOH to be difficult.Instead, the * HCOOH further reduces to * CHO, that is, PDS; the corresponding U L is 0.29 V, which is lower than that of most reported catalysts.−0.22, −1.33, −1.46, −0.49, −1.74, and −0.58 V (the corresponding free energy profiles are given in Figure S11b-g, Supporting Information).However, these metals are difficult to capture CO 2 molecules based on the HGY support.Therefore, to activate these metals for the CO 2 RR, rational design of catalytic support is of great importance.

Effect of Solvation on Catalytic Activity
Considering that the actual catalytic reaction is performed in aqueous solution, the solvation effect was also calculated by the implicit solvation model with the relative permittivity of 78.4 implemented in VASPsol. [33]Take M/GDY (M = V, Cr, Mn, Fe, Ru, Os) and M/HGY (M = Fe and Co) as an example, the detailed Gibbs free energy profiles of CO 2 RR along the most favorable reaction pathway are displayed in Figure S12, Supporting Information.The results indicate that the final products of these catalysts (except Mn) were still CH 4 and the U L values of M/GDY (M = V Cr, Fe, Ru, Os) and M/HGY (M = Fe and Co) are −0.74,−0.53, −0.28, −0.36, −0.57, −0.50 and −0.31 V, respectively.For Mn/GDY, HCOOH is the final product with the U L of −0.56 V.Although Mn/GDY has changed slightly in CO 2 RR after considering solvation effect, the U L of these eight catalysts is still lower than that of Cu( 211) and the overall trend has not changed significantly.

The CO 2 RR Activity Origin for M/GDY and M/HGY
To understand the activity trend of the M/GDY and M/HGY, it is important to find a descriptor to identify and describe the activity trend.To achieve that, we explored possible correlations between the U L and the adsorption energy of each reaction intermediates.6a-h), the adsorption strength of * OH intermediate was found to be able to govern the CO 2 RR activity of these catalysts with the R 2 value of 0.90.As shown in Figure 6i, the U L gradually decreases with the increase of the ΔE * OH , and the lowest U L value is achieved on the Fe/GDY with the less negative ΔE * OH of −2.81 eV, corresponding to the highest CO 2 RR activity.Note that though * COOH and * OH interact with metal centers via different atoms (i.e., carbon and oxygen), they displayed a good adsorption relationship mostly because of the strong correlation between the adsorption strengths of * C and * O on most transition metals. [34]his could be the reason why using a single descriptor of ΔE * OH can generally understand the catalytic activity of different SACs.
To reveal the internal origin of SACs on the CO 2 RR activity, the general correlation between ΔE * OH and intrinsic properties of the SACs was investigated by machine learning (ML) model. [35]In the model, a feature set with 12 intrinsic descriptors of these catalysts is selected, which consists of five structural descriptors (M atomic number (Z), M atomic radius (r), the shortest M─C bond length (d min ), average M─C bond length (d M−C ), and ICOHP between M and GDY/HGY), seven electronic descriptors (d-band center ( d ), the valence electrons (N a ), polarized charges (Q), electronegativity (EN), electron affinity (EA), first ionization energy (I 1 ), and magnetic moment (Mag) of the central metal atom in M/GDY and M/HGY).The specific data is listed in Table S4, Supporting Information.We first analyzed the Pearson correlation among these 12 descriptors to find their internal correlations and reduce the number of random variables.The calculated Pearson correlation heat map is shown in Figure S13, Supporting Information, where the feature pairs with the Pearson coefficient greater than 0.6 or less than −0.6 are regards to have strong correlation.Abandoning the strong correlation features, we finally determine a non-redundant feature set consisting of ICOHP, N a , Q, EA, and Mag for the further ML modeling.
To adopt a suitable algorithm for machine modeling, we explored the relationship between ΔE * OH and five relatively independent features (ICHOP, Na, Q, EA, Mag), as shown in Figure S14, Supporting Information.Notably, although a series of linear relationships between ΔE * OH and ΔE of some other reaction intermediates were found, there were no clear linear correlations between ΔE * OH and these independent features as indicated by the very small R 2 values (0.37, 0.05, 0.57, −0.04, and 0.26).Therefore, we used the random forest method for further modeling.7a,36] The random forest algorithm based on Scikit-learn was used to optimize the model.The dataset was randomly divided into two parts, including 20% test samples and 80% training samples (more details can be found in Tables S5 and S6, Supporting Information).Based on our previous work, the number of trees and the max depth were set at 200 and 5, respectively.In addition, random state = 9 was set.The random forest regression model was established after data standardization.Promisingly, the trained model exhibited satisfactory performance between predicted ΔE * OH and DFT calculated ΔE * OH (Figure 7a), where the train and test scores (R 2 ) were 0.94 and 0.97, respectively.The calculated rooted mean square error (RMSE) and mean absolute error (MAE) are also given in Table S7, Supporting Information.In order to reveal the internal influencing factors, we further analyzed the feature importance in the model.As shown in Figure 7b, ΔE * OH is mainly related to ICOHP, Q and Mag, of which polarized charges (Q) of the central M atom in M/GDY and M/HGY with the feature importance of 50.5% play a decisive role.Promisingly, we extended the ML model to the widely studied graphene-based SACs, that is, single iron atom coordinated with different number of carbon and boron (the corresponding feature values are given in Table S8, Supporting Information).The model still displays excellent performance with the high train and test scores (R 2 ) of 0.94 and 0.95 (Figure S15a and Table S7, Supporting Information).The feature importance in the model maintains a similar trend (Figure S15b, Supporting Information), in which Q (52.4%) dominates, followed by ICOHP (26.8%) and Mag (16.1%).These results imply a robust capability and universality of ML framework to design other types of SACs.
To clear the relationship between Q and CO 2 RR activity, we drew the violin plot through the probability density obtained by the statistical analysis between Q and ΔE * OH (Figure 7c).Obviously, central metals with smaller Q are more likely to induce high CO 2 RR activity of catalysts, where ΔE * OH ≥ −4.04 eV.In addition, charge transfer from M/GDY and M/HGY to CO 2 and polarized charges (Q) of the central M atom in M/GDY and M/HGY are also associated to further reveal the relationship between CO 2 activation and catalytic activity of catalysts.Figure 7d shows that there is an approximately linear dependence, where the larger Q of the central metal corresponds to the larger charge transfer to CO 2 , indicating that the empty d orbit is favorable to the capture and activation of CO 2 molecules.However, larger Q of the central metal is easy to cause excessive activation of CO 2 molecules, which makes the binding between the central metal and O atom excessively strong to produce more negative ΔE * OH (see Figure 6i), resulting in higher U L .These results indicate that changing the support can modulate the polarized electrons of the central metal d orbit, and the central metal atom with a smaller polarized charge is more likely to induce the moderate activation of CO 2 molecules to produce high CO 2 RR catalytic activity.

Conclusion
In summary, we comprehensively explored the CO 2 RR performance of GDY and HGY supported SACs using the DFT computations and unveiled their underlying catalytic origin in com-bination with the machine learning method.25 catalysts were identified as efficient CO 2 RR catalyst candidates due to effective activation of CO 2 molecules and significant inhibiting competitive HER reaction, of which eight catalysts, including V/GDY, Cr/GDY, Mn/GDY, Fe/GDY, Ru/GDY, Os/GDY, Fe/HGY, and Co/HGY, had high CO 2 RR activity for the CH 4 production with low limiting potentials (−0.67, −0.72, −0.22, −0.13, −0.49, −0.58, −0.50, and −0.29 V).The activity trends of these SACs showed a strong correlation with the binding energy of * OH.In particular, on the basis of ML model with random forest algorithm, three important properties were identified that determined the * OH adsorption as well as the CO 2 conversion, including the polarized charge (Q) of the central metal, the interaction (ICOHP) between the central metal and support, and the magnetic moment (Mag) of the central metal atom, where the polarized charge played a critical role.Promisingly, such ML model showed a robust capability to predict the catalytic activity of other types of SACs, offering a promising aspect to quickly screen the candidate materials for the CO 2 RR.Therefore, this work not only helped to understand the CO 2 RR mechanism of SACs based on GDY ang HGY supports but also to guide the design and synthesis for high-performance SACs toward CO 2 RR.

Figure 1 .
Figure 1.a) Optimized 2 × 2 GDY and HGY monolayer supercells and possible adsorption sites for single metal (M) atoms.b) Calculated binding energy and the energy difference between the binding energy and cohesive energy for 32 single metal (M) atom anchored GDY (M/GDY) and HGY (M/HGY) monolayer.For comparison, the E b values of Co/GDY, Mn/HGY, Y/HGY, and Zr/HGY are also given as shown in Figure 1b.c) AIMD energy profiles of Cr/GDY and Cr/HGY at 400 K, where the fluctuation of energy comes from the thermal disturbance of temperature.Optimized initial structure (IS) and final structure (FS) are as the insert.

Figure 2 .
Figure 2. a) The adsorption energy (ΔE) of CO 2 molecule and H atom on M/GDY and M/HGY.The relationship among b) average C─O bond length, c) the O─C─O bond angle, and the charge transfer amount of the adsorbed CO 2 molecule.d) The most stable adsorption structure, charge density difference and e) PDOS and COHP for CO 2 adsorbed on Ti/GDY.0.005 e Å −3 is set for the isosurface level and yellow and cyan denote the charge accumulation and depletion.0 eV is set for the fermi level.The bonded and anti-bonded states of COHP are yellow and cyan.f) The "donationbackdonation" mechanism between the adsorbed CO 2 molecule and M d orbital.
genation step of CO 2 RR ( * CO 2 + H + + e − → * COOH or * CO 2 + H + + e − → * OCHO) and HER (H + + e − → 1 2 H 2 (g)) on the 32 SACs, as shown in Figure 4. Figure 4 is divided into two areas, where the region below the dashed line suggests the high selectivity for CO 2 RR. 6 M/GDY (M = Ti, V, Y, Nb, Hf, Ta) and 6 M/HGY (M = Ti, V, Fe, Co, Nb, Hf) are below the dashed line, regardless of forming * COOH or * OCHO intermediate.Another 12 catalysts that preferentially adsorb CO 2 are easier to get intermediate * OCHO instead of hydrogen.Interestingly, although CO 2 and H have comparable adsorption capacity, more negative ΔG[ * OCHO] enables Re/GDY preference to catalyze CO 2 to yield * OCHO intermediate rather than HER.Therefore, these 25 catalysts are demonstrated to have high CO 2 RR selectivity; and thus, are further explored for the CO 2 reduction reaction.

Figure 3 .
Figure 3.The calculated PDOSs of Ti d orbital in Ti/GDY a) before and b-f) after CO 2 adsorption.0 eV is set for the fermi level.

Figure 4 .
Figure 4. Comparison of Gibbs free energy changes of CO 2 RR versus HER on M/GDY and M/HGY in the initial protonation step.
), indicating the generation of two-electron reduction product HCOOH is suppressed.Then, * H 2 COO intermediate is hydrogenated to form * H 2 COOH with the free energy changes of −0.05 and 0.70 eV for V/GDY and V/HGY, respectively.Once the * H 2 COOH is formed, it can be readily reduced to * O + CH 4 by the consecutive hydrogenation steps ( * H 2 COOH → * CH 2 O → * CH 3 O → * O + CH 4 ).No potential barriers are found in these steps, indicating these processes can easily occur.Subsequently, the formed * O intermediate is converted to H 2 O after two consecutive hydrogenation steps, of which * OH + H + + e − → * H 2 O has the largest potential barriers, which are regarded as the potential determining steps (PDS) of V/GDY and V/HGY, and the corresponding U L are estimated to be −0.67 and −0.95 V.
13 eV.Next, * CHO can be hydrogenated to * CHOH or * CH 2 O.As the formation barrier for * CHOH (0.15 eV) is larger than that of * CH 2 O (0.12 eV), the * CHO → * CH 2 O is more favorable.In the following steps, CH 4 can be formed via * CH 2 O → * CH 3 O → * O + CH 4 → * OH → * H 2 O with the continuous downhill free energy, and the large desorption energy (0.38 eV) of CH 2 O can effectively suppress the molecule desorption which makes the process feasible.Among all elementary steps, * HCOOH → * CHO shows the largest free energy change (0.13 eV); and thus, is the potential determining step.

Figure 5 .
Figure 5. a) Schematic illustration of the elementary reaction pathway during the CO 2 reduction.b) Summary of limiting potential (U L ) on 25 catalysts.The Gibbs free energy profiles of CO 2 RR on c) V/GDY, d) V/HGY, e) Fe/GDY, and f) Fe/HGY.

Figure 7 .
Figure 7. a) ML predicted ΔE * OH versus DFT calculated ΔE * OH on M/GDY and M/HGY and b) the feature importance in the random forest model.The inset shows the probability density distribution of R 2 prediction error for training set and test set.c) The violin plot for net charges of metal atoms (Q) versus ΔE * OH , where 0 and 1 present ΔE * OH < −4.04 eV and ΔE * OH ≥ −4.04 eV, respectively.d) Charge transfer from M/GDY and M/HGY to CO 2 as a function of polarized charges of metal atoms in M/GDY and M/HGY.