DFT investigation of the oxygen reduction reaction over nitrogen (N) doped graphdiyne as an electrocatalyst: the importance of pre-adsorbed OH* and the solvation effect

Abstract

Searching for novel electrocatalysts that can replace precious platinum for oxygen reduction reaction (ORR) is important for developing fuel cells. Recently, nitrogen (N) doped graphdiyne (GDY) has been synthesized and proved that the ORR electrocatalytic activity catalyzed by N-doped GDY is significantly improved, however, the roles of sp-N (including sp-N1 and sp-N2) and pyridinic (Pyri)-N dopants in mediating the ORR are still unclear. To clarify which sp-N or Pyri-N creates the active site for ORR, we systematically studied the ORR mechanism on sp-N1GDY and pyri-NGDY supported on graphene (G) with the solvation effect, which was performed using density functional theory (DFT) and ab initio molecular dynamics (AIMD) simulations. Firstly, we found that the dissociative mechanism is preferred on sp-N1GDY/G and the surface is easily terminated by the OH* intermediates, while the OH* pre-adsorbed surface (sp-N1GDY(OH)/G) prefers the associative mechanism. Pyri-NGDY/G also prefers the associative mechanism without any termination. Then, the solvation effect stabilizes all ORR intermediates in both cases. From the calculated free energy diagram, a model with water solvent gives a more appropriate estimation of the overpotential than the one without the water solvent, and sp-N1GDY/G with OH* pre-adsorbed has a lower overpotential (0.46 V), which is close to the experimental value (0.36 V), compared with Pyri-NGDY/G (0.75 V). Our study provides useful information for understanding the reaction mechanisms of ORR at the solid/liquid interface on the N-doped GDY surface.

---

1. Introduction

Fuel cells provide an innovative and efficient energy conversion technology and operate via electrochemical reactions. The use of clean and renewable fuels, such as hydrogen, methanol, and biomass, reduces dependence on fossil fuels and their associated environmental impacts.^1–3 The oxygen reduction reaction (ORR) is a key electrochemical reaction that takes place at the cathode of the fuel cells.⁴ Traditionally, platinum (Pt) and Pt-based materials have been widely used as electrocatalysts for ORR due to their excellent electrocatalytic activity, durability, and stability;⁵ however, the high cost and scarcity of Pt have led to significant research efforts to find alternative electrocatalysts that can perform as effectively as Pt-based materials while being more cost-effective and environmentally sustainable.^4,6 In 2009, Dai et al. reported that nitrogen-doped carbon nanotube arrays possessed high electrocatalytic activity for the ORR in alkaline media, even superior to Pt.⁷ This finding led to the development of an interesting field of carbon-based metal-free ORR electrocatalysts because of their advantages, including relatively abundant raw materials, economic feasibility, adjustable surface chemistry, easy processing, large specific surface area, high chemical stability, fast transfer kinetics, and a wide operating temperature range.⁶ Based on many findings,^8–11 N-doping is a particularly advantageous method to modify carbon materials for energy conversion, such as hydrogen evolution reaction (HER), oxygen evolution reaction (OER), ORR, nitrogen reduction reaction (N₂RR), and carbon dioxide reduction reaction (CO₂RR).

Recently, a new type of carbon material named graphdiyne (GDY), a rising star of 2D carbon allotropes with one-atom-thick planar layers, has achieved the co-existence of sp and sp² hybridized carbon atoms in a 2D planar structure.^12–14 GDYs with unique chemical and physical properties have attracted more attention in energy conversion.^15–17 However, pure GDY has low ORR electrocatalytic activity, and structural modification, including doping with metal or nonmetal atoms is a good way to improve the activity. Gao et al.¹⁸ reported that Fe-doped GDY (Fe@GDY) showed high catalytic activity towards ORR at even higher activity than the benchmark commercial Pt/C (20 wt%). Dai et al.¹⁹ predicted that Ni@GDY and Pt@GDY catalysts possessed comparable electrocatalytic activity for ORR and OER in alkaline media based on DFT studies. Li et al.²⁰ reported that N-doped GDY, especially sp-N doping showed much better ORR electrocatalytic performance than the commercial Pt-based catalyst in alkaline media and comparable activity in acidic media. On the other hand, Huang et al.²¹ reported that pyridinic (Pyri)-NGDY is mainly produced by N doping and Pyri-NDGY exhibits excellent ORR and OER catalytic activity. Despite experimental evidence demonstrating the excellent ORR electrocatalytic activity of N-doped GDY, the nature of the active sites and the dopants involved are still controversial. Recently, a few theoretical ORR studies on N-doped GDY with sp-N configurations have been reported.^22–24 Li et al.²² found that ORR proceeds on sp-N(II)GDY with the associative mechanism. The highest ORR electrocatalytic activity has a theoretical onset potential of 0.76 V with metastable adsorption of all the ORR intermediates. Lee et al.²³ found that double N-doped GDY has better ORR activity. In the above studies, it turns out that the metastable adsorption sites of ORR intermediates govern the ORR performance. However, the role of the most stable adsorption sites in the ORR mechanism on N-doped GDY remains unclear.

Density functional theory (DFT) calculation is a powerful tool to explore the ORR mechanism using the computational hydrogen electrode (CHE) model,²⁵ however, to accurately elucidate the ORR activity, it is important to incorporate the solvation effect. There are two approaches: (1) the explicit model, and (2) the implicit model. The former approach involves adding water molecules around the reactant or above the surface to represent the solvation effect,^26–29 and the latter approach replaces the liquid environment with a polarizable continuous medium (PCM).^30,31 Previous theoretical studies on metal-free catalysts mostly used the implicit model, which tends to underestimate the solvation effect during the ORR process. Boresch et al.³² reported that PCM cannot describe any direct solvent effects and might be inaccurate for specific hydrogen bonds. Wang et al.³³ reported that the explicit model can give a correct estimation of the ORR overpotential on MnN₄–graphene catalyst than the gas-phase model. Hansen et al.³⁴ investigated the solvation energy of ORR intermediates on N-doped graphene using the explicit and implicit models, and found that including explicit H₂O molecules is essential for a correct description of the adsorption energy of ORR intermediates on carbon materials, while continuum solvation models are unable to accurately describe the solvation energy.

