Origin of the N-coordinated single-atom Ni sites in heterogeneous electrocatalysts for CO₂ reduction reaction

Abstract

Heterogeneous Ni–N–C single-atom catalysts (SACs) have attracted great research interest regarding their capability in facilitating the CO₂ reduction reaction (CO₂RR), with CO accounting for the major product. However, the fundamental nature of their active Ni sites remains controversial, since the typically proposed pyridinic-type Ni configurations are inactive, display low selectivity, and/or possess an unfavorable formation energy. Herein, we present a constant-potential first-principles and microkinetic model to study the CO₂RR at a solid–water interface, which shows that the electrode potential is crucial for governing CO₂ activation. A formation energy analysis on several NiN_xC_4−x (x = 1–4) moieties indicates that the predominant Ni moieties of Ni–N–C SACs are expected to have a formula of NiN₄. After determining the potential-dependent thermodynamic and kinetic energy of these Ni moieties, we discover that the energetically favorable pyrrolic-type NiN₄ moiety displays high activity for facilitating the selective CO₂RR over the competing H₂ evolution. Moreover, model polarization curves and Tafel analysis results exhibit reasonable agreement with existing experimental data. This work highlights the intrinsic tetrapyrrolic coordination of Ni for facilitating the CO₂RR and offers practical guidance for the rational improvement of SACs, and this model can be expanded to explore mechanisms of other electrocatalysis in aqueous solutions.

---

Introduction

The high dependence of non-renewable fossil resources in modern economics has resulted in increasing concerns regarding energy and the environment. A promising solution to address this issue is to electrochemically convert waste CO₂ to fuels to achieve a carbon-neutral cycle.^1–3 Unfortunately, due to its sluggish kinetics and the concomitant hydrogen evolution reaction (HER), the feasibility of implementing such a CO₂ reduction reaction (CO₂RR) process on industrial scales suffers from the coupled problems of poor energy efficiency and low product selectivity. In particular, electrodes composed of earth-abundant metals generally exhibit low activity and/or selectivity for facilitating the CO₂RR. Research on developing effective CO₂RR catalysts has brought about the advent of a new class of heterogeneous electrocatalysts, namely, M–N–C single-atom catalysts (SACs), which feature atomically dispersed single-metal-atom sites.^4–8 These metal sites exhibit unique electronic and structural properties compared with their bulk counterparts, thereby endowing M–N–C SACs with outstanding catalytic performance.⁹

Ni–N–C, one of the most represented CO₂RR SACs, has been established to possess great potential for facilitating selective CO₂-to-CO conversion. Determining the role of coordination environments is beneficial for implementing strategic improvements to CO₂RR catalysts. While significant efforts have been expended to investigate the fundamental nature of the active Ni sites through employing physical characterizations and modern quantum mechanics methods, a common consensus pertaining to such an issue has not been reached.^3,9,16,17 Several pioneering studies have proposed a pyridinic-type NiN₄ moiety (referred to as “α-NiN₄” herein) in which a central Ni atom is coordinated by four pyridinic N atoms, and CO₂RR processes on such a moiety were explored through comprehensive density functional theory (DFT) calculations.^6–8 It should also be noted that pyridinic-type MN₄ is a prototypical moiety that has been successfully established for describing the oxygen reduction reaction (ORR) on various M–N–C SACs.^18–22 However, a detailed examination of the computational results of the CO₂RR on α-NiN₄ indicates a failure to reproduce the observed CO₂RR activities of Ni–N–C SACs since large theoretical overpotentials η^t (>1.6 V) are involved.^6–8

Recently, a variety of hybrid pyridinic N and C coordination environments with respect to the Ni sites have been proposed. For instance, Zhao et al. demonstrated that a pyridinic-type NiN₁C₃ moiety (referred to as “α-NiN₁C₃” herein) is considerably active for facilitating the CO₂RR with a small η^t value of 0.48 V and exhibits high selectivity toward CO over H₂.¹⁰ Such a hybrid moiety was also suggested in a recent grand canonical DFT study by Li and co-workers;²³ their calculations also revealed that the charge capacity and the formation of hydrogen bonds between the adsorbates and the surrounding water molecules play vital roles in facilitating the activation of CO₂. In contrast, Hossain et al. showed that pyridinic-type NiN₂C₂ and NiN₃C₁ moieties (referred to as “α-NiN₂C₂” and “α-NiN₃C₁” herein) are responsible for facilitating the CO₂RR at an electrode potential of −0.8 V.²⁴

While previous studies have predicted that Ni moieties containing fewer N atoms, i.e., α-NiN_xC_4−x (x = 1–3), exhibit excellent CO₂RR activity, a general situation in which the formation energy of an MN₄ moiety is much lower than those of MN_xC_4−x (x = 1–3) moieties has been overlooked.^25–27 From a thermodynamics perspective, the predominant moieties of M–N–C SACs are expected to have the formula of MN₄ instead of MN_xC_4−x (x = 1–3). Indeed, the research of Li et al. indicated that no Ni–C bond exists in the Ni–N–C SACs they synthesized on the basis of the C K-edge X-ray absorption spectroscopy (XAS) characterization of the Ni–N–C SACs and a control N–C system.⁶ Moreover, numerous experimental and computational studies on the ORR occurring on M–N–C SACs have demonstrated that the dispersed metal sites are mainly contributed by MN₄.^18–22 Therefore, none of the commonly studied α-NiN_xC_4−x moieties are responsible for the excellent CO₂RR performance of Ni–N–C SACs. This situation presents an ambiguity regarding the structure–activity of the Ni sites of heterogeneous Ni–N–C SACs, which requires the development of (operando) characterization techniques and modelling approaches with high sensitivity and resolution.

