Mechanistic insights and predictive screening of M@C₂N catalysts for urea electrosynthesis from N₂ and CO₂

Abstract

Electrocatalytic urea synthesis via the co-reduction of N₂ and CO₂ under ambient conditions offers a sustainable alternative to energy-intensive industrial processes. However, this process is hindered by several challenges, including the inertness of the N≡N bonds, sluggish C–N coupling kinetics, competing side reactions, and the lack of predictive models to guide catalyst development. In this study, we conduct a comprehensive density functional theory (DFT) screening of 26 transition metal single atoms anchored on graphitic C₂N (M@C₂N) to identify active and selective electrocatalysts for urea synthesis. Four mechanistic pathways, CO₂, OCOH, CO, and NCON, are systematically explored, revealing that the initial and final protonation steps of adsorbed N₂ are critical in determining catalytic performance. Among the candidates, Nb@C₂N, Mo@C₂N, and Re@C₂N exhibit the most favorable activity, achieving low limiting potentials of −0.50, −0.51, and −0.51 V, respectively. To accelerate catalyst discovery, we introduce a physically grounded descriptor, Φ, based on the d-electron count and electronegativity of the anchored metals, which accurately captures structure–activity relationships and enables rapid screening across materials. Our results establish a mechanistic framework and a descriptor-driven strategy for the rational design of single-atom electrocatalysts for ambient urea synthesis.

---

1. Introduction

The growing need for sustainable technologies has brought increasing attention to the electrocatalytic conversion of CO₂ and N₂ into valuable-added chemicals such as urea, a nitrogen-rich fertilizer.^1–3 However, its conventional synthesis remains highly energy-intensive, relying on the Haber–Bosch process for ammonia (NH₃) production, followed by CO₂ coupling at 150–200 °C and 150–250 bar.^4,5 Although N₂ is abundant, its activation remains a significant challenge due to the high dissociation energy of the N≡N triple bond (940.95 kJ mol⁻¹).^6,7 Approximately 80% of global ammonia is used for urea synthesis,⁸ and the Haber–Bosch process alone contributes nearly 2% of global CO₂ emissions, highlighting the urgent need for low-energy, ambient condition alternatives.

Electrochemical co-reduction of CO₂ and N₂ under mild conditions presents a more sustainable route. Catalysts such as PdCu nanoparticles on TiO₂ and Te-doped Pd have been shown to promote C–N coupling via a tower-like *NCON intermediate.^9,10 Heterostructures such as Bi–BiVO₄ and frustrated Lewis pair (FLP)-active InOOH materials have been developed to co-activate N₂ and CO₂ through nucleophilic and electrophilic interactions.^11,12 However, these strategies are often highly system-dependent, and the mechanistic landscape remains poorly understood, limiting rational catalyst design. The limited availability of efficient catalysts is partly due to the lack of a comprehensive understanding of reaction mechanisms and the absence of practical theoretical principles to guide catalyst discovery. Density functional theory (DFT) has proven valuable in understanding and exploring materials such as MBenes,¹³ CuB₁₂ monolayers,¹⁴ single-atom catalysts (SACs) on divacancy boron nitride (M/p-BN),¹⁵ and double-atom catalysts (DACs).^16–18 These studies highlight promising directions for catalyst development but often lack generalizable design principles to enable high-throughput screening.

Ideal electrocatalysts must activate both N₂ and CO₂ while favoring C–N bond formation over competing pathways.^19,20 SACs are of particular interest due to their tunable coordination and high metal utilization efficiency.^15,21–24 Our previous work on dual metal–boron doped g-C₂N introduced a “dual electron donating strategy” that enhanced urea synthesis. Scheme S1† provides a proposed route for synthesizing Fe–C₂N (or Nb@C₂N, Mo@C₂N, and Re@C₂N). Additionally, Co–C₂N materials have shown promise in lithium–sulfur batteries, demonstrating the broader versatility of C₂N as a support. However, the large design space of potential metal–support combinations and the lack of mechanistic clarity still hinder rational screening of SACs on g-C₂N.

In contrast to well-studied reactions such as the HER, oxygen evolution/reduction reaction (OER/ORR), and CO₂ reduction reaction (CO₂RR), where predictive descriptors such as p/d-band centers,^25–28 spin states,^29,30 charge transfer,³¹ and site densities³² enable efficient screening, urea synthesis lacks such intrinsic, physically grounded descriptors. As a result, computational discovery remains slow and system-specific. To address this gap, we introduce a descriptor-based screening strategy combining mechanistic DFT insights with a simple yet effective descriptor, Φ, which captures the d-electron count and metal–support electronegativity (Scheme 1a). We systematically evaluate 26 transition metals (3d, 4d, and 5d) anchored on g-C₂N and investigate four plausible C–N coupling pathways: CO₂, OCOH, CO, and NCON. Using a Four-Step Screening Strategy, we assess thermodynamic stability, adsorption characteristics, and establish structure–activity relationships using Φ and two additional descriptors, φ, and ε_d.

Scheme 1 (a) Workflow for descriptor-based screening of urea electrocatalysts. (b) Schematic diagram of four mechanistic pathways: CO₂, OCOH, CO, and NCON routes.

Our results identify Nb@C₂N, Mo@C₂N, and Re@C₂N as top-performing catalysts, outperforming benchmark PdCu.³³ A proposed synthetic route for the three catalysts is provided in Scheme S1 of the ESI.† Notably, Φ demonstrates excellent correlation with catalytic performance. This work not only advances mechanistic understanding of electrochemical urea synthesis but also introduces a practical descriptor framework to accelerate the discovery of high-performance SACs for sustainable urea production under ambient conditions.

