Design of Fe₂Mo@γ-GDY triatomic catalyst for electrocatalytic urea synthesis of N₂ and CO: a theoretical study

INTRODUCTION

With a notable high nitrogen content of 46%, urea [i.e., CO(NH₂)₂] is widely used worldwide as a nitrogen fertilizer, whose annual synthesis has exceeded 100 million tons^[1]. In addition, urea has various applications, such as cosmetic formulations, synthesis of engineered wood products, and monosodium glutamate^[2]. As a consequence, sustainable and efficient urea synthesis is of great importance. In the present industry, the first step of urea synthesis relies predominantly on the Haber-Bosch process for ammonia (NH₃) synthesis^[3,4], which requires a high temperature (~700 K) and an elevated pressure (~100 bar). Here as well, as a second step, the energy-intensive Bosch–Meiser process that combines NH₃ with CO₂ to synthesize urea^[5] necessitates elevated temperature (~423-473 K) and pressure (~150-200 bar)^[6]. Annually, the entire synthesis process not only accounts for approximately 1% to 2% of total energy consumption^[7], but also consumes 80% of global NH₃ synthesis^[8]. Therefore, it is imperative to develop a strategy for urea synthesis with clean, environmentally friendly, energy-efficient, and resource-conserving features.

Electrocatalysis is recognized as a promising alternative^[9] since it directly converts N₂ to urea under milder conditions, thereby bypassing the need for NH₃ synthesis and significantly conserving energy and resources. The electrocatalytic strategy involves the Langmuir-Hinshelwood (LH) and Eley-Rideal (ER) mechanisms [Figure 1A]. The difference between them is that N₂ and CO are co-adsorbed or not on the catalyst. Following the adsorption process, a key intermediate ^*NCON is formed through the C-N coupling. Then, hydrogenation of ^*NCON occurs gradually to form urea by the distal or the alternate mechanism. We note that the intermediate ^*NCON has been formed by Chen et al. through the LH mechanism and experimentally identified by analyzing the infrared spectrum vibration signal^[10]. Over there, the energy barrier (E_b) of the C-N coupling and the maximum free energy change in the whole reaction pathway (ΔG_max) serve as the critical indicators for evaluating the catalyst performance. The key C-N coupling process calls for the challenging destruction of the N≡N, the transfer of CO from its adsorption site to N atoms, and the creation of a C–N bond. However, the cleavage of the N≡N is highly challenging due to its great activation energy of ~941 kJ/mol^[11]. To effectively catalyze this step, the catalyst should strongly bind with N₂, aiming to fully activate the corresponding N₂. In this context, triatomic catalysts are able to act as excellent choices, where two of the constituent atoms can be designed to work in tandem to activate N₂ while the third one is optimized to weaken CO adsorption. To achieve the above goal, transition metals (TMs) can be selected due to their unfilled d-orbitals, which participate in a “donor-acceptor” mechanism with the adsorbate^[12]. As an ideal substrate, the two-dimensional (2D) γ-graphdiyne (γ-GDY) can be used to support TMs, because it, as the most stable graphdiyne isomer^[13,14], provides large pores to accommodate three TM atoms and the excellent stability and electrical conductivity^[15-17].

Figure 1. (A) Mechanisms of electrocatalytic synthesis of urea; (B) Curves for each element as an electrocatalyst for urea synthesis: U_L curve, E_ad*N2+*CO - E_b curve for the C-N coupling process, and E_b curve for the product; (C) Doping of Fe and Mo atoms into γ-GDY; (D and E) S-type and D-type structures after doping. γ-GDY: γ-graphdiyne.

Based on the considerations mentioned above, the selection of elements that compose TMs@γ-GDY is crucial. In conventional electrocatalyst design, researchers often use the volcano plot to select elements and design catalysts, where the adsorption energy is taken as the variant. However, this method usually considers the first and last steps of the catalytic reaction, neglecting the overall intermediate processes. With the development of electrocatalysis, we have gradually realized that the active sites in the catalytic reaction process are not fixed to one location but change during the reaction progress. Therefore, considering the intermediate steps in catalyst design is particularly important. Since intermediates in the reaction process always move towards the sites with the highest adsorption energy, we calculate the adsorption energy of the catalyst for intermediates in all the reaction steps, thereby optimizing each reaction step to promote the transfer of active sites, break the scaling relationship, and improve catalytic performance.

