Ammonia holds a significant role as both a vital industrial chemical and a promising sustainable fuel source due to its impressive gravimetric hydrogen density, reaching 17.75% (in mass)^[1⇓-3]. This prevailing industrial approach for ammonia production is the Haber-Bosch method, employing an iron-based catalyst and operating at high temperature (400-600 ℃) and high pressure (150-300 atm). Regrettably, this process accounts for a substantial 1%-2% of global energy consumption and contributes to 1.44% of worldwide CO₂ emission^{[4⇓⇓⇓-8]}. The urgent need to reduce energy consumption and mitigate environmental impact necessitates the exploration of alternative ammonia production methods. In theory, this could be realized by consecutive protonation steps under an applied potential, bypassing the direct cleavage of robust N≡N bonds. However, this requires the development of highly efficient catalysts capable of enhancing the nitrogen reduction reaction, a pursuit that has sparked tremendous enthusiasm.

Prior researches have highlighted functional graphene as a promising candidate for enhancing ammonia production efficiency^[9⇓-11]. For instance, the introduction of molybdenum (Mo) as a dopant has been explored by Zhao et al.^[12], who reported a limiting potential of 0.49 V for Mo/N₃-embedded graphene. Similarly, Ma et al.^[13] found a thermodynamic barrier of 0.67 eV for Mo/N₄-embedded graphene along the distal reaction mechanism. While, Choi and colleagues^[14] noted limiting potential of 0.89 V for Mo/C₃ and 0.69 V for Mo/C₄-embedded graphene, respectively. In our previous work^[15], it was identified that Mo/C₆-embedded graphene presented a thermodynamic barrier of 0.37 eV. These results indicate that Mo, characterized by its partially-filled d band, strikes a delicate balance between reactant activation and potential intermediate inhibition. Furthermore, the variations in the reported limiting potentials hint at the influence of the metal center’s reactivity modulation by its coordination environment. This concept gains further support from a study conducted by Zhou et al.^[16], who introduced selenium (Se) heteroatom into Mo-doped graphene and observed a reduced limiting potential of 0.41 V for Mo/C₂Se-embedded graphene. Given the distinct electronic configurations of p-block elements, it is reasonable to anticipate that Mo’s reactivity in ammonia synthesis depends on the specific p-block elements that it coordinated with. Our previous work has demonstrated that graphene’s three-vacancy defects create sufficient space for the inclusion of p-block heteroatoms around the metal site^[15]. Therefore, a crucial question arises regarding the impact of p-block heteroatoms on the catalytic performance of N₂-to-NH₃ conversion.

To address this question, our study investigated the influence of secondary p-block dopants (X = B, P, or S) on the electrocatalytic conversion of N₂ to NH₃ in TM/C₆-embedded graphene. The transition metals (TM) included Fe, Nb, Mo, Ru, and W. Herein, the selection of B, P, or S stemed from the previous reports^[17⇓-19]. Wherein, B and P are identified as the active sites to boost nitrogen reduction reaction (NRR), meanwhile S is a coordination element found in azotobacter, and the selection of Fe, Nb, Mo, Ru, and W is due to their reported good nitrogen-fixation activity^{[12,15,20⇓⇓⇓⇓-25]}. Density functional theory (DFT) calculations were employed to evaluate the activity and selectivity by analyzing free energy profiles. We also assessed structural stability through binding energy and structural rigidity at room temperature. Furthermore, Mulliken charges and partial density of states were calculated to analyze the electronic structure and determine the extent to which the reactivity of the metal center is modulated by coordination with p-block heteroatoms.