2. Computational details

All calculations were performed using density functional theory (DFT) as implemented in the Vienna ab initio simulation package (VASP).^34–36 The generalized gradient approximation (GGA)³⁷ in the Perdew–Burke–Ernzerhof (PBE) scheme was employed to describe exchange–correlation interactions, while the core and valence electrons were treated by the projector augmented wave (PAW) method.³⁸ The cutoff energy for plane wave expansions was set to 450 eV. To model a two-dimensional system, 15 Å of the vacuum layer was imposed to avoid the interlayer interaction from the periodic images, while Grimme's DFT-D3 (ref. 39) method was used to treat van der Waals interactions during the adsorption process. The reciprocal space was sampled with 3 × 3 × 1 k-points for structural optimization and 5 × 5 × 1 for the electronic structure analysis to guarantee convergence. The convergence threshold was set to 1 × 10⁻⁵ eV for electronic optimization and 0.02 eV Å⁻¹ for ionic optimization. Vibrational frequency calculations were carried out under the harmonic approximation, with the surface atoms fixed and only the adsorbate vibrational modes considered. To investigate the electronic properties, Bader charge analysis was performed using the algorithm developed by the Henkelman group,⁴⁰ and integrated crystal orbital Hamilton population (ICOHP) values were obtained using the LOBSTER package^41,42 based on VASP wavefunctions.

The thermodynamic and electrochemical stabilities were assessed through formation energy (E_f) and dissolution potential (U_diss), with the computational procedures summarized in Table S1.† The formation energy was calculated using:

E_f = (E_M@C2N −E_C2N − E_M) (1)

where E_M@C2N, E_C2N and E_M are the total energies of the doped systems, pristine substrate, and bulk metal, respectively. The dissolution potential⁴³ was determined using:where is the standard dissolution potential of bulk metal, n is the number of electrons involved in the dissolution, and e is the elementary charge. A catalyst with E_f < 0 and U_diss > 0 is considered thermodynamically and electrochemically stable.^44,45

Gibbs free energies (G) were computed using the computational hydrogen electrode (CHE) method proposed by Nørskov et al.,⁴⁶ which approximates the chemical potential of a proton–electron pair by half the free energy of a hydrogen molecule at U = 0 V vs. the standard hydrogen electrode (SHE) using the following equation:

G = E + E_ZPE − TS + ∫C_pdT (3)

where E, E_ZPE, T, S, and ∫C_pdT are the total electronic energy, the zero-point energy, the temperature (set to 298.15 K), the entropy, and the enthalpic temperature correction, respectively. The Gibbs free energy change (ΔG) between two intermediates is calculated using:

ΔG = ΔE + ΔE_ZPE − TΔS + ∫C_pdT + ΔG_pH + ΔG_U (4)

where ΔE is the electronic energy difference between the two intermediate states, ΔE_ZPE, TΔS, and ∫C_pdT refer to the change in zero-point energies, the change in entropy at temperature T, and the enthalpic temperature correction. In this work, the temperature was set to 298.15 K. ΔG_pH, and ΔG_U represent the free energy correction corresponding to the pH of the solution and the electrode potential of the electrochemical reaction, respectively. The ΔG_pH can be estimated using the following equation.

ΔG_pH = 2.303k_BTpH, (5)

where k_B is the Boltzmann constant. In our calculation, we assumed acidic conditions (pH = 0), so ΔG_pH = 0. Moreover, ΔG_U = −neU, where n, e, and U are the number of transferred electrons, the electron charge, and the applied potential to the reaction, respectively. When there is no applied potential in the reaction, ΔG_U = 0.

In addition, the limiting potential (U_L) can be obtained to evaluate the activity of a catalyst, following the equation:

U_L = −ΔG_max/|e| (6)

where ΔG_max is the free energy change of the potential-determining step (PDS), defined as the step with the highest ΔG along the reaction pathway.⁴⁷

3. Results and discussion

3.1 Proposed mechanism of electrocatalytic urea synthesis

To validate the proposed reaction mechanisms, we conducted systematic DFT investigations on porous g-C₂N anchored with transition metal atoms (M@C₂N). These calculations provided detailed mechanistic insights and supported the development of a structure–activity descriptor, Φ, capable of predicting the catalytic performance of M@C₂N systems. The overall screening workflow based on this descriptor is outlined in Scheme 1a and comprises three steps: (i) data preparation: we first evaluated the electrochemical and thermodynamic stability of each M@C₂N system, including formation energies and dissolution potentials. Adsorption energies for *N₂ were also calculated to confirm reactant binding affinity. (ii) Screening criteria: catalysts exhibiting favorable thermodynamics for the first protonation of N₂ (ΔG < 0.78 eV) were shortlisted for full pathway analysis. This threshold was adopted from benchmark studies (e.g., PdCu) to ensure practical activity under ambient conditions. (iii) Descriptor construction process: intrinsic catalyst properties were used to construct the Φ descriptor, which was then fitted to the computed U_L to establish a predictive structure–activity relationship. This approach enables efficient pre-screening of M@C₂N catalysts, accelerating the discovery of high-performing candidates for urea electrosynthesis.

