Machine learning-assisted Ru-N bond regulation for ammonia synthesis

Machine learning prediction of Ru-IMCs

The workflow for designing efficient Ru-based intermetallic catalysts for ammonia synthesis was outlined in Fig. 1, comprising three main components: materials theoretical prediction, experimental verification, and reaction mechanism analysis. Firstly, a foundational dataset containing 201 Ru-based binary intermetallics was curated from the Inorganic Crystal Structure Database (ICSD) for ML models, which are listed in Supplementary Table 1. Training of the models involved using eXtreme Gradient Boosting (XGB) algorithms for regression predictive modeling. The regression model demonstrating the best performance was then leveraged to predict unexplored materials with promoted activity for ammonia synthesis. The SHapley Additive exPlanations (SHAP) value analysis, an interpretable ML technique, was employed to identify the parameters for optimizing the nitrogen interaction, which was crucial as the adsorption energies of N₂ (E_N2) and N atoms (E_N) were pivotal activity indicators, playing a key role in the N₂ activation, and N atoms hydrogenation in the subsequent ammonia formation pathway (Supplementary Fig. 1). To validate the predictions, materials preparations, structural characterization, and catalytic activity evaluations were conducted on the predicted Ru-IMCs materials. Finally, in-situ characterization and DFT calculations were employed to uncover insights into the reaction mechanism enhancing catalytic activities in ammonia synthesis.

For ML of Ru-based binary intermetallics, each case was described using multidimensional vectors that considered both element-specific and chemical environment features of Ru active sites. Here, the composition of the IMCs is characterized by the normalized stoichiometric ratio of the elements present in the material (denoted as NSR). The chemical environment of the active sites is defined by the Miller indices (denoted as (h k l)), the total number of atoms in the sites (denoted as site atoms), and the Ru ratio of the sites (denoted as Ru ratio). Furthermore, the coordination specificity of the elements is determined by the variance in the first ionization energy in comparison to Ru (denoted as Δ(X-Ru)) and the radius of the coordinating elements (denoted as r(X)). Here, the selected features possessed high irrelevance in heat map depicted Fig. 2a, indicating they independently represent an aspect of the atomic and electronic structure of the chemical system while being correlated and independent. To augment the prediction accuracy of ML, we evaluated the performance of multiple regression algorithms, including XGB, Multi-Layer Perception (MLP), K Nearest Neighbors (KNN), Gradient Boosting Decision Tree (GBDT), and Support Vector Regression (SVR) models (Fig. 2b). The obtained adsorption energy data of E_N2 and E_N was randomly shuffled and divided, with 80% for the training model and the remaining 20% for test model (Fig. 2c). The difference between the predicted and calculated values was defined as the error to assess the models’ performance. As depicted in the violin plot (Fig. 2b), the XGB model displayed the most competitive prediction accuracy for error values and probability density for the adsorption energy of both N₂ and N atoms among the five models, with details of the predictions of different models displayed in Supplementary Figs. 2–5. Additional evaluation metrics for model accuracy were summarized in Supplementary Tables 2 and 3. Across all metrics, including MAE (Mean Absolute Error), RMSE (Root Mean Square Error), and R² (R-square), the XGB model consistently demonstrated the best performance and can be utilized for the ML algorithms.

