Composition Engineering Opens an Avenue Toward Efficient and Sustainable Nitrogen Fixation

In this work, we open an avenue toward rational design of potential efficient catalysts for sustainable ammonia synthesis through composition engineering strategy by exploiting the synergistic effects among the active sites as exemplified by diatomic metals anchored graphdiyne via the combination of hierarchical high‐throughput screening, first‐principles calculations, and molecular dynamics simulations. Totally 43 highly efficient catalysts feature ultralow onset potentials (|Uonset| ≤ 0.40 V) with Rh‐Hf and Rh‐Ta showing negligible onset potentials of 0 and −0.04 V, respectively. Extremely high catalytic activities of Rh‐Hf and Rh‐Ta can be ascribed to the synergistic effects. When forming heteronuclears, the combinations of relatively weak (such as Rh) and relatively strong (such as Hf or Ta) components usually lead to the optimal strengths of adsorption Gibbs free energies of reaction intermediates. The origin can be ascribed to the mediate d‐band centers of Rh‐Hf and Rh‐Ta, which lead to the optimal adsorption strengths of intermediates, thereby bringing the high catalytic activities. Our work provides a new and general strategy toward the architecture of highly efficient catalysts not only for electrocatalytic nitrogen reduction reaction (eNRR) but also for other important reactions. We expect that our work will boost both experimental and theoretical efforts in this direction.


