Tuning the Site-to-Site Interaction of Heteronuclear Diatom Catalysts MoTM/C2N (TM = 3d Transition Metal) for Electrochemical Ammonia Synthesis

Ammonia (NH3) synthesis is one of the most important catalytic reactions in energy and chemical fertilizer production, which is of great significance to the sustainable development of society and the economy. The electrochemical nitrogen reduction reaction (eNRR), especially when driven by renewable energy, is generally regarded as an energy-efficient and sustainable process to synthesize NH3 in ambient conditions. However, the performance of the electrocatalyst is far below expectations, with the lack of a high-efficiency catalyst being the main obstacle. Herein, by means of comprehensive spin-polarized density functional theory (DFT) computations, the catalytic performance of MoTM/C2N (TM = 3d transition metal) for use in eNRR was systematically evaluated. Among the results, MoFe/C2N can be considered the most promising catalyst due to its having the lowest limiting potential (−0.26 V) and high selectivity in the context of eNRR. Compared with its homonuclear counterparts, MoMo/C2N and FeFe/C2N, MoFe/C2N can balance the first protonation step and the sixth protonation step synergistically, showing outstanding activity regarding eNRR. Our work not only opens a new door to advancing sustainable NH3 production by tailoring the active sites of heteronuclear diatom catalysts but also promotes the design and production of novel low-cost and efficient nanocatalysts.


Introduction
Ammonia (NH 3 ) is one of the most productive inorganic compounds in the world, with an annual output of about 150 million tons. It is used for the manufacture of many important chemicals, particularly fertilizers [1]. At present, industrial NH 3 synthesis still adopts a high-temperature and high-pressure Haber-Bosch process invented more than a century ago [2]. The Haber-Bosch method of ammonia synthesis has successfully altered the history of food production, fed explosive population growth, and laid the foundation for heterogeneous catalysis and chemical engineering as well. However, this synthetic method consumes 1-2% of the world's energy supply and about 2% of the natural gas supply and also produces 1.44% of global carbon dioxide (CO 2 ) emissions [1,[3][4][5][6]. Therefore, it is necessary to find more sustainable and environmentally friendly alternatives as soon as possible to replace the energy-and carbon-intensive Haber-Bosch process [1,3,7,8]. Making full use of the abundant nitrogen (N 2 ) in the Earth's atmosphere and using water as the proton source, an electrochemical nitrogen reduction reaction (eNRR), especially when driven by renewable energy, is generally regarded as an energy-efficient and sustainable process to synthesize NH 3 in ambient conditions [9,10]. However, the performance of the electrocatalyst is far below expectations, and the lack of a high-efficiency catalyst is the main obstacle to large-scale NH 3 synthesis.

Geometric Structures, Stabilities, and Electronic Properties of the MoTM/C 2 N (TM = 3d Transition Metal)
First, the geometric structure of the C 2 N monolayer was optimized and the lattice parameters were calculated to be a = b = 8.32 Å, which matches well with the experimental data and previous DFT reports (Figure 1a) [30,34]. In the C 2 N monolayer, benzene rings are bridged by pyrazine rings to form six-membered nitrogen holes, in which the distance between two opposite nitrogen atoms is 5.51 Å [33].

