First-Principles Study of Bimetallic Pairs Embedded on Graphene Co-Doped with N and O for N2 Electroreduction

The electrocatalytic nitrogen reduction reaction (NRR) is considered a viable alternative to the Haber–Bosch process for ammonia synthesis, and the design of highly active and selective catalysts is crucial for the industrialization of the NRR. Dual-atom catalysts (DACs) with dual active sites offer flexible active sites and synergistic effects between atoms, providing more possibilities for the tuning of catalytic performance. In this study, we designed 48 graphene-based DACs with N4O2 coordination (MM′@N4O2-G) using density functional theory. Through a series of screening strategies, we explored the reaction mechanisms of the NRR for eight catalysts in depth and revealed the “acceptance–donation” mechanism between the active sites and the N2 molecules through electronic structure analysis. The study found that the limiting potential of the catalysts exhibited a volcano-shaped relationship with the d-band center of the active sites, indicating that the synergistic effect between the bimetallic components can regulate the d-band center position of the active metal M, thereby controlling the reaction activity. Furthermore, we investigated the selectivity of the eight DACs and identified five potential NRR catalysts. Among them, MoCo@N4O2-G showed the best NRR performance, with a limiting potential of −0.20 V. This study provides theoretical insights for the design and development of efficient NRR electrocatalysts.


Introduction
Ammonia (NH 3 ) is a key precursor in fertilizer synthesis and a carbon-free energy carrier that possesses advantages such as emission-free combustion, convenient storage, and high energy density.Its role in sustainable development is crucial [1][2][3].However, current industrial ammonia production heavily relies on the energy-intensive Haber-Bosch process, requiring high temperatures and pressure (300-500 • C, 150-300 atm) for NH 3 synthesis [4].This process not only consumes a notable amount of energy (approximately 1-2% of global energy consumption annually), but also releases substantial greenhouse gases.As a result, there is an urgent need to explore sustainable alternatives [5].Electrocatalytic nitrogen reduction (eNRR) has emerged as a promising alternative method due to its mild reaction conditions, sustainability, and environmental friendliness, potentially replacing the conventional Haber-Bosch process [6].However, the activation of N 2 molecules faces challenges due to the high bond energy of the N≡N bond (941 kJ•mol −1 ), low polarizability, and lack of a dipole moment [7].Consequently, NRR catalysis often requires a high limiting voltage to overcome these obstacles.Additionally, the simultaneous hydrogen evolution reaction (HER) during eNRR compromises the NH 3 selectivity and Faradaic efficiency (FE) [8].Thus, the focus lies in developing electrocatalysts with outstanding catalytic activity, high selectivity, and superior FE for large-scale industrial eNRR applications.
In recent years, atomically dispersed transition metal (TM) catalysts, including singleatom catalysts (SACs) and dual-atom catalysts (DACs) have garnered substantial interest in the field of catalysis due to their high atomic utilization, tunable electronic structures, and unique local chemical environments [9][10][11][12].Currently, there have been a large number of computational and experimental studies on SACs for eNRR.In atomically dispersed transition metal catalysts, choosing suitable support materials is crucial to prevent metal agglomeration and enhance the catalyst stability [13][14][15][16].Two-dimensional materials, such as graphene, C 2 N, g-C 3 N 4 , h-BN, and MXene, have been commonly used as substrates for SACs or DACs owing to their large surface area, ordered structure, and controllable electronic properties.Among these, graphene stands out as an ideal substrate for SACs and DACs due to its high charge carrier mobility, excellent conductivity, and chemical stability [17][18][19][20][21]. Recently, special attention has been directed towards the family of transition metal atoms anchored on N-doped graphene (TM-N x /G), such as TM-N 4 [22], due to its high reactivity and stability.However, they often exhibit lower overpotential in the hydrogen evolution reaction (HER), leading to lower Faradaic efficiency for the NRR.To address this issue, strategies involving the adjustment of active sites and/or coordination environments have been proposed to improve the catalytic performance and stability [23][24][25].
As an extension of SACs, DACs utilize the synergistic effect introduced by the second metal to maintain the low oxidation state of transition metals and effectively activate inert molecules [26].Experimental and theoretical evidence indicates that some DACs, such as Mn 2 @C 2 N [27], FeM-N 6 -C [28]; Fe 2 N 4 @graphene [29], FeRu@N 4 -P [30], and Mo 2 @PC 6 [31], exhibit superior catalytic performance in eNRR compared to Ru (0001) (−0.98 V) [32].Nevertheless, the stability and catalytic performance of these catalysts are highly sensitive to the coordination environment of the metal centers [33].Previous studies have demonstrated that adjusting the atoms in the first coordination layer around the transition metal atoms is an effective strategy to modulate the catalytic performance.In TM-N x /G, the introduction of additional non-metallic dopants (such as oxygen, boron, sulfur, etc.) in the first coordination layer around the transition metal atom significantly enhances the catalytic performance beyond N atom coordination alone [34][35][36].For instance, in the M 2 N 6 /G system, the Mn 2 diatomic catalyst with O and N co-coordination (Mn 2 ON 5 /G α ) demonstrates significantly improved NRR activity and selectivity [37].In addition, by employing different transition metals and ligands, DACs have the potential to break linear proportionality in certain reactions [38].DACs exhibit promising prospects in eNRR applications due to the ability to flexibly choose different metal combinations and coordination environments, and the number of active centers.
A recent experimental approach demonstrated a synthetic method to produce a library of DACs using macrocyclic precursors through an encapsulation-thermal decomposition synthesis strategy.This led to the successful synthesis of a series of nitrogen and oxygen co-doped porous carbon-supported DACs (MM ′ @N 4 O 2 -G) [39].These complexes allow for a wide range of metal center modulation, including 3d transition metals (Mn, Fe, Co, Ni, Cu, Zn) and noble metals (Pd, Rh, Pt).Furthermore, by independently controlling elements, the formation of both homonuclear and heteronuclear bimetallic centers, such as Fe 2 , Co 2 , FeCu, and CuCo, can be achieved.The study revealed that these DACs exhibit excellent catalytic activity and stability in the oxygen reduction reaction (ORR), particularly with FeCu-DAC outperforming corresponding SACs and other Fe-based DACs.This progress has inspired us to design dual-atom catalysts with diverse metal combinations and N 4 O 2 coordination for the NRR.
In this study, we employed density functional theory (DFT) to investigate the catalytic performance of a series of DACs with the M 2 N 4 O 2 motif embedded in graphene (MM ′ @N 4 O 2 -G) for the NRR to NH 3 .In terms of metal selection, considering the common metals found in the active centers of nitrogenase (such as Mo, V, and Fe) [36], we chose these three metals as M, while M ′ included 3d transition metals (Sc~Zn), 4d transition metals (Zr, Nb, and Mo), and 5d transition metals (Hf, Ta, W, and Re).The computational results indicate that these five DACs exhibit good catalytic performance for the NRR while inhibiting the side reaction (HER).Detailed analysis of the electronic structure properties of these catalysts revealed that the high activity of the DACs stemmed from the effective modulation of the intermediate stability and synergistic effects between the active sites.This study is anticipated to offer an important reference and guidance for the development of highly active and selective DACs.