In this study, we systemically investigated the ORR activity on sp-N1 and Pyri-N-doped GDY with graphene (G) support. DFT calculations and ab initio molecular dynamics (AIMD) simulations were carried out to investigate the ORR mechanism as well as the solvation effect for ORR in both cases. Free energy diagrams showed that the active sites of sp-N1-doped GDY surfaces are easily terminated by OH*, and the neighboring C sites of the –COH–N moiety become real active sites for ORR to proceed with low overpotentials compared with those in Pyri-N. Moreover, the explicit solvation model can give a more consistent overpotential (0.46 V) with the experimental value (0.36 V) compared with that under vacuum conditions (0.72 V).

2. Computational methods

Most of our DFT calculations were performed using the Simulation Tool from the Atom Technology (STATE) program package.^35–38 Ionic cores were described using the ultrasoft pseudopotentials,³⁹ and valence electron states were expanded using a plane wave basis set with the kinetic energy cutoffs of 36 and 400 Ry for wave functions and augmented charge density, respectively. A 2 × 4 × 1 Monkhorst–Pack k-point sampling was used for the N-doped GDY 2 × 1 × 1 supercell. The revised PBE (RPBE) functional⁴⁰ with the Grimme's van der Waals (vdW) correction (D2)⁴¹ was employed to describe the weak dispersion interactions between various intermediates and surfaces as well as to describe water–water interactions more accurately.²⁹

Spin-polarized calculations were performed whenever triplet O₂ was simulated. We considered the intersystem crossing (ISC) upon O₂ adsorption on the substrate. The isolated O₂ in the triplet state (³Σ_g⁻) as its ground electronic state is relatively unreactive due to its low energy. However, when O₂ approaches surfaces such as carbon nanotubes and diamond (100),^42,43 it can undergo the intersystem crossing (ISC) to a more reactive singlet state (¹Δ_g). In our calculations, the adsorbed O₂ in ³Σ_g⁻ was modeled by fixing the difference between the number of spin-up and spin-down electrons to two, while the adsorbed O₂ in ¹Δ_g was modeled by a spin-restricted calculation to obtain a closed-shell configuration, which was necessary to prevent spin contamination. Activation reaction barriers (E_a) were calculated using the climbing image nudge elastic band method (CI-NEB).^44,45

According to the experimental conditions,²⁰ the overall reaction of O₂ reduction to H₂O on N-doped GDY under alkaline conditions is O₂(g) + 2H₂O(l) + 4e⁻ → 4OH⁻. The adsorption Gibbs free energy changes (ΔG_ads) of all the elementary steps in ORR on this N-doped GDY/G using the CHE method²⁴ were computed as:

ΔG_ads(X*) = ΔE_X* + (ΔH_X* − TΔS_X*)−neU_RHE + ΔSol(X*)

where ΔE_X*, ΔH_X* and ΔS_X* are the changes in the adsorption energies, enthalpies, and entropies of each intermediate on the surface under vacuum conditions, respectively. T is the temperature (298.15 K). ne is the number of electrons transferred, and the number of electrons transferred in 2O*, O* + OH*, OOH*, O*, and OH* are 4e⁻, 3e⁻, 3e⁻, 2e⁻ and 1e⁻, respectively. U_RHE is the potential of the electrode relative to the reversible hydrogen electrode (RHE).⁴⁶ We describe this equilibrium using the RHE as the reference electrode, which equals the chemical potential of H⁺ + e⁻ to the chemical potential of 1/2H₂ at an arbitrary pH (pH₂ = 1 bar and T = 298.15 K).⁴⁶ ΔSol(X*) is the solvation energy of each intermediate. To calculate the ΔSol of ORR intermediates, we used the AIMD⁴⁷ calculations to find hydrogen bond (H bond) networks at the interface. Further details on the computations on ΔG_ads, reaction Gibbs free energy, and ΔSol are summarized in the ESI.†

3. Results and discussions

3.1 Atomic and electronic structures of the N-doped GDY catalysts

Based on the experimental results,^20,21 we constructed two types of N-doped GDY (NGDY), namely sp-NGDY and Pyri-NGDY. As shown in Fig. 1, the carbon atom in the diacetylene linkage (sp hybridized C) is replaced by N to generate two types of sp-N (Fig. 1(b) and (c)), i.e., sp-N1GDY and sp-N2GDY. Pyri-NGDY is the sp² hybridized N atom bonded with two sp² hybridized carbon (C) neighbors in the carbon ring and with H-terminated C atoms (Fig. 1(d)). In the experiment,²⁰ N-doped GDY was supported onto glassy carbon, herein, to simplify the model, graphene (G) was applied as a support underneath the N-doped GDY. The obtained optimized lattice constants of the N-doped GDY (9.46 Å), and graphene (2.46 Å) unit cells were consistent with those from previous results.^48,49 Considering the lattice matching between the two components, the p(8 × 4) supercell of G and the p(2 × 1) supercell of N-doped GDY were employed to create the interface model (Fig. 1(e)–(g)) by applying a mixed tensile/compressive strain (4%) to the graphene phase, resulting in corrugated graphene. The calculated minimum distances between sp-NGDY or Pyri-NGDY, and the top of graphene were 2.27 and 2.17 Å, respectively, which are close to those of other heterostructures.⁵⁰ The band gap (E_g) of the isolated G, NGDY in three N configurations (sp-N1, sp-N2, and Pyri-N) and their interface models (sp-N1GDY/G, sp-N2GDY/G, and pyri-NGDY/G) were calculated using the PBE⁵¹ and HSE06 functionals^52,53 using the VASP code,^54,55 as shown in Table S1 (ESI†). We find that Pyri-NGDY is a semiconductor with a band gap of 0.60 eV by PBE and 0.97 eV by HSE06, while sp-N1GDY and sp-N2GDY are metallic, as predicted by both functionals. After introducing the G substrate as the support, all sp-N1GDY, sp-N2GDY, and Pyri-NGDY showed metallicity.

Fig. 1 (a) Schematic representation of different types of N doping configurations (sp-N1, sp-N2, and Pyri-N) in GDY. Optimized structure of the sp-N1GDY (b), sp-N2GDY(c), Pyri-NGDY (d), sp-N1GDY/G (e), sp-N2GDY/G, (f) and pyri-NGDY/G (g). Red, white, gray, pink, and blue balls are O atoms, H atoms, C atoms in N-doped GDY, C atoms in G, and N atoms, respectively.