Geometric Structures, Stabilities, and Electronic Properties of the MoTM/C2N (TM = 3d Transition Metal)
First, the geometric structure of the C2N monolayer was optimized and the lattice parameters were calculated to be a = b = 8.32 Å, which matches well with the experimental data and previous DFT reports (Figure 1a) [30,34]. In the C2N monolayer, benzene rings are bridged by pyrazine rings to form six-membered nitrogen holes, in which the distance between two opposite nitrogen atoms is 5.51 Å [33]. Then, we investigated the geometric structures and the stability of MoTM/C2N (TM = 3d transition metal). The optimized structures of MoTM/C2N are presented in Figure S1 and Table S1 in the Supplementary Materials. Clearly, the framework structure of C2N can be maintained well with the decoration of Mo and TM atoms. Except for the larger atomic radius of Sc, which slightly protrudes out of the C2N plane, all the other metal atoms can be embedded into the central cavity of C2N holes and form almost coplanar structures. The bond lengths of Mo-TM are in the range of 1.93 Å (for MoV/C2N) and 2.48 Å (for MoSc/C2N). Compared with the length of the metal bond in the corresponding bulk phase, Then, we investigated the geometric structures and the stability of MoTM/C 2 N (TM = 3d transition metal). The optimized structures of MoTM/C 2 N are presented in Figure S1 and Table S1 in the Supplementary Materials. Clearly, the framework structure of C 2 N can be maintained well with the decoration of Mo and TM atoms. Except for the larger atomic radius of Sc, which slightly protrudes out of the C 2 N plane, all the other metal atoms can be embedded into the central cavity of C 2 N holes and form almost coplanar structures. The bond lengths of Mo-TM are in the range of 1.93 Å (for MoV/C 2 N) and 2.48 Å (for MoSc/C 2 N). Compared with the length of the metal bond in the corresponding bulk phase, the bond lengths of Mo-TM in the system of MoTM/C 2 N are shorter. Such interconnections enable the heteronuclear diatom to respond cooperatively to the adsorbate, with communicative structural self-adaption and electronic transformation, which induces a different catalytic performance from that of its single-atom counterparts [3,27,36]. The binding energy (E b ) of MoTM and the Bader charge of anchored Mo and TM were calculated to examine structural stability, as shown in Table S1 in the Supplementary Materials. The smallest E b of MoTM on C 2 N is −9.10 eV, suggesting strong coupling and good stability. Meanwhile, the Bader charge analysis showed that the atomic charge of the anchored TM atoms (from Sc to Ni) decreased gradually from +2.32 to +0.43, which is due to the increasing electronegativity of TM atoms (from Sc to Ni). Conversely, the atomic charge of anchored Mo atoms gradually increased.
In order to evaluate the thermal stability of MoFe/C 2 N, we simulated MoFe/C 2 N using the ab initio molecular dynamics (AIMD) method. As shown in Figure 1c, the energy fluctuation oscillates around the equilibrium state according to time evolution. Even at 800 K, the structure of MoFe/C 2 N remained stable after the 100 ps simulation, showing high thermal stability. This high stability stems from the strong bonding between Mo and Fe, as well as between the MoFe and N atoms of the C 2 N substrate. In addition, the AIMD results also show that the anchored Mo and Fe atoms cannot escape from the holes in the substrate, thus excluding the problems of diffusion and agglomeration.

Adsorption of N 2
It is well known that the adsorption of reactants on the catalyst surface is the key step to initializing a reaction. In order to comprehensively study the adsorption of N 2 on the catalyst's surface, we considered seven possible adsorption configurations, in which N 2 is adsorbed on the top site of a single metal or bimetallic bridge site via the endon or side-on mode, as shown in Figure S2 in the Supplementary Materials. After full structural relaxation, the most favorable adsorption configurations of N 2 on MoTM/C 2 N are presented in Figure S3 in the Supplementary Materials. The adsorption energies of N 2 were calculated and are displayed in Figure 2a.
the bond lengths of Mo-TM in the system of MoTM/C2N are shorter. Such interconnections enable the heteronuclear diatom to respond cooperatively to the adsorbate, with communicative structural self-adaption and electronic transformation, which induces a different catalytic performance from that of its single-atom counterparts [3,27,36]. The binding energy (Eb) of MoTM and the Bader charge of anchored Mo and TM were calculated to examine structural stability, as shown in Table S1 in the Supplementary Materials. The smallest Eb of MoTM on C2N is −9.10 eV, suggesting strong coupling and good stability. Meanwhile, the Bader charge analysis showed that the atomic charge of the anchored TM atoms (from Sc to Ni) decreased gradually from +2.32 to +0.43, which is due to the increasing electronegativity of TM atoms (from Sc to Ni). Conversely, the atomic charge of anchored Mo atoms gradually increased.
In order to evaluate the thermal stability of MoFe/C2N, we simulated MoFe/C2N using the ab initio molecular dynamics (AIMD) method. As shown in Figure 1c, the energy fluctuation oscillates around the equilibrium state according to time evolution. Even at 800 K, the structure of MoFe/C2N remained stable after the 100 ps simulation, showing high thermal stability. This high stability stems from the strong bonding between Mo and Fe, as well as between the MoFe and N atoms of the C2N substrate. In addition, the AIMD results also show that the anchored Mo and Fe atoms cannot escape from the holes in the substrate, thus excluding the problems of diffusion and agglomeration.