In this study, we addressed this ambiguity by employing first-principles calculations and microkinetic modelling to investigate pyrrolic-type NiN_xC_4−x (β-NiN_xC_4−x) moieties as well as α-NiN_xC_4−x moieties, in tandem with a constant-potential method and the explicit treatment of solvents. β-NiN_xC_4−x moieties were explored because pyrrolic-type moieties have been demonstrated to inherently exist (or coexist with pyridinic-type moieties) in heterogeneous M–N–C SAC,^28–32 but little is known of their CO₂RR mechanisms. Remarkably, we demonstrated that the electrode potential is crucial for governing CO₂ activation, and in contrast to other moieties that are inactive and/or exhibit low-selectivity, the β-NiN₄ moiety not only possesses a favourable formation energy but also exhibits considerable activity and selectivity for CO₂-to-CO conversion.

Results and discussion

Structural models of Ni–N–C SACs and their formation energy

Good knowledge of the active structure of a catalyst is a solid foundation for understanding its catalytic nature.^33,34 It has been established that the local coordination environment of a single-atom metal site is essential for determining its catalytic performance.^35,36 We began by constructing the potential configurations of single-atom Ni sites and evaluating their formation energy values. For the local structure of the Ni sites of Ni–N–C SACs, the central Ni atom adopts a (quasi-)planar tetracoordination configuration, and as mentioned earlier, the ligand species can be either pyridinic N, pyrrolic N, or C. Interestingly, the prototypical α-NiN_xC_4−x can be conveniently constructed by creating a double-vacancy in the basal plane of graphene and replacing certain C atoms with N atoms (Fig. 1a), which may explain why α-NiN_xC_4−x moieties are typically chosen as structural models for XAS fitting and DFT calculations in existing literature.^6–15,23,24

Fig. 1 Twelve different N-coordinated single-atom Ni configurations: (a) the commonly studied pyridinic-type Ni moieties in literature, α-NiN_xC_4−x (x = 1–4); (b) the pyrrolic-type Ni moieties, β-NiN_xC_4−x. The grey, blue, orange, and white spheres represent C, N, Ni, H atoms, respectively. Calculated formation energy of (c) α-NiN_xC_4−x and (d) β-NiN_xC_4−x. The formation energy values of α-NiN₄ and β-NiN₄ are set as zero.

In contrast, the construction of the β-NiN_xC_4−x moieties requires more skilful manipulation of the C atoms to create pyrrolic-type coordination environments (Fig. 1b), and according to previous studies,^28–30 the unsaturated C atoms have been terminated by H atoms. Such pyrrolic-type moieties have been demonstrated to typically exist within the micropores of graphene nanosheets. α-NiN₂C₂ exhibits three possible configurations that correspond to different combinations of its N and C ligands; each of these configurations have been explored in various studies, such as α-NiN₂C₂-a,¹⁰ α-NiN₂C₂-b,^12,23 and α-NiN₂C₂-c.²⁴ Similarly to α-NiN₂C₂, β-NiN₂C₂ also corresponds to three possible configurations. In this work, only the configurations with a relative low formation energy were adopted in the later CO₂RR calculations for NiN₂C₂.

The relative formation energy of the α-NiN_xC_4−x and β-NiN_xC_4−x moieties are presented in Fig. 1c and d, in which the formation energy values of α-NiN₄ and β-NiN₄ are set to be zero. For both the α-NiN_xC_4−x and β-NiN_xC_4−x moieties, the formation energy increases significantly as the number of N atoms decreases. In particular, the difference in formation energy values between NiN₁C₃ and NiN₄ reaches up to 2.90 eV for the pyridinic-type configuration and 4.44 eV for the pyrrolic-type configuration. These results suggest that compared with N, the C ligand is not a good support for anchoring single-atom Ni species.

Please note that such formation energy calculations were performed to provide a general demonstration of the single-atom Ni sites based on a thermodynamics criterion, and they do not indicate whether the sites are pyridinic or pyrrolic. Additionally, they do not exclude the presence of energetically unfavourable NiN_xC_4−x (x = 1–3) moieties existing in nature. Indeed, the local structure of the single metal site should be determined by a synthesis parameter, such as the precursor materials, synthetic atmosphere, and pyrolysis temperature. However, the favourable formation energy of NiN₄ over that of NiN_xC_4−x (x = 1–3) provides clear evidence that in principle, the predominant Ni moieties of Ni–N–C SACs are expected to correspond to the formula of NiN₄.

CO₂ reduction reaction on the pyrrolic-type NiN₄

We next assessed the intrinsic CO₂RR activities of the Ni moieties by employing a constant-potential/hybrid-solvent model. First, the constant-potential method considers the surface charge effect to determine potential-dependent energy values. Please note that the surface charge of an electrocatalyst is essential for assessing its chemical reactivity.^{23,24,37–39} Secondly, the explicit solvent model of the electrochemical interface allows the interactions between intermediates and the surrounding water molecules to be captured. Thirdly, the implicit solvation model accurately describes long-range electrostatic interactions and takes advantage of the fact that the electrostatic potential goes to zero in the electrolyte region. More details can be found in the Methods section.