To explore the charge transfer of the three N-doped GDY/G interfaces, we calculated the projected density of states (PDOS) (Fig. 2), charge density difference (CDD), and Bader charge of the three cases (Fig. S1, ESI†). The PDOS results showed that after the introduction of G, the Dirac point of G in sp-N1GDY/G, sp-N2GDY/G, and Pyri-NGDY/G is upshifted above the Fermi level, on the other hand, the C(p)s of sp-N1GDY, sp-N2GDY, and Pyri-NGDY in sp-N1GDY/G, sp-N2GDY/G, and Pyri-NGDY/G are downshifted below the Fermi level. The C (p) band centers in the three N-doped GDYs are shifted to a lower energy region upon the introduction of G (Fig. 2). The above results indicate that the use of G can increase the conductivity of N-doped GDY and induce the charge transfer from G to N-doped GDY. The CDD (Fig. S1, ESI†) results also revealed that the charge densities are redistributed by forming electron- and hole-rich regions within the three N-doped GDY/G interfaces. Charge depletion occurs on the G surfaces, resulting in the hole-rich sites. However, strong charge accumulation occurred on the three N-doped GDY surfaces, forming electron-rich sites. Therefore, the electrons are mainly transferred from G to the three N-doped GDY surfaces. To confirm this, the Bader charge analysis (Fig. S1, ESI†) and work function change (Table S2, ESI†) were also conducted. From the Bader charge analysis, 0.11 e⁻ are transferred from the G substrate to sp-N1GDY or sp-N2GDY, and 0.13 e⁻ are transferred from the G substrate to Pyri-NGDY. As shown in Table S2 (ESI†), the work functions of planar G (4.20 eV) and corrugated G (4.10 eV) are smaller than those of pure GDY (5.10 eV), sp-NGDY (4.53 eV), and Pyri-NGDY (4.93 eV), which is the fundamental cause for charge transfer from G to GDY and NGDY.

Fig. 2 (a) The projected density of states (PDOS) on atomic orbitals of sp-N1GDY, G, and sp-N1GDY/G; (b) PDOS of sp-N2GDY, G and sp-N2GDY/G; (c) Pyri-NGDY, G and Pyri-NGDY/G. The C(p) energy window of the integration is [−10.00, 0.00].

3.2 ORR mechanism and free energy analysis of sp-N1GDY/G, sp-N2GDY/G, and Pyri-NGDY/G under vacuum conditions

ORR mechanisms on clean surfaces. O₂ adsorption and dissociation are the two important steps that govern the ORR pathway. Thus, we first evaluated O₂ adsorption and dissociation on N-doped GDY models. We investigated all possible adsorption sites (C1–C5) of single O₂ on sp-N1GDY/G, sp-N2GDY/G, and Pyri-NGDY/G in two different geometries, namely end-on and side-on (Fig. S2, ESI†). We found that O₂ preferably adsorbed on the top sites of the C atoms near the N dopants with different strengths between sp-N and Pyri-N configurations. O₂ was chemisorbed at the C3C4 site (Table 1 and Fig. S2(b), ESI†) and the C2C4 site (Table S3 and Fig. S2(c), ESI†) on sp-N1GDY/G and sp-N2GDY/G, respectively. In contrast, O₂ is physisorbed on Pyri-NGDY/G (Table 1 and Fig. S2(e), ESI†).

Table 1 The adsorption energy of O₂ and OOH, and the activation energy (E_a) of O₂ dissociation and O₂ protonation to OOH on sp-N1GDY/G, sp-N1GDY(OH)/G, and Pyri-NGDY/G

Structure	ΔEads(O2)/eV	ΔEads(OOH)/eV	E a/eV
			O2 dissociation			O2 protonation
			O2 → O2* → 2O*	O2 → 2O*	O2 → 2O* with one H2O	O2 + H2O → OOH* + OH*
sp-N1GDY/G	−0.6	—	0.25	0.22	0.09	0.88
sp-N1GDY(OH)/G	−0.15	−1.23	—	1.86	1.64	1.40
Pyri-NGDY/G	−0.10	−0.93	—	1.52	1.35	1.17

---

Next, we examined the reaction barriers for O₂ dissociation to 2O* and O₂ protonation to OOH* using the CI-NEB method with the consideration of the ISC. We considered two probable pathways for O₂ dissociation cases, namely (1) gas-phase O₂ adsorption on the surface followed by O₂* dissociation to 2O* (O₂ → O₂* → 2O*); (2) direct dissociative adsorption of gas-phase O₂ to 2O* (O₂ → 2O*). The triplet and singlet potential energy surfaces (PES) are shown in Fig. S3 (ESI†). For all reaction processes, O₂ starts in a triplet state, reaches the ISC state, and changes to a singlet state. The activation reaction energies (E_a) of O₂ dissociation and protonation on the three N-doped GDY/G are shown in Table 1 and Fig. S3 (ESI†). We found that on two sp-NGDY/G systems (sp-N1GDY/G and sp-N2GDY/G), the direct O₂ → 2O* path is kinetically the most preferable with relatively low E_a (∼0.22 eV on sp-N1GDY/G, and ∼0.08 eV on sp-N2GDY/G). Moreover, water induces O₂ dissociation by lowering the activation barrier by ∼0.2 eV on all substrates due to the H bond effect between water and O₂.^56,57 In contrast, O₂ protonation to OOH* exhibits higher activation energy (∼0.88 eV on sp-N1GDY/G and ∼1.41 eV on sp-N2GDY/G) compared with that of O₂ dissociation. Consequently, we assume that the ORR mechanism on sp-N-doped catalysts will follow the dissociative mechanism (O₂ → 2O* → O* + OH* → O* → OH* → OH⁻). On the Pyri-NGDY/G surface, the O–O bond of O₂ does not activate upon the adsorption due to its physisorption state.

From Table 1, protonation to OOH* has a lower activation energy than its dissociation with one water (1.17 eV vs. 1.35 eV), however, this protonation barrier (1.17 eV) is still rather high; thus protonation process is kinetically unfavorable due to physisorption O₂.