Catalyst Structure and Stability
Based on the successful synthesis of DACs possessing N 4 O 2 coordination in experiments, we constructed 48 types of MM ′ @N 4 O 2 -G DACs on this substrate, including 3 homonuclear DACs, VV/MoMo/FeFe@N 4 O 2 -G, and 45 heteronuclear DACs, MM ′ @N 4 O 2 -G, where M = V/Fe/Mo and M ′ = Sc ~Zn, Zr, Nb, Mo, Hf, Ta, W, and Re.The optimized structures are shown in Figure 1a.To assess the thermodynamic stability of these catalysts, we calculated the formation energy (E f ) as follows: where E cat is the total energy of the optimized catalyst, E sub is the energy of the catalyst substrate, and E M is the energy of an individual metal atom in its most stable bulk phase structure.Subsequently, we calculated the dissolution potential (U diss ) to evaluate the electrochemical stability: /N e where U diss represents the standard dissolution potential of the bulk metal, denoted as ) /2 for the diatomic system; N e represents the number of electrons transferred during metal dissolution, expressed as N e = [N e (M)+ N e (M ′ )]/2 [28].According to the calculated results shown in Figure 1b, the E f values of these catalysts are all negative, indicating their thermodynamic stability.Moreover, among the 48 diatomic catalysts, 33 exhibit U diss values greater than 0, indicating good electrochemical stability while maintaining thermodynamic stability.Consequently, they are potential candidate catalysts for further investigation.Additionally, Bader charge analysis was used to calculate the charge transfer between the anchored metal atoms and the substrate.As shown in Table S1, a substantial number of electrons transfer from the anchored metal atoms to N 4 O 2 -G, indicating a strong interaction between the metal and the substrate.properties of these catalysts revealed that the high activity of the DACs stemmed from the effective modulation of the intermediate stability and synergistic effects between the active sites.This study is anticipated to offer an important reference and guidance for the development of highly active and selective DACs.