We conducted all DFT calculations for the geometric optimizations and electronic properties using the DMol³ software package^[26-27]. The Perdew-Burke- Ernzerh (PBE) functional with generalized gradient approximation (GGA) was employed to describe the exchange-correlation effect^[28]. To account for relative effects, the DFT Semi-core Pseudopots (DSPP) core treat method was used, which replaces core electrons with a single effective potential^[29]. The double numerical atomic orbital augmented by a polarization function (DNP) was chosen as the basis set^[26]. To achieve accurate electronic convergence, a smearing of 0.005 Ha (1 Ha = 27.21 eV) was applied to the orbital occupation. Structure optimization was performed until convergence critera were met, with an energy, maximum force and maximum displacement of 1.0×10^-5 Ha, 0.02 Ha/nm and 0.0005 nm, respectively. For all calculations, the spin-unrestricted method was used. Additionally, a conductor-like screening model (COSMO) was employed to simulate a solvent environment resembling H₂O^[30]. COSMO is a continuum model where the solute molecule creates a cavity within the dielectric continuum. We extended the DMol³/COSMO method to accommodate periodic boundary conditions. The dielectric constant of H₂O was set as 78.54. The ab initio molecular dynamic (AIMD) simulations were conducted within the NVE ensemble to check the structural stability at room temperature.

(5×5) supercell was adopted to model the graphene monolayer consisting of 50 C atoms and the functional graphene was created via replacing three C atoms by one TM atom and one C atom by one X atom. To prevent artificial interactions between the catalyst and its images, a 1.5 nm thick vacuum was introduced. The adsorption energy, E_ads, of the NRR intermediate was calculated as follow:

Here, E_system, E_catalyst, and E_m represent the total energies of the adsorption system, the catalyst, and the adsorbates, respectively.

To access the thermodynamic stability of functional graphene, we calculated binding energy, E_b, using the following formula, widely used in the literature^[31]:

Where E_XG is the energy of graphene doped with non-metal atoms, and E_TM is the total energy of the transition metal atoms.

The Gibbs free energy change (ΔG) of the elementary steps was determined based on the CHE

model developed by Nørskov et al.^[32], where the chemical potential of the (H⁺ + e^-) pair in solution equals half that of gas-phase H₂^[32-33]. The corresponding ΔG was calculated as follow:

Here, ∆E represents the energy difference, ∆E_ZPE is the change in zero point energy, and ∆S is the entropy change. The temperature, T, was set as 298.15 K. E_ZPE and S were determined from vibrational frequencies using standard methods. To simplify calculations, the substrates were fully constrained, following the suggestion of Wilcox et al.^[34]. The term ∆G_U = - eU accounts for the free energy contribution caused due to variation in electrode potential U, while ΔG_pH is the pH correction of the free energy, equal to kTln10 × pH, with k representing the Boltzmann constant. In this case, the solution is acidic with a pH of zero. ΔG <0 indicates exothermic adsorption and vice versa. These methods have previously proven successfully in analyzing the N₂ fixation process, demonstrating their robustness and feasibility^[14,35].

The limiting potential U_L is given by:

Where ΔG_PDS is the free energy change of the elementary step with the maximum endothermic character during the N₂-to-NH₃ conversion.