In this study, we propose a comprehensive mechanism for urea formation that enables the formation of desirable C–N bonds without complete cleavage of the inert N≡N bond. Previous studies have identified the tower-like *NCON intermediate as a key species for urea formation.^9,15,48,49 However, the direct insertion of *CO into molecular *N₂ is kinetically difficult due to N₂'s high stability. To address this challenge, we explored multiple mechanistic pathways for C–N bond formation from *N₂ and *CO₂-derived species.

As outlined in Scheme 1b, the full mechanism is complex and strongly dependent on the catalyst surface. A widely referenced pathway involves six sequential protonation steps. Starting from CO₂ adsorption, two protonation steps yield a *CO intermediate, which subsequently couples with *N₂ to form *NCON (tower-like intermediate), a key intermediate for urea production. This pathway, referred to as the NCON pathway (purple), continues with protonation steps that lead to final urea formation. In addition to the NCON pathway, alternative C–N bond formation pathways have been proposed, including coupling of *CO, *COOH, or *CO₂ with *N₂H_x intermediates. For instance, *CO–*NH₂ coupling was a favorable step on Te-doped Pd catalysts.¹⁰ More recently, Zhu et al. proposed three additional pathways (Scheme 1b).⁵⁰ (1) CO₂ pathway (red): *CO₂ remains adsorbed while *N₂ undergoes stepwise reduction to form *N₂H_x, enabling C–N coupling at the latter stage. (2) OCOH pathway (orange): *CO₂ is reduced to *OCOH, which couple with *N₂H_x intermediates. (3) CO pathway (blue): CO₂ is reduced to *CO, which then couples with hydrogenated *N₂ intermediates, such as *NH₂. These diverse pathways highlight the mechanistic complexity of urea electrosynthesis and emphasize the importance of catalyst-specific properties in guiding the preferred route.

3.2 Screening of M@C₂N for urea electrosynthesis

Electrocatalytic urea synthesis proceeds via a six-electron reaction: N₂ + CO₂ + 6H⁺ + 6e⁻ → CO(NH₂)₂ + H₂O. This process involves multiple competing reaction pathways and reactive intermediates. Therefore, effective descriptors are essential for rapidly screening promising electrocatalysts. First, we examined the geometric and electronic structures of 26 M@C₂N systems, as illustrated in Fig. S1–S3.†E_f = (E_M@C2N − E_C2N − E_M)U_diss. As shown in Fig. 1b, most M@C₂N systems meet this criterion, with the exception of Zn@C₂N, which was excluded due to its electrochemical instability.

Fig. 1 (a) Structural model of 2D M@C₂N. (b) Computed formation energies and dissolution potentials. (c) Two possible configurations of N₂ adsorption (side-on and end-on). (d) Adsorption energies of N₂ on 3d-, 4d-, and 5d-M@C₂N surfaces.

Activation of *N₂ is essential for enabling the co-reduction with *CO₂ in urea synthesis. Metal centers with unoccupied d-orbitals can accept lone-pair electrons from N₂, facilitating adsorption. To access this, we calculated the Gibbs free energies of *N₂ adsorption for both side-on and end-on configurations (Fig. 1d). All systems showed strong *N₂ binding, with end-on adsorption being energetically favored. As no catalysts were excluded at this stage, we focused on end-on adsorption in subsequent analyses. We also evaluated the co-adsorption of *CO₂ and *N₂ on the same metal site. Except for Pd@C₂N, most systems could co-adsorb both reactants, which is a required criterion for C–N coupling (Fig. 2a).

Fig. 2 (a) Optimized configurations and charge density difference (CDD) plots for *N₂ and *N₂ + *CO₂ on Nb@C₂N. Charge depletion and accumulation are depicted by cyan and yellow, respectively (isosurface = 0.003 e Å⁻³). (b) Molecular orbitals of free N₂, partial density of states (PDOS) and crystal orbital Hamilton population (COHP) analyses of *N₂ and *N₂+*CO₂ on Nb@C₂N. (c) Bonding and anti-bonding states in COHP are depicted by cyan and yellow, respectively.

Next, we analyzed *N₂ and *CO₂ + *N₂ co-adsorption on Nb@C₂N. Optimized structures and charge density difference (CDD) plots (Fig. 2b) revealed strong chemisorption with significant charge redistribution at the Nb center. Partial density of states (PDOS) showed strong orbital overlap between the 2s orbital of N and the 4s/3d_z² orbitals of Nb. According to the integrated crystal orbital Hamilton population (ICOHP) analysis (Fig. 2c), we evaluated the bonding and anti-bonding interactions. The COHP analysis offers valuable insights into the strength and nature of the chemical bonds. Negative projected COHP (pCOHP) values indicate bonding interactions, while positive values denote anti-bonding interactions. The ICOHP values for *N₂ were more negative compared to the co-adsorbed case, indicating a stronger *N₂ interaction that may hinder subsequent reaction steps. Additionally, CO₂ adsorption slightly weakens the Nb–N bond (lengthening from 1.97 Å to 2.06 Å), suggesting that *CO₂ moderates *N₂ interaction, potentially facilitating the reaction progress.

To ensure low energy input during urea formation, the ΔG values for the first hydrogenation step forming *NNH should be below 0.78 eV, based on the benchmark PdCu catalyst,³³ which was used as a general criterion for evaluating the catalytic activity.^16,48,51,52 We calculated ΔG for the initial *NNH formation across the four pathways: CO₂, OCOH, CO, and NCON, where this protonation step consistently served as the PDS.⁵³ Applying the ΔG < 0.78 eV cutoff,³³ we eliminated 11 catalysts (Sc, Fe, Co, Ni, Cu, Rh, Pd, Os, Ir, Pt, and Au). In contrast, ten early transition metals (Ti, V, Cr, Mn, Nb, Mo, Ru, Hf, Ta, W, and Re) met the criteria, outperforming PdCu.