Taking β-NiN₄, a Ni moiety with a favourable formation energy, as an example, the geometric structures of its CO₂RR reaction species are shown in Fig. 2a. For each species, five individual simulations with system charge values of −1.75e, −1.5e, −1.25e, −1e, and −0.75e were conducted. Based on these calculations, the electrode potential (U) of each species for a given charge can be determined (Fig. 2b). A wider range of charge values (nine groups from −0.25e to −2.25e) has also been sampled, and negligible difference was found between the results (Fig. S1†). Due to the non-identical electronic structures of these species, the same charge does not correspond to the same value of U. However, the total energy of each species can be described as a continuous function of U by fitting five potential–energy points (Fig. 2c and Table S1†). With these potential–energy curves, the energy of these species for a specific value of U can be obtained, and the corresponding free energy changes ΔG of the reaction steps can be determined by further considering contributions of the zero-point energy, enthalpy, and entropy (Fig. 2d).

Fig. 2 (a) Views of the geometric structures of adsorbate-free β-NiN₄ slab (*), COOH*, and CO* optimized using a charge value of −1.25e. The blue dashed lines represent hydrogen bonds. (b) Calculated electrode potentials of the studied species at charge values of −1.75e, −1.5e, −1.25e, −1e, and −0.75e. (c) Total energy of the studied species as a function of the electrode potential. (d) Free energy changes ΔG of the reaction steps as a function of the electrode potential. (e) Total energy of physiosorbed and chemisorbed CO₂. (f) The Ni–C length of as a function of the electrode potential. The insets are geometric structures that were optimized using charge values of −1.5e (left) and −0.25e (right).

It can be observed from Fig. 2d that CO₂ adsorption on β-NiN₄ ( formation step) is exothermic in nature below U = −0.26 V and is more energetically favourable than the competing H* formation above U = −1.20 V. This suggests that β-NiN₄ can efficiently activate CO₂ molecules at typical working potentials while maintaining high selectivity for the CO₂RR compared with the HER over a wide range of U. The conversion of to COOH* (COOH* formation step) is the potential-limiting step, which becomes exothermic at U = −0.79 V. Subsequently, COOH* is reduced to CO* (CO* formation step) with a negative ΔG value, and the adsorbed CO molecule can be easily released from the β-NiN₄ surface. A more appropriate comparison of the onset potentials and overpotentials with experimental data can be established based on microkinetic modelling, as discussed later.

Specifically, we determined that the consideration of U is crucial for reproducing the experimentally observed chemisorption of CO₂ on the Ni sites of Ni–N–C SACs. As shown in Fig. 2e and f, the adsorption of CO₂ on β-NiN₄ preferentially adopts a chemisorption configuration when the U is lower than −0.27 V, while a physisorption configuration is dominant above this potential value. Chemisorption configurations can be distinguished from physisorption configurations by the formation of a Ni–C bond. For example, a CO₂ molecule is bent when adsorbed onto β-NiN₄ while forming a Ni–C bond (1.97 Å) at U = −0.59 V (q = −1.5e), but it remains straight with a Ni–C distance of 3.09 Å at U = 0.03 V (q = −0.25e), which is within the range of typical vdW interactions. Please note that the different adsorption configurations were obtained by employing different initial structures (i.e., different Ni–C distances). The surrounding water molecules are also essential for the chemisorption of CO₂ since the hydrogen bonds between them stabilize the intermediates.²³ The removal of the surrounding water molecules can also result in the adsorption of CO₂ transiting from chemisorption to physisorption.

Comparison of the CO₂RR performance of different Ni moieties

The constant-potential results of other Ni moieties are shown in Fig. 3a, b, S2 and S3†, and their energy–potential relationships are summarized in Tables S2–S8.† For α-NiN₄, another Ni moiety exhibiting a favourable formation energy, its and COOH* formation steps are only exothermic at very negative potentials (−0.81 and −1.10 V, respectively), as shown in Fig. 3a. This suggests that α-NiN₄ is inactive for facilitating the CO₂RR.

Fig. 3 (a and b) Free energy changes of the reaction steps of (a) α-NiN₄ and (b) α-NiN₁C₃ as a function of the electrode potential. (c) Total energy of H* on the Ni top site and Ni–C bridge site of α-NiN₁C₃ as a function of the electrode potential U. The insets present the atomic moieties of H* adsorbed on the Ni top site (upper) and Ni–C bridge site (bottom). (d and e) Views of the atomic structures of H* adsorbed onto (d) α-NiN₄ and (e) β-NiN₄. (f) Free energy changes of formation, COOH* formation, and H* formation on the Ni moieties at U = −0.60 V. For NiN₂C₂, the configurations with a relatively lower formation energy, i.e., α-NiN₂C₂-a and β-NiN₂C₂-b, were adopted for energy calculations. (g) Difference in free energy between H* formation and formation at U = −0.60 V. A more positive value corresponds to higher CO₂RR selectivity, and vice versa.

