Introduction

The exploration of clean and renewable energy materials as well as their devices becomes an urgent demand when being confronted with the energy crisis and environmental problems1,2. Owing to high energy density and environmental friendliness, metal-air batteries and fuel cells have been become promising candidates for energy storage devices2,3,4. The slow oxygen reduction reaction (ORR) kinetics at the cathode, however, severely impedes their practical applications. Despite excellent ORR kinetic activity of Pt catalysts, the rare reserve and high cost heavily limit their practical use5,6. Therefore, developing high-efficiency, low-cost, and well durable alternatives to Pt catalysts still remains a persistent challenge for electrocatalysis.

Recently, dual-atom catalysts (DACs), especially N-doped graphene-based dual transition metal atoms systems (denoted as M1M2-N-Gs), have shown great application prospects7,8,9,10,11. Since the coordination structure of active site determines the catalytic activity11,12,13,14,15, more efforts focus on the improvement of the catalytic performances of M1M2-N-Gs through changing the local environment of active site. It is known that dispersing the active site at graphene edge is a potent approach for boosting the catalytic activity of graphene-based catalysts towards various chemical reactions. For example, the defective NiN3 active site at the edge achieved a carbon dioxide reduction reaction (CO2RR) performance and durability superior to the interior NiN3 active site16. Atomically dispersed FeN4 active sites near graphitic edges also exhibited higher ORR activity and durability compared with commercial Pt/C catalyst17,18,19. Density-functional theory (DFT) calculations further verified that the enhanced ORR activity originates from tilted-FeN4 configurations at edges20. The edges with vacancy defects lowered the energy barrier of ORR on CoZnN6- and FeCoN6-Gs, thus enhancing their catalytic activities21,22,23. Besides widely investigated Fe- and Co-based DACs11,24, NiNiN6-G25, NiCuN6-G26, and NiZnN6-G27 were screened out by combining DFT calculations and machine learning techniques, demonstrating their good ORR catalytic activities. Pd, as a congener element of Ni, is expected to have a synergistic effect with Ni atom, thus achieving the suitable electronic structures to adsorb ORR intermediates. Moreover, Ni-Pd nanoparticles impregnated on nitrogen-doped graphene and reduced graphene oxide graphene have been confirmed to have outstanding ORR catalytic performances28,29,30. Therefore, despite the unsatisfactory catalytic activities of NiN4 and PdN4 embedded graphenes31,32,33, NiPdN6-G may possess good activity for ORR and their catalytic performance could be further enhanced by edge terminal engineering. Furthermore, due to the multiple possibilities of coordination structures of edge terminated NiPdN6-G, there still lacks a fundamental understanding on how the structure of edge termination affects the catalytic activity.

Additionally, O doping may modulate the electronic structures of active sites and promote the ORR activity. For example, DFT calculations demonstrated that O-doped graphene supported Pt catalyst has excellent ORR catalytic activity among a series of B, N, and O doped structures34. O and N co-coordination can induce the shifts in the d-band center of Mo atom towards the Fermi-level, resulting in the prominent ORR activity of MoN2O2-G35. Recently, synthesized FeCuN4O2-G and FeCuN5O-G presented remarkably lowered activation energy and enhanced intrinsic ORR activity compared with FeCuN6-G12,36,37. Nevertheless, whether O doping at edge-terminated NiPdN6-G can further improve their catalytic performance towards ORR remains unknown. Hence, there is an urgent demand to comprehensively explore the modulation of ORR activity for NiPdN6-G through edge termination and O dopant.

Herein, in this work we performed DFT calculations to thoroughly investigate the electrochemical activity of NiPdN6-G with and without edge termination and O doping. Our results show that edge termination and O doping can effectively promote the catalytic activity for ORR, and the lowest overpotential is only 0.31 V. This work provides a theoretical guidance for rational design of low-cost and high-efficiency graphene-based DACs.

Computational method