Introduction
[7] Therefore, it is highly desirable and quite urgent to explore ecofriendly and sustainable alternatives for green and energysaving ammonia production.[18][19] To realize electrochemical Haber-Bosch process, the innovation of catalysts is the key.
The recent emerging atomistic-level catalysts (ALC) show great perspective in diverse important chemical processes with high activity, stability, and selectivity. [20]Among many supports of ALC, two-dimensional (2D) graphdiyne (GDY) monosheets show many fantastic physicochemical properties and diverse applications in interdisciplinary fields. [21,22][25][26][27][28] Due to the cooperativity between different active sites, [29,30] modulation of the components will greatly affect the catalytic performance, thus becoming an important way to optimize and improve the activities of target catalysts.Therefore, as a simple and typical synergistic system, the GDY-based diatomic catalysts (DACs) provide an ideal platform to construct efficient catalysts for eNRR via composition engineering.][33][34] This is the motivation and also the propellant for us to launch the present study.
In this work, we open an avenue toward efficient and sustainable nitrogen fixation for ammonia synthesis through composition engineering strategy.The details are as follows, homonuclear diatoms (TM 2 ) were introduced into the hole of GDY for screening to get viable compositions of catalysts and to further learn and guide the screening of heteronuclear diatoms (TM-TM') with great efficiency and unexpected successful rate.Surprisingly, this strategy leads to the discovery of highly efficient Rh-Hf@ and Rh-Ta@GDY with ultralow onset potentials of 0 and −0.04 V, respectively.The closely relevant experimental evidence and clue come from the recently synthesized GDY-based atomic catalysts showing ultralow potentials for eNRR. [23,31,32]n this work, we open an avenue toward rational design of potential efficient catalysts for sustainable ammonia synthesis through composition engineering strategy by exploiting the synergistic effects among the active sites as exemplified by diatomic metals anchored graphdiyne via the combination of hierarchical high-throughput screening, first-principles calculations, and molecular dynamics simulations.Totally 43 highly efficient catalysts feature ultralow onset potentials (|U onset | ≤ 0.40 V) with Rh-Hf and Rh-Ta showing negligible onset potentials of 0 and −0.04 V, respectively.Extremely high catalytic activities of Rh-Hf and Rh-Ta can be ascribed to the synergistic effects.When forming heteronuclears, the combinations of relatively weak (such as Rh) and relatively strong (such as Hf or Ta) components usually lead to the optimal strengths of adsorption Gibbs free energies of reaction intermediates.The origin can be ascribed to the mediate d-band centers of Rh-Hf and Rh-Ta, which lead to the optimal adsorption strengths of intermediates, thereby bringing the high catalytic activities.Our work provides a new and general strategy toward the architecture of highly efficient catalysts not only for electrocatalytic nitrogen reduction reaction (eNRR) but also for other important reactions.We expect that our work will boost both experimental and theoretical efforts in this direction.
To our knowledge, this is the first reported case of a catalyst material (Rh-Hf@GDY) without onset potential for eNRR.Through this prototypical example, we want to convey the valuable messages of the importance and universality of composition engineering strategy in architecture of highly active catalysts for eNRR (and other reactions) by exploiting the synergistic effects realized through the combination of the different components.

Computational Details
The spin-polarized density functional theory (DFT) methods, implemented in the Vienna ab initio simulation package (VASP), [35] were employed for geometry optimization and electronic property calculations.The ion-electron interactions and interactions between electrons were described by the projector-augmented wave (PAW) [36] and the Perdew-Burke-Ernzerhof (PBE) [37] exchange-correlation functional within generalized gradient approximation (GGA), respectively.The van der Waals (vdW) interactions were accounted by DFT-D3 method. [38]The kinetic energy cut-off for the plane wave basis was set to 450 eV.A 2 × 2 supercell of GDY monolayer containing 72 C atoms was used as the substrate to anchor the transition metal atoms.The k-points were sampled using the Γ-centered 3 × 3 × 1 and 7 × 7 × 1 Monkhorst-Pack meshes [39] for structural optimization and density of states (DOS) calculations, respectively.The convergence criteria were 1 × 10 −5 eV for the total energy and 0.02 eV/ Å for the force on each atom, respectively.To eliminate the interactions between adjacent layers, the thickness of vacuum layer along crystal axis c was set to 16 Å.
The computational hydrogen electrode (CHE) model [40] was applied to compute the free energy profiles for eNRR.The Gibbs free energy change (ΔG) for each elementary step of eNRR can be calculated as follows: where ΔE is the electronic energy difference, ΔZPE, ΔTS and ∫C p dT are changes in zero-point energies, entropy, and enthalpy corrections, respectively.ΔG U = neU, where U is the applied electrode potential, n is the number of transferred electrons, e is the charge of an electron.ΔG pH = k B T × ln10 × pH, where k B is the Boltzmann constant and T = 298.15K. pH is set to 0 for strong acidic medium.The elementary step involving proton-electron pair transfer with the most positive ΔG (ΔG max ) is the potentialdetermining step (PDS), and the onset potential is calculated as U onset = −ΔG max /e.Due to the heavy computational demanding and small effects on the results (solvation effects may lower the theoretical overpotentials for individual materials by up to 0.1 eV), [41] and the recent experimental results on GDY-based atomic catalyst show very good catalytic performance with ultralow potential (−0.045 V), [23] which provide experimental evidence and clues to support our current results to some extent, giving us the confidence on the reliability of our calculated results without solvation effects, thus the solvation effects have not been considered in the present study.The values of entropy, enthalpy corrections, and zero-point energies of gas phase molecules (H 2 , N 2 , NH 3 ) are obtained from the NIST database. [42]The climbing image nudged elastic band (CI-NEB) method [43] with six images was used to identify the transition state (TS) structures of different migration pathways of transition metal atoms on GDY substrate.The minimum energy path was optimized using the force-based conjugate-gradient method until the maximum force was less than 0.05 eV/ Å. Ab initio molecular dynamics (AIMD) simulations were carried out to evaluate the thermal stability of the most efficient catalysts with NVT ensembles under 500 K for 20 ps with a time step of 2 fs, and the Nosé-Hoover thermostat method was used to control the temperature. [44]A 2 × 1 × 1 supercell containing 148 atoms was adopted for the simulations without any symmetry constraints.The bond overlap population (BOP) values of adsorbed N 2 molecules were computed with on-the-fly pseudopotentials as implemented in the CASTEP module of Materials Studio 2018. [45]

Results and Discussion
For the convenience of discussion and analysis, here we assume the meaning of the labels throughout the manuscript, TM 2 @GDY and TM-TM'@GDY represent homonuclear and heteronuclear diatomic catalysts, respectively.For brevity and convenience, TM 2 @GDY and TM-TM'@GDY are abbreviated as TM 2 and TM-TM', respectively.Both full name and abbreviation are flexibly used throughout the paper for convenience under certain circumstance.
The architecture and workflow of the research are illustrated in Figure 1.For the convenience of discussion and analysis along a logical idea of hierarchical high-throughput screening, we label each turning point of screening as set-i, i = 1, 2, 3, 4, 5, 6, 7, 8 to represent the number of catalyst candidates left before or after each screening.We assume set-1 as the number of homonuclear diatomic catalysts TM 2 in Figure 1c, so, set-1 = 30.Following this line, set-2 = 21 is the number of left candidates TM 2 after the first criteria (Figure 1d); then set-3 = 13 is the number of left candidates after the second criteria (Figure 1e); next set-4 = 2 is the number of eligible TM 2 catalysts (Figure 1f).It should be pointed out that during the procedures of hierarchical high-throughput screening, some candidate materials will be excluded in each step if they do not meet the screening criteria.Only when the candidate materials meet all the screening criteria, can they become eligible catalysts.
Now the set of homonuclear candidates is used to construct the heteronuclears TM-TM'.Based on the permutation and combination of homonuclears in set-2, we obtained set-5 = 210 candidate materials for heteronuclear bimetallic catalysts in Figure 1g.After screening processes with three criteria, we get set-6, -7, -8 equal to 205, 148, 55 shown in Figure 1h-j, respectively.
To facilitate the analysis and discussion, this section is divided into the following six parts, that is, Section 3.1-3.6.In Section 3.1, candidate catalysts are screened out through hierarchical high-throughput screening.Followed by the reaction mechanisms of the selected catalysts are investigated in Section 3.2.Then, the selectivity and stability are discussed in Sections 3.3 and 3.4, respectively.After that, the intrinsic properties of efficient catalysts are analyzed in detail in Section 3.5.Finally, in Section 3.6, the synergistic effects of heteronuclear diatomic catalysts (Rh-Hf and Rh-Ta) are analyzed to unveil the origin for their high activities.

Architecture and Hierarchical High-Throughput Screening of Diatomic Catalysts
This part is divided into homonuclear and heteronuclear diatomic catalysts.The homonuclear part is the basis.The screening results of homonuclears are used to guide the combinations of different components for the architecture of highly efficient heteronuclears.
To better understand the essence and uniqueness of our highthroughput computational screening method, and the difference between our method and the reported method in literature, [28,46] here, we add some explanations and comments.Our screening method is an intelligent method featuring self-learning and self-improving.The specific process of our screening method is as follows: we first consider the interactions between the homonuclear diatoms and the N 2 molecules, two types of results can be obtained: valid interactions (vi) and invalid interactions (ii), corresponding to the Gibbs free energy of N 2 adsorption less than 0 and larger than 0, respectively, that is, vi ↔ ΔG (*N 2 ) < 0, ii ↔ ΔG(*N 2 ) > 0. Based on the screening results of the interactions between homonuclear diatoms and N 2 molecules, we carried out the screening of valid interactions between heteronuclear diatoms and N 2 molecules.That is, the interactions of the heteronuclear diatomic systems composed of two metals with valid interactions with the N 2 molecule is still valid with ΔG(*N 2 ) < 0, that is, vi & vi → vi.Whereas, the interactions of the heteronuclear diatomic systems contains invalid components, then their combinations lead to invalid results, that is, valid and invalid combinations lead to invalid results (vi & ii → ii), invalid and invalid combinations also lead to invalid results (ii & ii → ii), both combinations lead to the Gibbs free energy of N 2 adsorption is larger than 0 (i.e., ΔG(*N 2 ) > 0).Through such intelligent screening judgment, many invalid combinations can be excluded, thereby greatly improving the screening efficiency.Moreover, previous high-throughput screening method was mainly used for single-atom catalysts (relatively simple systems), while this method was used for diatomic catalysts (rather complex systems).From single-atom catalysts to diatomic catalysts, formally, it is just a simple increase in number from 1 to 2, actually, it is not just a simple increase in quantity, but a qualitative leap, from a relatively simple system to a rather complex system due to much more possibilities of combinations and interactions and synergic effects.
Due to the huge workload, hierarchical high-throughput screening was adopted to improve the screening efficiency and reduce the computational resources and keep the balance between efficiency and accuracy.On the basis of literature and our experience, [47][48][49][50][51][52][53] the spontaneous nitrogen adsorption, as well as the first (*N 2 + H + + e − → *NNH) and last (*NH 2 + H + + e − → *NH 3 ) protonation steps are very likely to be the PDS of eNRR.These three steps have profound effects on the results.Therefore, the screening criteria are set as follows: 1) ΔG(*N 2 ) < 0 eV; 2) ΔG(*N 2 → *NNH) ≤ 0.40 eV; 3) ΔG(*NH 2 → *NH 3 ) ≤ 0.40 eV.Only those materials that meet simultaneously three screening criteria will be selected out as eligible candidates.
The screening results for TM 2 @GDY are shown in Figure 2 and Table S2, Supporting Information.Six possible nitrogen adsorption configurations (Figure S1, Supporting Information) were considered.Only the most stable side-on and end-on configurations were adopted for the following study.From Figure 2a, one can see that only TM 2 with TM = Sc ~V, Mn ~Ni, Y ~Tc, Lu ~Re are valid for side-on adsorption.For end-on adsorption mode, TM 2 with TM = Sc ~V, Mn ~Ni, Y ~Rh, Lu ~Ir are effective.The points of end-on mode for TM 2 (TM = Cu, Zn, Pd, Pt) are missing in Figure 2b as the initial end-on mode is transferred to side-on mode after structural optimization.
Next, we turn to the second screening criterion, and the corresponding results are depicted in Figure 2c,d.One can see that TM 2 (TM = Sc ~V, Mn, Y ~Mo, Lu ~Ta) with side-on mode as well as TM = Sc ~V, Mn ~Co, Y, Nb, Lu ~Ta with end-on mode are screened out, respectively.
Finally, we turn to the third screening criterion, and the results are shown in Figure 2e,f.None of the materials with side-on mode (Figure 2e) and only Fe 2 and Co 2 with end-on mode (Figure 2f) are screened out.
In short, only Fe 2 @ and Co 2 @GDY are selected out for homonuclear diatomic catalysts after three-step high-throughput screening, and they will be further investigated below.