Adsorption of N2
It is well known that the adsorption of reactants on the catalyst surface is the key step to initializing a reaction. In order to comprehensively study the adsorption of N2 on the catalyst's surface, we considered seven possible adsorption configurations, in which N2 is adsorbed on the top site of a single metal or bimetallic bridge site via the end-on or sideon mode, as shown in Figure S2 in the Supplementary Materials. After full structural relaxation, the most favorable adsorption configurations of N2 on MoTM/C2N are presented in Figure S3 in the Supplementary Materials. The adsorption energies of N2 were calculated and are displayed in Figure 2a.  From the point of view of adsorption energy, N 2 tends to be adsorbed onto MoTM/C 2 N (except MoSc/C 2 N) via the end-on mode. Because of the geometric constraint that the Sc atom protrudes out of the catalyst surface, N 2 is preferably adsorbed on MoSc/C 2 N via the side-on mode instead of the end-on mode. It can be clearly seen from Figure  is a bimetallic bridge site. The side-on adsorption configuration of N 2 on MoFe/C 2 N as a representative is displayed in Figure 2b. According to previous studies [27], the moderate adsorption strengths and small differences in N 2 adsorption energies between side-on and end-on adsorption configurations suggest that MoFe/C 2 N may be a good eNRR catalyst.
Besides the adsorption energy, the Bader charge, charge density difference, and bond length of the adsorbed N 2 were also measured and are shown in Figure S3 in the Supplementary Materials. Charge density difference plots indicate that the charge transfer between the adsorbed N 2 and the anchored MoTM is bidirectional, and charge accumulation and depletion can be observed on both sides. This phenomenon is in perfect accordance with the "push-pull" hypothesis, in which the anchored MoTM can "push" the d electrons into the antibonding orbitals of N 2 and simultaneously "pull" the lone-pair electrons from N 2 . More importantly, we found that there is a good correlation between the N-N bond length and the Bader charge of N 2 ; that is, the more negative the Bader charge of N 2 is, the longer the N-N bond length is. For the side-on adsorption configurations of N 2 on the bimetallic bridge site of MoTM/C 2 N (TM = Sc, V, Fe and Co), the Bader charge of N 2 decreases to between −0.82 |e| and −1.05 |e|, with their N-N bond lengths obviously elongated to 1.213~1.293 Å (1.12 Å in the free gas phase). In the case of the side-on adsorption configurations of N 2 on the top site of the Mo atom of the MoTM/C 2 N (TM = Cr and Ni), the Bader charge of N 2 decreases to between −0.46 |e| and −0.49 |e|, with their N-N bond lengths elongated to 1.176-1.178 Å. In comparison, the end-on adsorption configurations of N 2 on the top site of the Mo-atom of MoTM/C 2 N, the Bader charge of N 2 decreases to between −0.34 |e| and 0.42 |e| with their N-N bond lengths slightly elongated to 1.138~1.143 Å. As is consistent with previous reports [27,37,38], this means that the activation of adsorbed N 2 with the side-on mode is stronger than that with the end-on mode.