Next, a volcano plot was generated using the XGB model algorithm to evaluate the catalytic activity of various Ru-IMCs catalysts (Fig. 2d). All predicted E_N/E_N2 values for the Ru-IMC catalysts are listed in Supplementary Table 1 and Supplementary Fig. 6. Considering the elementary steps of ammonia synthesis, the energy barriers of each step can be depicted by the adsorption energies of N₂ and N atoms. Subsequently, we calculated their adsorption energies on different surfaces of Ru-IMCs through high-throughput methods, represented by the gray square points in Fig. 2d. To validate the accuracy of these data points, a selection of them was prepared, and their catalytic activities for ammonia synthesis were measured under conditions of 0.1 MPa and 400 °C. The reaction rates were determined for YRu₂, LaRu₂, and Ga₃Ru, showing the order of Ga₃Ru (N. D.) 2 (0.25 mmol g⁻¹ h⁻¹) < YRu₂ (0.57 mmol g⁻¹ h⁻¹) (Supplementary Fig. 7), aligning well with the trend observed in the volcano plot (Fig. 2d). Additionally, in Fig. 2d, the positions corresponding to Ru-IMCs represented by the blue circles were identified through ML screening, with NdRu₂ emerging as a promising candidate for ammonia synthesis due to its proximity to the apex of the volcano plot. Next, we obtained NdRu₂ pure phase through arc-melting and further tested its catalytic performance (Supplementary Fig. 8). The ammonia production rate of NdRu₂ is as high as 0.86 mmol g⁻¹ h⁻¹, surpassing that of YRu₂, LaRu₂, and Ga₃Ru, which is in line with the predicted trend of the volcano plot, but it still does not quite reach the peak of the volcano plot (Fig. 2d). Referring to Fig. 2d, focusing on NdRu₂, enhancing N₂ adsorption without altering the interaction of N atoms too much could bring us closer to the volcano plot’s apex. Typically, element doping serves as an effective method to tailor the properties of the host material, potentially modifying the surface adsorption strength of the reactant molecules. Nevertheless, in the absence of clear guidance, the options for doping can become overwhelming and arbitrary.

In order to uncover the underlying chemical patterns within the model and facilitate optimization, we utilized the SHAP value method to delve into the complex relationship between features and adsorption energies. A negative SHAP value for a feature indicates an enhancement on E_N or E_N2 through the XGB model’s output. As illustrated in Fig. 2e, instances of robust N atom adsorption tendencies were typically linked to high feature values of NSR and atoms, but low feature values of r(X) and Ru ratio. Conversely, cases involving N₂ adsorption usually demonstrated elevated feature values of Δ(X-Ru) and Ru ratio, alongside lower feature values of r(X). By combining the feature importance values outlined in Fig. 2f and Supplementary Table 4, replacing an element with a lower r(X) at the Nd site could enhance N₂ adsorption while minimally impacting the strength of N atoms, thereby potentially enhancing the catalytic performance of NdRu₂.

Preparation and characterization of Ru-IMCs

By using the arc-melting method, we successfully synthesized the Laves phase NdRu₂. In this structure, Ru atoms occupy sites with a coordination number (CN) of 12, forming a Frank-Kasper polyhedron. The Nd atoms, on the other hand, have a CN of 16. This well-defined crystal structure and the combination of elements spanning the periodic table provide a unique opportunity to tune the properties of Ru by modifying the Nd site. Since Nd is a rare earth element, we use a light rare earth element, Sc as the doping element and designed a novel Sc-doped NdRu₂ intermetallic catalyst. Here we have optimized the molar ratio of Sc doping to be 12.5% (Supplementary Fig. 9), that can be denoted as Sc_1/8Nd_7/8Ru₂. Given a smaller radius of Sc (160 pm) compared to Nd (185 pm), the substitution of Nd atoms with Sc is expected to induce lattice contraction of NdRu₂. As shown in Fig. 3b, Rietveld refinement was conducted to precisely verify the expected lattice change observed in the X-ray diffraction (XRD) patterns. Upon doping with Sc, the lattice parameters of Sc_1/8Nd_7/8Ru₂ are smaller than those of NdRu₂ (contracting from 7.613 Å to 7.573 Å), consistent with the structural relaxation predicted by our DFT calculations (Fig. 3a). The Extended X-ray Absorption Fine Structure (EXAFS) of Ru K-edge as depicted in Fig. 3c reveals a slight shift of the first nearest-coordination peaks in R space of NdRu₂ and Sc_1/8Nd_7/8Ru₂. The curve fitting results of all three samples are listed in Supplementary Figs. 10 and 11 and Supplementary Tables 5 and 6, which demonstrates that the specific Ru−Nd bond length decreases from 2.55 Å (NdRu₂) to 2.48 Å (Sc_1/8Nd_7/8Ru₂), further confirming the lattice contraction induced by Sc incorporation.