Heteronuclear Diatomic Catalysts
From the above discussion, we know that only 21 types of homonuclear diatomic catalysts can adsorb N 2 spontaneously and favors the follow-up reduction reaction.They are Sc 2 ~V2 , Mn 2 ~Ni 2 , Y 2 ~Rh 2 , Lu 2 ~Ir 2 , as listed in set-2.We randomly pick out 2 metals from set-2, combine them and anchor into GDY to form 210 (C 2 21 ) heteronuclear bimetallic catalysts (Figure 1g).These TM-TM'@GDY will have high probabilities for the spontaneous nitrogen adsorption and be potentially good candidates for eNRR.
The same screening procedures and criteria as homonuclear diatomic catalysts are applied to those 210 heteronuclear diatomic catalysts.The obtained results (Figure S2 and Table S2 In short, totally 57 DACs (2 TM 2 and 55 TM-TM') are selected out from 240 possible candidates (30 TM 2 and 210 TM-TM') via highthroughput screening.In the following sections, we will focus on the selected 57 materials to explore the detailed reaction mechanisms.

Reaction Mechanisms for Selected Catalysts
It should be pointed out that the first and last protonation steps may not be the PDS in some cases.Under this circumstance, the full reaction pathway search is a must to explore various reaction channels and finally determine the reaction mechanism.So, we comprehensively considered various reaction pathways for the selected 57 DACs.To get highly efficient catalysts, we use ΔG max ≤ 0.40 eV as the criterion to screen the pathways, that is, if the maximum free energy change is larger than 0.4 eV for any elementary step of a specific channel in eNRR, then this pathway will be excluded.
Energy Environ.Mater.2024, 7, e12600 The results from full reaction pathway search are shown in Figure S3 and Table S3, Supporting Information.From there we can see that the ΔG max of excluded 9 catalysts (V-Re, Mn-Hf, Fe-Ta, Zr-Tc, Tc-Lu, Ta-Re, Ta-Ir, W-Os, W-Ir) are larger than 0. Obviously, the whole proton-electron transfer processes along this consecutive pathway are exothermic.Thus, Rh-Hf has no onset potential.This is unexpected and surprising.To our knowledge, this is the first example of eNRR catalyst without onset potential and has profound effects for the follow-up research.

Selectivity
[56] Therefore, the selectivity of eNRR vs HER on 48 screened catalysts are calculated and compared to reflect the catalytic performance.The onset potentials of eNRR and HER on 48 catalysts are plotted in Figure 4 and summarized in Table S5, Supporting Information.From Figure 4, one can see that the |U onset | of eNRR on most catalysts are smaller than that of HER, and only the |U onset | of HER on Sc-V, Sc-Fe, Mn-Mo, Mo-Ru, and Ru-Re are smaller or equal to that of eNRR.This indicates that eNRR on most catalysts are energetically more favorable and have higher probabilities than HER.43 diatomic catalysts are left for further studied.As the most promising catalysts, the required free energy input of Rh-Hf and Rh-Ta to initiate HER are 0.77 and 0.74 eV, respectively, which are much higher than the free energy change of PDS for eNRR, indicating their higher selectivity toward eNRR.In short, Rh-Hf and Rh-Ta not only possess excellent catalytic activity but also exhibit good selectivity.

Stability
The stability is an essential aspect to assess performance of the catalysts in practical applications.Therefore, we calculate the binding energy (E b ) between metal atoms and GDY and the cohesive energy (E c ) of bimetallic pairs for 43 screened catalysts to evaluate the stability by the equations: , where E catalyst , E GDY , E A , E B , A E bulk , and B E bulk are the energies of catalysts, substrate GDY, single metal A atom, single metal B atom, metal A crystal, and metal B crystal, respectively, m and n represent the number of metal atoms in the unit cells of A and B crystals, respectively.The results are shown in Figure 5 and Table S6, Supporting Information.The magnitudes of E b and E c of screened catalysts range from −6.98 to −12.87 and from −4.93 to −9.35 eV, respectively.As shown in Figure 5, E b below E c for all 43 catalysts, indicate that these catalysts are stable, and the bimetallic pairs can be firmly anchored on GDY rather than migration and aggregation.
Moreover, the CI-NEB method is used to identify the transition state structures of different migration pathways of transition metal atoms on GDY substrate, and the results are collected and presented in Figure 6.Let us first analyze the migration energy barrier of single-atom catalysts TM@GDY.We consider four TM (Mn, Rh, Hf, Ta) as our diatomic catalysts involves these 3 transition metal elements and also the recent experimental progress on Mn@GDY [23] for eNRR show promising catalytic performance.From Figure 6a, we can see that the migration barrier of Mn atom on GDY substrate of experimentally fabricated singleatom catalyst Mn@GDY is 1.33 eV, indicating the good stability.This is consistent with the experimental facts that the high durability of Mn@GDY for eNRR under experimental conditions.From Figure 6b-d, one can see that the migration barriers are 1.35, 2.24, 2.21 eV for Rh, Hf, Ta atoms on GDY substrate, respectively.All are larger than that of experimentally synthesized Mn@GDY.Compared with the experimentally made catalyst Mn@GDY, all the three catalysts TM@GDY Figure 3. Gibbs free energy diagrams and geometry structures of reaction intermediates along the most favorable reaction pathways a) on Rh-Hf@GDY and b) on Rh-Ta@GDY at zero (blue) and onset (red) potentials, respectively.
Energy Environ.Mater.2024, 7, e12600 (TM = Rh, Hf, Ta) are stable against metal atom migration on the substrate.Now, we turn to discuss and analyze the migration barriers of diatomic catalysts.For diatomic catalysts, there are more possible paths for atomic migration than that of single-atom catalysts, so we need to comprehensively consider all possible scenarios for a correct and unbiased evaluation of the stability.For TM-TM'@GDY, three possible migration pathways are considered, that is, (1) single-atom migration and (2) double-atom migration.Both (1) and ( 2) include two possibilities.For (1), the migration of TM or TM' in TM-TM'@GDY.For (2), two atoms migrate together (from one hole of GDY to another hole), two atoms migrate separately (one atom migrates to a certain hole and the other atom migrates to another hole).
Figure 6e-h show the migration barriers of a single atom in Rh-Hf@GDY and Rh-Ta@GDY.Figure 6e,f display the migration barriers of Hf and Rh in Rh-Hf@GDY, respectively.Figure 6g,h display the migration barriers of Ta and Rh in Rh-Ta@GDY, respectively.We can see the migration barriers of Hf and Rh in Rh-Hf@GDY are 3.10 and 1.76 eV, which are much higher than that of Hf (2.24 eV) in Hf@GDY and Rh (1.35 eV) in Rh@GDY, respectively.Similarly, the migration barriers of Rh and Ta in Rh-Ta@GDY are 3.45 and 3.63 eV, which are much higher than that of Rh (1.35 eV) in Rh@GDY and Ta (2.21 eV) in Ta@GDY, respectively.These indicate that the migration energy barrier of each metal atom in diatomic catalysts is increased (compare with the corresponding single-atom catalysts) due to the presence of another metal atom, indicating that the synergistic effect between the two metals is of great significance for increasing the migration energy barrier and stabilizing the catalyst.
Figure 6i-l show the migration barriers of two atoms in Rh-Hf@GDY and Rh-Ta@GDY.Figure 6i,j display the migration barriers of Hf and Rh migrate separately and together in Rh-Hf@GDY, respectively.Figure 6k,l display the migration barriers of Ta and Rh migrate separately and together in Rh-Ta@GDY, respectively.We can see that the energy barriers for two atoms to migrate separately are 3.31 and 3.70 eV in Rh-Hf@GDY and Rh-Ta@GDY, respectively.The energy barriers for two atoms to migrate together are 2.18 and 3.92 eV in Rh-Hf@GDY and Rh-Ta@GDY, respectively.In general, the diatomic migration energy barrier is higher than the single atomic migration energy barrier in the catalysts with the same component.
In summary, the calculated results of the energy barriers of atomic migration show that the migration energy barriers are sufficiently large to prevent the migration of metal atoms, so no atomic migration occurs, and transition metal atoms can be firmly anchored on the GDY substrate in the diatomic catalysts, and the catalysts show high stabilities.This is consistent with the available experimental observations as a series of transition metal atoms anchored GDY have been used as atomic catalysts for nitrogen reduction synthesis of ammonia and other chemical reactions (such as, CO 2 RR, OER, HER, ORR, etc).
Furthermore, the results from AIMD simulations with 20 ps duration are used to verify the thermal stabilities under 500 K.The evolution of the total energy and temperature versus the simulation time, as well as the geometry structures before and after simulations are displayed in Figure S6, Supporting Information.Although there are some slight oscillations in the total energy and temperature, and some wrinkles in the final geometry configurations, the framework and structural integrity are kept quite well without melting holes and dissociated atomic fragments.Bimetallic pairs are still firmly anchored on GDY substrate instead of migration and aggregation into large particles.All these  Energy Environ.Mater.2024, 7, e12600 demonstrate their good thermal stabilities.Based on the above analysis, we can conclude that Rh-Hf@GDY and Rh-Ta@GDY can be durable under practical reaction conditions.