It is noted that even though O₂ is weakly adsorbed and the activation barrier for the direct protonation of the adsorbed O₂ is high, ORR can still proceed with high activities on carbon materials via a different reaction pathway, namely, a process of long-range electron transfer (ET) to nonadsorbed O₂ in the outer Helmholtz plane (ET-OHP).^58,59 Choi et al.⁵⁹ recently revealed that the first proton-coupled electron transfer (PCET) step (O₂(aq) + (H⁺ + e⁻) → ˙OOH (aq)) can occur via ET-OHP, where O₂ directly forms ˙OOH (aq). Then, ˙OOH (aq) subsequently adsorbs on the catalytic site as OOH* (i.e., ET-OHP mechanism). This explains the problem of not finding suitable sites for O₂ binding on N-doped graphene and shows that the O₂ chemisorption is not essential for the occurrence of ORR on carbon-based catalysts. We found that O₂ could not chemisorb and that the O₂ protonation barrier (1.17 eV) was high on Pyri-NGDY/G. However, OOH* can chemisorb on the surface with ΔE_ads of −0.93 eV (Table 1). Therefore, we assume that the ORR on Pyri-NGDY/G follows ET-OHP associative mechanism (O₂ → OOH* → O* → OH* → OH⁻).

Then, the reaction and adsorption Gibbs free energies of each ORR intermediate on three N-doped models, i.e., sp-N1GDY/G, sp-N2GDY/G, and Pyri-NGDY/G, were calculated and are shown in Fig. 3, Fig. S4, and Tables S4–S6 (ESI†), respectively, for all the considered systems. We find that ORR proceeds with a rather high overpotential (η) on the three substrates. For Pyri-NGDY/G (Fig. 3e), η is 1.01 V, and the potential-determining step (PDS) is the O₂ → OOH* (ΔG₆). High η of ORR on sp-N1GDY/G (Fig. 3(a)), and sp-N2GDY/G (Fig. S4(a), ESI†) arises from the strong interactions of ORR intermediates with the substrate. The PDS for ORR on the N-doped catalyst in sp-N configurations is found to be the last step OH* → OH⁻ (ΔG₅). The third ORR step (O* + OH* → O*) is exergonic and becomes endergonic upon applying a potential of 1.23 V. The last steps (OH* → OH⁻) are endergonic reactions even at U_RHE = 0 V (0.13 eV on sp-N1GDY/G and 0.38 eV on sp-N2GDY/G), showing that OH* has a strong binding energy with the surface and is difficult to further hydrogenate. As a result, the first ORR could not be completed, and some ORR intermediates, such as O* or OH*, remain adsorbed on sp-NGDY/G.

Fig. 3 (a) Free energy diagram and structures of each ORR intermediate on sp-N1GDY/G. (b) Free energy of each ORR intermediate versus electrode potential (vs. RHE) on sp-N1GDY/G. (c) Free energy diagram and structures of each ORR on sp-N1GDY(OH)/G. (d) Free energy of each ORR intermediate versus electrode potential (vs. RHE) on sp-N1GDY(OH)/G. (e) Free energy diagram and structures of each ORR intermediate on Pyri-NGDY/G. (f) Free energy of each ORR intermediate versus electrode potential (vs. RHE) on Pyri-NGDY/G. Purple circles in (b), (d), and (f) represent the lowest lines crossed at different potentials. Pristine in (b), (d), and (f) represents sp-N1GDY/G, sp-N1GDY(OH), and Pyri-NGDY/G, respectively. Red, white, gray, pink, and blue balls are O atoms, H atoms, C atoms in N-doped GDY, C atoms in G, and N atoms, respectively. U_RHE is the potential of the electrode relative to the RHE.

To confirm this, the adsorption Gibbs free energy of each ORR intermediate as a function of the electrode potential (U_RHE) was evaluated. As shown in Fig. 3(b), sp-N1GDY/G was terminated by OH* at 0 V < U_RHE < 0.86 V, by O* at 0.86 V < U_RHE < 0.88 V, by O* + OH* at 0.88 V < U_RHE < 1.07 V, and by 2O* at 1.07 V < U_RHE < 1.23 V. Similarly, for sp-N2GDY/G (Fig. S4(b), ESI†), sp-N2GDY/G was terminated by OH* at 0 V< U_RHE < 0.49 V, by O* at 0.49 V < U_RHE < 0.73 V, and by 2O*at 0.73 V < U_RHE < 1.23 V. In contrast, on the Pyri-NGDY/G surface, the most stable structure is the Pyri-NGDY/G pristine surface under 0 V < U_RHE <1.11 V (Fig. 3(f)). It should be noted that in this work, we only considered single ORR intermediates as a function of potential, the coverage⁶⁰ of ORR intermediates may or may not affect the results, this subject will be investigated in the next project.

Consequently, upon applying a limiting potential of 0.22 V, the complete ORR process occurring on Pyri-NGDY/G, and ORR species will never be terminated on the surface. In contrast, OH* will be terminated on sp-N1GDY/G and sp-N2GDY/G surfaces without any applied potential (sp-N1GDY(OH)/G and sp-N2GDY/G hereafter). Therefore, we explored the ORR mechanism of sp-N1GDY(OH)/GDY and sp-N2GDY(OH)/GDY in the next subsection.

ORR mechanisms on OH pre-adsorbed surfaces. We find that the strength of O₂ adsorption is decreased upon adsorption on sp-N1GDY(OH)/GDY and sp-N2GDY(OH)/GDY compared to that on clean surfaces. As shown in Table 1 and Table S3 (ESI†), the vdW attractions dominate the interaction between O₂ and sp-N1GDY(OH)/GDY, resulting in ΔE_ads of O₂ and distance between O₂ and sp-N1GDY(OH)/GDY as −0.15 eV and 3.0 Å, respectively. While, on sp-N2GDY(OH)/GDY, O₂ is weakly chemisorbed with ΔE_ads of −0.58 eV and an adsorption distance of 1.426 Å. Then, the E_a of the O₂ dissociation and protonation on sp-N1GDY(OH)/GDY and sp-N2GDY(OH)/GDY calculated using the CI-NEB method is shown in Table 1 and Fig. S5 (ESI†). We obtained a high activation energy for O₂ dissociation and protonation of 1.40 eV on sp-N1GDY(OH)/GDY; thus the ET-OHP associative mechanism is assigned for ORR on this surface. On the other hand, for sp-N2GDY(OH)/GDY, O₂ dissociation is more favorable than O₂ protonation; thus, we assume ORR will proceed based on the dissociative mechanism.