We subsequently investigated the full reaction pathways for these 10 candidates (Fig. S17–S35†). End-on *N₂ adsorption was confirmed to be essential for *CO insertion and tower-like *NCON formation. However, consistent with previous work, the NCON pathway was not favored for any M@C₂N systems. To ensure product selectivity, we also evaluated the limiting potentials of competing pathways (e.g., CO₂RR and HER) and confirmed that favorable candidates showed higher selectivity toward urea formation (Fig. 3).

Fig. 3 (a) Maximum Gibbs free energy changes (ΔG) for urea production across four pathways on M@C₂N (M = 3d-, 4d-, and 5d-M): CO₂ (*CO₂ + *N₂ → *CO₂ + *NNH), OCOH (*OCOH + *N₂ → *OCOH + *NNH), CO (*CO + *N₂ → *CO + *NNH), and NCON (*CO + *N₂ → *NCON) (b) optimized structures of critical intermediates for each pathway on Nb@C₂N.

3.3 Catalytic activity of urea synthesis on screened M@C₂N catalysts

Based on detailed reaction pathway analyses, we identified distinct catalytic trends among the ten selected M@C₂N systems. The potential-determining steps (PDSs) and preferred pathways varied significantly across different metals, highlighting the diversity in their electronic structures and surface reactivity. For Ti, Mn, Mo, Hf, and W@C₂N, urea formation proceeds via the CO pathway, with the Gibbs free energy at PDS (ΔG_max) of 0.85, 0.73, 0.51, 0.80, and 0.56 eV, respectively. In most cases, the first protonation step (*CO + *N₂ → *CO + *NNH) was identified as the PDS. Interestingly, for W@C₂N, the PDS shifted to the final protonation step (*CO + *NHNH₂ → *CO + *NH₂NH₂), suggesting a change in the reaction-limiting feature based on surface electronic configuration. In the OCOH pathway, catalysts such as V, Cr, Nb, Ta, and Re@C₂N show feasible energetics, with PDS values of 0.71, 0.73, 0.57, 0.53, and 0.94 eV, respectively. For V, Cr, and Nb@C₂N, the first protonation step (*OCOH + *N₂ → *CO + *NNH) remains the bottleneck. Interestingly, for Ta@C₂N, the final protonation step (*OCOH + *NHNH₂ → *CO + *NH₂NH₂) was found to limit the reaction.

In the CO₂ pathway, the PDS is typically the last protonation step (*CO₂ + *NHNH₂ → *CO₂ + *NH₂NH₂), except for W and Re@C₂N, where the first *N₂ protonation dominates. These trends are summarized in Table S2,† highlighting the importance of protonation steps in determining catalytic performance. Among the ten candidates, we found that Mn, Mo, and W@C₂N preferentially follow the CO pathway, showing relatively low PDS values. For Cr and Ta@C₂N, they favor the OCOH pathway with PDS values of 0.73 and 0.53 eV, respectively. Lastly, we found that V, Zr, Nb, Hf, and Re@C₂N are likely to proceed through the CO₂ pathway, as summarized in Table 1. Interestingly, the distinction in the rate-limiting steps between pathways and catalysts demonstrates the diverse reactivity patterns, providing insights into tailoring catalysts for specific pathways to optimize urea electrosynthesis efficiency.

Table 1 Gibbs free energy (ΔG_max) of the potential determining step (PDS) of urea production on the screened M@C₂N surfaces (M = V, Cr, Mn, Nb, Mo, Ru, Hf, Ta, W, and Re)

Catalysts	Pathway	Potential determining step	ΔGmax (eV)
Mn@C2N	CO	*CO + *N2 → *CO + *NNH	0.73
Mo@C2N	CO	*CO + *N2 → *CO + *NNH	0.51
W@C2N	CO	*CO + *NHNH2 → *CO + *NH2NH2	0.56
Cr@C2N	OCOH	*OCOH + *N2 → *OCOH + *NNH	0.73
Ta@C2N	OCOH	*OCOH + *NHNH2 → *OCOH + *NH2NH2	0.53
V@C2N	CO2	*CO2 + *NHNH2 → *CO2 + *NH2NH2	0.60
Zr@C2N	CO2	*CO2 + *NHNH2 → *CO2 + *NH2NH2	0.64
Nb@C2N	CO2	*CO2 + *N2 → *CO2 + *NNH	0.50
Hf@C2N	CO2	*CO2 + *NHNH2 → *CO2 + *NH2NH2	0.78
Re@C2N	CO2	*CO2 + *N2 → *CO2 + *NNH	0.51

---

Among these, Nb@C₂N, Mo@C₂N, and Re@C₂N exhibit the lowest ΔG_max and were investigated in detail (Fig. 4 and Scheme 1b). Structural properties of the three catalysts are shown in Fig. S36.† For Nb@C₂N and Re@C₂N, the urea formation proceeds via the CO₂ pathway. After *CO₂ adsorption, *N₂ is subsequently adsorbed and hydrogenated, with the first *N₂ protonation (*CO₂ + *N₂ → *CO₂ + *NNH) serving as the PDS (ΔG_max = 0.51 and 0.50 eV), corresponding to low limiting potentials (U_L) of −0.50 V and −0.51 V, respectively. Subsequent hydrogenation of *NNH and coupling with *CO₂ proceed through thermodynamically downhill steps (ΔG = −1.22 eV for Nb and −0.23 eV for Re), facilitating product formation.