Intrinsic Properties of 43 Selected Catalysts
The above sections demonstrate the stability and selectivity of screened catalysts.In this section, we focus on the intrinsic electronic properties of selected catalysts.
The Bader charge [57] distributions for different fragments of catalysts are shown in Figure 7. Figure 7a-g demonstrate the Bader charges distributed on *N 2 , TM, TM', and TM-TM' dimer (include special cases TM = TM' for Fe and Co), C 8 , C 64 , C 72 , respectively.From Figure 7a we can see that the adsorbed N 2 obtains electrons (0.27-0.90 |e|) from catalysts and carries the negative charges to attract protons for the subsequent reactions.From Figure 7b-d one can see that bimetallic pairs lose electrons (0.88-2.62 |e|) when anchored onto GDY.After N 2 adsorption, bimetallic pairs further lose electron and finally bear positive charges (1.16-2.93|e|).After nitrogen adsorption, the charge change on the C 8 fragment is larger than that of the C 64 fragment since C 8 is subjected to the large direct influence due to the connection to the metal dimers, while C 64 is far away from the metal pairs and suffers from relatively small influence.Figure 7g demonstrates the charge evolution of the C 72 fragment, which is the combined effects of C 8 and C 64 .
To further understand how the charge transfer affects the activity, we calculated and plotted the charge density differences (Figure S7, Supporting Information) of N 2 adsorbed on these selected catalysts.
From there we can clearly see the charge accumulation and depletion regions on bimetallic pairs and adsorbed N 2 .The charge depletion on adsorbed N 2 is attributed to the lone-pair electrons of N 2 filling into the empty d orbitals of bimetallic pairs, and the charge accumulation on N 2 can be ascribed to the back donation of the d electrons from bimetallic pairs into antibonding π* orbitals of adsorbed N 2 .Both contribute to activate N 2 molecule and weaken and elongate the inert N≡N triple bond, which is vital for the nitrogen reduction reaction.
To get some insights into the intrinsic properties, the linear correlations among the Bader charge, general physical quantity (Ω) of transition metals, bond length, and BOP of the adsorbed N 2 are established in Figure 8.The general physical quantity ) represent the electronegativity (Pauling scale), the number of valence electron, and the atomic radii of transition metals, respectively.Here, Ω can be used a new descriptor to establish the linear relations with charges located on metal dimers in Figure 8a.The linear correlations between the charges on adsorbed N 2 and the N-N bond length, and between the BOP value and the N-N bond length are depicted in Figure 8b,c, respectively.From Figure 8b, c, we can see that the more charge distributed on N 2 , the longer N-N bond length, the lower BOP value, thus the weaker N-N bond interaction and the more activation of N 2 .All the linear relations in Figure 8 shed insights on the origin of catalytic activities.
As a new degree of freedom, the calculated magnetic moments of catalysts before and after N 2 adsorption are shown in Figure S9 and Table S10, Supporting Information.We can see that the metal atoms always bear large magnetic moments both before and after N 2 adsorption and dominate the whole systems with negligible contribution from carbon atoms of GDY.Interestingly, the magnetic moments of both the Figure 6.Diverse paths and their barriers of migration of transition metal atoms on GDY substrate in different atomic catalysts.a-d) are the migration of TM in TM@GDY (TM = Mn, Rh, Hf, Ta), e-h) are the migration of single TM in Rh-Hf@GDY and Rh-Ta@GDY, i-l) are the migration of two TMs in Rh-Hf@GDY and Rh-Ta@GDY.
Energy Environ.Mater.2024, 7, e12600 metal atoms and the catalysts decrease after N 2 adsorption, which may be ascribed to charge transfer between N 2 and metal atoms.Additionally, the magnetic moments distributed on the adsorbed N 2 are small for most catalysts.To vividly view the distribution of magnetic moments, the spin density plots for N 2 adsorption on screened catalysts are depicted in Figure S10, Supporting Information.From there, we can clearly see the spin density mainly distributed on the metal dimers, which is the origin for them to become the active sites for eNRR.
Electronic structure calculations are used to gain more insights into catalysts.The calculated band structures and density of states are shown in Figure S11, Supporting Information.When anchoring bimetallic pairs into GDY, the electronic conductivity is improved as the semiconducting GDY is transformed into the metallic TM-TM'@GDY from the band structure computation results.This will promote the charge transfer and boost eNRR.The calculated PDOS demonstrate the obvious hybridizations between the d states of metal atoms and the 2p states of C 8 -segment, indicating the substantial interactions between metal atoms and GDY.This contributes to the stability of catalyst materials.After N 2 adsorption, there are significant overlap of electronic state between the d states of metal atoms and the 2p states of N 2 , indicating the appreciable interactions between metal atoms and N 2 .This benefits to the adsorption and activation of N 2 and promotes the follow-up nitrogen reduction reaction.
After much effort and attempt, we established the relationship between catalytic activities (U onset ) and adsorption energies of intermediates (*N) for those highly efficient catalysts with U onset ≤ 0.2 V. From Figure 9 we can see that a nice linear correlation (R 2 = 0.88) is uncovered.For future perspective, this will shed some insights into further study on these complex systems.