The electronic structure of NdRu₂ and Sc_1/8Nd_7/8Ru₂ was analyzed using X-ray absorption near-edge structure (XANES) and X-ray photoelectron spectroscopy (XPS). Illustrated in Fig. 3d, the absorption pre-edge positions follow the sequence of Ru foil > Sc_1/8Nd_7/8Ru₂ > NdRu₂, indicating the negatively charged nature of the Ru species in both NdRu₂ and Sc_1/8Nd_7/8Ru₂, with the electron density of Ru being slightly lower in Sc_1/8Nd_7/8Ru₂ than that in NdRu₂. Figure 3e, f depict the XPS Ru 3d spectra with respect to the binding energy of Sc_1/8Nd_7/8Ru₂ and NdRu₂. The noticeable shifts of typical peaks to lower energy also indicate negatively charged feature of surface Ru sites, which is consistent with the XANES results. While the Nd 3 d peaks exhibit a significant shift to higher energy compared to Nd⁰ 3d (Supplementary Fig. 12). These findings suggest that Nd and Sc donate electrons while Ru accepts electrons in NdRu₂ and Sc_1/8Nd_7/8Ru₂ systems. The enrichment of electrons on the surface Ru highlights the potential for promoting the back-donation of electrons to the antibonding π-orbitals of adsorbed N₂.

Transmission electron microscopy (TEM) was employed to investigate the microstructure of NdRu₂ and Sc_1/8Nd_7/8Ru₂. The elemental distribution in Sc_1/8Nd_7/8Ru₂ was assessed using Energy Dispersive Spectroscopy (EDS) mapping, depicted in the middle of Fig. 3g, revealing the uniform dispersion of Nd, Sc, and Ru elements in Sc_1/8Nd_7/8Ru₂ particle. On the right side of Fig. 3g, a distinct lattice corresponding to the (3 1 1) index was evident, consistent with the Fourier transform pattern in the inset. Additional crystal plane orientations can be observed in Supplementary Fig. 13. The interplanar spacings of NdRu₂ and Sc_1/8Nd_7/8Ru₂ measured by TEM are summarized in Fig. 3h. It is apparent that following Sc doping, all interplanar spacings demonstrate a contraction of 4.3–7.4%, aligning with theoretical calculations (Fig. 3a), XRD patterns (Fig. 3b), and EXAFS results (Fig. 3c).

Ammonia synthesis performance of Ru-IMCs

We successfully doped NdRu₂ by Sc with a smaller ion radius to create a new Ru-IMCs, Sc_1/8Nd_7/8Ru₂, and sought to evaluate its performance in catalytic ammonia synthesis. In Fig. 4a, the catalytic activity is depicted as a function of temperature over Sc_1/8Nd_7/8Ru₂, NdRu₂, and Ru metal. The ammonia production rate of Sc_1/8Nd_7/8Ru₂ reached 2.77 mmol g⁻¹ h⁻¹, which is 3 times higher than that of NdRu₂ (0.86 mmol g⁻¹ h⁻¹) at 400 °C and 0.1 MPa, while pure Ru metal exhibited negligible activity. Note that the catalytic activity of Sc_1/8Nd_7/8Ru₂ surpasses that of other intermetallic catalysts investigated under similar reaction conditions, such as LaRuSi (1.76 mmol g⁻¹ h⁻¹), LaCoSi (1.25 mmol g⁻¹ h⁻¹), and YRu₂ (0.57 mmol g⁻¹ h⁻¹) (Fig. 4b). The activities of Ru/MgO and Ru/AC catalysts are also referenced as typical Ru-based catalysts; they exhibited lower activities for ammonia production, at 0.80 mmol g⁻¹ h⁻¹ and 0.77 mmol g⁻¹ h⁻¹ (Supplementary Fig. 14), respectively, significantly less than that of the Sc_1/8Nd_7/8Ru₂ catalyst. Despite having a relatively low surface area of only 0.34 m² g⁻¹, Sc_1/8Nd_7/8Ru₂ yields a remarkable reaction rate of 8.18 mmol m⁻² h⁻¹ when assessed by specific activity, positioning it as one of the most active binary Ru-based intermetallic catalysts (Fig. 4b and Supplementary Tables 7 and 8).