Fig. 4 Gibbs free energy diagrams and important intermediates along the most favorable urea formation pathways: (a) CO₂ pathway on Nb@C₂N, (b) CO pathway on Mo@C₂N, (c) CO₂ pathway on Re@C₂N, and (d) representative structures of key steps on Nb@C₂N.

Notably, Nb@C₂N achieves performance comparable to or exceeding benchmark catalysts such as PdCu (−0.64 V) and Fe-pBN (−0.63 V),¹⁵ (Table S9†), suggesting its superior potential for urea synthesis. Interestingly, when B is co-doped (Nb–B@C₂N),⁵¹ the pathway shifts to CO coupling, with the PDS being *CO + *N₂ → *CO + *NNH (ΔG = 0.56 eV), slightly less favorable than that of pristine Nb@C₂N. For Mo@C₂N, urea formation proceeds through the CO pathway. After *CO₂ is adsorbed and reduced to *CO, *N₂ adsorption and subsequent hydrogenation lead to the formation of *NNH. This first *N₂ reduction step is the PDS (ΔG = 0.51 eV). The final coupling of *CO and *NH₂NH₂ is thermodynamically favorable. In comparison, the Mo–B@C₂N catalyst⁵¹ favors the OCOH pathway, with a slightly higher ΔG (0.53 eV), indicating that Mo@C₂N alone is more efficient.

3.4 Descriptor for electrochemical urea production

The accurate prediction of electrocatalytic activity for urea synthesis typically requires computationally intensive modeling of reaction intermediates and transition states. This poses a major challenge for large-scale screening of potential catalysts. To overcome this limitation, identifying simple, intrinsic-property-based descriptors that capture the structure–activity relationship is essential for accelerating catalyst discovery. Previous studies have explored the use of adsorption energies, local charges, and intermediate adsorption energies in descriptor construction.^16,54–56 However, these quantities often require DFT-level calculations, limiting their scalability. A descriptor, Φ, was originally introduced in earlier work for the 2N₂ + CO₂ mechanism in the N₂ reduction reaction (NRR).^57–59 In our study, we extend the application of Φ to the urea synthesis reaction proceeding via distinct mechanistic pathways involving *CO₂, *CO, *OCOH, and *NCON intermediates. These mechanisms differ significantly from the previously studied route and include multiple possible C–N coupling steps. By evaluating Φ across these mechanistic variations, we show that it effectively correlates with catalytic activity (as reflected in the limiting potential, U_L) for a wide range of M@C₂N catalysts.

The descriptor, Φ, is defined as:

where N_d represents the number of d electrons of the metal atom, and c_M and c_sub are the electronegativities of the transition metal atom and substrate, respectively. Here, c_sub is defined as nc_N, where c_N represents the electronegativity of nitrogen, and n is the number of nitrogen atoms coordinating the M atom. The term quantifies the influence of metal–support interactions on the M's electron-attracting capability. Φ can be interpreted as the effective number of d-electrons and provides a direct relationship between the electronic structure and catalytic performance. This formulation allows for a rapid, physically interpretable screening of electronic structure effects across different SACs anchored on g-C₂N. Detailed parameters are summarized in Table S7.†

In addition to Φ, we evaluated two other descriptors: φ (ref. 60) and ε_d. The descriptor φ is an averaged electronegativity-based term incorporating the metal and its coordinated N atoms, and ε_d refers to the metal-centered d-band energy. While both φ and ε_d are established indicators in the SAC literature, we find that Φ exhibits the strongest correlation with U_L for urea formation within the M@C₂N catalyst set. This is reflected in the volcano plots and R² values shown in Fig. 5a–c, and summarized in Table S8.†

Fig. 5 Volcano plot showing the relationship between the limiting potential (U_L = −ΔG_max/|e|) and (a) Φ, (b) φ, and (c) ε_d for M@C₂N systems. (d) Volcano plot of U_Lvs. Φ including reported electrocatalysts (e.g., TM₂@C₄N₃, MBenes, PdCu, TM-B@C₂N, etc.). The grey dot indicates an outlier or an exception.

To further assess the generalizability of Φ, we collected data from various reported systems including TM₂@C₄N₃,⁵⁷ PdCu alloy,⁹ MBenes,⁴⁸ CuB₁₂,⁵⁰ M/p-BN,¹⁵ V₂N₆C,⁶¹ MoP(111),⁶² CuPc⁶³, and our previously reported TM–B@C₂N.⁵¹ We found that Φ is currently most effective within the g-C₂N support framework and closely related materials such as C₄N₃, as also shown in Fig. 5d and data are available in Table S9.† The g-C₂N framework provides a highly symmetric, nitrogen-rich coordination environment with six evenly distributed N atoms around the anchoring site. This uniformity ensures a consistent electronic interaction between the metal center and the support, minimizing geometric and electronic heterogeneity across different SACs. The planar conjugated structure also facilitates electron delocalization and effective charge redistribution upon adsorption of intermediates. These features contribute to the stability of the M–N_x motifs and enhance the reliability of structure–activity correlations, such as those captured by the descriptor Φ. In contrast, when applied to catalysts with other supports, coordination environments, or active site geometries, the descriptor's predictive power diminishes. Φ is considered a system-specific yet practical design tool within the context of M@C₂N-based electrocatalysts. Nonetheless, the demonstrated applicability of Φ across multiple mechanisms and 26 SACs highlights its utility for identifying promising catalysts efficiently within a large design space. This work lays a foundation for future development of more generalized, cross-platform descriptors that incorporate both electronic and geometric features for rational catalyst discovery in complex multistep electrosynthesis reactions such as urea production.