Catalyst Structure and Stability
Based on the successful synthesis of DACs possessing N4O2 coordination in experiments, we constructed 48 types of MM′@N4O2-G DACs on this substrate, including 3 homonuclear DACs, VV/MoMo/FeFe@N4O2-G, and 45 heteronuclear DACs, MM′@N4O2-G, where M = V/Fe/Mo and M′ = Sc ~ Zn, Zr, Nb, Mo, Hf, Ta, W, and Re.The optimized structures are shown in Figure 1a.To assess the thermodynamic stability of these catalysts, we calculated the formation energy (Ef) as follows: where Ecat is the total energy of the optimized catalyst, Esub is the energy of the catalyst substrate, and EM is the energy of an individual metal atom in its most stable bulk phase structure.Subsequently, we calculated the dissolution potential (Udiss) to evaluate the electrochemical stability: where Udiss represents the standard dissolution potential of the bulk metal, denoted as      /2 for the diatomic system; Ne represents the number of electrons transferred during metal dissolution, expressed as Ne = [Ne(M)+ Ne(M′)]/2 [28].
According to the calculated results shown in Figure 1b, the Ef values of these catalysts are all negative, indicating their thermodynamic stability.Moreover, among the 48 diatomic catalysts, 33 exhibit Udiss values greater than 0, indicating good electrochemical stability while maintaining thermodynamic stability.Consequently, they are potential candidate catalysts for further investigation.Additionally, Bader charge analysis was used to calculate the charge transfer between the anchored metal atoms and the substrate.As shown in Tables S1, a substantial number of electrons transfer from the anchored metal atoms to N4O2-G, indicating a strong interaction between the metal and the substrate.2.2.Performance of MM ′ @N 4 O 2 -G for Electrocatalytic NRR 2.2.1.Screening of Catalysts N 2 adsorption.The adsorption and activation of N 2 molecules represent the initial and crucial steps in the NRR process.To determine the most favorable adsorption structures, three types of adsorption configurations were considered for homonuclear DACs: one involving end-on adsorption, where one N atom of N 2 forms a bond with a single metal atom, and two types of side-on adsorption, one with two N atoms of N 2 adsorbed on a single metal atom and the other with two N atoms of N 2 bonded to two different metal atoms.Regarding heteronuclear DACs, given the diverse nature of the two metal active sites, five adsorption configurations were considered, comprising two end-on adsorptions and three side-on adsorptions, as shown in Figure 2a. Figure 2b presents the adsorption free energy (∆G *N2 ) of the optimized N 2 adsorption configurations on 33 catalysts.Notably, except for the physical adsorption of N 2 molecules on four catalysts (MM ′ @N 4 O 2 -G, MM ′ = FeTa, FeMn, FeCo, and FeNi), the other 29 catalysts demonstrate the chemical adsorption of N 2 molecules on the surface, with adsorption free energy values ranging from −0.06 to −1.07 eV.Among them, N 2 molecules exhibit side-on adsorption on four catalysts (MM ′ @N 4 O 2 -G, MM ′ = VV, VTi, MoTa, and MoTi), while end-on adsorption is mainly observed on 25 other catalysts, mainly at the M or M' active sites of MM ′ @N 4 O 2 -G.Thus, based on the calculated adsorption free energy of N 2 , only four catalysts (MM ′ @N 4 O 2 -G, MM ′ = FeTa, FeMn, FeCo, and FeNi) are excluded.Thermal stability of the catalyst.The stability of catalysts under operating conditions is crucial for their practical application.To further assess the thermal stability of the 8 catalysts, we conducted a 10 ps ab initio molecular dynamics (AIMD) simulation at a temperature of 500 K. Figures 4 and S1 illustrate the variations in temperature and energy over time for these catalysts.Clearly, during the 10 ps duration, there were negligible structural changes observed in the catalysts, indicating excellent thermal stability.Competitive adsorption at active sites.In the NRR, competitive adsorption among different species at active sites is an important concern, especially under experimental conditions where hydrogen protons and water molecules in the solution may compete for adsorption with N 2 molecules at the active sites on the catalyst surface.If the catalyst exhibits stronger adsorption of hydrogen atoms compared to N 2 molecules, this could potentially lead to the occupation of active sites by hydrogen atoms, favoring the promotion of the HER over the NRR.Such a case leads to the decreased Faradaic efficiency of the NRR, thereby impacting the selectivity and efficiency of the overall reaction.To evaluate this, we compared the adsorption free energies of *H (∆G *H ) and ∆G *N2 .Figure 2c illustrates that 17 catalysts exhibit stronger adsorption towards N 2 molecules than towards H atoms, indicating that these catalysts are more conducive to facilitating an efficient electrochemical NRR.Additionally, if solvent molecules cover metal atoms instead of N 2 molecules, it can affect the sustained NRR and lead to the oxidation of transition metals in the aqueous electrolyte, ultimately hindering the progress of the reaction [40].To address this issue, the adsorption free energy of H 2 O molecules (∆G *H2O ) on the remaining 17 catalysts was calculated, as shown in Figure 3d.Ten MM ′ @N 4 O 2 -G catalysts (MM ′ = MoFe, MoRe, VFe, FeFe, VV, MoMn, MoCo, VMn, MoCr, and MoNi) were identified for their favorable adsorption behavior towards N 2 molecules within the competitive adsorption environment.These catalysts show promising potential for an efficient NRR.The protonation of the first step and the last step.The eNRR involves six protoncoupled electron transfer (PCET) steps.Despite the specific mechanism of the NRR, the hydrogenation reactions of the first step (*N 2 + H + + e − → *NNH) and the sixth step (*NH 2 + H + + e − → *NH 3 ) are the most common elementary steps.Previous studies have indicated [41][42][43] that these two steps usually act as the potential determining steps (PDS) in the NRR, demonstrating the largest free energy change throughout the reaction process.To evaluate these crucial steps, we calculated the reaction free energy change for the protonation of the first step (∆G N2→NNH ) and the last step (∆G NH2→NH3 ).Using a standard value of 0.65 eV, we performed preliminary screening for the aforementioned ten candidate diatomic catalysts.As shown in Figure 3, eight catalysts met the criteria set in this study, including one homonuclear DAC (VV@N 4 O 2 -G) and seven heteronuclear DACs (MM ′ @N 4 O 2 -G, MM ′ = MoCo, MoFe, MoCr, MoMn, MoRe, VFe, and VMn).Additionally, it is notable that, for most diatomic catalysts, the free energy change for the protonation of N 2 to form NNH in the first step is greater, except for MoCo@N 4 O 2 -G, where the ∆G NH2→NH3 (0.20 eV) is greater compared to ∆G N2→NNH (0.07 eV).
Thermal stability of the catalyst.The stability of catalysts under operating conditions is crucial for their practical application.To further assess the thermal stability of the 8 catalysts, we conducted a 10 ps ab initio molecular dynamics (AIMD) simulation at a temperature of 500 K. Figure 4 and Figure S1 illustrate the variations in temperature and energy over time for these catalysts.Clearly, during the 10 ps duration, there were negligible structural changes observed in the catalysts, indicating excellent thermal stability.

NRR Reaction Mechanism
For the eight DACs selected through the aforementioned screening process, we conducted detailed calculations of the possible NRR pathways to assess their catalyti performance.Based on the N2 adsorption configurations and different hydrogenation sequences of the two N atoms, the electrocatalytic NRR on DAC surfaces typically in volves various pathways, as shown in Figure 5.For N2 adsorption through the end-on pattern, the NRR can proceed via either distal or alternating pathways for protonation reactions.In the distal pathway, the proton-electron pairs initially react with the N atom away from the adsorption site, resulting in the formation of the first NH3.Subsequently consecutive protonation steps lead to the formation of the second NH3.In the alternating pathway, the proton-electron pairs alternately attack the two N atoms, eventually lead