In this work, using density functional theory (DFT) calculations, we designed the Fe₂Mo@γ-GDY catalyst through the aforementioned design strategy. The catalyst exhibits a favorable limiting potential (U_L = -0.22 V) and an optimized C-N coupling energy barrier (E_b = 0.34 eV). Also, a detailed discussion is carried out on how the modulation of adsorption energy affects the reduction of ΔG_max through the regulation of adsorption configuration and sites. Additionally, the study explores the effects of N₂ and CO adsorption on the E_b of C-N coupling. Furthermore, enhanced selectivity for hydrogen evolution reaction (HER), carbon oxygen reduction reaction (CORR), and nitrogen reduction reaction (NRR) is observed, thereby offering a strategic approach for the design of efficient catalysts in the electrocatalytic synthesis of urea.

MATERIALS AND METHODS

All spin-polarized DFT calculations were performed within the DMol³ codes from the Materials Studio 2019 software^[18]. The exchange-correlation effects were described using the Perdew-Burke-Ernzerhof (PBE) functional within the generalized gradient approximation (GGA)^[19]. The Grimme method for DFT-D correction was employed^[20]. The core treatment utilized a DFT semi-core pseudopotential with a basis set of double numerical plus polarization (DNP)^[21]. Convergence was accelerated by adopting smearing with a value of 0.005 Ha, while electron convergence accuracy was ensured with a real-space global orbital cutoff radius of 4.9 Å. The convergence criteria for electronic relaxation, force, and displacement were set to 10^-5 eV, 0.002 eV/Å, and 0.005 Å, respectively. To analyze the bonding nature of the N≡N, the crystal orbital Hamilton population (COHP) analysis was performed using the Lobster software^[22,23].

The substrate was modeled by a 2 × 2 monolayer of γ-GDY (18.92 × 18.92 Ų). A vacuum of 20 Å was introduced along the z to avoid interactions from the periodic images in the direction normal to the substrate. Sampling of the Brillouin zone was performed using a 2 × 2 × 1 Monkhorst-Pack k-points mesh. The transition state (TS) analysis of the key intermediate (^*NCON) formed by C-N coupling was implemented using the linear synchronous transit (LST)/quadratic synchronous transit (QST) method^[24]. The basic parameter settings are the same as those used in other calculations, with the root mean square (RMS) convergence being set to 0.01 Ha/Å. The spin population and charge transfer were calculated using the Hirshfeld population analysis^[25]. The Gibbs free energy changes for each reaction step in the electrocatalytic synthesis of urea were evaluated based on the computational hydrogen electrode (CHE) model proposed by Nørskov et al.^[26]. Under standard reaction conditions (pH = 0, 298.15 K, and 1 atm) and at a potential of 0 V vs. reversible hydrogen electrode (RHE), the free energy associated with the transfer of a proton and electron pair is defined as half the value of gaseous hydrogen^[27],

The Gibbs free energy during the electrocatalytic synthesis process was calculated by

where ΔE represents the change of free energy obtained directly from DFT calculations, ΔZPE denotes the variation in zero-point energy, T is the Kelvin temperature (298.15 K), and ΔS indicates the change of entropy. For N₂, CO, and H₂, the ZPE and S values are obtained from the experiments^[28]. The limiting potential is computed as

where the variable e represents the number of electrons transferred during potential limiting step (PLS) and ΔG_max determines the PLS. The adsorption energy (E_ad) is calculated by

where E_total represents the total energy of the adsorption system, E_ads is the energy of the adsorbed species, and E^* denotes the total energy of the pristine catalyst. The mass loading of TMs atoms on the substrate of Fe₂Mo@γ-GDY is determined by

where m_TMs and m_catalyst denote the mass of Fe/Mo atoms loaded on the catalyst and the mass of the whole catalyst system, respectively. In our catalyst system, there are 72 C atoms and three TMs; therefore, m_TMs = n × N_Fe + (3 - n) × N_Mo,m_catalyst = 72 × N_C + m_TMs, where n represents the number of Fe atoms, and N_Fe, N_Mo and N_C are the corresponding molar mass.

RESULTS AND DISCUSSION

