Electrochemical carbon dioxide reduction reaction (eCO₂RR) provides an attractive way to convert CO₂ into valuable fuels and chemicals using renewable energy under mild and controllable conditions. On the one hand, it can realize the efficient storage of electric energy and alleviate the energy crisis. On the other hand, it can reduce the emission of CO₂ into the atmosphere and alleviate the greenhouse effect [1,2]. Therefore, eCO₂RR has broad application prospects. Despite a considerable amount of experimental and theoretical reports [3–8], there are still some critical challenges to the technological viability of eCO₂RR. In the process of eCO₂RR, high overpotential is required to activate the stable C=O bond [9]. At the same time, the competition between hydrogen evolution reaction (HER) and multiple reaction paths leads to low Faradaic efficiency and poor product selectivity [10]. Therefore, finding economical and efficient catalysts is the key to promote the development and application of electrocatalytic CO₂ reduction technology.

Natural enzymes play a vital role in living systems due to their excellent catalytic activities and selectivity. However, the high preparation cost, easy deactivation and limited reaction conditions hinder their practical application [11]. One of the improved strategies is to develop enzyme-mimic materials that can overcome the above drawbacks. In particular, single-atom nanozymes (SANs) have recently received considerable attention [12,13]. SAN is a special kind of single-atom catalysts (SACs). It is well known that SACs have the structural characteristics of high atom dispersion and ultra-high atom utilization, and they show high activity in many catalytic reactions [14,15]. In addition to these characteristics of SACs, SANs also mimic the coordination environment or geometric configuration of the enzymes and exhibit specific product selectivity for specific reactions. Moreover, the catalytic activity of SANs can be easily regulated by the microenvironment, and the required reaction conditions are relatively mild [16]. Due to the characteristics of the structure and performance, SANs have aroused increasing research interest in the electrocatalytic field [17–19]. For example, TM-N-C metal enzymes with active sites similar to peroxidase (POD), oxidase (OXD) and porphyrin have been extensively studied. Fe-N-C SAN with FeN₄ sites showed high peroxidase-like activity of 25.33 U/mg for the electroreduction of H₂O₂ [20]. SAN with nanoframe-confined FeN₅ active sites, which was designed by Huang et al. [21] with a similar structure to the axial ligand–coordinated heme of cytochrome P450, showed excellent catalytic activity for oxygen reduction reaction. In addition to the N-coordinated structure, the TM-S site is also the active center of many natural enzymes, such as nitrogenase, formate dehydrogenase (FADH), and CO-dehydrogenase (CODH). Wen et al. [22] theoretically designed a series of graphene-supported MoX₄ (X = C, N, B, S, P) SACs, and found that MoS₄/GR with the same coordination as nitrogenase can catalyze nitrogen reduction reaction (NRR) at a low limiting potential of −0.29 V. Recently, the active sites [NiS₄] of FADH and CODH were simulated experimentally [23], and the prepared MOFs incorporated with nickel bis(dithiolene) ligands significantly improved the conversion rate and Faradaic effciency (89%) of CO₂ electroreduction to HCOOH. In addition, carbon-supported CoS₄ SAC was also successfully synthesized by pyrolysis of Co-MOF-74 in a strongly polar molten salt system [24], and the counter electrode of this SAC showed higher photo-electric conversion effciency than Pt counter electrode, and has good structural stability.

Inspired by the active site structures of FADH and CODH and the latest experimental synthesis of CoS₄, in this study, we used density functional theory (DFT) to calculate and design a series of enzyme-like single-atom catalysts supported on graphene with TMS₄-C active center for electrochemical CO₂RR, where the TM atoms are the Group Ⅷ transition metals (TM = Co, Ni, Rh, Ir). Since eCO₂RR occurs in the electrolyte, in order to better consider the solvent effect, the hybrid solvation model was adopted in the calculation. To further clarify the influence of coordination environment on catalytic performance, we also studied the CO₂RR mechanism of typical TMN₄-C catalysts for the selected TM of above TMS₄-C SANs. DFT results show that CoS₄-C SAN has a theoretical limiting potential of −0.07 V for the highly selective electroreduction of CO₂RR to HCOOH. TMS₄-C SANs have better CO₂ catalytic performance than the corresponding TMN₄-C SACs. In addition, two descriptors of TMX₄-C (X = N, S) CO₂RR activity based on intrinsic properties were proposed. Our current results may provide new insights into the design of advanced CO₂RR electrocatalysts.