NRR Reaction Mechanism
For the eight DACs selected through the aforementioned screening process, we conducted detailed calculations of the possible NRR pathways to assess their catalytic performance.Based on the N 2 adsorption configurations and different hydrogenation sequences of the two N atoms, the electrocatalytic NRR on DAC surfaces typically involves various pathways, as shown in Figure 5.For N 2 adsorption through the end-on pattern, the NRR can proceed via either distal or alternating pathways for protonation reactions.In the distal pathway, the proton-electron pairs initially react with the N atom away from the adsorption site, resulting in the formation of the first NH 3 .Subsequently, consecutive protonation steps lead to the formation of the second NH 3 .In the alternating pathway, the proton-electron pairs alternately attack the two N atoms, eventually leading to the sequential production of two NH 3 molecules.Regarding the side-on adsorption pattern, the NRR occurs via two pathways: the enzymatic (red line) and the consecutive pathways (brown line).Additionally, the NRR can also take place through a mixed pathway, alternating between the distal and alternating pathways or between the enzymatic and consecutive pathways.
ing to the sequential production of two NH3 molecules.Regarding the side-on adsorption pattern, the NRR occurs via two pathways: the enzymatic (red line) and the consecutive pathways (brown line).Additionally, the NRR can also take place through a mixed pathway, alternating between the distal and alternating pathways or between the enzymatic and consecutive pathways.It is noteworthy that the desorption of the two NH3 molecules from VV@N4O2-G requires relatively high energies, at 1.34 eV and 0.82 eV, respectively.However, previous studies have demonstrated that NH3 generated in strong acid solutions can be easily reduced to NH4 + [44]; hence, the desorption of NH3 is not extensively considered here.Our calculations show that for VV@N4O2-G, the most probable reaction pathway is the enzymatic pathway, with a UL of −0.32 V. ; hence, the desorption of NH 3 is not extensively considered here.Our calculations show that for VV@N 4 O 2 -G, the most probable reaction pathway is the enzymatic pathway, with a U L of −0.32 V.
ing to the sequential production of two NH3 molecules.Regarding the side-on adsorption pattern, the NRR occurs via two pathways: the enzymatic (red line) and the consecutive pathways (brown line).Additionally, the NRR can also take place through a mixed pathway, alternating between the distal and alternating pathways or between the enzymatic and consecutive pathways.It is noteworthy that the desorption of the two NH3 molecules from VV@N4O2-G requires relatively high energies, at 1.34 eV and 0.82 eV, respectively.However, previous studies have demonstrated that NH3 generated in strong acid solutions can be easily reduced to NH4 + [44]; hence, the desorption of NH3 is not extensively considered here.Our calculations show that for VV@N4O2-G, the most probable reaction pathway is the enzymatic pathway, with a UL of −0.32 V. Due to the end-on adsorption of N2 on heteronuclear DACs, when a N2 molecule adsorbs on one metal atom and undergoes the NRR through a distal pathway, the other metal atom can also serve as a reactive site.We further investigated the NRR mechanism when two N2 molecules were simultaneously adsorbed on these seven heteronuclear DACs.Firstly, we studied the co-adsorption of two N2 molecules on the diatomic sites.It Due to the end-on adsorption of N 2 on heteronuclear DACs, when a N 2 molecule adsorbs on one metal atom and undergoes the NRR through a distal pathway, the other metal atom can also serve as a reactive site.We further investigated the NRR mechanism when two N 2 molecules were simultaneously adsorbed on these seven heteronuclear DACs.Firstly, we studied the co-adsorption of two N 2 molecules on the diatomic sites.S2).Conversely, on the MoRe@N 4 O 2 -G surface, the ∆G of this step is lower than the ∆G max when a single N 2 is adsorbed, suggesting that the MoRe@N 4 O 2 -G catalyst can simultaneously adsorb two N 2 molecules for the NRR, preferentially inducing the first hydrogenation of the N 2 adsorbed on the Mo atom.
In the subsequent reaction processes, because both N 2 molecules have the potential for hydrogenation, we compared the free energy changes of two hydrogenation elementary the rate-determining step's free energy obtained when a single N2 molecule is adsorb indicating that simultaneously adsorbing two N2 molecules is not feasible on these surfaces (Table S2).Conversely, on the MoRe@N4O2-G surface, the ΔG of this ste lower than the ΔGmax when a single N2 is adsorbed, suggesting that the MoRe@N4O catalyst can simultaneously adsorb two N2 molecules for the NRR, preferentially ind ing the first hydrogenation of the N2 adsorbed on the Mo atom.
In the subsequent reaction processes, because both N2 molecules have the poten for hydrogenation, we compared the free energy changes of two hydrogenation elem

Origin of NRR Catalytic Activity
To investigate the underlying factors influencing the activity of DACs in the N we conducted electronic structure calculations on these eight DACs.Firstly, we analy the charge transfer between the N2 molecule and the catalyst through charge den