In the periodic table, the elements from groups V, VI, and VIII are often favored in the field of electrocatalysis due to their unique outer electron arrangements. Therefore, we select V, Cr, Fe, Nb, Mo, Ru, Ta, W, and Os as candidate elements for the catalyst design. For electrocatalytic urea synthesis, U_L, E_b, and E_ad of the products are critical factors. Therefore, we first calculate the U_L value for each element as a catalyst in the electrocatalytic synthesis of urea, the E_ad*N2+*CO - E_b value for the C-N coupling process, and the E_ad value for the product [Figure 1B], and select the top three elements with the best performance. As seen in the results from Table 1, it is clear that Fe, Mo, and W present great potential. Considering that these three factors do not fully represent the entire reaction process, we also calculated the E_ad - ΔE curve for each reaction step in the electrocatalytic synthesis of urea using each element as a catalyst [Supplementary Figures 1-6]. Since each reaction step requires suitable ΔE and E_ad, we select four elements in the middle of each step’s curve; the element screen is shown in Table 2. Here as well, we found that Mo, W, and Fe exhibit good performance during the reaction process. According to our design, among the three atoms, there must be both strong and weak adsorptions to allow intermediates to change configurations on the catalyst and regulate adsorption energy. In addition, considering that the melting point of W is as high as 3,400 °C, which poses an obstacle to its synthesis, we finally selected Mo and Fe for the catalyst design (here, Fe was chosen as the weak adsorption atom). We note that Fe and Mo have been widely considered as active sites to weaken the N≡N, such as the applications of FeMo₃S₄ and Fe/Mo-N-C in NRR^[29,30]. The unoccupied d-orbitals of the FeMo cluster can accept lone pair electrons from N₂. Conversely, the occupied d-orbitals of the FeMo cluster can donate electrons to the antibonding orbitals of N₂, thereby weakening the N≡N^[31]. Additionally, the synergistic interaction between Fe and Mo can significantly enhance catalytic performance. Efficient electron transfer is facilitated by the strong orbital coupling between Fe-3d and Mo-4d. This ensures not only excellent electron conduction from the catalyst to the adsorbed intermediate, but also remarkable stability of the catalyst^[32]. Therefore, we select Fe and Mo loaded onto γ-GDY to form Fe_xMo_3-x (x = 0, 1, 2, 3) as the electrocatalyst.

Figure 1C builds a Fe_xMo_3-x@γ-GDY (x = 0, 1, 2, 3) model, which generates four structures: Mo₃@γ-GDY, Fe₃@γ-GDY, Fe₂Mo@γ-GDY, and FeMo₂@γ-GDY. Table 3 presents the calculated ΔG_max and E_b values for all the four cases. It is noteworthy that the electrocatalytic urea synthesis process on Fe₃@γ-GDY cannot be successfully completed. According to the findings of Zhou et al.^[33], a lower Hirshfeld charge on N points [Supplementary Table 1] indicates weaker nucleophilicity, making N₂ less sensitive to adsorption (or alternatively, allowing H⁺ to be more easily adsorbed by the underlying Fe). It is clear from Table 3 that Fe₂Mo@γ-GDY possesses excellent electrocatalytic thermodynamic and kinetic performances for urea synthesis (i.e., much lower ΔG_max and E_b values). Here, the excellent performance of Fe₂Mo@γ-GDY is in line with our design described in the Introduction section.

The optimal structure of Fe₂Mo@γ-GDY is shown in Figure 2A. The average bond length of the three longer edges of the C18 hexagonal hole is 6.78 Å, and the average bond length between the metal and C atoms is 2.07 Å, which is sufficient to accommodate the Fe₂Mo cluster. To assess the electronic conductivity of Fe₂Mo@γ-GDY, an essential feature for electrocatalysis, the density of states (DOS) was calculated. As shown in Figure 2B, the DOS crosses the Fermi level, which demonstrates the metallic characteristics of Fe₂Mo@γ-GDY. To explore the stability of Fe₂Mo@γ-GDY, its frequency and formation energy were evaluated. Figure 2C provides the frequency distribution for Fe₂Mo@γ-GDY, where the absence of imaginary frequencies indicates excellent stability. In addition, a negative formation energy of E_form = -0.49 eV again demonstrates the stability of Fe₂Mo@γ-GDY. Furthermore, the ab initio molecular dynamics (AIMD) simulations within the canonical ensemble (NVT) at a temperature of 300 K were performed to further evaluate the dynamic stability under aqueous conditions. The time step was set to 1.0 fs and the total simulation time was set to 10.0 ps. The thermostat was set to Massive GGM, with a chain length of 2, a relaxation time of 10.0, and a Yoshida parameter of 3. It is clear from the snapshots of Fe₂Mo@γ-GDY at 0, 2, 4, 6, 8 and 10 ps [Figure 2D] that the structural deformation is insignificant.

Figure 2. (A) Optimized structures of Fe₂Mo@γ-GDY; (B) DOS of Fe₂Mo@γ-GDY with the Fermi level set at 0 eV, where the black and red lines represent spin up and spin down, respectively; (C) Frequency of Fe₂Mo@γ-GDY; (D) AIMD simulation of Fe₂Mo-D@γ-GDY at the temperature of 300 K under aqueous conditions. γ-GDY: γ-graphdiyne; DOS: density of states.