Catalytic Activity for eNRR
Next, we move to the assessment of the catalytic performance of MoTM/C 2 N for the reduction of N 2 into NH 3 . Four possible associative NRR catalytic mechanisms, namely, distal, alternative, enzymatic, and mixed mechanisms, are possible with MoTM/C 2 N catalysts, as is schematically illustrated in Figure 3a. In the case of end-on adsorption, the conversion of N 2 to NH 3 follows a distal or alternative or mixed mechanism. In the case of side-on adsorption, the conversion of N 2 to NH 3 follows either an enzymatic or mixed mechanism.
In principle, the intrinsic activity of electrocatalysts can be evaluated by limiting potential (U L ) or overpotential. Considering that overpotential can be affected by electrolyte conditions, we choose the limiting potential (U L ) as the evaluation standard in this paper. Among all the protonation steps, the protonation step with the largest positive Gibbs free energy difference is defined as the potential-determination step (PDS), while the lowest negative potential that promotes PDS to become exothermic is defined as the limiting potential (U L = −∆G PDS /e). As shown in Figure 3b, the U L of all catalysts MoTM/C 2 N (TM = 3d transition metal) were marked, and the catalytic activity of the widely reported Ru (0001) stepped surface (−0.98 V) was taken as a benchmark. As can be seen, the catalytic activity of all candidate catalysts of MoTM/C 2 N is superior to that of the Ru (0001) stepped surface.
Detailed Gibbs free energy diagrams of eNRR on MoTM/C 2 N are given in Figures  S4-S11. In the case of end-on adsorption, the PDS of eNRR is generally the first protonation step, that is, the generation of *NNH (where * stands for end-on adsorption), owing to the relatively low charge transfer from the substrate to N 2 . In the case of side-on adsorption, the PDS of the eNRR is the sixth protonation step, that is, the generation of the second ammonia, owing to the effective activation of N 2 . Among them, the limiting potential of eNRR on MoFe/C 2 N with side-on adsorption is the lowest (−0.26 V). Therefore, we now focus on analyzing the reaction path of eNRR on MoFe/C 2 N. S4-S11. In the case of end-on adsorption, the PDS of eNRR is generally the first protonation step, that is, the generation of *NNH (where * stands for end-on adsorption), owing to the relatively low charge transfer from the substrate to N2. In the case of side-on adsorption, the PDS of the eNRR is the sixth protonation step, that is, the generation of the second ammonia, owing to the effective activation of N2. Among them, the limiting potential of eNRR on MoFe/C2N with side-on adsorption is the lowest (−0.26 V). Therefore, we now focus on analyzing the reaction path of eNRR on MoFe/C2N.  The Gibbs free energy diagrams of eNRR on the MoFe/C2N catalyst with side-on adsorption are presented in Figure 3c. In the reduction process, ** NNH (where ** stands for side-on adsorption) is formed by the interaction between the **N2 and H + /e − pair, and the first protonation step (** N2 + H + + e − → ** NNH) is 0.22 eV uphill in the free energy profile. When ** NNH is protonated further, the H + /e − pair can continue to attack the same N atom to form ** NNH2 or attack another N atom to form ** NHNH. Comparing the Gibbs free energy changes of these two elementary steps, we found that it prefers to go through a step (** NNH + H + + e − → ** NHNH) with a Gibbs free energy uphill of 0.06 eV. Subsequently, the ** NHNH can be protonated to ** NH2NH easily, and this step (** NHNH + H + + e − → ** NH2NH) is exothermic at 0.02 eV. Once the ** NH2NH is formed, the first NH3 molecule can be readily desorbed from the catalyst surface, leaving * NH (the symbol * here only represents the adsorption state) on the catalyst surface and releasing 1.24 eV of heat. The fifth protonation step (* NH + H + + e − → * NH2) is exothermic at 0.52 eV. The sixth protonation step (* NH2 + H + + e − → * NH3) is endothermic at 0.26 eV. Among all these elementary steps, the sixth protonation step has the largest positive Gibbs free energy difference of 0.26 eV; therefore, it is considered the PDS. To sum up, the eNRR reaction on MoFe/C2N follows the mixed mechanism, and the UL is −0.26 V.
For comparison, we also investigated the catalytic performance of the homonuclear diatomic catalysts, MoMo/C2N and FeFe/C2N, with N2 side-on adsorption for eNRR, as shown in Figure 4. Note that the sixth protonation step (* NH2 + H + + e − → * NH3) on The Gibbs free energy diagrams of eNRR on the MoFe/C 2 N catalyst with side-on adsorption are presented in Figure 3c. In the reduction process, **NNH (where ** stands for side-on adsorption) is formed by the interaction between the **N 2 and H + /e − pair, and the first protonation step (**N 2 + H + + e − → **NNH) is 0.22 eV uphill in the free energy profile. When **NNH is protonated further, the H + /e − pair can continue to attack the same N atom to form **NNH 2 or attack another N atom to form **NHNH. Comparing the Gibbs free energy changes of these two elementary steps, we found that it prefers to go through a step (**NNH + H + + e − → **NHNH) with a Gibbs free energy uphill of 0.06 eV. Subsequently, the **NHNH can be protonated to **NH 2 NH easily, and this step (**NHNH + H + + e − → **NH 2 NH) is exothermic at 0.02 eV. Once the **NH 2 NH is formed, the first NH 3 molecule can be readily desorbed from the catalyst surface, leaving *NH (the symbol * here only represents the adsorption state) on the catalyst surface and releasing 1.24 eV of heat. The fifth protonation step (*NH + H + + e − → *NH 2 ) is exothermic at 0.52 eV. The sixth protonation step (*NH 2 + H + + e − → *NH 3 ) is endothermic at 0.26 eV. Among all these elementary steps, the sixth protonation step has the largest positive Gibbs free energy difference of 0.26 eV; therefore, it is considered the PDS. To sum up, the eNRR reaction on MoFe/C 2 N follows the mixed mechanism, and the U L is −0.26 V.
For comparison, we also investigated the catalytic performance of the homonuclear diatomic catalysts, MoMo/C 2 N and FeFe/C 2 N, with N 2 side-on adsorption for eNRR, as shown in Figure 4. Note that the sixth protonation step (*NH 2 + H + + e − → *NH 3 ) on MoMo/C 2 N is the PDS with a maximum free energy change of 0.77 eV, and the first protonation step (**N 2 + H + + e − → **NNH) on FeFe/C 2 N is the PDS with a maximum free energy change of 0.71 eV. As mentioned in the literature [38], the hydrogenation of *NH 2 species to *NH 3 is identified as the PDS on Mo-based catalysts, while the formation step of *NNH species is identified as the PDS on Fe-based catalysts. However, the energy changes of the first protonation step and the sixth protonation step on the heteronuclear diatomic catalyst MoFe/C 2 N are both small (0.22 eV and 0.26 eV, respectively). Notably, compared with MoMo/C 2 N and FeFe/C 2 N, MoFe/C 2 N can balance the first protonation step and the sixth protonation step synergistically, showing outstanding activity for eNRR. From this perspective, the heteronuclear diatomic catalyst is more effective than its homonuclear counterparts. Moreover, except for their high catalytic activity, MoFe-based catalysts also possess other merits, such as nontoxicity and a low price [39][40][41][42].
MoMo/C2N is the PDS with a maximum free energy change of 0.77 eV, and the first protonation step (** N2 + H + + e − → ** NNH) on FeFe/C2N is the PDS with a maximum free energy change of 0.71 eV. As mentioned in the literature [38], the hydrogenation of * NH2 species to * NH3 is identified as the PDS on Mo-based catalysts, while the formation step of * NNH species is identified as the PDS on Fe-based catalysts. However, the energy changes of the first protonation step and the sixth protonation step on the heteronuclear diatomic catalyst MoFe/C2N are both small (0.22 eV and 0.26 eV, respectively). Notably, compared with MoMo/C2N and FeFe/C2N, MoFe/C2N can balance the first protonation step and the sixth protonation step synergistically, showing outstanding activity for eNRR. From this perspective, the heteronuclear diatomic catalyst is more effective than its homonuclear counterparts. Moreover, except for their high catalytic activity, MoFebased catalysts also possess other merits, such as nontoxicity and a low price [39][40][41][42].