In the earlier section on ML, we noted that the small radius of the coordinating atom in Ru-IMCs would notably improve their N₂ adsorption, thereby enhancing the activity of ammonia synthesis. Here, if we introduce the element Y with a slightly larger atomic radius (180 pm) than Sc and incorporate it into NdRu₂ at the same ratio, the resulting Y_1/8Nd_7/8Ru₂ catalyst displays a catalytic activity of 0.94 mmol g⁻¹ h⁻¹ for ammonia production (Supplementary Fig. 15). This places it between Sc_1/8Nd_7/8Ru₂ (2.77 mmol g⁻¹ h⁻¹) and NdRu₂ (0.86 mmol g⁻¹ h⁻¹), underscoring the significant impact of the doping atomic radius on enhancing the activity of NdRu₂. The outstanding catalytic performance observed with small atomic radius element doping in NdRu₂ for ammonia synthesis aligns with the trend depicted in the volcano plot in Fig. 2d, which emphasizes our methodology of utilizing the XGB algorithm for screening and optimizing Ru-IMCs catalysts for ammonia synthesis through ML as a reliable and effective approach.

The kinetics of ammonia synthesis over NdRu₂ and Sc_1/8Nd_7/8Ru₂ catalysts were analyzed using the ammonia production rates (r) at different reaction temperatures. A series of Arrhenius plots were generated by plotting lnr against 1/T. The slope of the Arrhenius plot allowed for the calculation of the apparent activation energy (E_a) for Sc_1/8Nd_7/8Ru₂, resulting in a value of 38.7 kJ·mol⁻¹ (Fig. 4c). This represents a reduction of approximately 30% compared to the E_a for the undoped NdRu₂ catalysts (58.9 kJ·mol⁻¹), providing evidence that Sc doping significantly enhances the ammonia synthesis reaction over NdRu₂. The suppressed E_a is attributed to Sc doping, which improves the adsorption of N₂, thereby reducing the energy barrier for N₂ dissociation on Sc_1/8Nd_7/8Ru₂. This hypothesis is initially supported by the N₂ reaction orders, with the N₂ reaction order for undoped NdRu₂ being close to unity (α = 0.94). While, the N₂ reaction order was measured to be 0.77 over Sc_1/8Nd_7/8Ru₂ (Fig. 4d), indicating higher efficiency of N₂ dissociation on Sc doped NdRu₂. Compared to NdRu₂, Sc_1/8Nd_7/8Ru₂ not only exhibits a lower activation energy but also demonstrates a fundamental improvement in resistance to hydrogen poisoning. Figure 4d illustrates the H₂ reaction orders, with the undoped NdRu₂ catalyst showing a negative value of −0.22, which demonstrates a strong hydrogen adsorption on its surface, named hydrogen poisoning effect. Conversely, the H₂ reaction order of Sc_1/8Nd_7/8Ru₂ is positive, reaching 0.79 (Fig. 4d), indicating that hydrogen poisoning has been mitigated by Sc doping. This substantial difference results in markedly different catalytic performance under high-pressure conditions, where the reaction rate over Sc_1/8Nd_7/8Ru₂ shows an approximately linear response to pressure increase. This led to a fourfold increase to 9.03 mmol g⁻¹ h⁻¹ when the reaction pressure rose from 0.1 to 0.9 MPa (Fig. 4e). Such enhancement of catalytic activity with pressure is not observed for NdRu₂, as the catalytic performance remains independent of total pressure, probably attributed to the hydrogen poisoning effect over NdRu₂. We further investigated the high-pressure responsiveness of Sc_1/8Nd_7/8Ru₂. As shown in Supplementary Fig. 16, the reaction rate increased monotonically to 15.85 mmol g⁻¹ h⁻¹ at 2.0 MPa, indicating a positive pressure effect for Sc_1/8Nd_7/8Ru₂ in ammonia synthesis. The exceptional stability of Sc_1/8Nd_7/8Ru₂ is evidenced by its sustained activity over 105 h at 400 °C under pressures of 0.1 MPa and 0.9 MPa, respectively (Fig. 4f). XRD patterns show the lattice parameters of the used Sc_1/8Nd_7/8Ru₂ sample are slightly larger than those of the fresh sample (increasing from 7.573 to 7.584 Å, Supplementary Figs. 17 and 18), likely due to the incorporation of hydrogen, which will be discussed below. Repeatable ammonia synthesis activity was observed for different batches of Sc_1/8Nd_7/8Ru₂ over an 80-h period (Supplementary Fig. 19).

