The Screening of Homo‐ and Hetero‐Dual Atoms Anchored Graphdiyne for Boosting Electrochemical CO2 Reduction

Developing electrocatalysts with high catalytic performance and selectivity is crucial for electrochemical CO2 reduction reaction (CRR). There are many catalyst studies of transition metal (TM) atom doping to sp2 carbon material, such as graphene or carbon nanotubes. On the other hand, graphdiyne (GDY) has both sp and sp2 hybridization and stable pores, so we can tune its interaction with TM. Following the successful experimental synthesis of Ni atom doping to GDY monolayer, the CRR activity of double–atom catalysts was evaluated, including homo and hetero metal‐Ni doped on the GDY monolayer (MNi@GDY where M is Ti, V, Cr, Mn, Fe, Co, Ni, and Cu) using the density functional theory calculations. The valence‐electron number of the catalytic center shows a strong positive correlation to the limiting potentials in the volcano plot. NiNi@GDY is the most promising candidate for converting CO2 to produce CH4 with a remarkable low limiting potential of −0.28 V, which is better than Ni@GDY and Ni3@GDY counterparts. NiNi@GDY shows excellent thermal stability and ability to suppress the competing hydrogen evolution reaction, showing its high selectivity to CH4.


Introduction
Overusing fossil fuels hastens the depletion of energy resources and the excessive release of CO 2 and contributes to the energy crisis, pollution, and global warming. [1] One of the most effective approaches to address these difficulties is to convert CO 2 into valuable chemical fuels. [2,3]  For example, through the "precursor-preselected" strategy, Li et al. have synthesized the Fe dimers doped C 3 N 4 surface. This DAC shows outstanding catalytic performance for the epoxidation of trans-stilbene to trans-stilbene oxide. [20] Several types of DACs have also gotten much attention from theoretical calculations. [30][31][32][33][34][35][36][37] Surprisingly, it has been discovered that in some reports, DACs outperform their SAC counterparts in terms of catalytic performance. Double transition metal atoms doped on g-C 3 N 4 are good electron donors and show a superior catalytic performance than single atom counterparts for O 2 activation. [38] Sun et al. studied phthalocyanine sheets with Mn dimer, which exhibited greater activity for CRR than single Mn embedded phthalocyanine. [30] Chen et al. investigated Cu dimers doped on a C 2 N surface, which illustrated improved catalytic efficacy for CRR with a low onset potential of −0.23 V. [31] By density functional theory (DFT), Jiang et al. studied non-noble metal atom dimers supported on C 2 N surfaces as electrocatalysts for NRR. Compared to the performance of corresponding SACs, these studies show that all DACs had lower onset potentials. [33] Moreover, in addition to homonuclear DACs, heteronuclear DACs have also been synthesized recently. Zhang et al. synthesized the PtRu dimer supported on the N-doped carbon nanotube (NCNT), using atomic layer deposition (ALD) technology. This PtRu-NCNT shows greater HER activity and good stability than the commercial Pt/C catalysts. [21] FeCo [22,23] and NiFe [24] bimetallic anchored on the N-doped carbon substrate also have been used for ORR [22,23] and CRR, [24] respectively. Heteronuclear NiMn doped on g-C 3 N 4 was successfully synthesized by Ding et al. [39] and showed superior catalytic performance for CRR than mono-metal counterparts. Our previous report demonstrated that NiMn anchored g-C 3 N 4 substrates are promising candidate materials in CRR by booting the CC coupling processes. [40] Based on the literature search above, DACs have considerable potential in electrocatalysis.
One of the most promising carbon allotropes with a singleatom thickness is the graphdiyne (GDY) monolayer. [41] Because of its distinct structural, physical, and chemical features, it has numerous potential use in catalysis, rechargeable batteries, electronics, and other fields. [33,[42][43][44][45][46] The pores of GDY are consistently distributed in sub nanoscale and are surrounding sp-hybridized carbon atoms. This makes GDY an ideal substrate for firmly anchoring atomic metal catalysts containing one or more transition metal atoms. [46][47][48][49][50][51][52][53][54][55][56] Single atomic catalysts Ni@GDY, [49] Fe@GDY, [49] Mo@GDY, [46] and Pd@GDY [53] have been synthesized for some crucial electrochemical reactions, including HER (Pd@GDY [53] and Fe@GDY [49] ), NRR (Mo@GDY), [46] and ORR (Fe@GDY). [50] The potentials of GDYbased SACs for CO oxidation, [49] water splitting, ORR, [48] and NRR [55] have all been proposed theoretically. Furthermore, Xing et al. investigated Ru x Os 3−x trimers embedded in GDY sheets for semi-hydrogenation of acetylene to ethylene and found that all of them have good catalytic activity and selectivity. [52] He et al. presented a theoretical study of the active-site dependent activity and selectivity for CRR on Fe n clusters with a few atoms (n = 1-4) doped on the GDY surface. Fe dimer and trimer show the highest catalytic activity and selectivity. [57] Ma et al. have investigated single atom, homonuclear double, and triple atoms of Mn, Fe, Co, and Ni anchored GDY surface for NRR, and found that the DACs exhibit greater NRR activity than the single or triple metal catalysts. [56] Feng et al. constructed single-, double-, and triple-Cu doped on GDY monolayer (Cu 1-3 @GDY) and discovered that Cu 2 @GDY catalyst could efficiently convert CO 2 to CH 4 . [58] As previously stated, the DACs with homo-and heteronuclear metal dimers have been synthesized. [21][22][23][24] Moreover, Ni@GDY has also been synthesized; [49] here, we theoretically modeled the doping of M-Ni-pairs (M = Ti, V, Cr, Mn, Fe, Co, Ni, and Cu) embedded GDY surface (denoted as MNi@GDY) as the electrocatalysts for CRR. We propose to achieve three major goals: 1) examine the structural stability and electronic structure of these DACs; 2) investigate energetically favorable reaction routes and provide the most efficient electrocatalyst for CRR; 3) examine the original activity of these DACs in CRR reaction to understand and develop a best CRR candidate. Our computational findings suggest that NiNi@GDY is the most promising. Its computed limiting potential (U L ) is lower than its heteronuclear counterparts. The HER is suppressed by NiNi@GDY, which is proposed to have high catalytic activity and product selectivity. Furthermore, NiNi@GDY shows much superior activity and selectivity toward CRR against HER than its single or triple atomic catalysts, which were reported to have a high CRR activity. [26,58] This work shows that precisely doping NiNi dimer in GDY can boot the CRR with high selectivity. The study not only predicts the cost-effective electrocatalysts for CRR, but also suggests ways to improve the catalytic activity and selectivity of the atomic metal catalysts supported on a suitable substrate.