Selectivity of MoFe/C2N
As we all know, the hydrogen evolution reaction (HER) is the major competing reaction in eNRR and significantly affects the Faradaic efficiency of eNRR. It is necessary for an excellent eNRR catalyst to effectively suppress HER. In order to test the selectivity of MoFe/C2N for eNRR, we calculated the adsorption free energy of H on MoFe/C2N and used it to estimate the overpotential of HER. The *H structure of the adsorption site and

Selectivity of MoFe/C 2 N
As we all know, the hydrogen evolution reaction (HER) is the major competing reaction in eNRR and significantly affects the Faradaic efficiency of eNRR. It is necessary for an excellent eNRR catalyst to effectively suppress HER. In order to test the selectivity of MoFe/C 2 N for eNRR, we calculated the adsorption free energy of H on MoFe/C 2 N and used it to estimate the overpotential of HER. The *H structure of the adsorption site and corresponding free energy diagram are depicted in Figure 5a. The free energy barrier (0.45 eV) is considerably larger than the PDS barrier for eNRR (0.26 eV) on heteronuclear MoFe/C 2 N. Thus, MoFe/C 2 N exhibits an excellent suppression effect on HER during eNRR. In summary, as we expected, the heteronuclear diatom catalyst MoFe/C 2 N exhibits superior catalytic activity for eNRR, with a limiting potential of −0.26 V, and possesses high selectivity toward eNRR. corresponding free energy diagram are depicted in Figure 5a. The free energy barrier (0.45 eV) is considerably larger than the PDS barrier for eNRR (0.26 eV) on heteronuclear MoFe/C2N. Thus, MoFe/C2N exhibits an excellent suppression effect on HER during eNRR. In summary, as we expected, the heteronuclear diatom catalyst MoFe/C2N exhibits superior catalytic activity for eNRR, with a limiting potential of −0.26 V, and possesses high selectivity toward eNRR. Because the catalyst MoFe/C2N itself contains the N element, another problem to be solved is to exclude the nitrogen source of produced NH3 from the nitrogen-containing catalyst, that is, to judge whether the well-known Mars-van Krevelen (MvK) pathway [43,44] will occur on MoFe/C2N. It is worth noting that in the MvK pathway, a lattice N atom of the catalyst is reduced to NH3; the catalyst is subsequently regenerated with gaseous N2, in which the adsorption of *H species on the N site of the catalyst is the first and most critical step. In order to solve this problem, we examined the *H adsorption on the N site of the MoFe/C2N (Figure 5b). The result showed that the Ead of *H on the N site of MoFe/C2N is almost zero. Thus, the N site of MoFe/C2N is unlikely to be the active site for eNRR, implying that the reactant N2 molecule is the only N source for the NH3 product.