While the and COOH* formation steps of α-NiN₁C₃ are more energetically favourable than its H* formation step, the latter surpasses the formation step below U = −0.45 V (Fig. 3b). On α-NiN₁C₃, the location of H* with the lowest energy is the Ni–C bridge site instead of the Ni top site (Fig. 3c), which differs from both α-NiN₄ and β-NiN₄ that preferentially adopt the Ni-top H* configuration (Fig. 3d and e). Although it appears that H* could be weakly adsorbed onto the Ni top sites of α-NiN₁C₃,²³ the strong binding of H* on the Ni–C bridge sites prevents the reduction of CO₂, thereby resulting in poor CO₂RR selectivity. Meanwhile, strong H* adsorption on the bridge sites has also been identified in other Ni moieties that involve the hybrid C and N coordination (Fig. S4†).

Fig. 3f shows the values of ΔG for the three key reaction steps, i.e., formation, COOH* formation, and H* formation, of the studied Ni moieties at U = −0.60 V, at which a partial current density of CO can be detected experimentally for Ni–N–C SACs.^10–15 The results for the CO* formation and CO desorption steps are excluded since the two steps are typically exothermic. At U = −0.60 V, formation on α-NiN₄ remains endothermic by 0.20 eV, and COOH* formation is endothermic by 0.41 eV. These results indicate that α-NiN₄ displays a low activity for facilitating the CO₂RR, which is consistent with previous studies.

Although the formation steps of the α- and β-NiN_xC_4−x (x = 1–3) moieties correspond to negative ΔG values, their H* formation steps are highly exothermic. Since H* formation is more energetically favourable than formation, the α- and β-NiN_xC_4−x (x = 1–3) moieties are expected to facilitate the HER preferentially instead of the CO₂RR at U = −0.60 V. Evidently, such a situation directly contradicts the experimental observation of Ni–N–C SACs exhibiting high selectivity for the CO₂RR. Hence, apart from α-NiN₄, Ni moieties with hybrid C and N coordination environments, i.e., the α- and β-NiN_xC_4−x (x = 1–3) moieties, cannot be responsible for the excellent CO₂RR performance of Ni–N–C SACs.

Remarkably, β-NiN₄ exhibits a negative ΔG value for formation (−0.19 eV) and a positive ΔG value for H* formation (0.40 eV), as shown in Fig. 3f, which indicates that at U = −0.60 V, CO₂ molecules can be activated effectively over β-NiN₄, and the competing HER process is suppressed. Please note that the difference in energy between formation and H* formation is as large as 0.59 eV (Fig. 3g); such a value implies that β-NiN₄ possesses high selectivity toward the CO₂RR. Moreover, for COOH* formation, it can be shown from Fig. 3f that aside from α-NiN₁C₃ (−0.04 eV), the other Ni moieties exhibit positive ΔG values at U = −0.60 V; however, the ΔG values of α-NiN₂C₂, β-NiN₄, and β-NiN₁C₃ are only slightly positive, indicating that COOH* formation on these three moieties is still feasible. Overall, the above results have established that among these Ni moieties, only β-NiN₄ has a favourable formation energy value and is active and selective toward the CO₂RR. Such characteristics of β-NiN₄ may potentially account for the experimentally observed high CO₂RR performance of Ni–N–C SACs.

Kinetics and free energy diagrams

To obtain more insight into the CO₂RR occurring on the active β-NiN₄ moiety, we further explored the kinetics of the key reaction steps involved, and the thermodynamic and kinetic energy values were compiled in a free energy diagram (Fig. 4a). The potential dependence of the kinetic barrier ΔG^‡(U) was determined through multiple barrier calculations at varying potentials (Fig. 4b and c). According to our calculations, at U = −0.60 V, β-NiN₄ exhibits a small barrier of 0.05 eV for formation and a moderate barrier of 0.69 eV for COOH* formation, while it displays a quite high barrier of 1.31 eV for H* formation (Volmer reaction). After the formation of the COOH* or H* intermediates, the subsequent protonation steps are kinetically feasible, e.g., ΔG^‡(−0.6 V) = 0.42 eV for the reduction of COOH* to CO* and ΔG^‡(−0.6 V) = 0.10 eV for the reduction of H* to H₂ (Heyrovsky reaction). By extrapolating these thermodynamic and kinetic energy values, we thereby conclude that at U = −0.60 V, the efficient conversion of CO₂ to CO on β-NiN₄ can be achieved.

Fig. 4 (a) Free energy diagram of the CO₂RR and competing HER on β-NiN₄ at U = −0.60 V. The insets are snapshots of the transition states (TS). The arrow symbols indicate the directions of the HER and CO₂RR processes. (b and c) Calculated kinetic barriers ΔG^‡ for the (b) CO₂RR and (c) HER on β-NiN₄ as a function of the electrode potential U.

The kinetic barriers of a prototypical Ni moiety, α-NiN₁C₃, were also investigated. The geometric structures of its transition states are shown in Fig. S5.† As expected, at U = −0.60 V, α-NiN₁C₃ was predicted to exhibit a relatively favourable ΔG^‡ value for COOH* formation (0.39 eV) compared with β-NiN₄. However, upon considering the lowest-energy location of the H* intermediate (i.e., H* adsorption on the Ni–C bridge sites), the value of ΔG^‡ for H* formation (0.36 eV) was found to be lower than that for COOH* formation. Therefore, the kinetic results have also established the poor CO₂RR selectivity of α-NiN₁C₃, which has already been identified by the above thermodynamic calculations.