We performed the spin-polarized DFT calculations with the Vienna ab initio simulation package (VASP)38,39. The Perdew-Burke-Ernzerhof (PBE) functional method within the generalized gradient approximation (GGA) was used to treat electronic exchange–correlation interactions40. M1M2N6 respectively located in the interior and edge sites of graphene were constructed by using 6 × 6 graphene supercell containing 72 C atoms and graphene nanoribbon containing 72 C and 12 H atoms (see Fig. S1), respectively. The kinetic energy cutoff was set as 500 eV to keep the energy convergence. All structures were allowed to relax until the electronic structure iteration and atomic force were respectively smaller than 10−5 eV and 0.02 eV/Å. The sufficient vacuums of 15 Å in the c direction for interior Ni-Pd catalyst and about 15 Å in the b direction for edge-hosted Ni-Pd catalysts were used to avoid the interactions between the models and their mirror images. The k-meshes of 3 × 3 × 1 and 3 × 1 × 1 were used for the structural optimizations of interior and edge-hosted structures, respectively. The van der Waals interactions between catalysts and reaction intermediates were described using the DFT + D3 dispersion correction method41.

The formation energy (Ef) of all catalysts was calculated according to the equation as follow

$${E}_{\text{f}}={E}_{\text{t}}-{\text{E}}_{\text{Ni}}-{E}_{\text{Pd}}-{n\mu }_{\text{X}},$$
(1)

where Et, ENi, and EPd are the energies of catalysts, isolated Ni and Pd atoms, respectively. μX represents the chemical potentials of C, N, O, or H atom, which is defined as the energy of C atom in perfect graphene, N atom in N2 molecule, O atom in O2 molecule, or H atom in H2 molecule, respectively. n denotes the numbers of C, N, O, or H atoms.

The adsorption energy (Eads) of reaction intermediate on catalyst was defined as

$${E}_{\text{ads}}={E}_{\text{catal}+\text{inter}}-{E}_{\text{inter}}-{E}_{\text{catal}},$$
(2)

where \({E}_{\text{catal}+\text{inter}}\), \({E}_{\text{catal}}\), and \({E}_{\text{inter}}\) respectively represent the energies of catalyst with adsorbed reaction intermediate, catalyst, and isolated reaction intermediate.

The free energy change (ΔG) of each reaction step was obtained according to the equation as follow

$$\Delta G = \Delta E + \Delta ZPE - T\Delta S + \Delta G_{U} + \Delta G_{{pH}} ,$$
(3)

where ΔE, ΔZPE, and ΔS correspond to the changes in total energy, zero-point energy, and entropy of each reaction step, respectively. T is set to be 298.15 K, corresponding to the room temperature. ΔGU =  − meU, in which m denotes the number of electrons transferred during the reaction process, and U denotes the electrode potential. ΔGpH is equal to kBT × ln10 × pH, where kB is the Boltzmann constant, and pH is set to be 0 in an acidic medium.

For O2 → *OOH → *O + H2O → *OH + H2O → 2H2O along the four-electron reaction pathway, the free energy change (ΔG) for each reaction step was respectively calculated according to the equations as follows

$$\Delta G_{{1}} = \Delta G_{{*{\text{OOH}}}} {-}{ 4}.{92 } + eU,$$
(4)
$$\Delta G_{{2}} = \Delta G_{{*{\text{O}}}} {-} \Delta G_{{*{\text{OOH}}}} + eU,$$
(5)
$$\Delta G_{{3}} = \Delta G_{{*{\text{OH}}}} {-}\Delta G_{{*{\text{O}}}} + eU,$$
(6)
$$\Delta G_{{4}} = \, {-}\Delta G_{{*{\text{OH}}}} + eU.$$
(7)

Finally, the overpotential (η) was obtained by the following equation

$$\eta = { 1}.{23 } + {\text{ max}}(\Delta G_{{1}} ,\Delta G_{{2}} ,\Delta G_{{3}} ,\Delta G_{{4}} )/{\text{e}}.$$
(8)

Results and discussion