The free energy diagram, reaction, and adsorption Gibbs free energy of ORR intermediates on sp-N1GDY(OH)/GDY and sp-N2GDY(OH)/GDY are shown in Fig. 3(c), Fig. S4(c), and Tables S4–S6 (ESI†). For sp-N1GDY(OH)/GDY, the stabilities of the three ORR intermediates, i.e., OOH*, O*, and OH* were estimated. We find that all ORR intermediates preferably adsorb on top sites of the neighboring C atom of the –COH–N moiety. Interactions between sp-N1GDY(OH) and intermediates become weak after involving the pre-adsorbed OH*. The ΔG_OH* of sp-N1GDY(OH)/G is 1.15 eV larger than that of sp-N1GDY/G (−0.13 eV), indicating that the pre-adsorbed OH* is improving rather than poisoning. Moreover, the PDS is O₂ → OOH* with a limiting potential of 0.51 eV and η of 0.72 V (Table 2). As shown in Fig. 3(b) and (d), we can also prove that under U_RHE = 0.51 V, sp-N1GDY(OH) is the most stable structure. For sp-N2GDY(OH)/G, we found 2O* → O* + OH* is a strongly endergonic reaction with ΔG₂ of 1.06 eV (Fig. S4(c) and Table S5, ESI†), which may cause possible poisoning of active sites by 2O*. Therefore, sp-N2GDY/G cannot be used as an ORR electrocatalyst.

Table 2 Overpotential (η) and potential-determining step (PDS) of ORR on sp-N1GDY(OH), sp-N1GDY(OH)/G, Pyri-NGDY, and Pyri-NGDY/G with/without water using RPBE + D2 and PBE + D2, respectively

	RPBE + D2		PBE + D2
	η	PDS	η	PDS
sp-N1GDY(OH)	0.75	O2 → OOH*	—	—
sp-N1GDY(OH)/G	0.72	O2 → OOH*	0.54	O2 → OOH*
sp-N1GDY(OH)/G with water	0.46	O2 → OOH*	0.54	OH* → OH^−
Pyri-NDGY	1.17	O2 → OOH*	—	—
Pyri-NDGY/G	1.01	O2 → OOH*	0.89	O2 → OOH*
Pyri-NDGY/G with water	0.75	O2 → OOH*	0.65	O* → OH*
Experiment^20		η = 0.36

---

Under vacuum conditions, we found that sp-N1GDY(OH)/G can create a conceivable active site for ORR with low overpotential. The clean sp-N1GDY/G surfaces are easily terminated by OH* intermediate and the neighboring C site of -the COH-N moiety of sp-N1GDY(OH)/G is the real active site for ORR to proceed with low η compared with Pyri-NGDY/G.

The electronic structure of the active site relates to O₂ activation. In general, the interactions between C 2p of the active sites and O₂ π* orbitals govern the adsorption strength of O₂ and the elongation of O₂ upon adsorption on the substrate. Hybridizations between C 2p and O₂ π* orbitals facilitate the back donation to π* orbitals, thus increasing the occupations of these orbitals upon adsorption (Fig. S6, ESI†). As a result, the PDOS of C 2p near the Fermi level is quite important for determining the O₂ adsorption and O₂ dissociation barriers.

We plotted the PDOS of C2p of the active site in sp-N1GDY/G, sp-N2GDY/G, sp-N1GDY(OH)/G, sp-N2GDY(OH)/G, and Pyri-NGDY/G. In Fig. 4(a), at the Fermi level, there is almost no partially occupied state of C 2p in sp-N1GDY(OH)/G and Pyri-NGDY/G, while for sp-N1GDY/G, sp-N2GDY/G, and sp-N2GDY(OH)/G, there exists partially occupied state of C 2p. From Fig. 4(b), we also found that there is a linear relationship between the PDOS height value at the Fermi level and the adsorption energy of O₂. The results can also prove that the presence of a partially occupied state of the C 2p of the active site can induce O₂ adsorption. Upon adsorption of O₂, there is almost no hybridization in sp-N1GDY(OH)/G and Pyri-NGDY/G, thus the adsorption energies are weak (−0.15 eV and −0.10 eV) and O₂ dissociation barrier is high (1.9 eV and 1.57 eV). On the other hand, O₂ has strong hybridizations with C on sp-N1GDY/G, sp-N2GDY/G, and sp-N2GDY(OH)/G due to the presence of partially occupied states. The adsorption energy (−0.6 eV, −1.37 eV, and −0.58 eV) and NEB barrier (0.22 eV, 0.08 eV, and 0.81 eV) results show that O₂ is chemisorbed and easily dissociates on sp-N1GDY/G, sp-N2GDY/G, and sp-N2GDY(OH)/G, respectively.

Fig. 4 (a) PDOS of C 2p of the active site in sp-N1GDY/G, sp-N1GDY(OH)/G, sp-N2GDY/G, sp-N2GDY(OH)/G, and Pyri-NGDY/G. (b) The relationship between PDOS height value at the Fermi level and the adsorption energy of O₂.

3.3 ORR mechanism and free energy analysis on sp-N1GDY/G and Pyri-NGDY/G under water conditions

Experimentally, electrochemical reactions occur at solid–liquid interfaces, thus it is necessary to incorporate the solvation effects in an explicit water environment when investigating the ORR using the CHE model. The solvation energy of the ORR intermediate can be estimated from the AIMD simulation of the adsorbed systems with water solvents explicitly. However, this approach incurs a huge computational cost. To this end, we first simulated the interfaces of NGDY's, namely sp-N1GDY/G and Pyri-NGDY/G with water to elucidate the H bond networks. Next, we only kept the water configurations in the contact region with the surface of the five most stable AIMD snapshots of the clean systems and replaced one nearest water molecule above the active site by each ORR intermediate to construct the H bond networks of ORR intermediate and water. The ΔSol of each ORR intermediate is then elucidated from an average of five AIMD snapshots. Moreover, we also constructed H bond networks using the ice-like bilayer model. The details of the AIMD analysis are summarized in the ESI,† Fig. S7–S9, and Table S7.