Electronic Structures
To reveal the origin of MoFe/C2N as an efficient catalyst for eNRR, the electronic structure was investigated. As can be seen from Figure 6a, before N2 adsorption on MoFe/C2N, the d orbital hybridized peaks between Mo and Fe are observed, which means that the bond of Mo and Fe is formed in the C2N substrate. After N2 adsorption with the side-on mode, the energy levels of MoFe d orbitals and N2 π* orbitals are well matched, leading to an anti-bonding orbital d-π*(unocc) above the Fermi level (located at 2.88 eV and 4.01 eV) and a bonding orbital d-π*(occ) below the Fermi level (located at −0.98 eV and −0.58 eV). This strong interaction between N2 and MoFe metal atoms leads to the formation of N-Mo and N-Fe, as well as the weakening of the N-N bond, which facilitates the subsequent reduction reaction of N2. Meanwhile, this is also clearly confirmed by the charge density difference diagram (Figure 6b). The charge accumulation (shown in yellow) Because the catalyst MoFe/C 2 N itself contains the N element, another problem to be solved is to exclude the nitrogen source of produced NH 3 from the nitrogen-containing catalyst, that is, to judge whether the well-known Mars-van Krevelen (MvK) pathway [43,44] will occur on MoFe/C 2 N. It is worth noting that in the MvK pathway, a lattice N atom of the catalyst is reduced to NH 3 ; the catalyst is subsequently regenerated with gaseous N 2 , in which the adsorption of *H species on the N site of the catalyst is the first and most critical step. In order to solve this problem, we examined the *H adsorption on the N site of the MoFe/C 2 N (Figure 5b). The result showed that the E ad of *H on the N site of MoFe/C 2 N is almost zero. Thus, the N site of MoFe/C 2 N is unlikely to be the active site for eNRR, implying that the reactant N 2 molecule is the only N source for the NH 3 product.