3.5 Selectivity of promising M@C₂N

In the electrocatalytic urea production, the HER is the major competing process in aqueous solutions, thereby reducing the faradaic efficiency. To assess the selectivity of Nb@C₂N, Mo@C₂N, and Re@C₂N, we compared the adsorption energies (ΔE) of H⁺ (Volmer step) and H₂O with N₂ adsorption. Fig. S37† shows that the adsorption energies for N₂ are −1.26 eV (Nb@C₂N), −0.90 eV (Mo@C₂N), and −1.27 eV (Re@C₂N), which are significantly lower (more favorable) than those for H⁺ and H₂O. These results suggest that the co-reduction of N₂ and CO₂ is energetically favored over the HER under acidic to neutral conditions (pH ≤ 7), supporting the high selectivity of these catalysts for urea production.

To further evaluate competing CO₂RR, particularly the formation of C1 products such as CO and CH₄ on Nb@C₂N, we analyzed the relevant reaction pathways. As depicted in Fig. 6a, CO desorption is energetically unfavorable (ΔG = 1.80 eV), while CH₄ desorption is thermodynamically favorable (ΔG = −0.15 eV), suggesting that once formed, CH₄ can readily desorb. The CO₂RR pathway proceeds via *CO₂ + H⁺ + e⁻ → *OCOH → *CO. Subsequently, *CO may either undergo further reduction or participate in C–N coupling via an Eley–Rideal mechanism with *N₂H_y intermediates, favoring urea production. Alternatively, *CO can be hydrogenated to *CHO or *COH, initiating the multi-step reduction to CH₄. The formation of *CHO and *COH requires 0.62 eV and 1.06 eV, respectively, with the PDS being *CO → *CHO (ΔG = 0.62 eV). This value exceeds the PDS for urea formation (0.50 eV), indicating that CH₄ formation is less favorable. Similar energetic profiles were observed for Mo@C₂N and Re@C₂N, confirming the suppression of CH₄ formation on these catalysts.

Fig. 6 (a) Schematic diagram of CO₂ reduction reaction (CRR) and N₂ reduction reaction (NRR) pathways. Gibbs free energy diagrams on Nb@C₂N for (b) the CRR to C1 products (CO and CH₄) and (c) the NRR to NH₃.

Given that nitrogen reduction to NH₃ is thermodynamically more favorable than urea-formation, evaluating the selectivity against the N₂ reduction reaction (NRR) is also crucial. We explored the NRR via three pathways: alternating (A), distal (D), and mixed (DA) mechanisms (Fig. 6c, S39, S41, and S43†). The most favorable route for NH₃ formation is: N₂ → *NNH → *NNH₂ → *N → *NH → *NH₂ → *NH₃ → NH₃. This mechanism is consistent with our previous findings.⁵¹ As shown in Fig. 6c, the first protonation step (*N₂ → *NNH) is the PDS for NH₃ production (ΔG = 0.40 eV), the same step as that for urea formation. Thus, NH₃ is likely co-produced alongside urea on these catalysts. Interestingly, CO₂ co-adsorption can alter the reaction pathway. As demonstrated in Fig. 4a, the presence of CO₂ modifies the NRR profile, especially on Nb@C₂N. This suggests that tuning the relative concentrations of CO₂ and N₂ could be an effective strategy to steer selectivity toward urea over NH₃.

In summary, Nb@C₂N, Mo@C₂N, and Re@C₂N exhibit high selectivity for urea synthesis under electrochemical conditions. Competing pathways such as the HER and CO₂RR (CH₄ formation) are suppressed due to higher energy barriers, while NH₃ may emerge as a co-product due to shared initial steps. These catalysts strike a promising balance between activity and selectivity, making them strong candidates for high-performance urea electrosynthesis.

4 Conclusions

In this work, we systematically investigated the mechanistic pathways of electrocatalytic urea synthesis from N₂ and CO₂ using high-throughput density functional theory (DFT) calculations. By exploring a broad range of 3d, 4d, and 5d transition metals anchored on graphitic C₂N (M@C₂N), we applied a Four-Step Screening Strategy to identify catalysts that are simultaneously thermodynamically stable, catalytically active, and highly selective. Among the 25 evaluated M@C₂N systems, Nb@C₂N, Mo@C₂N, and Re@C₂N emerged as the most promising candidates, exhibiting low limiting potentials of −0.50 V, −0.51 V, and −0.51 V, respectively, surpassing the benchmark PdCu catalyst (−0.78 V). Their outstanding performance arises from the strong and selective adsorption of N₂ and favorable energetics of key protonation steps, particularly the first and final hydrogenation events of *N₂ intermediates.

To accelerate catalyst discovery, we proposed a simple yet predictive descriptor, Φ, derived from the number of d-electrons and the electronegativity of the metal and support. This descriptor captures the intrinsic electronic structure–reactivity relationship and successfully reproduces activity trends across both the screened M@C₂N systems and reported electrocatalysts. Our findings validate Φ as a transferable and computationally inexpensive tool for rapid catalyst screening. Altogether, this study provides critical mechanistic insights into the co-reduction of N₂ and CO₂, identifies top-performing SACs on C₂N supports, and introduces a robust descriptor-driven framework for rational catalyst design. These findings pave the way for the development of next-generation electrocatalysts for sustainable urea synthesis under ambient conditions.