Microkinetic modelling

Finally, we constructed a microkinetic model to further examine the CO₂RR performance of the active β-NiN₄ moiety and to provide a qualitative comparison with existing experimental data of Ni–N–C SACs. Fig. 5a shows the reaction steps of the CO₂RR (i = 1–5) and HER (i = 6, 7) processes considered in our model. The terms CO₂(aq) and CO₂(dl) represent CO₂ molecules within the electrolyte and at the catalyst–electrolyte interface, respectively. The term k_i is the rate constant of the step i, which depends directly on the potential-based kinetic barrier ΔG^‡(U) in accordance with the transition state theory, and k_−i is the rate constant for the reverse reaction. The CO partial current density can be estimated by solving these rate equations at a steady state. Further details are summarized in the Methods section.

Fig. 5 (a) Schematic plot of reaction steps considered in the microkinetic modelling. (b) Simulated partial CO current density of β-NiN₄ compared with existing experimental data (catalyst loading: 1.0 mg cm⁻² in ref. 11 and 0.75 mg cm⁻² in ref. 12). For the polarization curve simulation, eρ_Ni = 10.5 μC cm⁻² was adopted; the term ρ_Ni is the density of active single-Ni-atom sites. (c) Tafel plot of β-NiN₄.

Fig. 5b displays the model partial CO current density (j_CO) as a function of the electrode potential. The onset-potential of β-NiN₄ was calculated to be −0.55 V, beyond which the value of j_CO increases continuously with U and reaches 10 mA cm⁻² at U = −0.77 V. Such a curve is consistent with its experimentally measured counterparts.^11,12 Please note that the current density is highly dependent on the catalyst morphology and electrochemical measurement parameters (e.g., electrode material, catalyst mass loading, and CO₂ diffusion rate), and the present calculation was performed to provide a qualitative demonstration of CO₂RR activity as the effective density of active Ni sites (ρ_Ni) is unknown experimentally, even though the Ni content and loading can be determined.

A Tafel analysis was further conducted (Fig. 5c). According to our calculations, the Tafel slope of β-NiN₄ was estimated to be ∼101 mV dec⁻¹, which is in reasonable agreement with the typical values for Ni–N–C SACs in experiments (∼108 mV dec⁻¹).^6,8,10,11 The polarization curve for the competing HER is shown in Fig. S6,† in which the poor HER activity of β-NiN₄ is revealed (U = −1.57 V at 10 mA cm⁻²). This also implies that the experimentally observed concomitant H₂ production at U < −0.7 V may be attributed to the HER occurring on the N-neighbouring C sites of Ni–N–C SACs, considering that the N-doped graphene nanosheets are active for facilitating the HER.¹¹

Conclusions

In summary, we have presented a fundamental understanding of the CO₂RR occurring on heterogeneous Ni–N–C SACs by considering twelve possible single-atom Ni moieties. We discovered that an intrinsic tetrapyrrolic coordination environment with respect to Ni, i.e., β-NiN₄, can boost the CO₂RR activity, which rationalizes existing experimental observations. Such an understanding is based on the first-principles and microkinetic modelling of reaction processes at a solid–water interface, in which a constant-potential/hybrid-solvent model was adopted. According to an analysis of the formation energy, the predominant Ni moieties in Ni–N–C SACs are expected to correspond to the formula of NiN₄. After determining the potential-dependent thermodynamic and kinetic energy values of the studied Ni moieties, we found that among them, only β-NiN₄ exhibits considerable activity for the selective CO₂-to-CO conversion compared with the competing H₂ production. Specifically, its theoretical polarization curve and Tafel plot show reasonable agreement with existing experimental data. Our work discloses the activity origin of an important model SAC, offers practical guidance for the rational design of electrocatalysts, and serves as a prime example for future studies exploring the underlying mechanisms of other heterogeneous electrocatalysis in aqueous solutions.

Methods

First-principles computations

The first-principles calculations were performed with Vienna Ab initio Simulation Package (VASP),^40,41 and the projector-augmented wave approach was employed to describe the ion–electron interactions.⁴² The electronic exchange–correlation energy was modelled using the BEEF-vdW functional,⁴³ which accurately describes vdW interactions while retaining a good description of chemisorption of adsorbates on surfaces. All calculations included spin polarization and were performed using a cutoff energy of 400 eV with a Gaussian smearing of 0.1 eV for slabs and 0.01 eV for gas-phase species. The k-space samplings were set as 3 × 3 × 1. The optimized lattice parameters of the α-NiN_xC_4−x slabs were a = 8.39 Å and b = 12.43 Å, and those of the β-NiN_xC_4−x slabs were a = 8.12 Å and b = 12.02 Å. For all slabs, the distance of the vacuum space was set as ∼25 Å to minimize interactions between periodic images. The lowest energy location of each intermediate was reported.

Formation energy calculations

The formation energy of a NiN_xC_4−x (x = 1, 2, 3, or 4) moiety is defined as:

E_form = E_NiNxC4−x − N_Cμ_C − N_Nμ_N − N_Hμ_H − μ_Ni (1)