Integrating N6-coordinated M1-M2 active site into the interior, zigzag-edge and armchair-edge of graphene can form one interior (M1M2N6-I), six zigzag-edge (M1M2N6-A1 ~ A6), and four armchair-edge (M1M2N6-Z1 ~ Z4) configurations, as shown in Fig. 1. However, due to the differences in the local N coordination environments of M1 and M2 atoms located in the edge, the exchange of both atoms can form different active site. Therefore, there exist twenty-one possible Ni-Pd co-doped N6-coordinated graphene structures in total, namely NiPdN6-I, NiPdN6-A1 ~ A6, PdNiN6-A1 ~ A6, NiPdN6-Z1 ~ Z4, and PdNiN6-Z1 ~ Z4. To identify their thermodynamic stability, the formation energy (Ef) was illustrated in Fig. 2 and Table S1. One can find that for all structures, the Ef are negative, implying that they are structurally stable and plausible to exist on the graphene. It is noted that, except NiPdN6-A6 and PdNiN6-A6 configurations, other NiPdN6 and PdNiN6 structures in armchair-edge have smaller Ef than those in zigzag-edge, so the probability of these structures in armchair-edge should be higher than that in zigzag-edge. Moreover, the Ef of these structures in armchair-edge were close to or less than that of the standard NiPdN6-I structure, suggesting that the NiPdN6 and PdNiN6 structures may be prone to be formed in the armchair-edge. In addition, the exchange of Ni and Pd atoms led to certain differences of Ef (0.02 ~ 0.41 eV) for each edge structure. Among which, the energy differences between NiPdN6-A2/A4/A6/Z2/Z4 and PdNiN6-A2/A4/A6/Z2/Z4 are smaller than 0.1 eV, while NiPdN6-A1/A3/Z3 and PdNiN6-A5/Z1 have markedly smaller Ef than that of PdNiN6-A1/A3/Z3 and NiPdN6-A5/Z1. Therefore, except PdNiN6-A1/A3/Z3 and NiPdN6-A5/Z1, all the NiPdN6 and PdNiN6 structures were considered in the following discussion to identify their ORR catalytic activities.

Fig. 1
figure 1

Schematic diagram of M1M2N6 located in the interior (-I), armchair-edge (− A1 ~ A6), and zigzag-edge (− Z1 ~ Z4) of graphene. The pink, green, blue, yellow, and balls denote H, C, N, M1, and M2 atoms, respectively.

Fig. 2
figure 2

The formation energy of M1M2N6 with different active sites.

Since the adsorption of reaction intermediate has a great effect on the ORR catalytic activity, the adsorption characteristics of all ORR intermediates (*O2, *OOH, *O, and *OH) were firstly investigated. By optimizing their adsorptions on all catalysts, the most stable adsorption structures were illustrated in Fig. 3a. One can find that the *O2 and *OOH are only stably adsorbed on Ni atom in all catalysts, while *O and O in *OH may bond with single Ni or both Ni and Pd atoms. Therefore, Ni atom is the main active site for the adsorption of ORR intermediates, which was confirmed by the fact that Ni atom loses more charges than Pd atom in all NiPdN6 and PdNiN6 structures (see Table S2). Considering that the adsorption strength to reactants directly affects the activity of catalysts, Eads for these intermediates were further calculated by Eq. (2) and listed in Table S3. By comparison with NiPdN6-I, the Eads for all intermediates were changed at the edge active sites to different extent. According to available studies, M1M2-N-Gs with the Eads for O2 in the range of − 0.24 to − 0.80 eV may possess good ORR activity25,42,43,44,45,46. Evidently, all Eads of NiPdN6 and PdNiN6 for O2 in this work (− 0.38 ~ − 0.72 eV) fall in the above range, thus some excellent ORR catalysts are likely to exist in these structures.

Fig. 3
figure 3

(a) Schematic diagram for the intermediate adsorption structures. (b) The overpotentials of ORR on NiPdN6 and PdNiN6 structures. The yellow, blue, red, and pink balls denote Ni, Pd, O, and H atoms, respectively.

To identify the ORR activity of all NiPdN6 structures, the free energy change for each reaction step, the free energy diagram, and the overpotential were calculated by Eqs. (3)(8), as shown in Table S3 and Figs. S2 and 3b. For NiPdN6-I, the rate-determining step (RDS) of ORR is the hydrogenation of *O to *OH, and the corresponding overpotential is 0.49 V. This value is slightly higher than those in Pt (111) (0.44 V)47 and Pt (100) (0.47 V)48, indicating the relatively lower ORR catalytic activity of NiPdN6-I than Pt catalyst. After introducing the edge active sites, the RDS of ORR was transformed into the formation of the second H2O molecule for NiPdN6-A2, PdNiN6-A6, PdNiN6-Z1, NiPdN6-Z4, and PdNiN6-Z4 as well as the hydrogenation of O2 to *OOH for NiPdN6-A4, PdNiN6-A4, NiPdN6-Z2, and NiPdN6-Z3. The overpotentials of ORR on NiPdN6-A2, PdNiN6-A6, NiPdN6-Z2, PdNiN6-Z2, and PdNiN6-A4 were respectively lowered to 0.41, 0.42, 0.45, 0.45, and 0.46 V. However, the overpotential was increased in different extent for other structures, which should be attributed to the changed Eads for the ORR intermediates at edge active sites, as indicated by Table S3. Therefore, NiPdN6-A2 and PdNiN6-A6 possess excellent ORR catalytic activity superior to Pt catalyst, while the ORR catalytic activity of PdNiN6-A4, NiPdN6-Z2, and PdNiN6-Z2 is comparable to that of Pt catalyst.