H bond networks of water on sp-N1GDY/G and Pri-NGDY/G surfaces are shown in Fig. S10 (ESI†). On both substrates, the water layer of the H bond networks at the contact region with N-doped GDY is mainly composed of six-membered ring components after optimization, i.e., one H₂O has three H bonds with three neighboring H₂O. This is similar to the ice-like bilayer that is often observed for water/flat metal interfaces.^28,61 The ΔSol of 2O*, O* + OH*, O*, and OH* on sp-N1GDY/G and the ΔSol of OOH*, O*, and OH* on sp-N1GDY(OH)/G or Pyri-NGDY/G are shown in Tables S8 and S9 (ESI†). We found that water stabilizes all the ORR intermediates, which arise from H bonds with water. Moreover, we also calculated the ΔSol of each ORR intermediate using an ice-like bilayer (Tables S8 and S9, ESI†). The results showed that the difference in the ΔSol‘s between the H bond networks model from AIMD and the ice-like bilayer model varies by only ∼0.06 eV. We found that the ΔSol of these surfaces is independent of the water model, which arises from a similar H bond between adsorbates and water in the ice-like bilayer and AIMD H bond networks.

The free energy diagrams of ORR with water on sp-N1GDY/G, sp-N1GDY(OH)/G, and Pyri-NGDY/G are shown in Fig. 5. Because of the stabilization driven by the solvation effect, the reaction Gibbs free energies of all ORR intermediates are more stable than those under vacuum conditions. Thus, OH* is more easily terminated on the sp-N1GDY/G surface (Fig. 5(a) and (b)), and further ORR steps are considered on sp-N1GDY(OH)/G (Fig. 5(c) and (d)). The PDS of ORR on sp-N1GDY(OH)/G in water remains the same as that under vacuum conditions (O₂ → OOH*), while the limiting potential is changed from 0.51 V under vacuum condition to 0.77 V under water conditions. At U_RHE = 0.77 V (Fig. 5(b) and 5(d)), sp-N1GDY(OH) is the most stable structure under water conditions. Therefore, introducing ΔSol, the η of sp-N1GDY(OH)/G is 0.46 V (Table 2), which is close to the experimentally reported value (0.36 V). Similarly, on Pyri-NGDY/G, even though the solvation effect is included, the PDS of ORR on Pyri-NGDY/G in water conditions does not change compared with that under vacuum conditions (O₂ → OOH*). We obtained a limiting-potential increase to 0.48 eV (Fig. 5(e) and (f)) and an η decrease to 0.75 V (Table 2), which is higher than that of sp-N1GDY(OH)/G. Finally, we conclude that ORR easily proceeds on sp-N1 doping with pre-adsorbed OH*.

Fig. 5 (a) Free energy diagram and structures of each ORR intermediate with water on sp-N1GDY/G. (b) Free energy of each ORR intermediate with water versus electrode potential (vs. RHE) on sp-N1GDY/G. (c) Free energy diagram and structures of each ORR intermediate with water on sp-N1GDY(OH)/G. (d) Free energy of each ORR intermediate with water versus electrode potential (vs. RHE) on sp-N1GDY(OH)/G. (e) Free energy diagram and structures of each ORR intermediate with water on Pyri-NGDY/G. (f) Free energy of each ORR intermediate versus electrode potential (vs. RHE) on Pyri-NGDY/G. Pristine in (b), (d), and (f) represents sp-N1GDY/G, sp-N1GDY(OH), and Pyri-NGDY/G, respectively. U_RHE is the potential of the electrode relative to the RHE. Red, green, white, gray, pink, and blue balls are O atoms in ORR intermediates, O atoms in water, H atoms, C atoms in N-doped GDY, and C atoms in G, and N atoms, respectively. The purple dashed lines in structures of (a), (c), and (e) represent the H bonds. Purple circles in (b), (d), and (f) represent the lowest lines crossed at different potentials. We regard an H bond as being formed when the O–O distance of adjacent water molecules is smaller than 3.5 Å and the angle between the O–H vector of one molecule and the O–O vector is smaller than 30° in ref. 62.

3.4 Effects of exchange–correlation energy functionals and graphene support on ORR

Here, we discuss the effect of exchange–correlation energy functionals on the adsorption Gibbs free energy of ORR. We employed the PBE + D2 energy functional to perform the adsorption Gibbs free energy of ORR on sp-N1GDY/G, sp-N1GDY(OH)/G, and Pyri-NGDY/G with and without water. As shown in Table S10 (ESI†), we found that PBE tends to overestimate the binding energy of chemisorption species. This is due to the overestimation of attractive interactions in molecular systems. Under vacuum conditions, the adsorption Gibbs free energies of ORR intermediates on sp-N1GDY/G, sp-N1GDY(OH)/G, and Pyri-NGDY/G are more stable using PBE + D2 compared with RPBE + D2. Thus, the η of ORR on sp-N1GDY(OH)/G and Pyri-NGDY/G is slightly decreased using PBE + D2 compared with RPBE + D2. However, the trend of ORR activity remains the same, which is sp-N1GDY(OH)/G > Pyri-NGDY/G. Under water conditions, we only used the ice-like bilayer water to simulate the ΔSol because we already proved that the ΔSol using the H bond networks from AIMD is not much different from that using the ice-like bilayer water based on RPBE + D2 energy functional. Table S11 (ESI†) shows that the ΔSol of each ORR intermediate is more negative using the PBE + D2 energy functional compared with that using the RPBE + D2 energy functional. For Pyri-NGDY/G with water, η is decreased to 0.65 V (Table 2) and PDS is changed from O₂ → OOH* (vacuum conditions) to O* → OH* (water conditions). For sp-N1GDY(OH)/G, PDS is changed from O₂ → OOH* (vacuum conditions) to OH* → OH⁻ (water conditions) but η is coincidentally not changed (0.54 V) (Table 2). Although the trend of ORR activity in water is unchanged, the η of sp-N1GDY(OH)/G using RPBE + D2 (0.46 V) is closer to the experimental result (0.36 V) compared with that using PBE + D2 (0.54 V). Therefore, the RPBE + D2 energy functional gives more reasonable binding energies and solvation energies of adsorbates.

In Section 3.1, we discussed that the G support on NGDY enhances the metallic property and causes charge transfer from G to NGDY. To investigate the effect of G support on ORR electrocatalytic activity, calculated adsorption Gibbs free energies of ORR with G support and without G support are shown in Tables S5 and S6 (ESI†), respectively. The results revealed that the G support stabilizes all ORR intermediates in these cases, especially on Pyri-NGDY, which stabilizes ∼0.2 eV on OOH*, O*, and OH* (Table S6, ESI†). Higher charge transfer from G to the substrate as indicated by work function change upon graphene support (0.31 eV for Pyri-NGDY/G vs. 0.13 eV for sp-N1GDY/G) results in a strong interaction of ORR intermediate with N-doped GDY. Briefly, the introduction of the G support enhances the ORR electrocatalytic activity of N-doped GDY.