Data availability

The authors confirm that the data supporting the findings of this study are available within the article and its ESI.†

Conflicts of interest

There are no conflicts to declare.

Acknowledgements

This work was supported by the National Research Council of Thailand, Thailand (NRCT): N42A670101 and Program Management Unit for Human Resources & Institutional Development, Research and Innovation (grant number: B41G680026). We would like to thank the NSTDA Supercomputer Center (ThaiSC) for computational resources.

References

1. J. W. Erisman, M. A. Sutton, J. Galloway, Z. Klimont and W. Winiwarter, Nat. Geosci., 2008, 1, 636–639.
2. M. Aresta, A. Dibenedetto and A. Angelini, Chem. Rev., 2014, 114, 1709–1742.
3. B. M. Comer, P. Fuentes, C. O. Dimkpa, Y.-H. Liu, C. A. Fernandez, P. Arora, M. Realff, U. Singh, M. C. Hatzell and A. J. Medford, Joule, 2019, 3, 1578–1605.
4. M. Pérez-Fortes, A. Bocin-Dumitriu and E. Tzimas, Energy Procedia, 2014, 63, 7968–7975.
5. F. Barzagli, F. Mani and M. Peruzzini, Green Chem., 2011, 13, 1267–1274.
6. R. F. Service, Science, 2014, 345, 610.
7. D. Bao, Q. Zhang, F. L. Meng, H. X. Zhong, M. M. Shi, Y. Zhang, J. M. Yan, Q. Jiang and X. B. Zhang, Adv. Mater., 2017, 29, 1604799.
8. S. Giddey, S. Badwal and A. Kulkarni, Int. J. Hydrogen Energy, 2013, 38, 14576–14594.
9. C. Chen, X. Zhu, X. Wen, Y. Zhou, L. Zhou, H. Li, L. Tao, Q. Li, S. Du and T. Liu, Nat. Chem., 2020, 12, 717–724.
10. Y. Feng, H. Yang, Y. Zhang, X. Huang, L. Li, T. Cheng and Q. Shao, Nano Lett., 2020, 20, 8282–8289.
11. M. Yuan, J. Chen, Y. Bai, Z. Liu, J. Zhang, T. Zhao, Q. Wang, S. Li, H. He and G. Zhang, Angew. Chem., 2021, 133, 11005–11013.
12. M. Yuan, H. Zhang, Y. Xu, R. Liu, R. Wang, T. Zhao, J. Zhang, Z. Liu, H. He and C. Yang, Chem Catal., 2022, 2, 309–320.
13. Y. Xiao, C. Shen, Z. Xiong, J. Li and W. Zhang, Mater. Today Phys., 2022, 26, 100726.
14. C. Zhu, C. Wen, M. Wang, M. Zhang, Y. Geng and Z. Su, Chem. Sci., 2022, 13, 1342–1354.
15. L. Kong, D. Jiao, Z. Wang, Y. Liu, Y. Shang, L. Yin, Q. Cai and J. Zhao, Chem. Eng. J., 2023, 451, 138885.
16. C. Zhu, M. Wang, C. Wen, M. Zhang, Y. Geng, G. Zhu and Z. Su, Adv. Sci., 2022, 9, 2105697.
17. P. Roy, A. Pramanik and P. Sarkar, J. Phys. Chem. Lett., 2021, 12, 10837–10844.
18. Z. Zhang and L. Guo, Dalton Trans., 2021, 50, 11158–11166.
19. H.-R. Pan, T. Tang and J.-S. Hu, Fundam. Res., 2021, 1, 461–465.
20. J. Fu, Y. Yang and J.-S. Hu, ACS Mater. Lett., 2021, 3, 1468–1476.
21. X. Li, L. Liu, X. Ren, J. Gao, Y. Huang and B. Liu, Sci. Adv., 2020, 6, eabb6833.
22. C.-X. Zhao, B.-Q. Li, J.-N. Liu and Q. Zhang, Angew. Chem., Int. Ed., 2021, 60, 4448–4463.
23. Z. Xiong, Y. Xiao and C. Shen, Chem. Mater., 2022, 34, 9402–9413.
24. L. Cheng, Z. Cheng, M. Fan, Y. Ren, Y. Liu, T. Yuan and Z. Shen, Electrochim. Acta, 2023, 470, 143283.
25. J. K. Nørskov, F. Abild-Pedersen, F. Studt and T. Bligaard, Proc. Natl. Acad. Sci. U. S. A., 2011, 108, 937–943.
26. J. K. Nørskov, T. Bligaard, J. Rossmeisl and C. H. Christensen, Nat. Chem., 2009, 1, 37–46.
27. S. Liu, Z. Li, C. Wang, W. Tao, M. Huang, M. Zuo, Y. Yang, K. Yang, L. Zhang and S. Chen, Nat. Commun., 2020, 11, 938.
28. W. T. Hong, M. Risch, K. A. Stoerzinger, A. Grimaud, J. Suntivich and Y. Shao-Horn, Energy Environ. Sci., 2015, 8, 1404–1427.
29. W. Zhong, Y. Qiu, H. Shen, X. Wang, J. Yuan, C. Jia, S. Bi and J. Jiang, J. Am. Chem. Soc., 2021, 143, 4405–4413.
30. K. Liu, J. Fu, Y. Lin, T. Luo, G. Ni, H. Li, Z. Lin and M. Liu, Nat. Commun., 2022, 13, 2075.
31. S. Ringe, Nat. Commun., 2023, 14, 2598.
32. X. Zhang, L. Wang, M. Li, W. Meng, Y. Liu, X. Dai, G. Liu, Y. Gu, J. Liu and L. Kou, Mater. Today, 2023, 67, 23–32.
33. C. Chen, X. Zhu, X. Wen, Y. Zhou, L. Zhou, H. Li, L. Tao, Q. Li, S. Du, T. Liu, D. Yan, C. Xie, Y. Zou, Y. Wang, R. Chen, J. Huo, Y. Li, J. Cheng, H. Su, X. Zhao, W. Cheng, Q. Liu, H. Lin, J. Luo, J. Chen, M. Dong, K. Cheng, C. Li and S. Wang, Nat. Chem., 2020, 12, 717–724.
34. G. Kresse and J. Hafner, Physical Review B:Condensed Matter, 1993, 47, 558.
35. G. Kresse and J. Furthmüller, Comput. Mater. Sci., 1996, 6, 15–50.
36. G. Kresse and J. Furthmüller, Physical Review B:Condensed Matter, 1996, 54, 11169.
37. J. P. Perdew, K. Burke and M. Ernzerhof, Phys. Rev. Lett., 1996, 77, 3865.
38. G. Kresse and D. Joubert, Physical Review B:Condensed Matter, 1999, 59, 1758.
39. S. Grimme, J. Antony, S. Ehrlich and H. Krieg, J. Chem. Phys., 2010, 132, 154104.
40. W. Tang, E. Sanville and G. Henkelman, J. Phys.: Condens. Matter, 2009, 21, 084204.
41. R. Nelson, C. Ertural, J. George, V. L. Deringer, G. Hautier and R. Dronskowski, J. Comput. Chem., 2020, 41, 1931–1940.
42. P. C. Müller, C. Ertural, J. Hempelmann and R. Dronskowski, J. Phys. Chem. C, 2021, 125, 7959–7970.
43. D. Ren, H. Li and X. Cheng, Solid State Commun., 2015, 223, 54–59.
44. X. Lv, W. Wei, B. Huang, Y. Dai and T. Frauenheim, Nano Lett., 2021, 21, 1871–1878.
45. X. Guo, J. Gu, S. Lin, S. Zhang, Z. Chen and S. Huang, J. Am. Chem. Soc., 2020, 142, 5709–5721.
46. J. K. Nørskov, J. Rossmeisl, A. Logadottir, L. Lindqvist, J. R. Kitchin, T. Bligaard and H. Jónsson, J. Phys. Chem. B, 2004, 108, 17886–17892.
47. J.-H. Liu, L.-M. Yang and E. Ganz, J. Mater. Chem. A, 2019, 7, 11944–11952.
48. X. Zhu, X. Zhou, Y. Jing and Y. Li, Nat. Commun., 2021, 12, 4080.
49. D. Jiao, Z. Wang, Y. Liu, Q. Cai, J. Zhao, C. R. Cabrera and Z. Chen, Energy Environ. Mater., 2024, 7, e12496.
50. C. Zhu, C. Wen, M. Wang, M. Zhang, Y. Geng and Z. Su, Chem. Sci., 2022, 13, 1342–1354.
51. N. Yodsin, P. Mano, K. Takahashi and s. namuangruk, ChemCatChem, 2024, 16(21), e202400304.
52. Y. Xiao, C. Shen, Z. Xiong, J. Li and W. Zhang, Mater. Today Phys., 2022, 26, 100726.
53. N. Yodsin, P. Mano, K. Takahashi and S. Namuangruk, ChemCatChem, 2024, 16, e202400304.
54. C. Ling, Y. Ouyang, Q. Li, X. Bai, X. Mao, A. Du and J. Wang, Small Methods, 2019, 3, 1800376.
55. H. Li, Y. Liu, K. Chen, J. T. Margraf, Y. Li and K. Reuter, ACS Catal., 2021, 11, 7906–7914.
56. Z. Ren, X. Wang, S. Wang, H. Zhang, B. Huang, Y. Dai and W. Wei, J. Mater. Chem. A, 2023, 11, 11507–11516.
57. J. Liu, X. Lv, Y. Ma, S. C. Smith, Y. Gu and L. Kou, ACS Nano, 2023, 17, 25667–25678.
58. Z. Fu, M. Wu, Q. Li, C. Ling and J. Wang, Mater. Horiz., 2023, 10, 852–858.
59. H. Niu, X. Wang, C. Shao, Z. Zhang and Y. Guo, ACS Sustain. Chem. Eng., 2020, 8, 13749–13758.
60. M. Luo, X. Cai, Y. Ni, Y. Chen, C. Guo and H. Wang, ACS Appl. Mater. Interfaces, 2023, 16, 807–818.
61. Z. Zhang and L. Guo, Dalton Trans., 2021, 50, 11158–11166.
62. D. Jiao, Y. Dong, X. Cui, Q. Cai, C. R. Cabrera, J. Zhao and Z. Chen, J. Mater. Chem. A, 2023, 11, 232–240.
63. J. Mukherjee, S. Paul, A. Adalder, S. Kapse, R. Thapa, S. Mandal, B. Ghorai, S. Sarkar and U. K. Ghorai, Adv. Funct. Mater., 2022, 32, 2200882.

---

Footnote

† Electronic supplementary information (ESI) available. See DOI: https://doi.org/10.1039/d5ta03783b

---