Fig. 1(a) illustrates the structural model of TMX doped graphene. Fig. 1(b) shows the bond lengths, l_TMX, of TM-X bonds after geometrical optimization. Wherein, l_TMB is the shortest while l_TMP and l_TMS are similar in magnitude. For example, the l_TMX values are 2.108, 2.359 and 2.331 Å (1 Å = 10^-10 m) for Mo-B, Mo-P and Mo-S bonds, respectively. To further characterize the TM-X interaction, Fig. 1(c) shows the binding energies, E_b(TMX), between TM and its coordination atom. It reveals that E_b(TMS) is the weakest while E_b(TMB) is similar to E_b(TMP). For example, E_b(TMX) are -6.71, -6.95 and -5.80 eV for the MoB, MoP and MoS doped graphene, respectively. The disparity in E_b(TMX) stems from valence electron configurations (2s²2p¹ for B, 3s²3p³ for P and 3s²3p⁴ for S). Due to lone pair, the S prefers to form two bonds with its coordination environment, for example S at the edge of graphene; however, it is different for B and P, which prefer to form three bonds with its coordination environment^[36]. The presence of lone pair in the p band of S leads to the weakened E_b(TMS), regardless of TM selection. As known, a high E_b represents strong chemical interaction, indicating that the metallic atom is firmly embedded into the graphene^[20]. To further investigate the origin of the E_b differences, the partial density of states (PDOS) between the d band of Mo and the p band of X is shown in Fig. 1(d). The existence of pd coupling below the Fermi energy level confirms the formation of the Mo-X covalent bond. Furthermore, the reduced pd hybridization of the Mo-S bond is consistent with the weakened E_b(MoS) value. Furthermore, the introduction of X into graphene also changes the d band distribution of Mo. The d band center, ɛ_d, is used to quantify this variation. Herein, the corresponding values are -1.54, -1.93 and -1.82 eV for MoB, MoP and MoS, respectively. A similar trend is observed for the rest TMX systems. The ɛ_d data are listed in Table S1 in Supporting Materials. With exception of FeX system, it is general tendency that ɛ_d of TMS is in the medium. However, E_b(TMS) is always the weakest. Clearly, ɛ_d and E_b do not have the same trend. To shed the lights for the mentioned observations, we further analyzed the charge distribution and the Mulliken charge Q is also tabulated in Table S1 in Supporting Materials. As we known, one assumption of d band theory is that there is no electron change or electron filling of one element is fixed^[37]. However, it is different from our case that Mulliken charge of TM is varied. Therefore, it is one reason that d band center ɛ_d could not describe the tendency of E_b(TMX). Another possible reason is from the structural flexibility. The binding energy E_b is a sum of interaction energy E_int and deformation energy E_def. As shown in Table S2 in Supporting Materials, the E_def changes from 1.23 eV to 3.37 eV, indicating the structural deformation, which may raise the different trend. Herein, it is clear that the secondary dopant X is expected to change the TM reactivity.

Next, the adsorption behavior of the nitrogen molecule on graphene is considered. Fig. 2(a) illustrates the different adsorption structures, with aN₂ and bN₂ referring to the end-on and side-on N₂ adsorption, respectively, on the side of graphene with TM protrusion, while cN₂ and dN₂ stand for the corresponding adsorptions on the side of graphene with X protrusion. Fig. 2(b) shows the difference in adsorption free energy between N₂ molecules and H atoms, i.e. ΔΔG_ads = G_ads(*N₂) - G_ads(*H), as a function of the free energy of nitrogen adsorption, G_ads(*N₂), for the different structures. Possessing a good affinity toward nitrogen and a relatively weak hydrogen adsorption are prerequisite for improving the selectivity of NRR vs hydrogen evolution reaction (HER). Therefore, a good catalyst should possess the properties of G_ads(*N₂) < 0 eV and ΔΔG_ads < 0 eV^[38]. Based on the G_ads(*N₂), aN₂ is identified as the most stable configuration, and hence is the only initial state considered in the subsequent activity investigation. Then, the H adsorption on top site of TM atom was considered, which is the same as aN₂ adsorption, in order to reveal the competition of reactant adsorption. After optimization, H is still located at the top site of TM atom with only exception of FeB wherein H adsorption is on bridge site between Fe atom and B atom. The corresponding adsorption structures are shown in Fig. S1 in Supporting Materials. From ΔΔG_ads values, H adsorption is energetically favorable on WP and WS while N₂ adsorption is favored on the rest systems. To further investigate the N₂ adsorption, Fig. 2(c) shows the PDOS between the p band of adsorbed N₂ and the d band of Mo as a representative system. It is evident that there is orbital hybridization around the Fermi energy level, which is consistent with the negative G_ads(*N₂). Additionally, electron occupation in the Fermi energy level is beneficial for electron communication upon adsorption. To quantitatively describe the activation of the N₂ adsorbate, its Mulliken charge, Q_N2, and bond length, l_N2, are shown in Fig. 2(d). The results show that the N₂ adsorbate is negatively charged on the MoX, NbX, and WX systems. Due to electron accumulation, the value of l_N2 increases from its original (unadsorbed) value of 1.108 Å to above 1.120 Å, indicating weakened N≡N bonds. Hence, this data validates the activation of N₂ on these systems. Furthermore, electron accumulation on the N₂ adsorbate is crucial for triggering the subsequent protonation progress due to electrostatic attraction between the negatively charged N₂ adsorbate and positively charged H⁺. However, the N₂ adsorbate is positively charged on FeX and RuX systems, implying insufficient N₂ activation.