Origin of NRR Catalytic Activity
To investigate the underlying factors influencing the activity of DACs in the NRR, we conducted electronic structure calculations on these eight DACs.Firstly, we analyzed the charge transfer between the N 2 molecule and the catalyst through charge density difference (CDD) and Bader charge analysis.Taking MoCo@N 4 O 2 as an example, as shown in Figure 9a, evident charge transfer between the active site and N 2 is observed, with a tendency for charge accumulation near the proximal N atom of N 2 , reducing the charge density between the two N atoms, thereby weakening the chemical bond and facilitating N 2 activation.Additionally, the increased N-N bond length after N 2 adsorption also reflects the activation of N 2 .Compared to free N 2 molecules (d N-N = 1.114Å), the adsorbed N 2 exhibits a significantly increased N-N bond length, ranging from 1.148 Å to 1.269 Å, indicating that the electron transfer between the catalyst and N 2 effectively activates the N 2 molecule.
Molecules 2024, 29, x FOR PEER REVIEW Compared to the free N2 molecular orbitals, the 2π and 3σ orbitals of the adsor shift upwards and exhibit significant hybridization with the Mo 3d orbitals bel Fermi level.This indicates that the unoccupied 3d orbitals of the Mo atom acce trons from the 2π and 3σ orbitals of the N2 molecule, forming bonding states th mote nitrogen adsorption.On the other hand, the unoccupied 2π* orbitals of N towards the Fermi level after adsorption, forming partially occupied 2π* orbita suggests that the occupied 3d orbitals of the Mo atom donate electrons to the antib orbitals of N2, thereby weakening the strength of the N-N bond and facilitating quent hydrogenation reactions.Similar situations are observed in the PDOS o DACs (Figure S4), indicating that N2 activation on these catalysts follows an "accep donation" mechanism.Moreover, to investigate the mechanism of the synergistic effects between d sites, we explored the influence of the d-band center (ɛd) of the active sites on the r activity.Specifically, we evaluated the d-band center of the Mo atom in the case Co@N4O2-G and assessed the d-band center of the two V atoms in VV@N4O2-G.As in Figure 10, a distinct volcano-shaped relationship exists between the limiting po of these eight DACs and the d-band centers of the active metal atoms.Notably, the efficient MoCo@N4O2-G is located near the peak of the volcano plot.The volcan suggests that the superior NRR performance of DACs is attributed to the appr position of the d-band center.Additionally, it is evident that although the metal M case of heteronuclear DAC) does not directly participate in the hydrogenation pro N2, the synergistic effects between the M′ and M sites effectively regulate the pos the d-band center of the active site, thereby impacting the reaction activity.In order to gain deeper insights into the fundamental electron transfer mechanism during N 2 activation, using MoCo@N 4 O 2 -G as an example, the partial density of states (PDOS) of the DACs before and after N 2 adsorption was studied, as shown in Figure 9b,c.Compared to the free N 2 molecular orbitals, the 2π and 3σ orbitals of the adsorbed N 2 shift upwards and exhibit significant hybridization with the Mo 3d orbitals below the Fermi level.This indicates that the unoccupied 3d orbitals of the Mo atom accept electrons from the 2π and 3σ orbitals of the N 2 molecule, forming bonding states that promote nitrogen adsorption.On the other hand, the unoccupied 2π* orbitals of N 2 move towards the Fermi level after adsorption, forming partially occupied 2π* orbitals.This suggests that the occupied 3d orbitals of the Mo atom donate electrons to the antibonding orbitals of N 2 , thereby weakening the strength of the N-N bond and facilitating subsequent hydrogenation reactions.Similar situations are observed in the PDOS of other DACs (Figure S4), indicating that N 2 activation on these catalysts follows an "acceptance-donation" mechanism.
Moreover, to investigate the mechanism of the synergistic effects between diatomic sites, we explored the influence of the d-band center (E d ) of the active sites on the reac-tion activity.Specifically, we evaluated the d-band center of the Mo atom in the case of MoCo@N 4 O 2 -G and assessed the d-band center of the two V atoms in VV@N 4 O 2 -G.As shown in Figure 10, a distinct volcano-shaped relationship exists between the limiting potentials of these eight DACs and the d-band centers of the active metal atoms.Notably, the highly efficient MoCo@N 4 O 2 -G is located near the peak of the volcano plot.The volcano curve suggests that the superior NRR performance of DACs is attributed to the appropriate position of the d-band center.Additionally, it is evident that although the metal M ′ (in the case of heteronuclear DAC) does not directly participate in the hydrogenation process of N 2 , the synergistic effects between the M ′ and M sites effectively regulate the position of the d-band center of the active site, thereby impacting the reaction activity.

NRR Selectivity of the MM'@N4O2-G Catalysts
The main competing reaction during the NRR is the HER.This competition greatly influences the selectivity of the catalyst.The difference between UL(NRR) and UL(HER) is commonly used to assess the selectivity of a catalyst.S3).The results show that all four catalysts are more likely to form *NNH. Therefore, these five catalysts (MM′@N4O2-G, MM′ = MoFe, MoMn, MoCo, MoCr, and VV) exhibit good selectivity and hold potential as catalysts for the NRR.

Computational Methods
All computations in this study, based on spin-polarized density functional theory (DFT) [45,46], were conducted using the Vienna Ab Initio Simulation Package (VASP 5.4.4) [47,48].The projector augmented wave (PAW) [49] method was employed to deal with the ion-electron interactions.The cutoff energy for the plane-wave basis set was set to 450 eV.The Perdew-Burke-Ernzerhof (PBE) functional [50] within the general gradient approximation (GGA) was used to describe the electronic exchange-correlation interactions.A (5 × 5) graphene supercell was adopted as the catalyst substrate, with a vacuum layer of 20 Å introduced along the z-axis to eliminate the interaction between periodic images.For structure optimization and electronic structure calculations, Monkhorst-Pack k-point grids of 3 × 3 × 1 and 11 × 11 × 1 were utilized to sample the Brillouin zone.The convergence criteria for energy and forces were set to 10 −5 eV and 0.02 eV/Å, respectively.To account for van der Waals (vdW) interactions, the DFT-D3 method proposed by Grimme et al. [51] was employed in all calculations.The implicit solvent model implemented in the VASPsol software package (VASPsol 5.4.1) was used to treat the solvation effects [52,53].To investigate the thermal stability of the DACs, ab initio molecular dynamics (AIMD) simulations [54] were performed for 10 ps at 500 K with a time step of 2 fs.
The Gibbs free energy change (ΔG) for each elementary step of the NRR was calculated using the computational hydrogen electrode (CHE) model proposed by Nørskov et al. [55].The formula for calculating ΔG is as follows: where ΔE is the reaction energy of each step calculated by DFT.ΔEZPE and ΔS are the changes in zero-point energy and entropy at 298.15 K, respectively, obtained by calculating the vibrational frequencies.The vibrational frequencies and entropy of gas molecules (N2, H2, NH3) are obtained from the NIST database [56].U represents the electrode potential, n represents the number of transferred electrons, and ΔGpH represents the free energy correction value at pH, defined as ΔGpH = 2.303 × kBT × pH.In this work, the pH is set to 0. The highest positive ΔG value (ΔGmax) throughout the process was employed to derive the limiting potential (UL), i.e., UL = −ΔGmax/e.MoCo@N 4 O 2 -G exhibits the most favorable catalytic performance among the five DACs, with a limiting potential of −0.20 V. Its catalytic activity surpasses that of several other catalysts, including TiV-CG (−0.30V) [23], FeMo-N 6 -C (−0.63 V) [28], Fe 2 N 4 @graphene (−0.32 V) [29], and Mn 2 ON 5 /G α (−0.27V) [37].