Electronic Structures
To reveal the origin of MoFe/C 2 N as an efficient catalyst for eNRR, the electronic structure was investigated. As can be seen from Figure 6a, before N 2 adsorption on MoFe/C 2 N, the d orbital hybridized peaks between Mo and Fe are observed, which means that the bond of Mo and Fe is formed in the C 2 N substrate. After N 2 adsorption with the side-on mode, the energy levels of MoFe d orbitals and N 2 π* orbitals are well matched, leading to an anti-bonding orbital d-π*(unocc) above the Fermi level (located at 2.88 eV and 4.01 eV) and a bonding orbital d-π*(occ) below the Fermi level (located at −0.98 eV and −0.58 eV). This strong interaction between N 2 and MoFe metal atoms leads to the formation of N-Mo and N-Fe, as well as the weakening of the N-N bond, which facilitates the subsequent reduction reaction of N 2 . Meanwhile, this is also clearly confirmed by the charge density difference diagram (Figure 6b). The charge accumulation (shown in yellow) is mainly distributed between N and Mo, as well as between N and Fe, indicating the formation of N-Mo bonds and N-Fe bonds. In contrast, charge depletion (shown in cyan) is mainly distributed in the N-N bond, indicating that the N-N bond is weakened.
Molecules 2023, 28, x FOR PEER REVIEW 10 of 14 is mainly distributed between N and Mo, as well as between N and Fe, indicating the formation of N-Mo bonds and N-Fe bonds. In contrast, charge depletion (shown in cyan) is mainly distributed in the N-N bond, indicating that the N-N bond is weakened.

Computational Methods
In this work, all the computations were carried out using the spin-polarized DFT method, which includes van der Waals (vdW) corrections, as implemented in the Vienna ab-initio simulation package (VASP) [45][46][47]. The projector-augmented-wave (PAW) method and the generalized gradient approximation (GGA) with a Perdew-Burke-Ernzerhof (PBE) exchange-correlation function were used [48][49][50]. Grimme's semiempirical DFT-D3 scheme of dispersion correction was employed to account for the weak interactions [51,52]. Lattice parameters were fully relaxed, and the crystal geometry was fully optimized before electron structure and total energy calculations. A plane-wave cutoff energy of 600 eV was set to evaluate the core-electron interaction, and the total energy and force converged within 10 −6 eV and 0.02 eV·Å −1 , respectively. The convergence standard for the calculation of the vibrational frequency was set at 10 −7 eV. The integration of the Brillouin zones was sampled with 2 × 2 × 1 Monkhorst-Pack meshes [53]. In order to prevent any artificial interaction between periodically repeated images, a vacuum space of 15 Å was inserted along the z direction to the monolayers. The atomic charge was computed using Bader's charge population analysis [54]. To evaluate the stability of the catalysts, ab initio molecular dynamics (AIMD) simulations were performed in an NVT ensemble. The AIMD simulations lasted 100 ps, with a time step of 2.0 fs, and the temperature was controlled using the Nosé-Hoover method [55].
The binding energy (Eb) of the MoTM atoms was computed from the following equation:

Computational Methods
In this work, all the computations were carried out using the spin-polarized DFT method, which includes van der Waals (vdW) corrections, as implemented in the Vienna abinitio simulation package (VASP) [45][46][47]. The projector-augmented-wave (PAW) method and the generalized gradient approximation (GGA) with a Perdew-Burke-Ernzerhof (PBE) exchange-correlation function were used [48][49][50]. Grimme's semiempirical DFT-D3 scheme of dispersion correction was employed to account for the weak interactions [51,52]. Lattice parameters were fully relaxed, and the crystal geometry was fully optimized before electron structure and total energy calculations. A plane-wave cutoff energy of 600 eV was set to evaluate the core-electron interaction, and the total energy and force converged within 10 −6 eV and 0.02 eV·Å −1 , respectively. The convergence standard for the calculation of the vibrational frequency was set at 10 −7 eV. The integration of the Brillouin zones was sampled with 2 × 2 × 1 Monkhorst-Pack meshes [53]. In order to prevent any artificial interaction between periodically repeated images, a vacuum space of 15 Å was inserted along the z direction to the monolayers. The atomic charge was computed using Bader's charge population analysis [54]. To evaluate the stability of the catalysts, ab initio molecular dynamics (AIMD) simulations were performed in an NVT ensemble. The AIMD simulations lasted 100 ps, with a time step of 2.0 fs, and the temperature was controlled using the Nosé-Hoover method [55].
The binding energy (E b ) of the MoTM atoms was computed from the following equation: where E MoTM/C2N , E C2N , E Mo , and E TM represent the total energy of MoTM/C 2 N, the energy of the substrate C 2 N, the energy of the gas-phase Mo atom, and the energy of the gas-phase TM atom, respectively. The adsorption energy (E ad ) of the adsorbates was computed according to the following equation: where E adsorbate + E MoTM/C2N , E MoTM/C2N , and E adsorbate represent the total energy of the catalyst MoTM/C 2 N and the adsorbate, the energy of the catalyst MoTM/C 2 N, and the energy of free adsorbate, respectively. According to this definition, the larger the negative value of E ad is, the more stable the adsorption state is. The computational hydrogen electrode (CHE) model proposed by Nørskov et al. was used to compute the free energy change of each elementary step for the N 2 reduction reaction [56]. According to the CHE model, the chemical potential of the H + /e − pair is equal to half that of H 2 with the standard conditions of pH = 0 and p(H 2 ) = 1 bar: 2 µ H 2 . The free energy change (∆G) of each protonation step was computed according to the following equation: where ∆E is the reaction energy obtained from DFT computations, ∆E ZPE is the change of zero-point energy, T is the temperature (298.15 K), and ∆S is the entropy change. The entropies of the gas phase were taken from the NIST database. The entropies and zero-point energy of the adsorbates were calculated based on their vibrational frequencies; the entropy can be determined as follows [57,58]: where R is the universal gas constant, h is Planck's constant, k B represents Boltzmann's constant, N A indicates Avogadro's number, ν i represents the frequency of the normal mode, and N is the number of adsorbed atoms. ∆G U could be computed using ∆G U = −eU, where e represents the charge transferred and U represents the potential at the electrode. ∆G pH is the correction of the H + free energy by the concentration and can be expressed as ∆G pH = k B T ln 10 × pH.