After establishing the catalyst structure and adsorption properties, a complete study of the NRR was conducted. Due to the energetic preference of N₂ adsorption in an end-on manner, the feasible reactions are via the distal and alternative pathways, not the enzymatic pathway^[17,39-40]. This means that the proton-electron pairs will either continuously attack one N atom until NH₃ is formed or alternately attack two N atoms and thereby produce two NH₃ molecules simultaneously. What is more, we did not consider the direct splitting pathway, which is generally infeasible for single-active site. Taken MoX as an example, optimization calculations on the adsorption structures consisting of two N atoms co-adsorption on Mo top site were performed. After optimization, two N atoms adsorbed at MoB and MoS systems would spontaneously combine together to form N₂, and the co-adsorption of two N atoms is only possible for MoP system. Then the calculation of TS search was conducted to identify the kinetic barrier of N₂ splitting catalyzed by MoP moiety. As shown in Fig. S2 in Supporting Materials, the reaction barrier is about 3.91 eV. The extra-high barrier in combination with the endothermic reaction energy means the impossibility for direct N₂ splitting. Therefore, we considered the distal and alternative pathways in the following discussions.

Taking MoS as an example, Fig. 3(a, b) describe the free energy profiles under the distal and alternative pathways, respectively. For the former, ΔG of the elementary (consecutive) protonation steps are 0.47, -0.22, -0.54, -0.42, -0.66, -0.33, and 1.49 eV, respectively. Therefore, the potential-determining step (PDS) along the distal pathway is the first protonation step, i.e. N₂ + H⁺ + e^- → *NNH, and its corresponding thermodynamic barrier, ΔG_PDS, is 0.47 eV. For the alternating pathway, the corresponding ΔG are 0.47, 0.51, -0.91, 0.15, -1.59, -0.33, and 1.49 eV, respectively. In this case, the PDS is the second protonation step, i.e. *NNH + H⁺ + e^-→ *NHNH, and ΔG_PDS is 0.51 eV. Hence, the distal pathway is slightly preferred under the limiting potential, U_L(NRR), of 0.47 V. Following a similar analysis, ΔG and U_L(NRR) for all TMX systems are summarized in Table S3 in Supporting Materials. The lower U_L(NRR), the less energy is required to catalyze the reaction, indicating higher activity. Generally, a U_L(NRR) less than 0.55 V indicates good reactivity toward nitrogen fixation. Using this criterion, WS, WP, MoS, MoP, and NbB are expected to be active, with U_L(NRR) of 0.22, 0.32, 0.47, 0.49, and 0.54 V, respectively, and their activities are in the order of WS > WP > MoS ≈ MoP ≈ NbB. The rest systems considered are not very reactive for ammonia synthesis due to their relatively high U_L(NRR), which is consistent with the Q_N2 analysis.