Computational Methods
All computations in this study, based on spin-polarized density functional theory (DFT) [45,46], were conducted using the Vienna Ab Initio Simulation Package (VASP 5.4.4) [47,48].The projector augmented wave (PAW) [49] method was employed to deal with the ion-electron interactions.The cutoff energy for the plane-wave basis set was set to 450 eV.The Perdew-Burke-Ernzerhof (PBE) functional [50] within the general gradient approximation (GGA) was used to describe the electronic exchange-correlation interactions.A (5 × 5) graphene supercell was adopted as the catalyst substrate, with a vacuum layer of 20 Å introduced along the z-axis to eliminate the interaction between periodic images.For structure optimization and electronic structure calculations, Monkhorst-Pack k-point grids of 3 × 3 × 1 and 11 × 11 × 1 were utilized to sample the Brillouin zone.The convergence criteria for energy and forces were set to 10 −5 eV and 0.02 eV/Å, respectively.To account for van der Waals (vdW) interactions, the DFT-D3 method proposed by Grimme et al. [51] was employed in all calculations.The implicit solvent model implemented in the VASPsol software package (VASPsol 5.4.1) was used to treat the solvation effects [52,53].To investigate the thermal stability of the DACs, ab initio molecular dynamics (AIMD) simulations [54] were performed for 10 ps at 500 K with a time step of 2 fs.
The Gibbs free energy change (∆G) for each elementary step of the NRR was calculated using the computational hydrogen electrode (CHE) model proposed by Nørskov et al. [55].The formula for calculating ∆G is as follows: where ∆E is the reaction energy of each step calculated by DFT.∆E ZPE and ∆S are the changes in zero-point energy and entropy at 298.15 K, respectively, obtained by calculating the vibrational frequencies.The vibrational frequencies and entropy of gas molecules (N 2 , H 2 , NH 3 ) are obtained from the NIST database [56].U represents the electrode potential,

Figure 1 .
Figure 1.(a) Schematic of the optimized catalyst structure.(b) Corresponding formation energy and dissolution potential of the catalyst.Figure 1.(a) Schematic of the optimized catalyst structure.(b) Corresponding formation energy and dissolution potential of the catalyst.

Figure 1 .
Figure 1.(a) Schematic of the optimized catalyst structure.(b) Corresponding formation energy and dissolution potential of the catalyst.Figure 1.(a) Schematic of the optimized catalyst structure.(b) Corresponding formation energy and dissolution potential of the catalyst.

Molecules 2024 ,
29,  x FOR PEER REVIEW 5 of 16 of N2 to form NNH in the first step is greater, except for MoCo@N4O2-G, where the ΔGNH2→NH3 (0.20 eV) is greater compared to ΔGN2→NNH (0.07 eV).

Figure 2 .
Figure 2. (a) Schematic structure of N2 adsorption on catalyst.(b) The adsorption free energy of N2 on the catalyst.(c) Comparison between ΔG*H and ΔG*N2.(d) Comparison between ΔG*N2 and ΔG*H2O, where the red coordinates represent the catalysts that meet the requirements.

Figure 2 .
Figure 2. (a) Schematic structure of N 2 adsorption on catalyst.(b) The adsorption free energy of N 2 on the catalyst.(c) Comparison between ∆G *H and ∆G *N2 .(d) Comparison between ∆G *N2 and ∆G *H2O , where the red coordinates represent the catalysts that meet the requirements.

Figure 3 .
Figure 3. Free energy changes (∆G N2→NNH and ∆G NH2→NH3 , eV) for the first and last protonation steps on MM ′ @N 4 O 2 -G.The black dashed line indicates the screening criteria (∆G = 0.65 eV).

Figure 4 .
Figure 4. Total energy variation of MoFe@N4O2-G for AIMD simulation at 500 K for 10 ps.The C, N O, Mo, and Fe atoms are labeled as gray, blue, red, green, and lavender balls, respectively.

Figure 4 .
Figure 4. Total energy variation of MoFe@N 4 O 2 -G for AIMD simulation at 500 K for 10 ps.The C, N, O, Mo, and Fe atoms are labeled as gray, blue, red, green, and lavender balls, respectively.

Figure 5 .
Figure 5. Diagram of possible reaction mechanisms for NRR.VV@N4O2-G was the only homonuclear DAC that remained after the screening process.N2 molecules exhibit side-on adsorption on the catalyst surface, where two N atoms bond to two V atoms, with ΔG*N2 of −0.76 eV.The Gibbs free energy diagram for the NRR on VV@N4O2-G and corresponding intermediate structures are shown in Figure 6.As shown in the figure, for the VV@N4O2-G catalyst, the potential-determining step (PDS) in the consecutive pathway is the second step of the protonation reaction (*NNH + H + + e − → *NNH2), with a ΔG value of 0.49 eV.The PDS of both the enzymatic and mixed pathways is the first protonation step (*N2 + H + + e − → *NNH), with a ΔG value of 0.32 eV.The first four protonation steps in the two pathways lead to the *NHNH2 intermediate.In the fifth protonation step, a proton-electron pair attacks one N atom in the *NHNH2 intermediate, forming *NH2NH2 or *NHNH3, with ΔG values of 1.84 and −1.35 eV, respectively, indicating that the former is more feasible in thermodynamical terms.Subsequently, *NH2 + *NH2 undergoes two hydrogenation steps to produce two adsorbed NH3 molecules, with ΔG values of 0.02 and −0.02 eV.It is noteworthy that the desorption of the two NH3 molecules from VV@N4O2-G requires relatively high energies, at 1.34 eV and 0.82 eV, respectively.However, previous studies have demonstrated that NH3 generated in strong acid solutions can be easily reduced to NH4 +[44]; hence, the desorption of NH3 is not extensively considered here.Our calculations show that for VV@N4O2-G, the most probable reaction pathway is the enzymatic pathway, with a UL of −0.32 V.