The theoretical feasibility of triatomic doping in graphene has been proposed^[34]. The distinctive sp-hybridized triple bond structure of γ-GDY suggests significant advantages in anchoring surface metal atoms^[35]. Moreover, triatomic doping in similar 2D materials has been realized. For instance, Wang et al. successfully synthesized the CNT@C₃N₄-Fe₂Cu catalyst, where Fe and Cu were predominantly loaded as triatomic clusters^[36]. In particular, CNT@C₃N₄-Fe₂Cu exhibits an 18.3% mass loading rate of metal atoms, while the Fe₂Mo@γ-GDY model presents 19.4%, which denotes a reasonable design of metal atom loading in our study. Based on the experimental results of CNT@C₃N₄-Fe₂Cu, we design a feasible experimental method for the fabrication of Fe₂Mo@γ-GDY: forming intermediate products between γ-GDY and precursor containing Mo and Fe by hydrothermal method and annealing them to produce Fe₂Mo@γ-GDY in the atmosphere of H₂ and Ar.

The adsorption of N₂ on Fe₂Mo@γ-GDY has two possible configurations: the S-type (Fe₂Mo-S@γ-GDY, Figure 1D) with N₂ on two identical atoms and the D-type with N₂ on two distinct atoms (Fe₂Mo-D@γ-GDY,Figure 1E). In order to reveal optimal adsorption configuration, the E_coad values of Fe₂Mo-S@γ-GDY (-0.35 eV) and Fe₂Mo-D@γ-GDY (-2.21 eV) were calculated. Obviously, N₂ always preferentially adsorbs on Fe₂Mo-D@γ-GDY (hereafter, Fe₂Mo@γ-GDY will refer to Fe₂Mo-D@γ-GDY). The mechanism of electrocatalytic urea synthesis on Fe₂Mo@γ-GDY involves four reaction pathways: LH and ER mechanisms during the C-N coupling and distal and alternative mechanisms during the proton-coupled electron transfer (PCET) steps. Given that the reaction pathways initially involve LH and ER mechanisms, the corresponding E_b values calculated are 0.34 and 0.64 eV, respectively, indicating the more favorable LH mechanism for the C-N coupling. Namely, N₂ and CO are firstly co-adsorbed on the Fe₂Mo@γ-GDY catalyst, followed by the formation of ^*NCON via the C-N coupling. Then, the ^*NCON transforms gradually to urea through the PCET steps via either the distal or the alternating mechanism. The ΔG_max value for the alternative mechanism is 0.22 eV, lower than 0.30 eV for the distal mechanism. Overall, the electrocatalytic urea synthesis process on Fe₂Mo@γ-GDY will proceed with the LH and alternative mechanisms.

To investigate the activation of N₂, we conducted a comprehensive analysis by calculating the partial DOS (PDOS) before and after N₂ adsorption [Figure 3A]. The findings reveal a remarkable alignment between the 2s and 2p orbitals of N₂ and the d-orbitals of Fe₂Mo, indicating a pronounced hybridization between N₂ and Fe₂Mo. Subsequent to this hybridization, a portion of the unoccupied 2π^* orbitals of N₂ falls below the Fermi level, signifying their partial occupation by the d-orbitals of Fe₂Mo. This observation supports the notion of a substantial electron transfer between Fe₂Mo and N₂, ultimately leading to the activation of the N≡N. The activated N≡N species demonstrate enhanced propensity to form a 2N^* adsorption state on the catalyst surface, thereby facilitating the C-N coupling process and reducing the E_b. Furthermore, the COHP for ^*N₂ on Fe₂Mo@γ-GDY and pure N₂ molecule are analyzed [Figure 3B], where red and blue areas represent the antibonding and bonding contributions, respectively. Here, the integrated COHP (ICOHP) can quantitatively describe the strength of N≡N, in which a less negative ICOHP value corresponds to a weaker N≡N^[37]. It is obvious that Fe₂Mo@γ-GDY weakens the N≡N.

Figure 3. (A) PDOS of Fe₂Mo, N₂ adsorbed at hollow site on Fe₂Mo@γ-GDY, and free N₂ molecule, surface. The black, blue, and red lines represent the outermost d orbitals, p orbitals, and s orbitals of the corresponding structures, respectively; (B) COHP for N≡N bond of free N₂ molecule, adsorbed N₂ on Fe₂Mo@γ-GDY. PDOS: Partial density of states; γ-GDY: γ-graphdiyne; COHP: crystal orbital Hamilton population.