Computational Details
Spin-polarized and periodic DFT calculations were performed using the VASP code [59,60] with projector augmented wave method, [61] and Perdew-Burke-Ernzerhof functional of generalized gradient approximation. [62] Grimme's van der Waals correction, often known as the DFT-D3 approach, was used to account for the dispersion interactions between the adsorbates and the substrate. [63] The plane-wave cutoff energy was set to 450 eV. The conjugate gradient algorithm was employed in ionic optimization, and the convergence criteria were set to the residual energy and force less than 10 −5 eV and 0.02 eV Å −1 , respectively. A 2 × 2 unit cell of GDY surface with a vacuum region in the z-direction of 15 Å was employed as the substrate. The simulated supercell of GDY is shown in Figure S1, Supporting Information. k-points were sampled using 3 × 3 × 1 and 11 × 11 × 1 Monkhorst-Pack mesh [64] for the geometry optimization and electronic structure calculations, respectively. The Bader charge analysis was used to evaluate the charge transfer between the substrate and the adsorbate molecules. The climbing image nudged elastic band (CI-NEB) [65] was adopted in the transition state search to determine the kinetic barrier of the reactions. Ab initio molecular dynamic (AIMD) simulations [66] were employed to assess the thermal stability of the catalysts. Details of the calculations, including the free energy calculations and the limiting potentials of the CRR and HER, are given in the Supporting Information.
where E MNi@GDY , E M@GDY , E Ni@GDY , E M , and E GDY are the total energy of MNi@GDY, M@GDY, Ni@GDY, the isolated metal atom in vacuum, and the GDY substrate, respectively. Figure 1c shows that all the M atoms have a large negative binding energy value, indicating a strong interaction between the M atom and the substrate. One important aspect is that the binding energy of M in the Ni@GDY is stronger than in the GDY. This indicates that the presence of the Ni atom helps in stabilizing an adjacent M at H1 sites. Such cooperative effect is obviously seen from the binding of Ni on Ni@GDY, −4.55 eV, which is more negative than the binding on the bare GDY substrate, −4.32 eV (Table S1, Supporting Information). Furthermore, a recent report has shown that on a bare GDY surface, the first row of transition metal atoms from Ti to Cu can be firmly trapped in the form of single-atom catalysts. [48] Thus, the presence of the Ni atom adjacent to the M atom has positive effects in stabilizing the M adsorption and preventing metal clustering. To gain an insightful understanding of the interaction between the embedded transition metal atoms and the Ni@GDY monolayer, we investigate the partial density of state (PDOS) curves of MNi@GDY systems, as displayed in Figure S3, Supporting Information. Most of the MNi@GDY systems show magneticity, as seen in asymmetric DOS peaks between the up and down spins. One crucial point is that this asymmetry mainly comes from the hetero metal atoms, the M 3d orbitals in MNi@ GDY. Since the asymmetrical DOS peaks are dominated by the metal 3d orbital, we expect that spin densities are localized on the M atom. The results in Table S2, Supporting Information, confirmed that the magnetic moment for the whole system is mainly attributed to the hetero metal atom M in MNi@GDY.