Spin-polarized density functional theory (DFT) calculations were performed by the Vienna Ab initio Simulation Package (VASP) using the generalized gradient approximation (GGA) with Perdew−Burke−Ernzerhof (PBE) functional, with a cutoff energy of 400 eV [25–28]. A Monkhorst−Pack mesh of 3 × 3 × 1 k-points was used for geometric optimization and 9 × 9 × 1 for electronic structure analysis. The long-range van der Waal's interactions were described using the DFT-D3 empirical dispersion correction [29]. The energy and force convergence criteria were set to 10⁻⁵ eV and 0.02 eV/Å, respectively. The vacuum was set as 20 Å to avoid interactions between periodic images in the z direction. Referring to the previous reports [30–32], the graphene monolayer was modeled using a 4 × 4 × 1 supercell. Moreover, the hybrid solvation was adopted to describe the solvation effect. A hexagonal water bilayer containing 16 H₂O molecules was added to simulate surrounding water molecules, and Poisson–Boltzmann implicit solvation model was utilized to describe long-range interactions, as implemented in the VASPsol [33].

To evaluate the thermodynamic and electrochemical stability of SANs, we calculated the formation energy (E_f) and dissolution potential (U_diss, vs. SHE), which are defined as

where E_total and E_sub are the energies of SAN and substrate, respectively, and E_TM is the energy of the single metal atom in the bulk; U_diss^o(metal, bulk) is the standard dissolution potential of bulk metals, and N_e is the number of electrons involved in the dissolution [34]. It can be seen that Ru-, Os- and PtS₄-C with positive E_f cannot maintain stability due to the thermodynamically favorable diffusion and aggregation of metal atoms, and FeS₄-C is easily deactivated due to the dissolution of Fe atoms under electrochemical conditions (Fig. 1a). The stability of the other 5 SANs is further confirmed through 20 ps AIMD simulations with a time step of 1 fs in the NVT ensemble at 400 K. It can be seen that Co-, Ni-, Rh- and IrS₄-C maintain stable TM-S₄ structures after 20 ps AIMD simulations (Fig. 1b and Fig. S1 in Supporting information). In addition, although PdS₄-C SAN with S-Pd-S coordination could hold stability through the AIMD simulation, it undergoes a great distortion when the intermediates are adsorbed on the surface (Fig. S2 in Supporting information). Thus, it was not considered in the following mechanism study.

Due to the single-atom characteristic of SANs, only C₁ products were considered in the current work. As shown in Scheme S1 (Supporting information), proton-electron pairs attack CO₂ to form different initial configurations (*COOH and *HCOO), resulting in different pathways of eCO₂RR, and finally generate C₁ products, HCOOH, CO, CH₃OH and CH₄ via multi-step proton coupled electron transfer (PCET). Fig. 2 (left) shows the free energy diagrams of CO₂RR on TMS₄-C. Reduction of intermediates (such as * H₂COO, * COH) that are too high in the free energy diagrams was not considered further. To facilitate comparison, the corresponding CO₂RR free energy curves of TMN₄-C catalyst are shown in Fig. 2 (right). The free energy change of each elementary step was calculated according to the computational hydrogen electrode (CHE) model proposed by Nørskov et al. [35]. The limiting potential of CO₂RR is defined as U_L = -ΔG_PDS/e, where the ΔG_PDS is the free energy change of the potential-determined step (PDS).