4. Conclusions

In this study, we performed DFT calculations and AIMD simulations to investigate the ORR mechanism on the sp-N1GDY/G and Pyri-NGDY/G with and without the solvation effect. We obtained the following important results:

(1) Under both vacuum and water conditions, ORR firstly proceeds on sp-N1GDY/G via a dissociative mechanism because O₂ can be chemisorbed on a clean surface and easily dissociated rather than protonated to OOH*. However, OH* is strongly adsorbed on the sp-N1GDY/G surface, resulting in the weakening of the second O₂ adsorption, and ORR occurs via the ET-OHP associative mechanism. Pyri-NGDY/G also prefers the ET-OHP associative mechanism.

(2) From the AIMD simulation, the H bond networks at the contact region of water and NGDY are mainly composed of six-membered ring components. H bond with water stabilizes each ORR intermediate and the free energy diagram with solvation effect agrees well with the experimental data. Moreover, we found that using the ice-like bilayer model to construct the H bond networks can also give a reasonable estimation of ΔSol. Therefore, the ΔSol of these surfaces is independent of the water model, which arises from a similar H bond between adsorbates and water in the ice-like bilayer and AIMD H bond networks.

(3) sp-N1GDY/G with OH* pre-adsorbed surface has the highest ORR electrocatalytic activity and the neighboring C site of –COH–N moiety is the active site for ORR. Incorporation of the solvation effect is of importance because η with the solvation effect (0.46 V) is much closer to the experimental one (0.36 V). Our work highlights the importance of considering solvation effects in designing and optimizing catalysts for ORR and other chemical reactions.

Conflicts of interest

There are no conflicts to declare.

Acknowledgements

This work was supported by a Grant-in-Aid for Transformative Research Areas (A) “Hyper-Ordered Structure Science” (grant no. JP20H05883) and a Grant-in-Aid for Scientific Research (B) (grant no. JP20H02569) and JSPS Fellows (22J10937) from the Japan Society for the Promotion of Science (JSPS) and Elements Strategy Initiative for Catalysts and Batteries (ESICB) (grant no. JPMXP0112101003) from the Ministry of Education, Culture, Sports, Science, and Technology, Japan (MEXT). The computations were performed using computer resources at the Research Center for Computational Science, Okazaki; the Cyberscience Center, Tohoku University; the Center for Computational Materials Science, Institute for Materials Research, Tohoku University; the Institute for Solid State Physics, University of Tokyo. This work was also supported by MEXT as “Program for Promoting Researches on the Supercomputer Fugaku” (Fugaku battery & Fuel Cell Project) and used the computational resources of the supercomputer Fugaku provided by the RIKEN Center for Computational Science.

Notes and references