Activation of CO 2 on MNi@GDY
The initial state of the CRR is the adsorption of CO 2 onto the active site of the catalysts. The CO 2 is activated by electron donation into the antibonding 2π u orbitals, which invariably reduces the OCO bond angle. As displayed in Table 1, the binding strength between the adsorbed CO 2 and a MNi@GDY is very strong, with adsorption energies decreasing in the order of TiNi (−1. , respectively, which correspond well to the bent degree of *CO 2 (∠OCO: 121.8° to 179.9°). As seen in Table 1, there is a significant charge transfer between the metal dimers and the CO 2 , and the stronger CO 2 binding, the more charge transfer from the surfaces. We selected the TiNi@GDY and CuNi@GDY systems before and after CO 2 adsorption to illustrate the difference between the strongest and weakest interactions, as shown in Figure 2. Those for other systems are given in Figures S4 and S5, Supporting Information. From the charge density plot shown in Figure 2a, we confirmed that charge transfers from TiNi@GDY catalyst to CO 2 , while in the CuNi@GDY, we hardly observed chargedensity difference between adsorbate and catalyst as depicted in Figure 2b. Moreover, the partial density of states or PDOS plot of TiNi@GDY and CuNi@GDY surface before and after CO 2 adsorption conditions displayed in Figure 2c and d clearly show a strong hybridization between Ti d-orbitals and p-orbitals of CO 2 www.advmatinterfaces.de molecule at the energy level of 0 to −4 eV. On the other hand, on the CuNi@GDY surface, we see a minor change in the PDOS before and after interaction with the CO 2 . Additionally, we found a peak near the Fermi level in the PDOS of TiNi@GDY surface, as displayed in Figure 2c, which mainly comes from the up spin Ti 3d-orbital. We confirmed that the spin density is mainly localized over the Ti atom, as shown in Table S2, Supporting Information, which means that the active unpair electrons in the Ti atom can readily couple with the approaching CO 2 molecule. For the CuNi@GDY surface given in Figure 2d, the PDOS is symmetric, the electrons form stable pairs, and the total magnetic moment is nearly zero, thus resulting in a non-reactive surface. These findings show that the interaction between the metal dimers in GDY and CO 2 molecule can be tuned. We noted that the binding of CO 2 side-on to the VNi@GDY surface results in a dissociation of the CO bond where the dissociated O atom binds strongly with the V atom while the CO binds with the bridge VNi atoms, as shown in Figure S6, Supporting Information.

Electrocatalytic CO 2 Reduction on MNi@GDY
Next, we evaluated the Gibbs free energy of all possible reaction intermediates to find the favored pathways during CRR by comparing the potential-determining step (PDS) of each pathway. Figure 3 depicts the various CRR reaction pathways considered on the MNi@GDY sheets. Interestingly, for the early transition metals (M = Ti, V, Cr, and Mn), the last hydrogenation step (*OH + H + + e − → *H 2 O), identified as the PDS with ∆G PDS , is greater than 0.75 eV (see Figure S7, Supporting Information). This rules out their potential to reduce CO 2 because the *OH would block the catalytic center, leading to the deactivation of the catalysts. Thus, TiNi@GDY, VNi@GDY, CrNi@GDY, and MnNi@GDY surfaces are screened out from the study. Following the screening mentioned above, the GDY doped with FeNi, CoNi, NiNi, and CuNi configurations is further studied for CRR. We have investigated the possible CRR pathways on these four MNi@GDY catalysts and the product distribution. The calculated free energy profiles of CRR together   Figure 4, while their corresponding atomic structures along the most preferable path are presented in Figure 5.    Table 2. This suggests the superior CRR catalytic activity of these Ni@GDY doped with Fe, Co, Ni, and Cu. Furthermore, we have selected NiNi@GDY as an example to investigate the C 2 productions from the CC coupling reaction shown in Figure S8, Supporting Information, and found that the free energy barrier for the first protonation step of the most preferable route (*CO+*CO → *CO+*COH) is extremely high about 1.00 eV. The calculation results suggest that the MNi@GDY catalyst does not favor the production of C 2 products.

Scaling Relationships and Intrinsic Descriptor for CO 2 Reduction
To understand the activity trend of the CRR on the various MNi@GDYs, we calculated the adsorption free energies of all CRR intermediate species, as illustrated in Table S4, Supporting Information. We also plotted these adsorption free energies of all CRR species versus the adsorption free energies of *OH as shown in Figure S9, Supporting Information, the adsorption free energies of CRR intermediates are well correlated with ∆G *OH with the R 2 near 1.00. This is consistent with the results in Figure 5, showing that the intermediates favor O binding to MNi@GDY. Considering that the PDS is the *OH + H + + e − → *H 2 O step, we expect the U L values to correlate to the adsorption free energies of *OH on the catalysts, confirmed in Figure 6a. Among all the MNi@GDY surfaces, NiNi@GDY has the lowest adsorption free energy for the *OH (−0.72 eV). Correspondingly, the CRR pathway has the smallest U L (−0.28 V). On the contrary, TiNi@GDY has the highest adsorption free energy for the *OH (−2.04 eV) and the largest U L (−1.11 V). We note that for M@GDY, the CRR PDS was at different steps, [68] which makes it hard to control. By co-doping with Ni in the form of MNi@GDY, we find that we only need to control the *OH + H + + e − → *H 2 O reaction. As demonstrated in Figure 6a, the limiting potential of CRR for this catalyst is related to the *OH adsorption on MNi@GDY. To better understand its activity, we calculated the electronic structures of *OH adsorbed on TiNi@GDY and NiNi@GDY, as shown in Figure 7a,b. As mentioned above, on the TiNi@GDY surface, we have a peak in the up spin around the Fermi level that mainly comes from the Ti-3d state. However, for the NiNi@GDY surface, the up and down spin contributions are symmetric in both Ni atoms, which means it is a closed-shell system. As a result, after *OH binding, TiNi becomes a closed shell, while NiNi becomes an open shell. We can also see that the adsorbed *OH shows strong hybridization between the OH-2p orbital and the NiNi-3d orbital near the Fermi level on NiNi@GDY, indicating that OH is still reactive. In contrast, the strong adsorption of *OH on TiNi@GDY causes its orbitals to shift to a lower energy level, lowering the reactivity of *OH and making it difficult for further protonation. Thus, the strong adsorption of OH should cause a poisoning phenomenon on  TiNi@GDY. This might suggest that the binding mechanism of OH* is mainly regulated by the binding with *OH radical electron, which necessitates a radical spin at the MNi site. Looking at the spin moments in Table S2, Supporting Information, we can understand that the binding of *OH is not preferred on the NiNi site. In addition, we also performed the projected crystal orbital Hamilton population (pCOHP) analysis. The integrated COHP (ICOHP) per bond was computed. We obtained the ICOHP value for the TiO bond in TiNi@GDY system as −0.24 eV and for the Ni 2 O bond in the NiNi@GDY system as −0.09 eV. Here a more negative value of ICOHP implies a stronger metalOH bond, which accounts for the origin of CRR activity on NiNi@GDY.
Although adsorption free energies ∆G OH* can be utilized to estimate the reaction activity, a prediction based on simpler descriptors is preferable. To construct bimetallic atomic catalysts for CO 2 reduction, an effective simple descriptor (ψ) based on the valence electron number and electronegativity of the metal at the active centers was recently presented. [76] This descriptor effectively captures the local-environment impacts of active centers, and the related parameters in the descriptor are conveniently accessible. The descriptor, which depends on valence-electron number and electronegativity, is determined as follows: where N is the number of atoms at the active centers, S vi and χ i are the valence-electron number, and electronegativity of the i th atom at the active centers. The data of S vi , χ i , and ψ for MNi@ GDY are listed in Table S5, Supporting Information. Figure 6c shows the adsorption free energy of *OH as a function of the descriptor ψ. We found that the descriptor ψ exhibits a volcano relationship when plotted against the adsorption free energy of *OH. For M from the IVB group to the VIIIB group (Ti to Ni), ∆G *OH increases linearly with ψ, while ∆G *OH decreases with ψ for Cu. Moreover, we also plot the sum of valence d-electron of bi-metal atoms (S v1 + S v2 ) against the adsorption free energy of *OH as shown in Figure 6b, and we found that it shows a similar trend with the descriptor ψ against ∆G *OH . NiNi@GDY is located on the peak of the broken line in each system for the adsorption of *OH. This indicates that the adsorption free energy of *OH species and the CRR activity in this catalyst is correlated to the number of valance electrons in the metal atom bound to Ni.

Competition between CRR and HER on MNi@GDY
Hydrogen evolution reaction (HER) is a major competitive reaction in the CRR process. We estimated the free energy barriers of HER using the Volmer-Heyrovsky mechanism, as shown in Figure 8a. The difference between the limiting potentials (U L ) of HER and CRR is displayed as a function of CRR U L in Figure 8b. Having a higher U L (CRR) − U L (HER) value suggests selectivity toward CRR over HER, whereas a less negative U L (CRR) value indicates more active CRR. As a reference, we use a metal-based benchmark of −0.74 V from Cu (211), [77] which is one of the finest cathodes for CRR. We see that FeNi@GDY, CoNi@GDY, NiNi@GDY, and CuNi@GDY are in the upper right-hand corner, indicating that they can combine better catalytic activity and selectivity than the pure copper surface.

Stability and Coverage Dependence of NiNi@GDY
Here, we evaluate the kinetic stability of NiNi-GDY. The E a of NiNi breakdown into two separate Ni atoms on the GDY sheet was computed. The results show that separating NiNi into two Ni atoms on GDY sheets needs an E a of more than 2.0 eV ( Figure S11, Supporting Information), which is too high to proceed under ambient conditions. In addition, the thermal stability of NiNi@GDY is investigated using Ab initio molecular dynamics (AIMD) simulations with a time step of 2 fs at 500 K. The geometric structure of NiNi@GDY is well retained for 10 ps ( Figure S12, Supporting Information), demonstrating that NiNi@GDY has excellent thermal stability. We also calculated the dissolution potential, which shows positive value suggesting that this surface has electrochemical stability, see Table S6, Supporting Information. Although the stability of NiNi@GDY in the real CRR process needs to be tested experimentally, our results show that NiNi@GDY is potentially synthesizable with high kinetics and thermal stabilities. The www.advmatinterfaces.de charge transfer at each hydrogenation step was investigated to analyze the catalytic activity of NiNi@GDY, as shown in Figure S10, Supporting Information. The result shows that each adsorbate always gains an electron from the NiNi@GDY. We also observe that during the CRR process, there is a change in positive and negative charges carried on Ni 2 metal and GDY sheet, suggesting that the cooperations between the Ni 2 and GDY substrates help catalyze the CRR reaction. Lastly, since it is hard to experimentally control which hole the Ni dimer will fill up in the GDY, we also evaluated the CRR activities at various Ni dimer concentrations. We focused on how it may affect the PDS, *OH + H + + e − → *H 2 O. We added Ni 2 in a neighbor hole or increased the supercell from 2 × 2 to 3 × 3 GDY monolayer, as illustrated in Figure S13, Supporting Information. In Table S7, Supporting Information, we show that having a more dense or dilute Ni dimer does not affect the binding energy of the two species related to the PDS, *OH, and *H 2 O. Therefore, owing to the large separation of the holes in GDY, the activity of the local NiNi site will not be affected by the NiNi concentration.

A Comparison of CRR on Ni 1-3 @GDY Monolayer
Previously, Zhao et al. [78] reported that Ni 1-3 @GDYs were ruled out from their CRR study since they showed low CO 2 activation. Here, we present some differences from this previous study supporting the distinct performance for CRR of the Ni dimer embedded on GDY. To do this, the complete CRR mechanism for Ni 1-3 @GDYs was investigated. We also tested the case of four Ni atoms anchored on a GDY monolayer. Unfortunately, the Ni 4 cluster in one GDY hole is unfavorable in adsorption energy (Table S8, Supporting Information). It is more favorable to bind as Ni 3 and Ni clusters in two different GDY holes. Thus, we only consider the Ni 1-3 @GDY systems in the following study. We label Ni 1 , Ni 2, and Ni 3 @GDY for one, two, and three Ni atoms embedded in GDY, respectively. Interestingly, we found that Ni 2 @GDY can activate CO 2 resulting in an anionic bent CO 2 structure. As depicted in Figure 9, the physisorption of linear CO 2 molecules was seen for Ni 1 @GDY and Ni 3 @GDY. This point differs from the previous study by Zhao et al., [78] which only found physisorption for Ni 1-3 @GDYs. The Gibbs free energy changes along the reaction pathway are shown in Figure 9. We found that the Ni 1 @GDY tends to produce HCOOH as the final product. In contrast, Ni 3 @GDY prefers to generate CO. The desorption of the two products is easier than further protonation, and the PDS is the first protonation step for both Ni 1 @GDY and Ni 3 @GDY surfaces. Our result is consistent with the previous report for the dominant pathway and limiting potential of CRR on Ni 1 @GDY [68] and Ni 3 @GDY. [79] The results indicate that Ni dimer shows the highest catalytic activity for CRR compared with the corresponding Ni atom and Ni cluster embedding on GDY.
From Figure 9a-c, we can see that the first protonation step is the potential determining step (PDS) for both Ni 1 and Ni 3 @GDY surfaces, while the Ni 2 @GDY surface is not. To gain more insights into the origin of the superior activity of the Ni 2 @GDY, we evaluated the electronic interaction of the most preferable first protonation intermediate over Ni 1-3 @GDYs. We calculated the projected DOS (PDOS) of Ni atoms and *OCHO for Ni 1-2 @GDY and *COOH for Ni 3 @GDY catalysts in Figure 10. Comparing Figure 10a,b versus Figure 10c, we see that the interactions between Ni-3d and O-2p or C-2p orbitals of the first protonation intermediate are clearly different. The O-2p orbitals of *OCHO overlap more with Ni-3d orbitals around the Fermi level than the C-2p orbitals of *COOH. The computed adsorption energy of the most preferable first protonation intermediate over Ni 1-3 @GDY surface is −0.20 eV for Ni 1 @GDY, −0.58 eV for Ni 2 @GDY, and 0.42 eV for Ni 3 @GDY, respectively. The DOS analysis is consistent with the binding energy. As illustrated in Figure 10a,b, the bonding distance between Ni atoms and O atoms of *OCHO on Ni 1-2 @GDY is 2.06 Å for Ni 1 @GDY and 1.94 Å for Ni 2 @GDY. The synergistic effect of Ni 2 @GDY shows stronger binding with *OCHO than Ni 1 @GDY, and as a consequence, the free energy barrier of the first protonation step of Ni 2 @GDY is smaller.

Conclusion
In conclusion, by means of DFT calculations combined with the CHE model, we have systematically studied the potential of the homo-and heteronuclear DACs MNi@GDY (M = Ti, V, Cr, Mn, Fe, Co, Ni, Cu, and Zn) as efficient CRR catalysts. We found that all the 8 DACs considered in this study are kinetic and thermodynamic stable and can effectively bind and activate Figure 9. Gibbs free energy diagrams for CRR toward C1 product with the most stable structures on a) Ni 1 @GDY, b) Ni 2 @GDY, and c) Ni 3 @GDY at 0 V (vs RHE) with corresponding atomic structures along the most preferable path are presented in the right. The brown, yale blue, and red balls represent the C, Ni, and O atoms, respectively.