As shown in Fig. 2a, the formation of intermediates *HCOO and *COOH on CoS₄-C SAN requires the free energy of 0.07 and 0.37 eV, respectively. The hydrogenation of these two intermediates directly forms HCOOH(dl) at the catalysts-electrolyte interface instead of adsorbed *HCOOH, which prevents the continuous hydrogenation of *HCOOH. In addition, *COOH can be further converted into *CO. However, the release of CO or further hydrogenation of *CO to *CHO needs to overcome a large free energy barrier, which is 1.00 and 2.01 eV, respectively, indicating the formation of other C₁ products via *CO, such as CH₃OH and CH₄, is unfavorable. Therefore, HCOOH is the predominant product of CO₂ reduction on CoS₄-C, and the PDS is * + CO₂ + H⁺ + e⁻ → *HCOO, with the U_L of −0.07 V. The reaction mechanism on the surface of CoN₄-C is quite different from that of CoS₄-C. It can be seen from Fig. 1b that the intermediate *HCOO is hydrogenated to form *HCOOH, and *COOH is hydrogenated to form *CO. *HCOOH is formed via the hydrogenation of *HCOO, and *CO is formed via the hydrogenation of *COOH. *HCOOH and *CO either desorb from the surface to generate 2e⁻ products HCOOH and CO due to the weak adsorption of *HCOOH and *CO, or occur subsequent hydrogenation to generate 6e⁻ product CH₃OH and 8e⁻ product CH₄. The PDS of the formation of HCOOH, CH₃OH and CH₄ is the step of *HCOO + H⁺ + e⁻ → *HCOOH, with a limiting potential of −0.30 V, and the PDS of the formation of CO is the step of * + CO₂ + H⁺ + e⁻ → *COOH with U_L= −0.46 V. These results show that CoN₄-C can catalyze the reduction of CO₂ to multiple C₁ products, but suffers poor product selectivity. Obviously, CoS₄-C shows better activity and product selectivity than CoN₄-C. As shown in Fig. 2b, both NiS₄-C and NiN₄-C tend to generate HCOOH through the HCOO pathway, and the first hydrogenation step * + CO₂ + H⁺ + e⁻ → *HCOO is the PDS, corresponding to the U_L of −0.42 and −0.66 V, respectively. So NiS₄-C exhibits better performance for HCOOH production than NiN₄-C. On RhS₄-C, as shown in Fig. 2c, CO₂RR to HCOOH via the *HCOO intermediate is more favorable than CO generation via *COOH. The first dehydrogenation step is the PDS for HCOOH formation, and the U_L is −0.84 V. For RhN₄-C, in addition to generating HCOOH, since *CO is facile to be hydrogenated to *CHO, which makes the subsequent PCET steps feasible, various products may be generated in the CO₂RR process, including CO, HCOOH, CH₃OH, CH₄, with the same PDS (* + CO₂ + H⁺ + e⁻ → *COOH). The limiting potentials of RhS₄-C and RhN₄-C are −0.84 and −1.08 V, respectively. IrS₄-C mainly generates HCOOH and CO via *COOH at the limiting potential of −0.93 V, while for IrN₄-C, four C₁ products CO, HCOOH, CH₃OH, and CH₄ will be produced at the limiting potential of −1.16 V. And the PDS of these two catalysts is the step of * + CO₂ + H⁺ + e⁻ → *COOH (Fig. S3 in Supporting information). Thus, Rh/IrS₄-C also shows better activity and product selectivity than Rh/IrN₄-C.

Hydrogen evolution reaction (HER) is the main competing side reaction involved in electrocatalytic CO₂ reduction process, resulting in low Faradaic efficiency. The free energy diagrams of HER are given in Fig. 3a. It can be seen that Co-, Rh- and IrS₄-C suppress the HER in comparison with corresponding TMN₄-C SACs, while NiS₄-C shows better HER performance than NiN₄-C. Furthermore, in order to evaluate the selectivity preference, we compared the difference of the limiting potentials between CO₂RR and HER. As shown in Fig. 3b, the U_L(CO₂RR) - U_L(HER) value of Co/NiS₄-C and NiN₄-C was positive, which means that they exhibit better selectivity to CO₂RR than to HER. In addition, although U_L(CO₂RR) - U_L(HER) values of RhS₄-C and IrS₄-C are negative, they are more positive than the corresponding values of TMN₄-C catalysts (−0.08 V vs. −0.67 V; −0.43 V vs. −0.88 V), indicating that compared with TMN₄-C coordination, TMS₄-C improves the selectivity of the catalyst towards CO₂RR.

The above results indicate that these TMS₄-C SANs show better activity and selectivity than widely studied TMN₄-C catalysts. Compared with Cu(211) (U_L = −0.61 V), CoS₄-C and NiS₄-C SANs exhibit better catalytic performance. Specially, the CoS₄-C SAN that has been prepared in experiment owns the FADH-like activity and selectivity for the production of HCOOH, and can generate HCOOH at a low U_L of −0.07 V, showing its great practical application potential. To explore the nature that TMS₄-C SANs have better catalytic activity than TMN₄-C, the electronic structure analysis of these catalysts was performed. The partial charge transfer of these two kinds of SACs was analyzed by Bader charge analysis [36]. The charge density difference (CDD) of TMS₄-C SANs and TMN₄-C SACs is plotted in Fig. 4. As can be seen from the figure, since the electronegativity of S is lower than that of N, the TM atom donates less electrons to S atoms than to N atoms, resulting in lower valence state. For the same metal, the greater the electron density of TM atom caused by the low electronegativity of S, the higher the reduction ability.