1. B. C. H. Steele and H. Angelika, Nature, 2001, 414, 345–352.
2. J. Kui, J. Xuan, Q. Du, Z. Bao, B. Xie, B. Wang, Y. Zhao, L. Fan, H. Wang, Z. Hou, S. Hu, N. P. Brandon, Y. Yin and M. D. Guiver, Nature, 2021, 595, 361–369.
3. D. Larcher and J. M. Tarascon, Nat. Chem., 2015, 7, 19–29.
4. R. Ma, G. Lin, Y. Zhou, Q. Liu, T. Zhang, G. Shan, M. Yang and J. Wang, npj Comput. Mater., 2019, 5, 78.
5. S. Sui, X. Wang, X. Zhou, Y. Su, S. Riffat and C. J. Liu, J. Mater. Chem. A, 2017, 5, 1808–1825.
6. L. Dai, Y. Xue, L. Qu, H. J. Choi and J. B. Baek, Chem. Rev., 2015, 115, 4823–4892.
7. K. Gong, F. Du, Z. Xia, M. Durstock and L. Dai, Science, 2009, 323, 760–764.
8. D. Guo, R. Shibuya, C. Akiba, S. Saji, T. Kondo and J. Nakamura, Science, 2016, 351, 361–365.
9. S. Tang, X. Zhou, T. Liu, S. Zhang, T. Yang, Y. Luo, E. Sharman and J. Jiang, J. Mater. Chem. A, 2019, 7, 26261–26265.
10. S. Tang, Q. Dang, T. Liu, S. Zhang, Z. Zhou, X. Li, X. Wang, E. Sharman, Y. Luo and J. Jiang, J. Am. Chem. Soc., 2020, 142, 19308–19315.
11. C. Z. Yuan, H. B. Li, Y. F. Jiang, K. Liang, S. J. Zhao, X. X. Fang, L. B. Ma, T. Zhao, C. Lin and A. W. Xu, J. Mater. Chem. A, 2019, 7, 6894–6900.
12. L. Jian, X. Gao, B. Liu, Q. Feng, X. B. Li, M. Y. Huang, Z. Liu, J. Zhang, C. H. Tung and L. Z. Wu, J. Am. Chem. Soc., 2016, 138, 3954–3957.
13. G. Xin, H. Liu, D. Wang and J. Zhang, Chem. Soc. Rev., 2019, 48, 908–936.
14. Y. Fang, Y. Liu, L. Qi, Y. Xue and Y. Li, Chem. Soc. Rev., 2022, 51, 2681–2709.
15. Y. Li, L. Xu, H. Liu and Y. Li, Chem. Soc. Rev., 2014, 43, 2572–2586.
16. C. Huang, Y. Li, N. Wang, Y. Xue, Z. Zuo, H. Liu and Y. Li, Chem. Rev., 2018, 118, 7744–7803.
17. Z. Zuo and Y. Li, Joule, 2019, 3, 899–903.
18. Y. Gao, Z. Cai, X. Wu, Z. Lv, P. Wu and C. Cai, ACS Catal., 2018, 8, 10364–10374.
19. Z. Feng, R. Li, Y. Ma, Y. Li, D. Wei, Y. Tang and X. Dai, Phys. Chem. Chem. Phys., 2019, 21, 19651–19659.
20. Y. Zhao, J. Wan, H. Yao, L. Zhang, K. Lin, L. Wang, N. Yang, D. Liu, L. Song, J. Zhu, L. Gu, L. Liu, H. Zhao, Y. Li and D. Wang, Nat. Chem., 2018, 10, 924–931.
21. T. Lu, X. Hu, J. He, R. Li, J. Gao, Q. Lv, Z. Yang, S. Cui and C. Huang, Nano Energy, 2021, 85, 106024.
22. X. Chen, W. J. Ong, Z. Kong, X. Zhao and N. Li, Sci. Bull., 2020, 65, 45–54.
23. B. Kang, S. Wu, J. Ma, H. Ai and J. Y. Lee, Nanoscale, 2019, 11, 16599–16605.
24. J. Gu, S. Magagula, J. Zhao and Z. Chen, Small Methods, 2019, 1800550.
25. J. K. Nørskov, J. Rossmeisl, A. Logadottir, L. Lindqvist, J. R. Kitchin, T. Bligaar and H. Jónsson, J. Phys. Chem. B, 2004, 108, 17886–17892.
26. L. Yu, X. Pan, X. Cao, P. Hu and X. Bao, J. Catal., 2011, 282, 183–190.
27. P. S. Rice, Y. Mao, C. Guo and P. Hu, Phys. Chem. Chem. Phys., 2019, 21, 5932–5940.
28. S. Schnur and A. Groß, New J. Phys., 2009, 11, 125003.
29. A. Groß and S. Sung, Chem. Rev., 2022, 122, 10746–10776.
30. S. N. Steinmann, P. Sautet and C. Michel, Phys. Chem. Chem. Phys., 2016, 18, 31850–31861.
31. J. Ho and M. Z. Ertem, J. Phys. Chem. B, 2016, 120, 1319–1329.
32. G. König and B. Stefan, J. Phys. Chem. B, 2009, 113, 8967–8974.
33. H. Cao, G. J. Xia, J. W. Chen, H. M. Yan, Z. Huang and Y. G. Wang, J. Phys. Chem. C, 2020, 124, 7287–7294.
34. M. Reda, H. A. Hansen and T. Vegge, Catal. Today, 2018, 312, 118–125.
35. Y. Morikawa, Phys. Rev. B: Condens. Matter Mater. Phys., 1995, 51, 14802–14805.
36. Y. Hamamoto, I. Hamada, K. Inagaki and Y. Morikawa, Phys. Rev. B: Condens. Matter Mater. Phys., 2016, 93, 245440.
37. Y. Morikawa, K. Iwata and K. Terakura, Appl. Surf. Sci., 2000, 11, 169–170.
38. S. A. Wella, Y. Hamamoto, F. Iskandar, Suprijadi, Y. Morikawa and I. Hamada, J. Chem. Phys., 2020, 152, 104707.
39. D. Vanderbilt, Phys. Rev. B: Condens. Matter Mater. Phys., 1990, 41, 7892–7895.
40. B. Hammer, L. B. Hansen and J. K. Nørskov, Phys. Rev. B: Condens. Matter Mater. Phys., 1999, 59, 7413.
41. S. Grimme, J. Comput. Chem., 2006, 27, 1787–1799.
42. G. Henkelman and H. Jonsson, J. Chem. Phys., 2000, 113, 9978–9985.
43. G. Henkelman, B. P. Uberuaga and H. Jónsson, J. Chem. Phys., 2000, 113, 9901–9904.
44. S. P. Chan, G. Chen, X. G. Gong and Z. F. Liu, Phys. Rev. Lett., 2003, 90, 086403.
45. J. I. Enriquez, F. Muttaqien, M. Michiuchi, K. Inagaki, M. Geshi, I. Hamada and Y. Morikawa, Carbon, 2021, 174, 36–51.
46. J. K. Nørskov, F. Studt, F. Abild-Pedersen and T. Bligaard, Fundamental Concepts in Heterogeneous Catalysis, Wiley, New York, NY, 2014.
47. R. Iftimie, P. Minary and M. E. Tuckerman, Proc. Natl. Acad. Sci. U. S. A., 2005, 102, 6654–6659.
48. X. Chen, Y. Zhang, Y. Ren, D. Wang and J. Yun, Mater. Res. Express, 2019, 6, 095610.
49. V. M. Karpan, G. Giovannetti, P. A. Khomyakov, M. Talanana, A. A. Starikov, M. Zwierzycki, J. V. D. Brink, G. Brocks and P. J. Kelly, Phys. Rev. Lett., 2007, 99, 176602.
50. A. S. Haile, H. A. Hansen, W. Yohannes and Y. S. Mekonnen, J. Phys. Chem. Lett., 2021, 12, 3552–3559.
51. J. P. Perdew, K. Burke and M. Ernzerhof, Phys. Rev. Lett., 1996, 77, 3865–3868.
52. J. Heyd and G. E. Scuseria, J. Chem. Phys., 2004, 121, 1187–1192.
53. J. Heyd, G. E. Scuseria and M. Ernzerhof, J. Chem. Phys., 2003, 118, 8207–8215.
54. G. Kresse and J. Furthmüller, Comput. Mater. Sci., 1996, 6, 15–50.
55. G. Kresse and J. Furthmüller, Phys. Rev. B: Condens. Matter Mater. Phys., 1996, 54, 11169.
56. M. Yan, Z. Q. Huang, Y. Zhang and C. R. Chang, Phys. Chem. Chem. Phys., 2017, 19, 2364–2371.
57. T. N. Pham, M. Sugiyama, F. Muttaqien, S. E. M. Putra, K. Inagaki, D. N. Son, Y. Hamamoto, I. Hamada and Y. Morikawa, J. Phys. Chem. C, 2018, 122, 11814–11824.
58. N. Ramaswamy and S. Mukerjee, J. Phys. Chem. C, 2011, 115, 18015–18026.
59. C. H. Choi, H. K. Lim, M. W. Chung, J. C. Park, H. Shin, H. Kim and S. I. Woo, J. Am. Chem. Soc., 2014, 136, 9070–9077.
60. Y. Wang, Y. J. Tang and K. Zhou, J. Am. Chem. Soc., 2019, 141, 14115–14119.
61. A. Groß, F. Gossenberger, X. Lin, M. Naderian, S. Sakong and T. Roman, J. Electrochem. Soc., 2014, 161, E3015.
62. A. Luzar and D. Chandler, Nature, 1996, 379, 55–57.

---

Footnote

† Electronic supplementary information (ESI) available. See DOI: https://doi.org/10.1039/d3ma00502j

---