As HER is the main side reaction, which directly consumes protons and electrons meant for NRR^[14,41-42], the selectivity of NRR vs HER was studied further via a comparison of limiting potentials. Fig. 3(c) shows the HER free energy profiles of WS, WP, MoS, MoP, and NbB systems. For WS, WP, MoS, and MoP, the first protonation step is exothermic while the second protonation step is endothermic, indicating that hydrogen evolution suffers from too strong H adsorption. As shown in Fig. 3(c), the decoration of X heteroatoms influences the HER performance, resulting in U_L(HER) values of 0.95, 0.84, 0.42 and 0.51 V for WS, WP, MoP and MoS, respectively. By contrast, the proton capture step on NbB is endothermic, and the HER occurs under a U_L(HER) of 0.15 V. Plotting U_L (NRR) vs U_L(HER) in Fig. 3(d) shows that WS, WP, and MoS are active for ammonia synthesis; however, the efficiencies of N₂-to-NH₃ conversion on MoP and NbB are limited by the undesirable HER. Furthermore, it should be noted that, as mentioned earlier, the positive ΔΔG_ads for WS and WP demonstrate that the H adsorption on these catalysts is stronger than N₂ adsorption, which implies that the active sites are inaccessible for N₂ molecules due to being covered by protons. Therefore, MoS appears to be the only promising candidate for nitrogen fixation electrocatalysis.

From the mentioned results, it is clear that the introduction of p-block heteroatoms influences the reactivity of the metal center in the application of N₂-to-NH₃ conversion. In order to further reveal the physical origin of this effect, Mulliken charge population analysis of MoX systems during the consecutive protonation progress was performed. According to previous report^[17], each intermediate is divided into three moieties: N_xH_y adsorbates (moiety 1), MoX (moiety 2) and the graphene support (moiety 3). The results presented in Fig. 4(a-c) show electron communication between the graphene support and the intermediates, indicating effective electron transfer during protonation steps, which is consistent with previous studies^[43-44]. Furthermore, the MoX moiety donates electrons to N_xH_y species, with MoS less positively charged compared with MoB and MoP. It demonstrates that the introduction of p-block heteroatoms influences the electrocatalysis of nitrogen reduction by affecting the electron transfer. To support the idea, we further plotted the fitting between the limiting potential U_L(NRR) and the Mulliken charge Q(N₂*) of the adsorbed N₂, as shown in Fig. S3 in Supporting Materials. Consistently, more charge accumulated in N₂ reactant, smaller U_L(NRR) required for nitrogen reduction. Our results are in line with the reported electron acceptance-back donation mechanism^[21,45]. Last but not least, AIMD simulation is preformed to check the structural robustness of MoS doped graphene. As shown in Fig. 4(d), its atomic structure is well-maintained at room temperature. Therefore, Mo/S co-doped graphene endowed with low limiting potential and good structural stability, holds feasibility in the application of ammonia synthesis electrocatalysis.

This study explores the potential of transition metal/non-metal co-doped graphene as electrocatalysts for ammonia synthesis through DFT calculations. Our findings indicate that W/S, W/P, Mo/S, Mo/P, and Nb/B exhibit promising activity with relatively low limiting potentials of 0.22, 0.32, 0.47, 0.49, and 0.54 V, respectively. Nevertheless, W/S, W/P, Mo/P, and Nb/B also enhance the side hydrogen evolution reaction, compromising ammonia synthesis selectivity. Consequently, Mo/S emerges as the most compelling candidate for experimental verification due to its remarkable efficiency and stability. Furthermore, our DFT calculations reveal that the metal center’s reactivity can be finely tuned by coordinating with p-block heteroatoms, presenting novel prospects for ammonia synthesis.

The supporting materials related to this article can be found at https://doi.org/10.15541/jim20230433.

Mo/S Co-doped Graphene for Ammonia Synthesis: a Density Functional Theory Study

LI Honglan¹, ZHANG Junmiao¹, SONG Erhong², YANG Xinglin¹

(1. School of Energy and Power Engineering, Jiangsu University of Science and Technology, Zhenjiang 212003, China; 2. State Key Laboratory of High Performance Ceramics and Superfine Microstructure, Shanghai Institute of Ceramics, Chinese Academy of Sciences, Shanghai 200050, China)