Mechanism studies of Ru-IMCs for ammonia synthesis

To explore the origin of high performance of Sc_1/8Nd_7/8Ru₂, we conducted computational studies using first-principles calculations. Bader charge analysis indicates that the valence states of Ru in NdRu₂ and Sc_1/8Nd_7/8Ru₂ are −0.66 and −0.61, respectively, signifying electron transfer from Nd and/or Sc to Ru in the lattice of these intermetallic compounds (Supplementary Table 9), consistent with the experimental results in XANES and XPS (Fig. 3d–f). It is important to note that the electronic structure analysis from the XANES and XPS results demonstrates that the electron density of Ru in Sc_1/8Nd_7/8Ru₂ is lower than that in NdRu₂ (Fig. 3d–f), yet N₂ activation is significantly enhanced over Sc_1/8Nd_7/8Ru₂ (Fig. 4d). This observation contradicts the well-known assumption that increasing the electron density of Ru active sites would enhance electron back donation for N₂ dissociation and subsequently ammonia production.

Motivated by the aforementioned discrepancy, we conducted DFT calculations to elucidate the electronic structure of these catalysts. Figures 5a–c displays the band structures of the Ru element in Ru, NdRu₂, and Sc_1/8Nd_7/8Ru₂, in which the electron states near the Fermi level of both NdRu₂ and Sc_1/8Nd_7/8Ru₂ are significantly higher than that of pure Ru. To discern the enhanced N₂ activation over NdRu₂ and Sc_1/8Nd_7/8Ru₂, we calculated the projected density of states (pDOS), and the 3σ and 2π orbitals of free N₂ contribute to the energy range of −4 to 0 eV, aligning well with the energy level of Ru d orbitals (Fig. 5d), and their overlap is most prominent over Sc_1/8Nd_7/8Ru₂ (Fig. 5e), indicating a strong tendency for Sc_1/8Nd_7/8Ru₂ to interact with N₂. Additionally, the pDOS results of adsorbed N₂ reveal that the p-band of N splits into a bonding state and a small antibonding state (below and above the Fermi level, respectively), with the bonding orbitals of N₂ hybridizing with Ru d orbitals at around −8 to −5 eV (Supplementary Fig. 20). Notably, the electronic state at the bonding orbitals of N₂-adsorbed Sc_1/8Nd_7/8Ru₂ exhibits a more negatively shifted energy compared with NdRu₂ and Ru, indicating a more stabilized N₂. The significant N₂ adsorption results in the weakening of the N ≡ N triple bond and a more pronounced elongation of the bond length over Sc_1/8Nd_7/8Ru₂ (Fig. 5e), which is consistent with the prediction made by the SHAP value method in ML, suggesting that doping small-radius elements in Ru-IMCs would enhance N₂ adsorption.