Figure 5 .
Figure 5. Diagram of possible reaction mechanisms for NRR.VV@N 4 O 2 -G was the only homonuclear DAC that remained after the screening process.N 2 molecules exhibit side-on adsorption on the catalyst surface, where two N atoms bond to two V atoms, with ∆G *N2 of −0.76 eV.The Gibbs free energy diagram for the NRR on VV@N 4 O 2 -G and corresponding intermediate structures are shown in Figure 6.As shown in the figure, for the VV@N 4 O 2 -G catalyst, the potential-determining step (PDS) in the consecutive pathway is the second step of the protonation reaction (*NNH + H + + e − → *NNH 2 ), with a ∆G value of 0.49 eV.The PDS of both the enzymatic and mixed pathways is the first protonation step (*N 2 + H + + e − → *NNH), with a ∆G value of 0.32 eV.The first four protonation steps in the two pathways lead to the *NHNH 2 intermediate.In the fifth protonation step, a proton-electron pair attacks one N atom in the *NHNH 2 intermediate, forming *NH 2 NH 2 or *NHNH 3 , with ∆G values of 1.84 and −1.35 eV, respectively, indicating that the former is more feasible in thermodynamical terms.Subsequently, *NH 2 + *NH 2 undergoes two hydrogenation steps to produce two adsorbed NH 3 molecules, with ∆G values of 0.02 and −0.02 eV.It is noteworthy that the desorption of the two NH 3 molecules from VV@N 4 O 2 -G requires relatively high energies, at 1.34 eV and 0.82 eV, respectively.However, previous studies have demonstrated that NH 3 generated in strong acid solutions can be easily reduced to NH 4 + [44]; hence, the desorption of NH 3 is not extensively considered here.Our calculations show that for VV@N 4 O 2 -G, the most probable reaction pathway is the enzymatic pathway, with a U L of −0.32 V.

Figure 5 .
Figure 5. Diagram of possible reaction mechanisms for NRR.VV@N4O2-G was the only homonuclear DAC that remained after the screening process.N2 molecules exhibit side-on adsorption on the catalyst surface, where two N atoms bond to two V atoms, with ΔG*N2 of −0.76 eV.The Gibbs free energy diagram for the NRR on VV@N4O2-G and corresponding intermediate structures are shown in Figure 6.As shown in the figure, for the VV@N4O2-G catalyst, the potential-determining step (PDS) in the consecutive pathway is the second step of the protonation reaction (*NNH + H + + e − → *NNH2), with a ΔG value of 0.49 eV.The PDS of both the enzymatic and mixed pathways is the first protonation step (*N2 + H + + e − → *NNH), with a ΔG value of 0.32 eV.The first four protonation steps in the two pathways lead to the *NHNH2 intermediate.In the fifth protonation step, a proton-electron pair attacks one N atom in the *NHNH2 intermediate, forming *NH2NH2 or *NHNH3, with ΔG values of 1.84 and −1.35 eV, respectively, indicating that the former is more feasible in thermodynamical terms.Subsequently, *NH2 + *NH2 undergoes two hydrogenation steps to produce two adsorbed NH3 molecules, with ΔG values of 0.02 and −0.02 eV.It is noteworthy that the desorption of the two NH3 molecules from VV@N4O2-G requires relatively high energies, at 1.34 eV and 0.82 eV, respectively.However, previous studies have demonstrated that NH3 generated in strong acid solutions can be easily reduced to NH4 +[44]; hence, the desorption of NH3 is not extensively considered here.Our calculations show that for VV@N4O2-G, the most probable reaction pathway is the enzymatic pathway, with a UL of −0.32 V.