where E_NiNxC4−x is the total energy of a NiN_xC_4−x slab, and N_C, N_N, and N_H are the number of C, N, and H atoms in the NiN_xC_4−x slab, respectively. The term μ denotes the chemical potential of a species, and herein μ_C, μ_N, μ_H, and μ_Ni were taken as the total energy of one C atom of pristine graphene, one N atom of a N₂ molecule, one H atom of a H₂ molecule, and one Ni atom of bulk Ni, respectively. According to this definition, a NiN_xC_4−x moiety with a lower formation energy corresponds to a higher likelihood of being present in Ni–N–C SACs.

Constant-potential/hybrid-solvent model

A constant-potential/hybrid-solvent model was employed to determine the potential-dependent thermodynamic and kinetic energy values of reactions at the catalyst–solvent interface. A water bilayer containing 24 H₂O molecules was constructed to explicitly describe the interactions between intermediates and the surrounding water molecules. The Poisson–Boltzmann implicit solvation model was utilized to establish the relationship between the system charge and electrode potential U, as implemented in the VASPsol.^44,45 This model takes advantage of the fact that the electrostatic potential goes to zero in the electrolyte region. The effective surface tension parameter was set as zero as the cavitation energy is negligible.⁴⁵ For an electrochemical system, the electrostatic potential is analogous to the workfunction. We added additional electrons (a net charge) and subsequently calculated the workfunction of the charged slab (Φ_slab), which is reflected by the Fermi level with respect to the electrostatic potential in the electrolyte region. After determining the value of Φ_slab, the electrode potential of the slab on the standard hydrogen electrode (SHE) scale, U_SHE, can be calculated by

U_SHE = Φ_slab/e − U⁰ (2)

where U⁰ is the absolute potential of the SHE. The value of U⁰ was predicted to be 4.6 V from the VASPsol simulation,⁴⁵ which is in good agreement with experimental estimates (4.4–4.8 V).⁴⁶ To facilitate comparisons between our results and existing experimental data, the values of U_SHE were regulated with respect to the reversible hydrogen electrode (RHE) scale by a potential shift of k_BT × pH × ln10. In this work, the pH was taken as 7, as the measurements of the CO₂RR on Ni–N–C SACs are typically preformed in a near neutral environment (e.g., CO₂-saturated 0.5 M KHCO₃). The total energy of a charged slab was further determined by considering energy corrections regarding the background charge and the electrostatic potential shift of the charged slab with respect to an uncharged slab.^47–52

Potential-dependent kinetic barrier

For adsorbates, the zero-point energy (E_ZPE), enthalpy (∫C_pdT), and entropy (−TS) contributions to free energy were calculated by vibrational frequencies calculations, in which all 3N degrees of freedom are treated as harmonic vibrational motions. For molecules, these contributions were obtained from the NIST database. Electronic energy corrections were applied to the calculations of E_DFT of CO₂ (+0.33 eV) and H₂ (+0.09 eV) molecules because of the errors residing in the description of the OCO backbone and H₂ from the BEEF-vdW functional.⁵³ Based on the above constant potential simulations and free energy contributions, the free energy change (ΔG) for each reaction step can be determined by the following equation:

ΔG = ΔE + ΔE_ZPE + Δ∫C_pdT − TΔS (3)

where ΔE, ΔE_ZPE, Δ∫C_pdT, and TΔS are the electronic energy difference, zero point energy difference, enthalpy difference and entropy difference, respectively.

The potential dependence of the kinetic barrier, ΔG^‡(U), was obtained via multiple barrier calculations at varying potentials (different system charges). The transition state was determined using the climbing-image nudged elastic band (CI-NEB) approach⁵⁴ with a given charge and has been further verified by vibrational frequency calculations (only one imaginary frequency). The values of ΔG^‡ have been corrected by considering the zero-point energy, but no configurational entropies were included according to the transition state theory.

Microkinetic model

The rate equation of each species in the CO₂RR process is summarized in eqn (4)–(9), and they are presented as follows:where θ and x are the coverage and mole fraction of a species, respectively. In this work, x_CO2(aq) is set to be 5.79 × 10⁻⁴, corresponding to CO₂ gas in equilibrium with aqueous CO₂; p_CO denotes the effective CO pressure, which is taken to be 1 bar and has been identified to exert a negligible influence on the theoretical polarization curve of β-NiN₄ due to its weak CO binding (Fig. S7†). According to the transition state theory, the rate constant k_i is calculated bywhere A′ is an effective prefactor. For CO₂ adsorption and CO desorption steps, A′ is taken as 1 × 10⁸ s⁻¹; for electrochemical steps, A′ is set to The current density is calculated by

j = 2eρ_NiTOF (10)

where ρ_Ni is the density of active Ni sites. In this simulation, the concentration of the active Ni sites in the Ni–N–C catalysts was assumed to be ∼1.5 at%, which corresponds to ρ_Ni = ∼6.6 × 10¹³ cm⁻² based on a graphene slab model. The value of the turnover frequency (TOF) was determined by numerically solving the above rate equations of reaction species at a steady state. The effective CO₂ species that can participate in the CO₂RR over the catalyst surface has been treated as CO₂ diffusing from the bulk electrolyte CO₂(aq) to the catalyst–electrolyte interface CO₂(dl); the reaction rate k₁ was set as 1 × 10⁶ s⁻¹, and the resultant polarization curve well matches existing experimental data. Overall, we expected to give a qualitative description for the CO₂RR performance instead of a quantitative demonstration.