Exploring the scaling relationship between the Gibbs adsorption free energy of reaction intermediates is significant for understanding the catalytic characteristics and developing more effective ORR catalysts. As illustrated in Fig. 4a, ΔG*OH and ΔG*OOH for all catalysts show a clear linear scaling relationship (ΔG*OOH = 0.74ΔG*OH + 3.28 eV, and the square of the correlation coefficient R2 is 0.93), indicating their strong relationships. However, the correlations between ΔG*O and ΔG*OHG*OOH are not impressive, and the R2 respectively are 0.84 and 0.74, as shown in Fig. 4b and c. According to the relationships between these intermediates, an approximately classic volcano curve between ΔG*OH and overpotential was presented in Fig. 4d. Overall, as ΔG*OH increases, the overpotential first reduces and subsequently rises. When ΔG*OH increases to 0.94 eV, the overpotential reaches the minimum (0.29 V), corresponding to the optimal ORR catalytic activity. According to the adsorption free energy of reaction intermediates and the identified RDS, all catalysts can be divided into three regimes under this volcano-like trend. There is a strictly linear left boundary of the volcano curve, and the overpotential is purely dictated by the ΔG*OH. The catalysts in this boundary have strong binding strength to *OH. However, the catalyst sites are more scattered on the right boundary. Among which, the RDS of NiPdN6-A4, PdNiN6-A4, NiPdN6-Z2, and NiPdN6-Z3 is the hydrogenation of O2 to *OOH. Although the hydrogenation of *O to *OH is the RDS of NiPdN6-I, NiPdN6-A3, and PdNiN6-Z2, the differences between ΔG1 and ΔG3 for these catalysts are only 0.07, 0.02, and 0.03 eV (see Table S3), respectively. Thus, the catalysts around the right leg are approximately governed by the weak *OOH binding. In addition, NiPdN6-A1, PdNiN6-A2, PdNiN6-A5, and NiPdN6-A6 with the RDS of the hydrogenation of *O to *OH are located in the middle regime of the above two boundaries, corresponding to the adsorption strength to *O. The less-defined nature reflects the weak correlations between ΔG*O and ΔG*OHG*OOH. Therefore, although the ΔG*OH of these catalysts are close to the optimal value of ΔG*OH, their ORR overpotentials are not optimal. In contrast, NiPdN6-A2 with ΔG*OH of 0.82 V is closest to the apex of the volcano plot, thus possesses the best catalytic activity. However, there are still some distances to reach the apex of the volcano plot for this structure. Therefore, it can be inferred that there exists a possibility to further improve the ORR catalytic activity.

Fig. 4
figure 4

The scaling relationships for the Gibbs adsorption free energy between (a) ΔG*OH and ΔG*OOH, (b) ΔG*OH and ΔG*O, and (c) ΔG*O and ΔG*OOH. (d) The volcano plot of overpotential with ΔG*OH.

Available works have demonstrated a great potential to enhance the ORR catalytic activity of graphene-based metal atom catalysts by doping heterotopic O atom15,34,35,36,37. Inspired by these works, we try to replace one N atom in NiPdN6-A2 (with the lowest overpotential) with O atom for further enhancing its ORR catalytic activity. In NiPdN6-A2, each N atom can be replaced by O atom (see Fig. 5a), thus there exist six possible Ni-Pd embedded O and N co-coordinated graphene structures (refer to Fig. S3), and the corresponding Ef were calculated by Eq. (1) and illustrated in Fig. 5b. The Ef for these structures are -6.74, -5.77, 7.39, -6.12, -5,75, and -5.64 eV, respectively. In general, the lower value of Ef represents the better structural stability. Therefore, the N1 atom is most likely to be replaced by O atom to form Ni-Pd embedded O and N co-coordinated graphene structure (denoted as NiPdN5O-A2). In addition, it should be noted that 7.39 eV is significantly larger than other values because the NiPdN6-A2 structure undergoes severe deformation after replacing the N3 by O atom (see Fig. S3).

Fig. 5
figure 5

(a) Schematic diagram of the N atom in NiPdN6-A2 that can be replaced by O atom. The pink and green balls denote H and C atoms, respectively. (b) Formation energy for the O doped NiPdN6-A2 structure.