Synergistic Effects for Heteronuclear Diatomic Catalysts
The intrinsic electronic properties of selected catalysts were elucidated in the above section.Here, we focus on unveiling the origin of synergistic effects for Rh-Hf and Rh-Ta with the highest activities for eNRR.The relevant homonuclear diatomic catalysts (Rh 2 , Hf 2 and Ta 2 ) are listed for comparison.The full reaction pathway search is performed for eNRR on Rh-Hf, Rh-Ta, Hf-Ta, Rh 2 , Hf 2, and Ta 2 .Note that Hf-Ta combination was excluded during the screening process.The important and representative paths of (Rh-Hf, Rh-Ta) and (Rh 2 , Hf 2 , Ta 2 ) are shown in Figures 3  and 10, respectively.The remaining paths are displayed in Figure S14, Supporting Information.From Figures 3 and 10, we can see that the ΔG max values are 0, 0.04, 0.73, 0.44, and 0.59 eV for Rh-Hf, Rh-Ta, Rh 2 , Hf 2 , and Ta 2 , respectively.These indicates that the catalytic activity of heteronuclear combination (Rh-Hf and Rh-Ta) are much higher than that of homonuclears (Rh 2 , Hf 2 , Ta 2 ).We ascribe the sharp difference of activity between heteronuclear and homonuclears to synergistic effects.
The calculated d-band centers of five selected catalysts (Figures 3 and  10) are shown in Figure 11 with detailed values listed in Table S12, Supporting Information.The d-band center in this paper was obtained by the VASPKIT software processing. [58]For the convenience of discussion and analysis, we classify the five catalysts into three groups with each group has an identical component, that is, (Rh-Hf, Rh-Ta, Rh 2 ), (Rh-Hf, Hf 2 ), (Rh-Ta, Ta 2 ) shown in Figure 11a/a', b/b', c/c', respectively.Let us first analyze the first group (Rh-Hf, Rh-Ta, Rh 2 ).In this group, three catalysts have the same Rh.From Figure 11a, we can see that the order of d-band centers of three different components follows Rh < Ta < Hf from low (negative) to high (positive) energy with the d-band centers of the same components (Rh) locating in comparable region.According to d-band theory, an upshift of d-band center results in a stronger adsorption due to up-shifting the antibonding state, [59][60][61] which is consistent with the observed adsorption strength of N 2 .The corresponding binding strength (ΔG(*N 2 )) of nitrogen adsorption shown in Figure 11a' follow the order of Rh 2 < Rh-Ta < Rh-Hf.Now, we turn to discuss the second group (Rh-Hf, Hf 2 ).The calculated results are displayed in Figure 11b/b'.From Figure 11b, we can see that the d-band centers of the same components (Hf) in both Rh-Hf and Hf 2 are similar, whereas the d-band centers of different components Energy Environ.Mater.2024, 7, e12600 9 of 13 (Rh vs Hf) in (Rh-Hf, Hf 2 ) follows the order of Rh < Hf, indicating the weaker binding strength between adsorbates and Rh-Hf compared with that of Hf 2 .The relevant free energies of nitrogen adsorption are shown in Figure 11b' with ΔG (*N 2 ) of Hf 2 is more negative than that of Rh-Hf.
Finally, we turn to analyze the third group (Rh-Ta, Ta 2 ) in Figure 11c/c'.The situation of the third group is quite similar to that of the second group.The order of d-band centers of different component (Rh vs Ta) follows Rh < Ta, and the free energies of nitrogen adsorption of Ta 2 is more negative than that of Rh-Ta.
From the above discussion and analysis, we can get further insights for synergistic effects.For the homonuclears, both dband centers (ε d ) and free energy changes of nitrogen adsorption (ΔG(*N 2 )) follow the order of Rh 2 < Ta 2 < Hf 2 .The PDS for Rh 2 is the first protonation step, whereas the last protonation step becomes the PDS for both Hf 2 and Ta 2 , which may be attributed to the relatively weak binding strength of adsorbates on Rh 2 and the relatively strong binding strength of intermediates on Hf 2 and Ta 2 .The combinations of relatively weak (such as Rh) and relatively strong (such as Hf or Ta) components can lead to the optimal strengths.When forming heteronuclears, both Rh-Hf and Rh-Ta show the mediate d-band centers, which lead to the optimal adsorption strengths of intermediates, thus bringing the high catalytic activities.According to Sabatier principle, [62] a moderate binding strength between reaction intermediates and catalyst surface will give an optimal activity since they can balance the adsorption strength of the multiple reaction intermediates.This is the reason why Rh-Hf and Rh-Ta are the most efficient catalysts.
In order to intuitively understand the connections between catalytic activities and d-band centers, the relations of U onset versus the region of ε d are vividly displayed in Figure 12a.From Figure 12a, we can clearly see that the ε d of Rh-Rh locates at relatively low (negative) energy region, whereas the ε d of Ta-Ta and Hf-Hf locate at relatively high (positive) energy region.The combinations of relatively low and relatively high ε d lead to the mediate ε d for Rh-Hf and Rh-Ta.
To further deepen understanding the underlying mechanism of catalytic activities of selected catalysts, here, we define ξ as the unified scaling of the general energies for the discussion of the free energy change (ΔG) of intermediates.ξ min , ξ max , and ξ avg = (ξ min + ξ max )/2 represent the minimum, the maximum and the arithmetic mean value of general energy, respectively.For specific application of ξ for N   Energy Environ.Mater.2024, 7, e12600 can be obtained in Figure 12b.These observations agree quite well with the Sabatier principle.
In summary, the synergistic effects for heteronuclear diatomic catalysts have been comprehensively analyzed.The combinations of a relatively weak and a relatively strong components for binding intermediates will have high probabilities to achieve the optimal catalytic activities.More generally, the present study can be viewed as a special case of universal composition engineering, which will open a new avenue toward efficient and sustainable nitrogen fixation under ambient conditions.