During the protonation process, ^*NCON and ^*NCONH are adsorbed by Fe₂Mo@γ-GDY at the hollow site. Then, ^*HNCONH and ^*H₂NCONH are located by Fe₂Mo@γ-GDY at the bridge site. Finally, the adsorbed intermediate (^*H₂NCONH₂) undergoes a transition from the bridge site to the top site. Analysis of intermediate adsorption energies reveals a dynamic shift during the reaction process, facilitating a transition from a weaker to a stronger adsorption site. During the reaction, active sites relocate throughout the process, disrupting the scaling relationship^[38] and leading to a reduction in reaction-free energy^[39]. From the evolution of adsorption energy for the entire process [Figure 4A], it can be seen that the initial adsorption energy of the protonation step is the strongest. This robust initial adsorption favors H⁺ binding, promoting the reaction and enhancing adsorption capacity for later H⁺ binding events. Subsequently, the adsorption energy gradually diminishes, culminating in adsorption at top site. This shift is conducive to reaction termination, facilitating the easy desorption of formed urea molecules from the reaction site. This enables the catalyst to progress seamlessly to the subsequent synthesis steps. The ΔG value at each step is shown in Figure 4B. It is worth stressing that the electrocatalytic performance of Fe₂Mo@γ-GDY is superior to other catalysts reported in the DFT-based studies on urea synthesis, as compared in Figure 4C^[40-46].

Figure 4. (A) Evolution of adsorption energy; (B) Free energy step on Fe₂Mo@γ-GDY; (C) Performance comparison. γ-GDY: γ-graphdiyne.

At this stage, the selectivity of Fe₂Mo@γ-GDY is considered. The competitive reactions of the electrocatalytic synthesis of urea mainly include the HER, CORR, and NRR; their structural diagrams are shown in Supplementary Figures 7-9. For HER, we analyzed the H⁺ competition in each reaction step prior to the protonation step. For the first and second steps, which relate to adsorption processes, the adsorption energies were used for comparison. For the C-N coupling process, the changes in free energy were used for comparison. For ^* - ^*N₂, the E_ad of the catalyst for N₂ is -0.93 eV, while for H⁺, it is only -0.23 eV; for ^*N₂ - ^*N₂ + ^*CO, the E_ad for CO is -1.28 eV, while for H⁺, it is a positive value of 0.06 eV; for ^*N₂ + ^*CO - ^*NCON, ΔG for the C-N coupling is -0.18 eV. The formation of ^*N₂H + ^*CO has a ΔG value of 0.17 eV, and the formation of ^*N₂ + ^*COH has a ΔG value of -0.06 eV. Therefore, HER is difficult to occur. There are two paths for the first hydrogenation of CORR^[47]. One is the COH formation by hydrogenation on O, and the other is the CHO formation by hydrogenation on C. The free energy values of the two pathways from ^*N₂ + ^*CO are -0.06 and 1.17 eV, respectively. Both are higher than the C-N coupling process (-0.18 eV). As to NRR, the free energy for forming ^*NNH + ^*CO is 0.00 eV, which is higher than the free energy of C-N coupling (-0.18 eV), indicating that retrograde C-N coupling is more preferred. The reaction free energy profile for NRR is depicted in Supplementary Figure 9, where ΔG_max equals 0.63 eV, significantly higher than that of urea synthesis (0.22 eV). Thus, Fe₂Mo@γ-GDY has a good selectivity for electrocatalytic urea synthesis.

CONCLUSIONS

We have designed a triatomic Fe₂Mo@γ-GDY catalyst for the electrochemical synthesis of urea by regulating the adsorption energy and adsorption configuration for each step of the reaction process. Through DFT calculations, we have demonstrated that this catalyst possesses good stability and electronic conductivity. The catalyst not only gives an excellent catalytic performance in theoretical calculations (U_L = -0.22 V, E_b = 0.34 eV), but also shows a good selectivity through various evaluation metrics. Additionally, the catalyst exhibits strong antioxidative performance under operational potentials. Compared to the current methods of urea synthesis under high-temperature and high-pressure conditions, the Fe₂Mo@γ-GDY catalyst offers a viable route for synthesizing urea under ambient conditions. Overall, the Fe₂Mo@γ-GDY catalyst provides a new and more sustainable method for the electrocatalytic synthesis of urea, and it is expected to replace traditional energy-intensive synthesis methods in future industrial applications.