The work function of Sc_1/8Nd_7/8Ru₂ on the crystal facet (001), the most stable surface, was calculated to be 2.95 eV (Supplementary Fig. 21), indicating its high electron donation ability. With the high electronic states near the Fermi level, Sc_1/8Nd_7/8Ru₂ is poised to facilitate electron injection into the adsorbed N₂. Hence, the 3D electron density isosurface map illustrates a pronounced electron transfer from Sc_1/8Nd_7/8Ru₂ to the adsorbed N₂ (Fig. 5d inset, Supplementary Figs. 22–24). The activated N₂ reaction pathway was further examined on the surfaces of Ru, NdRu₂, and Sc_1/8Nd_7/8Ru₂. NdRu₂ demonstrates stronger N₂ adsorption compared to Ru, with adsorption energies of −0.94 eV versus −0.62 eV (Fig. 5f). The addition of Sc enhances N₂ adsorption further, increasing it from −0.94 eV to −1.17 eV. This finding is in line with both the previously discussed reaction kinetics findings (Fig. 4d) and theoretical calculation results (Fig. 2d). The N₂ dissociation barrier over Ru, NdRu₂, and Sc_1/8Nd_7/8Ru₂ was determined to be 1.79 eV (for Ru flat site), 1.32 eV (for Ru B5 site see in Supplementary Fig. 25),1.15 eV, and 0.75 eV, respectively, which indicates that the adsorbed N₂ is more readily cleaved over Sc_1/8Nd_7/8Ru₂.

Experimental studies were conducted using the calculation data provided to assess the validity of the proposed N₂ activation mechanism. As depicted in Fig. 5g, N₂ desorption from Sc_1/8Nd_7/8Ru₂ occurs at 501 °C, a higher temperature range than that from NdRu₂ (379 °C), while there is almost no desorption from pure Ru, which suggests that the nitrogen interaction strength over Sc_1/8Nd_7/8Ru₂ is significantly stronger, consistent with the calculated adsorption energies shown in Fig. 5f. Quantification revealed 0.005 mmol of desorbed nitrogen, a value similar to the number of exposed Ru atoms (0.004 mmol) on the surface of Sc_1/8Nd_7/8Ru₂ catalyst (Supplementary Table 10). This, along with the XPS observation of a strong N 1s peak (~398.4 eV) on the used catalyst that almost disappeared after Ar plasma etching (Supplementary Fig. 26), supports the conclusion that the desorbed nitrogen originates primarily from surface species rather than from the bulk material of Sc_1/8Nd_7/8Ru₂. FT-IR results reveal broad peaks around 2300–2100 cm⁻¹, attributed to the stretching vibrational mode of chemisorbed N₂ on Ru sites (Fig. 5h). Due to N ≡ N stretching (ν_N2), the ¹⁴N₂-adsorbed NdRu₂ exhibits a broad band at 2222.6 cm⁻¹. While, ¹⁴N₂-adsorbed Sc_1/8Nd_7/8Ru₂ shows a broad band assignable to ν_N2 at around 2191.9 cm⁻¹ (Fig. 5h). The red-shift of the N₂ signal is clear evidence of a stronger electronic interaction between N₂ and Sc_1/8Nd_7/8Ru₂, in which N₂ accepts more electrons and is more dissociative than that on NdRu₂. For ¹⁵N₂ adsorbed on Sc_1/8Nd_7/8Ru₂, the adsorption band appears around 2118.4 cm⁻¹ (Fig. 5h), consistent with the value estimated based on the isotope effect (2191.9 cm⁻¹ × \(\sqrt{14/15}\) = 2117.6 cm⁻¹), indicating the chemisorbed N₂.

To determine the origin of the resistance to hydrogen poisoning over Sc_1/8Nd_7/8Ru₂, we conducted an H₂-TPD experiment to examine the thermal desorption of the used samples. As depicted in Fig. 5i, the mass signal of H₂ is scarcely detected on NdRu₂, whereas a prominent peak is observed at 171 °C on Sc_1/8Nd_7/8Ru₂, indicating that the ability to incorporate hydrogen has been achieved through the substitution of Sc. The mole ratio of the hydrogen can be estimated to be 0.52 through mass quantification, denoted as Sc_1/8Nd_7/8Ru₂H_0.52. DFT calculation results also reveal that the hydrogen-involving Sc_1/8Nd_7/8Ru₂ structure is much more energetically favored than that of NdRu₂, further confirming that H⁻ can be trapped in Sc_1/8Nd_7/8Ru₂, which accounts for the high resistance to hydrogen poisoning issue.