After introducing the O atom into NiPdN6-A2 to form NiPdN5O-A2 structure, the Eads for *O2 was decreased from -0.59 to − 0.61 eV, while the Eads for *OOH, *O, and *OH were increased from − 1.22, − 4.11, and − 2.57 eV to − 0.96, − 3.91, and − 2.39 eV, respectively. Therefore, the ORR catalytic activity of NiPdN5O-A2 may be different from that of NiPdN6-A2. Furthermore, the free energy change and overpotential of ORR on NiPdN5O-A2 were further explored. As displayed in Fig. 6a and Table S3, the downhill free energy for ORR on NiPdN5O-A2 at zero electrode potential (U = 0 V) demonstrates that the ORR process should be exothermic and thermodynamically feasible. With the increasing of electrode potential, NiPdN5O-A2 presents the working potential of 0.94 V during the exothermic reaction process. This value is larger than that (0.82 V) for NiPdN6-A2, indicating that ORR is more prone to proceed on the NiPdN5O-A2. Under the standard ORR potential (U = 1.23 V), the hydrogenation of *O to *OH and the hydrogenation of *OH to H2O became uphill, indicating these steps are thermodynamically unfavorable. In particular, the maximum free energy barrier of 0.31 eV needs to be overcome in the RDS of the hydrogenation of *OH to H2O. Hence, the overpotential of ORR on the NiPdN5O-A2 was only 0.31 V. Since the ORR usually takes place in an aqueous electrolyte solution, the solvent effect was evaluated using an implicit solvation model based on the VASPsol49, and the corresponding free energy diagram was presented in Fig. S4. By comparison with the free change in the gas phase, no remarkable difference in each reaction step was observed. The RDS still occurred in the hydrogenation of *OH to H2O, and the overpotential was slightly increased by 0.06 V, indicating the negligible solvent effect. In comparison to NiPdN6-A2 with larger overpotential (0.41 V), NiPdN5O-A2 should possess higher ORR catalytic activity. Moreover, the overpotential for the NiPdN5O-A2 is remarkably smaller than that for Pt (100) (0.47 V)47 and Pt (111) (0.44 V)48, indicating the ORR catalytic activity of the former superior to the latter.

Fig. 6
figure 6

(a) Free energy diagram of ORR on NiPdN5O-A2. (b) The PDOS and the d-band centers of Ni and Pd atoms in NiPdN6-A2 and NiPdN5O-A2.

To concretely explore the origin of the superior catalytic activity of NiPdN5O-A2, the projected density of states (PDOS) of d orbitals and the d-band centers of Ni and Pd atoms in NiPdN6-A2 and NiPdN5O-A2 were calculated, as illustrated in Fig. 6b. One can find that the distribution of d orbital and d-band center of Ni atom were closer to the Fermi-level than that of Pd atom in NiPdN6-A2, implying that Ni atom should be the main active site for ORR, in agreement with the adsorption structures of intermediates on NiPdN6-A2. When replacing N1 with O atom, the d orbital of Ni atom moved towards the Fermi-level, and the d-band center were increased from − 1.03 to − 0.88 eV for Ni atom while decreased from − 2.05 to − 2.23 eV for Pd atom. This led to a rising reaction activity of Ni atom and a lowering reaction activity of Pd atom, and the adsorption configuration for *OH was transformed from bonding with Ni and Pd atoms into only binding Ni atom. Eventually, NiPdN5O-A2 presented reduced adsorption capability for *OH (see Table S3) and enhanced ORR catalytic activity.

Conclusions

In summary, our DFT calculations demonstrated that armchair-edge termination was energetically more favorable than zigzag-edge termination. Edge-terminated NiPdN6-Gs, especially armchair edge-terminated NiPdN6-G, presented significantly enhanced ORR catalytic activity as compared with NiPdN6-G without edge termination. After doping O atom into armchair edge-terminated NiPdN6-G, the overpotential was further lowered to 0.31 V, indicating excellent ORR catalytic activity. This remarkable improvement in activity was attributed to the shift of the d-band center of Ni atom towards the Fermi-level while the shift of the d-band center of Pd atom away from the Fermi-level after O doping, leading to the transformation of adsorption configuration of *OH from binding with Ni and Pd atoms into only binding Ni atom and thus the decreased adsorption strength to *OH. This work sheds light on development of high-performance graphene-based DACs toward ORR by edge and doping engineering.