Figure 6 . 16 Figure 7 .
Figure 6.Gibbs free energy diagrams of NRR on VV@N 4 O 2 -G.The C, N, O, H, and V atoms are labeled as gray, blue, red, white, and light blue balls, respectively.In the other seven heteronuclear DACs (MM ′ @N 4 O 2 -G, MM ′ = MoCo, MoCr, MoFe, MoMn, MoRe, VFe, and VMn), N 2 is adsorbed in an end-on configuration on Mo or V atoms.As depicted in Figure7a, for the MoCo@N 4 O 2 -G catalyst, the protonation reactions follow
It was found that the second N 2 molecule only physisorbed on the catalyst surface in the cases of MoCr@N 4 O 2 -G, VMn@N 4 O 2 -G, VFe@N 4 O 2 -G, and MoMn@N 4 O 2 -G, while the adsorption free energies of the second N 2 molecule on MoFe@N 4 O 2 -G, MoRe@N 4 O 2 -G, and MoCo@N 4 O 2 -G were −0.34, −0.69, and −0.41 eV, respectively.Subsequently, we calculated the free energy changes from *N 2 + *N 2 to *NNH + *N 2 on the latter three catalysts.On MoFe, MoRe, and MoCo, the free energy changes from *N 2 + *N 2 to *NNH(Mo) + *N 2 are 0.32, 0.35, and 0.32 eV, respectively, while the free energy changes from *N 2 + *N 2 to *N 2 + *NNH(M') are 1.18, 1.76, and 0.62 eV, respectively.The results indicate that on the surfaces of MoFe@N 4 O 2 -G and MoCo@N 4 O 2 -G, the ∆G of this step is larger than the rate-determining step's free energy obtained when a single N 2 molecule is adsorbed, indicating that simultaneously adsorbing two N 2 molecules is not feasible on these two surfaces (Table steps starting from the intermediate *NNH: *NNH + *N 2 → *NNH + *NNH and *NNH + *N 2 → *NNH 2 + *N 2 .The former displays a much higher free energy change (1.11 eV) compared to the latter (0.05 eV).Similarly, for the two elementary steps starting from the intermediate *NNH 2 , the ∆G of *NNH 2 + *N 2 → *NNH 2 + *NNH (1.14 eV) is much higher that of *NNH + *N 2 → *NNH 2 + *N 2 .Therefore, we infer that on MoRe@N 4 O 2 -G, the NRR continuously hydrogenates one N 2 molecule while suppressing the hydrogenation of another N 2 molecule.The corresponding reaction free energy diagram and optimized intermediate structures are depicted in Figure 8. Computational results indicate that the mixed pathway is the most feasible route and the PDS remains as *N + *N 2 + H + + e − → *NNH + *N 2 , with a U L of −0.35 V, significantly lower than the U L (−0.56 V) corresponding to the case when a single N 2 adsorbs on the surface.Consequently, it can be inferred that on MoRe@N 4 O 2 -G, the NRR is more inclined towards the adsorbing two N 2 molecules and follows a mixed mechanism.Molecules 2024, 29, x FOR PEER REVIEW 10 o tary steps starting from the intermediate *NNH: *NNH + *N2 → *NNH + *NNH *NNH + *N2 → *NNH2 + *N2.The former displays a much higher free energy change ( eV) compared to the latter (0.05 eV).Similarly, for the two elementary steps starting f the intermediate *NNH2, the ΔG of *NNH2 + *N2 → *NNH2 + *NNH (1.14 eV) is m higher that of *NNH + *N2 → *NNH2 + *N2.Therefore, we infer that on MoRe@N4O the NRR continuously hydrogenates one N2 molecule while suppressing the hydroge tion of another N2 molecule.The corresponding reaction free energy diagram and o mized intermediate structures are depicted in Figure 8. Computational results indi that the mixed pathway is the most feasible route and the PDS remains as *N + *N2 + H e -→ *NNH + *N2, with a UL of −0.35 V, significantly lower than the UL (−0.56 V) co sponding to the case when a single N2 adsorbs on the surface.Consequently, it can inferred that on MoRe@N4O2-G, the NRR is more inclined towards the adsorbing two molecules and follows a mixed mechanism.

Figure 8 .
Figure 8. Gibbs free energy diagram of NRR on MoRe@N4O2-G adsorption of two N2 molecu The C, N, O, H, Mo, and Re atoms are labeled as gray, blue, red, white, green, and dark green b respectively.

Figure 8 .
Figure 8. Gibbs free energy diagram of NRR on MoRe@N 4 O 2 -G adsorption of two N 2 molecules.The C, N, O, H, Mo, and Re atoms are labeled as gray, blue, red, white, green, and dark green balls, respectively.

Figure 9 .
Figure 9. (a) Charge density differences after N2 adsorption on MoCo@N4O2-G and Bade (Q*N2, |e|) after N2 adsorption.The C, N, O, Mo, and Co atoms are labeled as gray, blue, red and pink balls, respectively.The electrons accumulation and loss are represented by yel cyan areas.(b) PDOS before and (c) after N2 adsorption on MoCo@N4O2-G.The black das represents the Fermi energy level.

Figure 9 .
Figure 9. (a) Charge density differences after N 2 adsorption on MoCo@N 4 O 2 -G and Bader charge (Q *N2 , |e|) after N 2 adsorption.The C, N, O, Mo, and Co atoms are labeled as gray, blue, red, green, and pink balls, respectively.The electrons accumulation and loss are represented by yellow and cyan areas.(b) PDOS before and (c) after N 2 adsorption on MoCo@N 4 O 2 -G.The black dashed line represents the Fermi energy level.
A positive value of UL(NRR)−UL(HER) indicates that the catalyst favors the NRR over the HER, while a negative value indicates the opposite.The results of UL(NRR)−UL(HER) of eight DACs are shown in Figure 11.It can be seen that MoRe@N4O2-G, VMn@N4O2-G, and VFe@N4O2-G exhibit negative UL(NRR)−UL(HER) values, indicating poor selectivity for the NRR.Conversely, the remaining five catalysts (MM′@N4O2-G, MM′ = MoFe, MoMn, MoCo, MoCr, and VV) all have positive values for UL(NRR)−UL(HER).Notably, because four (MM′@N4O2-G, MM′ = MoFe, MoMn, MoCo, and MoCr) of these five catalysts adsorb N2 onto their surfaces via an end-on mode, the initial hydrogenation reaction between *N2 and H + can either generate the *NNH intermediate as discussed above or form the *N2 + *H intermediate through H + directly adsorbing onto another active site.To determine the more feasible intermediate, we compared the free energy changes of these two intermediates formed on these four catalysts (Table

Figure 10 .
Figure 10.Relationship between U L and E d on MM ′ @N 4 O 2 -G.

2. 4 .
NRR Selectivity of the MM'@N 4 O 2 -G Catalysts The main competing reaction during the NRR is the HER.This competition greatly influences the selectivity of the catalyst.The difference between U L (NRR) and U L (HER) is commonly used to assess the selectivity of a catalyst.A positive value of U L (NRR)−U L (HER) indicates that the catalyst favors the NRR over the HER, while a negative value indicates the opposite.The results of U L (NRR)−U L (HER) of eight DACs are shown in Figure 11.It can be seen that MoRe@N 4 O 2 -G, VMn@N 4 O 2 -G, and VFe@N 4 O 2 -G exhibit negative U L (NRR)−U L (HER) values, indicating poor selectivity for the NRR.Conversely, the remaining five catalysts (MM ′ @N 4 O 2 -G, MM ′ = MoFe, MoMn, MoCo, MoCr, and VV) all have positive values for U L (NRR)−U L (HER).Notably, because four (MM ′ @N 4 O 2 -G, MM ′ = MoFe, MoMn, MoCo, and MoCr) of these five catalysts adsorb N 2 onto their surfaces via an end-on mode, the initial hydrogenation reaction between *N 2 and H + can either generate the *NNH intermediate as discussed above or form the *N 2 + *H intermediate through H + directly adsorbing onto another active site.To determine the more feasible intermediate, we compared the free energy changes of these two intermediates formed on these four catalysts (Table

Figure 11 .
Figure 11.Difference between the limiting potential of NRR (U L (NRR)) and the limiting potential of HER (U L (HER)) on the 8 DACs (U L (NRR)−U L (HER)).