Experimental Section
Detailed information related to the synthesis of active electrodes, physicochemical characterization, and electrochemical evaluation of bifunctional electrodes towards UOR and supercapacitor application is provided in Supporting Information.

Conclusion
In this work, we systematically investigated 240 combinations of bimetallic pairs anchored on GDY as electrocatalysts (TM 2 @GDY and TM-TM'@GDY) for eNRR by composition engineering, hierarchical highthroughput screening, and first-principles calculations.43 combinations are selected out as promising candidates with high stability, activity (|U onset | ≤ 0.40 V) and selectivity (suppress HER).Electronic structure analyses indicate that the bond lengths of adsorbed N 2 show good linear correlation with the BOP values and the charges distributed on N 2 .Surprisingly, Rh-Hf@GDY and Rh-Ta@GDY display ultralow onset potentials of 0 and −0.04 V, respectively.The underlying mechanism of excellent performance can be attributed to the synergistic effects of bimetallic pairs.The combinations of a relatively weak (d-band centers in low energy region) and a relatively strong (d-band centers in high energy region) components can balance the adsorption strength of multiple reaction intermediates, thereby having high probabilities to  Energy Environ.Mater.2024, 7, e12600 achieve the optimal activities.More generally, the present study can be viewed as a special case of universal composition engineering strategy, which will open a new avenue toward efficient and sustainable ammonia production under ambient conditions and inspire future experimental efforts in this direction.
Figure 3a.From Figure 3a, we can see that the free energy change of nitrogen adsorption on Rh-Hf through side-on mode is −0.27 eV.Then, a proton-electron pair attacks one N atom of *N 2 to form *NNH, which is exothermic (−0.01 eV).When another proton-electron pair approaches to the same N atom, *NNH 2 is generated with the free energy change of −0.02 eV.Subsequently, *NNH 2 is reduced to *N and releases NH 3 (g) with the downhill free energy of −0.01 eV.After that, *N is gradually protonated to form *NH, *NH 2 , and *NH 3 with free energy changes of −1.72, −0.20, and −0.15 eV, respectively.Obviously, the whole proton-electron transfer processes along this consecutive pathway are exothermic.Thus, Rh-Hf has no onset potential.This is unexpected and surprising.To our knowledge, this is the first example of eNRR catalyst without onset potential and has profound effects for the follow-up research.

Figure 2 .
Figure 2. Schematic representation of the processes (right) and results (left) of high-throughput screening for TM 2 @GDY; a, b) represent the results from ΔG(*N 2 ) < 0 for both side-on and end-on adsorption modes, c, d) represent the results from ΔG(*N 2 → *NNH) ≤ 0.4 eV for both side-on and end-on configurations; e, f) represent the results from ΔG(*NH 2 → *NH 3 ) ≤ 0.4 eV for both side-on and end-on modes.All the points below the dashed lines in the figure are valid results.

Figure 5 .
Figure 5. Binding energy and cohesive energy of 43 screened diatomic catalysts.

Figure 4 .
Figure 4. Onset potentials for eNRR and HER on selected 48 diatomic catalysts.

Figure 7 .
Figure 7. Bader charge distributions a) on the adsorbed N 2 , and different segments b) TM, c) TM', d) TM-TM' dimer, e) C 8 , f) C 64 , g) C 72 before and after N 2 adsorption.C 8 , C 64 , and C 72 stand for the key eight carbon atoms that coordinate to the two metal atoms, the rest 64 carbon atoms, and total carbon atoms of GDY, respectively.

Figure 8 .
Figure 8. Linear relations a) between the charges of TM-TM' and general physical quantity Ω b) between the charge of adsorbed N 2 and the N-N bond length and c) between the BOP of adsorbed N 2 and the N-N bond length.

Figure 9 .
Figure 9. Linear correlation between adsorption energy of N atom and onset potential for those highly efficient catalysts (U onset ≤ 0.2 V).