Conclusions
In order to obtain more efficient electrocatalysts with higher Faradaic efficiency, we have constructed a series of heteronuclear diatomic catalysts of MoTM/C 2 N (TM = 3d transition metal) for eNRR by tuning the interaction between Mo and TM. Then, the stability, activity, and selectivity of MoTM/C 2 N catalysts have been comprehensively evaluated using the spin-polarized DFT method. Among them, MoFe/C 2 N has been identified as the most promising catalyst, due to its low value of U L (−0.26 V) and high selectivity toward the eNRR. In order to reach a deeper understanding, we have compared MoFe/C 2 N with its homonuclear counterparts, MoMo/C 2 N and FeFe/C 2 N. The results show that MoFe/C 2 N can balance the first protonation step and the sixth protonation step synergistically, which is the origin of its excellent activity. We hope that our work will inspire more experimental and theoretical studies to further explore the potential of heteronuclear DACs in other chemical reactions.
Supplementary Materials: The following supporting information can be downloaded at: https: //www.mdpi.com/article/10.3390/molecules28104003/s1, Figure S1: Optimized structures of the most stable MoTM/C 2 N (TM = 3d transition metal). The C, N and Mo atoms are labeled as gray, blue and teal balls, respectively. Figure S2: Possible adsorption configurations diagram of N 2 on MoTM/C 2 N catalysts (take MoTi/C 2 N as an example). The C, N, Mo and Ti atoms are labeled as gray, blue, teal and red balls, respectively. Figure S3: Optimized adsorption configurations and charge density differences of N 2 chemisorbed on MoTM/C 2 N (TM = 3d transition metal). The charge accumulation and depletion were depicted by yellow and cyan, respectively. The isosurface value is 0.003 e/Å 3. Figure S4: Gibbs free energy diagrams for NRR on MoSc/C 2 N. Figure S5: Gibbs free energy diagram for NRR on MoTi/C 2 N. Figure S6: Gibbs free energy diagrams for NRR on MoV/C 2 N. Figure S7: Gibbs free energy diagrams for NRR on MoCr/C 2 N. Figure S8: Gibbs free energy diagram for NRR on MoMn/C 2 N. Figure S9: Gibbs free energy diagrams for NRR on MoFe/C 2 N. Figure S10: Gibbs free energy diagrams for NRR on MoCo/C 2 N. Figure S11: Gibbs free energy diagrams for NRR on MoNi/C 2 N. Figure S12: Gibbs free energy diagrams for NRR on MoMo/C 2 N and FeFe/C 2 N, respectively. Figure S13: Gibbs free energy diagram of HER on MoFe/C 2 N. Figure S14: The *H adsorption on the N site of MoFe/C 2 N. Table S1: The bond length of Mo-TM (d Mo-