Author contributions

Y. W. performed the calculations. Y. W., L. Y., and K. Z. analysed the data and wrote the paper. K. Z. supervised the study.

Conflicts of interest

There are no conflicts to declare.

Acknowledgements

The authors thank Yu-Jia Tang for many fruitful discussions. The authors acknowledge financial support from the Nanyang Environment and Water Research Institute (Core Fund), Nanyang Technological University, Singapore. The calculations were performed on computer clusters at High Performance Computing Centre at Nanyang Technological University.

References

1. Y. Y. Birdja, E. Pérez-Gallent, M. C. Figueiredo, A. J. Göttle, F. Calle-Vallejo and M. T. Koper, Nat. Energy, 2019, 4, 732–745.
2. S. Zhang, Q. Fan, R. Xia and T. J. Meyer, Acc. Chem. Res., 2020, 53, 255–264.
3. M. Li, H. Wang, W. Luo, P. C. Sherrell, J. Chen and J. Yang, Adv. Mater., 2020, 32, 2001848.
4. A. S. Varela, N. Ranjbar Sahraie, J. Steinberg, W. Ju, H.-S. Oh and P. Strasser, Angew. Chem., Int. Ed., 2015, 54, 10758–10762.
5. J. Gu, C.-S. Hsu, L. Bai, H. M. Chen and X. Hu, Science, 2019, 364, 1091–1094.
6. X. Li, W. Bi, M. Chen, Y. Sun, H. Ju, W. Yan, J. Zhu, X. Wu, W. Chu and C. Wu, J. Am. Chem. Soc., 2017, 139, 14889–14892.
7. W. Ju, A. Bagger, G.-P. Hao, A. S. Varela, I. Sinev, V. Bon, B. R. Cuenya, S. Kaskel, J. Rossmeisl and P. Strasser, Nat. Commun., 2017, 8, 944.
8. H. B. Yang, S.-F. Hung, S. Liu, K. Yuan, S. Miao, L. Zhang, X. Huang, H.-Y. Wang, W. Cai and R. Chen, Nat. Energy, 2018, 3, 140–147.
9. Q. Qu, S. Ji, Y. Chen, D. Wang and Y. Li, Chem. Sci., 2021, 12, 4201–4215.
10. C. Zhao, Y. Wang, Z. Li, W. Chen, Q. Xu, D. He, D. Xi, Q. Zhang, T. Yuan and Y. Qu, Joule, 2019, 3, 584–594.
11. K. Jiang, S. Siahrostami, T. Zheng, Y. Hu, S. Hwang, E. Stavitski, Y. Peng, J. Dynes, M. Gangisetty and D. Su, Energy Environ. Sci., 2018, 11, 893–903.
12. T. Möller, W. Ju, A. Bagger, X. Wang, F. Luo, T. N. Thanh, A. S. Varela, J. Rossmeisl and P. Strasser, Energy Environ. Sci., 2019, 12, 640–647.
13. Y. N. Gong, L. Jiao, Y. Qian, C. Y. Pan, L. Zheng, X. Cai, B. Liu, S. H. Yu and H. L. Jiang, Angew. Chem., Int. Ed., 2020, 59, 2705–2709.
14. H. Yang, Q. Lin, C. Zhang, X. Yu, Z. Cheng, G. Li, Q. Hu, X. Ren, Q. Zhang, J. Liu and C. He, Nat. Commun., 2020, 11, 593.
15. Z. Chen, X. Zhang, W. Liu, M. Jiao, K. Mou, X. Zhang and L. Liu, Energy Environ. Sci., 2021, 14, 2349–2356.
16. Y. Li and A. I. Frenkel, Acc. Chem. Res., 2021, 8738–8748.
17. T. N. Nguyen, M. Salehi, Q. V. Le, A. Seifitokaldani and C. T. Dinh, ACS Catal., 2020, 10, 10068–10095.
18. H. T. Chung, D. A. Cullen, D. Higgins, B. T. Sneed, E. F. Holby, K. L. More and P. Zelenay, Science, 2017, 357, 479–484.
19. J. Li, M. Chen, D. A. Cullen, S. Hwang, M. Wang, B. Li, K. Liu, S. Karakalos, M. Lucero and H. Zhang, Nat. Catal., 2018, 1, 935–945.
20. S. Kattel and G. Wang, J. Phys. Chem. Lett., 2014, 5, 452–456.
21. J. Sun, Y.-H. Fang and Z.-P. Liu, Phys. Chem. Chem. Phys., 2014, 16, 13733–13740.
22. Y. Wang, Y.-J. Tang and K. Zhou, J. Am. Chem. Soc., 2019, 141, 14115–14119.
23. X. Zhao and Y. Liu, J. Am. Chem. Soc., 2020, 142, 5773–5777.
24. M. D. Hossain, Y. Huang, H. Y. Ted, W. A. Goddard III and Z. Luo, Nat. Commun., 2020, 11, 1–14.
25. J. M. Ziegelbauer, T. S. Olson, S. Pylypenko, F. Alamgir, C. Jaye, P. Atanassov and S. Mukerjee, J. Phys. Chem. C, 2008, 112, 8839–8849.
26. S. Kattel, P. Atanassov and B. Kiefer, J. Phys. Chem. C, 2012, 116, 8161–8166.
27. Y. Li, X. Liu, L. Zheng, J. Shang, X. Wan, R. Hu, X. Guo, S. Hong and J. Shui, J. Mater. Chem. A, 2019, 7, 26147–26153.
28. A. Zitolo, V. Goellner, V. Armel, M.-T. Sougrati, T. Mineva, L. Stievano, E. Fonda and F. Jaouen, Nat. Mater., 2015, 14, 937–942.
29. A. Zitolo, N. Ranjbar-Sahraie, T. Mineva, J. Li, Q. Jia, S. Stamatin, G. F. Harrington, S. M. Lyth, P. Krtil and S. Mukerjee, Nat. Commun., 2017, 8, 1–11.
30. J. Li, M. T. Sougrati, A. Zitolo, J. M. Ablett, I. C. Oğuz, T. Mineva, I. Matanovic, P. Atanassov, Y. Huang and I. Zenyuk, Nat. Catal., 2021, 4, 10–19.
31. T. Mineva, I. Matanovic, P. Atanassov, M.-T. Sougrati, L. Stievano, M. Clémancey, A. Kochem, J.-M. Latour and F. Jaouen, ACS Catal., 2019, 9, 9359–9371.
32. S. Wagner, H. Auerbach, C. E. Tait, I. Martinaiou, S. C. Kumar, C. Kübel, I. Sergeev, H. C. Wille, J. Behrends and J. A. Wolny, Angew. Chem., Int. Ed., 2019, 58, 10486–10492.
33. Q. Wang, C. Cai, M. Dai, J. Fu, X. Zhang, H. Li, H. Zhang, K. Chen, Y. Lin, H. Li, J. Hu, M. Miyauchi and M. Liu, Small Sci., 2021, 1, 2000028.
34. M. Li, Y. Ma, J. Chen, R. Lawrence, W. Luo, M. Sacchi, W. Jiang and J. Yang, Angew. Chem., Int. Ed., 2021, 60, 11487–11493.
35. A. Chen, X. Zhang, L. Chen, S. Yao and Z. Zhou, J. Phys. Chem. C, 2020, 124, 22471–22478.
36. D. Gao, T. Liu, G. Wang and X. Bao, ACS Energy Lett., 2021, 6, 713–727.
37. Y. Huang, R. J. Nielsen and W. A. Goddard III, J. Am. Chem. Soc., 2018, 140, 16773–16782.
38. D. Kim, J. Shi and Y. Liu, J. Am. Chem. Soc., 2018, 140, 9127–9131.
39. C. Choi, G. H. Gu, J. Noh, H. S. Park and Y. Jung, Nat. Commun., 2021, 12, 4353.
40. G. Kresse and J. Hafner, Phys. Rev. B: Condens. Matter Mater. Phys., 1993, 47, 558.
41. G. Kresse and J. Furthmüller, Phys. Rev. B: Condens. Matter Mater. Phys., 1996, 54, 11169.
42. P. E. Blöchl, Phys. Rev. B: Condens. Matter Mater. Phys., 1994, 50, 17953.
43. J. Wellendorff, K. T. Lundgaard, A. Møgelhøj, V. Petzold, D. D. Landis, J. K. Nørskov, T. Bligaard and K. W. Jacobsen, Phys. Rev. B: Condens. Matter Mater. Phys., 2012, 85, 235149.
44. K. Mathew, R. Sundararaman, K. Letchworth-Weaver, T. Arias and R. G. Hennig, J. Chem. Phys., 2014, 140, 084106.
45. K. Mathew, V. C. Kolluru, S. Mula, S. N. Steinmann and R. G. Hennig, J. Chem. Phys., 2019, 151, 234101.
46. H. Reiss and A. Heller, J. Phys. Chem., 1985, 89, 4207–4213.
47. C. D. Taylor, S. A. Wasileski, J.-S. Filhol and M. Neurock, Phys. Rev. B: Condens. Matter Mater. Phys., 2006, 73, 165402.
48. J.-S. Filhol and M. Neurock, Angew. Chem., Int. Ed., 2006, 45, 402–406.
49. K. Chan and J. K. Nørskov, J. Phys. Chem. Lett., 2015, 6, 2663–2668.
50. N. G. Hörmann, O. Andreussi and N. Marzari, J. Chem. Phys., 2019, 150, 041730.
51. Z. Duan and G. Henkelman, ACS Catal., 2019, 9, 5567–5573.
52. C. Guo, X. Fu and J. Xiao, J. Phys. Chem. C, 2020, 124, 25796–25804.
53. F. Studt, F. Abild-Pedersen, J. B. Varley and J. K. Nørskov, Catal. Lett., 2013, 143, 71–73.
54. G. Henkelman, B. P. Uberuaga and H. Jónsson, J. Chem. Phys., 2000, 113, 9901–9904.

---

Footnotes

† Electronic supplementary information (ESI) available: Computational methods, supporting figures, and supporting tables. See DOI: 10.1039/d1sc04094d

‡ Current address: Jiangsu Collaborative Innovation Centre of Biomedical Functional Materials, Jiangsu Key Laboratory of New Power Batteries, School of Chemistry and Materials Science, Nanjing Normal University, 1 Wenyuan Road, Nanjing 210023, China.

---