The distribution of electrons in the d_z² orbitals has a great influence on the catalytic activity of SACs, because intermediates are usually axially adsorbed on the single-atom site, which has been reported in previous studies [37,38]. Here, the d-band centers of axial d_z² orbitals (ε(d_z²)) of transition metals in TMS₄-C and TMN₄-C were calculated and shown in Fig. 5. It is seen that the ε(d_z²) of TMS₄-C is closer to the Fermi level than the ε(d_z²) of the corresponding TMN₄-C SACs, indicating that compared with TMN₄-C, the S-coordination causes the d-band centers of d_z² orbitals to move upwards, thus enhancing the interaction between TMS₄-C and intermediates *HCOO (Co-, Ni-, RhS₄-C) and *COOH (IrS₄-C). As the first hydrogenation step of CO₂ to *HCOO/*COOH is the PDS of the electrocatalytic CO₂RR on TMS₄-C and TMN₄-C, the stronger interaction between TMS₄-C surface and *HCOO/*COOH makes it show better activity than TMN₄-C.

It is meaningful to explore the descriptors highly correlated to the catalytic activity in order to guide the design of novel CO₂RR electrocatalysts. Based on the intrinsic properties of the catalysts, we proposed a descriptor φ₁, which is defined as φ₁ = χ_TM/χ_NM - N_d – N_es, where χ_TM is the electronegativity of TM atom, χ_NM is the electronegativity of nonmetal atom (S or N) coordinated with the TM atom, and N_d and N_es are the number of d electrons and the number of electron shell layers, respectively. Fig. 6a shows that the U_L of TMS₄-C and TMN₄-C presents a linear relationship with φ₁, indicating that φ₁ can serve as an activity descriptor for CO₂RR, and the electronegativity of active center atoms and coordination atoms has a great influence on the catalytic activity. In addition, for TMS₄-C and TMN₄-C catalysts, considering the coordination environments and TM-support interaction, we plotted the variation of the limiting potentials of CO₂RR on all these SACs with φ₂, where φ₂ represents the difference between the d-band center of TM atom (ε_d) and p-band center of coordinated S or N atoms (ε_p), that is, φ₂ = ε_d – ε_p. Interestingly, it can be seen from Fig. 6b that the dependence of U_L as a function of φ₂ displays a volcano-shaped curve, in which the optimal CoS₄-C SAN with the best performance for CO₂RR lies around the peak of the volcano, with φ₂ of 2.31. Thus, φ₂ can be used as a good descriptor to exhibit the CO₂RR activity. The finding of these two descriptors with reliable R² will be beneficial to give a reliable prediction for the catalytic ability of more similar systems.

In summary, we designed carbon-supported Group Ⅷ TMS₄-C single-atom nanozymes (SANs) with FADH-mimic TMS₄ sites for electrochemical CO₂ reduction reaction (eCO₂RR) through DFT calculations combined with the hybrid solvation model. After the evaluation of thermodynamic and electrochemical stability, four transition metal Co, Ni, Rh and Ir SANs were selected to study the CO₂RR mechanism, and their catalytic performance was compared with that of TMN₄-C. Gibbs free energy profiles show that TMS₄-C selectively reduces CO₂ to HCOOH, while TMN₄-C generates multiple C1 products. Moreover, TMS₄-C SANs show better catalytic ability toward CO₂RR than corresponding TMN₄-C SACs, which can be attributed to the more electrons on TM atoms and the upshift of d-band center of d_z² orbitals to the Fermi level. More importantly, the origin of the CO₂RR activity of TMS₄-C and TMN₄-C is revealed by the two descriptors (φ₁ and φ₂), which are highly related to the intrinsic parameters and electronic structure properties. The limiting potentials of these catalysts exhibit a volcano-shaped curve with φ₂, and CoS₄-C stands near the top of the volcano. Its U_L (−0.07 V) is far lower than the results reported previously, indicating that CoS₄-C is a potential SAN catalyst for CO₂ reduction. These findings are helpful to predict the catalytic activity of SACs and can provide new clues and ideas for the design of novel electrocatalysts for CO₂ reduction.

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

This work was supported by the National Key Research and Development Program of China (No. 2021YFA1500403), and the National Natural Science Foundation of China (No. 21773083).

Supplementary material associated with this article can be found, in the online version, at https://doi.org/10.1016/j.cclet.2022.108018.