Comparative density functional theory study for predicting oxygen reduction activity of single-atom catalyst

All the calculations were performed using the projector augmented wave (PAW)[31] method as implemented in the Quantum-ESPRESSO code [32]. The PAW potentials were adopted from the PSLIBRARY [33] version 1.0.0. Wave functions and augmentation charge density were expanded in terms of a plane-wave basis set with the kinetic energy cutoffs of 80 and 800 Ry, respectively. The Marzari-Vanderbilt cold smearing [34] with a smearing width of 0.02 Ry was used to treat the Fermi level. In this work, we compare the performance of different XC functionals within GGA, specifically PBE [35] and revised PBE (RPBE) [36] functionals. Grimme’s semiempirical dispersion correction (DFT-D3) [37] is also include for both functionals namely, PBE+D3 and RPBE+D3 to further investigate the importance of the dispersion force in the systems considered in this work. Finally, we include the Bayesian error estimations functional with nonlocal correlation (BEEF-vdW) [38] to further assess the accuracy of the GGA functionals. The BEEF-vdW is designed specifically to address vdW forces reasonably well while maintaining an accurate description of chemisorption energies of molecule on surface. Moreover, the BEEF-vdW itself has an error estimate built in, to not only yield an accurate energy difference in condensed matter studies, but also estimate the errors on computed quantities. Therefore, BEEF-vdW serves as a benchmark and a baseline for comparison in this study. To employ the error estimation capabilities in the BEEF-vdW, a probability distribution for model parameters for the exchange and correlation is randomly sampled through an ensemble of density functionals generated around the BEEF-vdW [38, 39]. In this work, we generated an ensemble of 2000 functionals for each adsorption system, isolated surface and adsorbate considered, to obtain a distribution of the adsorption energies ($E_{\mathrm{ads}}$). Then, the standard deviation of the ensemble for $E_{\mathrm{ads}}$ around the BEEF-vdW value was used to estimate the error for each adsorption system. In addition to $E_{\mathrm{ads}}$, the standard deviation was also used to estimate the uncertainty in the limiting potential and overpotential. We note that we did not scale the error in such a way that it reproduces that for the benchmark adsorption energies, unlike Ref. [39]. Thus, our estimated error for both $E_{\mathrm{ads}}$ and limiting potential/overpotential can be overestimated.

Throughout this study, we used PAW potentials generated using the PBE functional. Detailed results of the optimized lattice constant and the adsorption energies of the ORR intermediates obtained using RPBE and RPBE+D3 functionals with PAW potentials generated using the RPBE functionals are included in Supplemental Materials [40] for further reference. We note that it was suggested [24] that hybrid functions such as PBE0 [41, 42] and HSE06 [43] give more accurate results for ORR on nitrogen-doped graphene than GGA or SCAN meta-GGA[44], as compared with those obtained using the highly accurate coupled cluster theory. Liu et al. [45] used various XC functionals including HSE06 and SCAN and reported that the absolute value of the magnetic moment, which is suggested to be a good descriptor for ORR on Fe-N-C catalysts, varies depending on the functional, but the trend is the same (the authors do not report the reaction free energies with different functionals). In this study, we limit ourselves to the GGA level of theory, in order to systematically and comprehensively study the roles of the XC functional and dispersion correction on not only the energetics, but also the vibrational contributions to the reaction free energies.

We considered Fe-N₄ and Co-N₄ moieties embedded in graphene. The lattice constant for graphene was optimized in the primitive hexagonal unit cell using a 16$\times$16 $\mathbf{k}$-point grid for each functional. The DFT-optimized lattice constants are given in the Table  1. For the calculation of the ORR intermediates, the Fe-$\mathrm{N_{4}}$-C and Co-$\mathrm{N_{4}}$-C models was constructed using a (5$\times$5) surface unit cell of graphene, in which we removed two carbon atoms to anchor a single Fe/Co atom and replaced four-coordinated carbon atoms surrounding it with nitrogen atom [Fig. 1(a)]. To eliminate the spurious electrostatic interaction with the periodic images, the effective screening medium method [47, 48] was employed, along with the vacuum spacing of 20 Å. The geometries of all the reaction intermediate ($\mathrm{{}^{\ast}O_{2}}$, $\mathrm{{}^{\ast}OOH}$, $\mathrm{{}^{\ast}O}$, and $\mathrm{{}^{\ast}OH}$) on both Fe-$\mathrm{N_{4}}$-C and Co-$\mathrm{N_{4}}$-C moeties were fully optimized until the residual forces on the constituent atoms became smaller than $10^{-4}$ Ry/Bohr ($2.57\times 10^{-3}$ eV/Å). A 4$\times$4 $\mathbf{k}$-point grid was found to give reasonably accurate results (see Supplemental Materials [40]) and 8$\times$8 $\mathbf{k}$-point grid was used for the electronic structure analyses.

The adsorption energy of the intermediate is defined by

where $E_{\mathrm{S/A}}$, $E_{\mathrm{S}}$, and $E_{\mathrm{A}}$, are the total energies of combined system, isolated surface structure, and isolated adsorbate, respectively. All energies above were obtained under the same parameter settings.

The change in adsorption energy with an applied electrode potential was calculated based on the CHE model proposed by Nørskov and co-workers [18]. In this study, we focus on the associative mechanism of the ORR process as it has been reported that it is energetically more favorable and that the activation barrier for dissociative O₂ adsorption is considerably high (1.16 eV)[23] in the cases of the Fe-N₄-C and Co-N₄-C active sites. The elementary steps for ORR along the four-electron associative mechanisms can be written as follows: where ^∗ denotes the adsorption site. The DFT reaction energies of the ${}^{\ast}\mathrm{OOH}$, ${}^{\ast}\mathrm{O}$ and ${}^{\ast}\mathrm{OH}$ are then defined relative to $\mathrm{H_{2}{O}}$ (liquid phase) and $\mathrm{H_{2}}$ (gas phase), to avoid the calculation of the $\mathrm{O_{2}}$ molecule or the radical of $\mathrm{OOH}$ and $\mathrm{OH}$, which is notoriously difficult to describe within the semilocal approximation to DFT [18] as: where $E_{\ast}$ is the total energy of the surface without adsorbates, $E_{\ast{\mathrm{OOH}}}$, $E_{\ast{\mathrm{O}}}$ and $E_{\ast{\mathrm{OH}}}$ are total energies of the surface with adsorbates (${}^{*}\mathrm{OOH}$, ${}^{*}\mathrm{O}$ and ${}^{*}\mathrm{OH}$, respectively) bound to the surface, $E_{\mathrm{H_{2}O(g)}}$ and $E_{\mathrm{H_{2}(g)}}$ are the total energies of $\mathrm{H_{2}O}$ and $\mathrm{H_{2}}$ molecules in the gas phase, respectively. The energy of H₂O in the liquid phase ($E_{\mathrm{H2O(l)}}$) is calculated at 0.035 bar (i.e. $E_{\mathrm{H_{2}O(l)}}$ = $E_{\mathrm{H_{2}O(g)}}$ + $k_{\mathrm{B}}T\ln 0.035$), which correspond to equilibrium pressure of $\mathrm{H_{2}O}$ at 298.15 K, and therefore this state corresponds to that of liquid water [49].

The DFT formation energies for the ORR intermediates are converted to the adsorption free energy by including corrections for the change in zero-point energy $\Delta\mathrm{ZPE}$ and entropy $T\Delta S$ at 298.15 K estimated using finite displacement method to calculate the vibrational frequencies of the adsorbates through the use of the atomic simulation environment (ASE) [50] and $\Delta G_{\mathrm{pH}}$ are the pH values (pH was defined as 0 for acidic medium) as

This allows us to define the adsorption free energies of ${}^{*}\mathrm{OOH}$, ${}^{*}\mathrm{O}$ and ${}^{*}\mathrm{OH}$ (i.e. along associative mechanism) relative to $\mathrm{H_{2}O}$ and $\mathrm{H_{2}}$ In the case of the adsorption free energy of O₂ (^∗ + $\mathrm{O_{2}}$ $\rightarrow$ ${}^{\ast}{\mathrm{O_{2}}}$), the adsorption energy is defined relative to an $\mathrm{O}_{2}$ molecule in the gas phase ($\mathrm{O_{2}}$(g)).

By setting the reversible hydrogen electrode (RHE) as the reference electrode, the effect of the applied potential ($U$ $\mathrm{V_{RHE}}$) on the Gibbs free energy can be approximated by

where $n$ is the number of electrons transferred in each consecutive step and $e$ is the elementary charge of an electron. The Gibbs free energy of the four elementary reactions steps in the associative mechanism (at $U$ = 0 $\mathrm{V_{RHE}}$) can be written as

where $\Delta{G_{\mathrm{ORR}}}$ is the overall formation free energy value of ORR through the relation (2$\mathrm{H_{2}}$ + $\mathrm{O_{2}}$ $\rightarrow$ 2$\mathrm{H_{2}{O}}$), defined experimentally ($\Delta{G_{\mathrm{ORR}}}$ = $-$4.92 eV) at 298.15 $\mathrm{K}$ with pressure of $\mathrm{O_{2}}$ and $\mathrm{H_{2}}$ of 1 bar [51]. With this value, the thermodynamic equilibrium potential will be 1.23 $\mathrm{V_{RHE}}$.

Based on the calculations described above, the free energy diagram along the reaction path considered can be constructed. The limiting-potential (${U_{\mathrm{L}}}$) is estimated from Eq. (8) by determining the minimum free energy change ($\Delta{G_{\mathrm{min}}}$) among all the electrochemical steps.

In the CHE model, the ${U_{\mathrm{L}}}$ is also known as potential at which all the reaction steps are exothermic and by this definition, it can be used to compare with the experimental half-wave potential ($\mathrm{E_{\frac{1}{2}}}$) value [52]. Accordingly, the overpotential ($\eta_{ORR}$) can be further derived from the (${U_{\mathrm{L}}}$) as Eq. (9), in which a lower overpotential implies a better ORR activity.

We first optimized the ORR intermediates ($\mathrm{{}^{\ast}O_{2}}$, $\mathrm{{}^{\ast}OOH}$, $\mathrm{{}^{\ast}O}$, and $\mathrm{{}^{\ast}OH}$) on Fe-$\mathrm{N_{4}}$-C and Co-$\mathrm{N_{4}}$-C active sites, and calculated their adsorption energies with PBE, PBE+D3, RPBE, RPBE+D3 and BEEF-vdW functionals (Fig. 2). The calculated energies are summarizes in Table  2.

The present PBE results on both active sites are in a good agreement with the previous works by Kattel et al. [17, 16]. In general, the adsorption energies obtained using PBE are larger in magnitude, and by the inclusion of dispersion correction, PBE+D3 gives adsorption energies larger in magnitude than those obtained using other functionals, and the energy differences are more significant in the case of Fe-$\mathrm{N_{4}}$-C active site. In contrast, RPBE predicts more repulsive interactions and resulted in more positive interaction energies on both active sites. However, RPBE+D3 improves the adsorption energies, which are larger in magnitude than those by RPBE, and the values are in good agreement with those obtained by using BEEF-vdW. The results indicate that the contribution of dispersion interaction is significant within these two systems, and the inclusion of the dispersion correction in RPBE+D3 leads to the adsorption energies of all the ORR intermediates with similar accuracy to the BEEF-vdW. We note that the estimated error in the BEEF-vdW adsorption energies is relatively large (0.2 - 0.3 eV for Fe-$\mathrm{N_{4}}$-C and 0.1 - 0.2 eV for Co-$\mathrm{N_{4}}$-C), but agrees well with the previously reported one for N-doped graphene [24]. We also note that PBE tends to overestimate the adsorption energy of chemisorption species, partly because the exchange in the PBE functional tends to lead too attractive interaction in molecular systems [36]. In addition, PBE can show spurious attractive interaction for weakly interacting systems from the exchange only. The detail discussion can be found in Ref. 53, 54, 30.

We next evaluated the catalytic activities of the Fe-$\mathrm{N_{4}}$-C and Co-$\mathrm{N_{4}}$-C for the ORR with different functionals. The free energy diagrams for ORR along the four-electron associative mechanism on the Fe-$\mathrm{N_{4}}$-C and the Co-$\mathrm{N_{4}}$-C catalysts are calculated and displayed in Fig. 3. The calculated adsorption free energies and details theoretical limiting potentials and overpotentials on the Fe-$\mathrm{N_{4}}$-C and the Co-$\mathrm{N_{4}}$-C catalysts are summarized in Tables  3 and 4, respectively.

Our calculated ORR free energy diagrams for the associative mechanism (Fig. 3) show that RPBE and RPBE+D3 functionals predict a similar potential determining step (PDS) with BEEF-vdW, which is the O protonation to form OH, while both PBE and PBE+D3 predict the formation of H₂O (last reaction steps) on the Fe-$\mathrm{N_{4}}$-C. On the other hand, all the functionals predict the OOH formation (initial reaction steps) as the PDS on the Co-$\mathrm{N_{4}}$-C catalyst. The difference in the estimated PDSs in the Fe-$\mathrm{N_{4}}$-C catalyst originates from the significant difference in the ORR intermediate adsorption energies predicted by the PBE and PBE+D3 functionals. The ${U_{\mathrm{L}}}$ values for Fe-$\mathrm{N_{4}}$-C predicted by PBE+D3 and RPBE are in good agreement with that by BEEF-vdW with a difference of $\sim$ 0.02 V, while PBE and RPBE+D3 is slightly overestimate it by $\sim$ 0.08 V. It should be stressed that all the functionals predicts ${U_{\mathrm{L}}}$’s for Fe-$\mathrm{N_{4}}$-C, which seemingly agree well with that from BEEF-vdW. However, PBE and PBE+D3 values are obtained for a different PDS, i.e., different mechanisms. On the other hand, in the case of Co-$\mathrm{N_{4}}$-C, only RPBE+D3 predicts ${U_{\mathrm{L}}}$, which is comparable to that from BEEF-vdW, while PBE, PBE+D3, and RPBE underestimate it by $\sim$ 0.25, $\sim$ 0.10, and $\sim$ 0.26 V, respectively.

We also estimated the $\eta_{\mathrm{ORR}}$ by using Eq. 9 (Table 4). The error bars for BEEF-vdW indicate a relatively large uncertainty (0.4 V) for the $\eta_{\mathrm{ORR}}$ on both active sites. Nevertheless, only RPBE+D3 consistently estimated a similar PDS and comparable $\eta_{\mathrm{ORR}}$ to BEEF-vdW. The inconsistent results on both catalysts obtained using PBE and PBE+D3 may be due to the too attractive exchange in PBE, as discussed above. With respect to the $\eta_{\mathrm{ORR}}$ values (Table 4), PBE+D3, RPBE+D3, and BEEF-vdW predict that the Co-$\mathrm{N_{4}}$-C catalyst is more reactive (lower $\eta_{\mathrm{ORR}}$ value), while PBE and RPBE estimate the catalytic performance of the Fe-$\mathrm{N_{4}}$-C and the Co-$\mathrm{N_{4}}$-C are comparable. It is noted that our results in Table 4 are limited to the vacuum condition. It is also noted that we estimated the adsorption free energies using RPBE+D3 with the on-site Coulomb interaction ($U$) and report in Supplemental Materials [40], which are consistent with those reported in Ref. [23]. However, calculated adsorption energies and potential determining steps can depend on the choice of the $U$ value, which calls for further investigation.

To gain insight into their catalytic activities, we calculated the densities of states for clean and $\mathrm{{{}^{\ast}}O_{2}}$ adsorbed Fe-$\mathrm{N_{4}}$-C and Co-$\mathrm{N_{4}}$-C sites using RPBE+D3 functional. We found that a distinct spin polarization at the Fe-$\mathrm{N_{4}}$-C site compared to the Co-$\mathrm{N_{4}}$-C sites, which can be also confirmed by a larger local magnetic moment (1.61 $\mu_{\mathrm{B}}$ and 0.68 $\mu_{\mathrm{B}}$ for the former and the latter, respectively). These results further explain stronger binding strength for all ORR intermediates that we obtained for the Fe-$\mathrm{N_{4}}$-C as compared to Co-$\mathrm{N_{4}}$-C site (Table 2). Besides, the spin down channels of $3d_{z^{2}}$ and $3d_{zy}$ in the Fe-$\mathrm{N_{4}}$-C are unoccupied and located above the Fermi level, while the Co-$\mathrm{N_{4}}$-C has the half-filled $3d_{z^{2}}$ spin down channel cross the Fermi level. This indicates that the Co-N₄-C active site has a half-metallic nature, and has less opportunities to form a bond with O. Therefore, the Fe-$\mathrm{N_{4}}$-C site can form a stronger Fe$-$O $\sigma$ bond compare to Co$-$O [55, 56, 57]. This behavior can be related to the Sabatier principle in which the strong bonding between active site and reactant effectively inhibits further reaction and results in decreasing catalytic activity [58, 59]. As expected, the moderate spin polarization and half-metallic properties shown by the Co-$\mathrm{N_{4}}$-C site ensures the moderate binding strength of the main ORR intermediates, and is responsible for the higher ORR catalytic activity [21].

In a recent theoretical work by Wang et al. [23], it was proposed that the intrinsic intermediate $\mathrm{{{}^{\ast}}OH}$ is co-adsorbed on the single-metal-atom active site, and is the origin of the improved activity of the ORR through the associative mechanism on the TM-$\mathrm{N_{4}}$-C catalyst. The presence of the OH species was also confirmed by the experiment [60]. Following the aforementioned studies, we further investigate the effect of OH-termination on the Fe-$\mathrm{N_{4}}$-C and Co-$\mathrm{N_{4}}$-C active sites without the solvent effect (Fig. 5). Here we employed RPBE+D3 functional, which predicts consistent results with BEEF-vdW and compares with the results obtained using PBE+D3. We also include PBE in this comparison, because it was used in Ref. 23. The calculated energies are summarizes in Table 5

Here, with the OH-termination of both Fe-N₄-C and Co-N₄-C active sites, PBE+ D3 and RPBE+D3 predict the formation of $\mathrm{H_{2}O}$ as the PDS, while PBE, the protonation of $\mathrm{O_{2}}$ to form $\mathrm{OOH}$ on the Fe-N₄-C. In the case of Co-N₄-C, all the functionals predicts the $\mathrm{O_{2}}$ protonation to form $\mathrm{OOH}$ as the PDS. We further estimated the $\eta_{\mathrm{ORR}}$ for the detailed comparison (Table 6). The $\eta_{\mathrm{ORR}}$ values for Fe-N₄-C predicted by PBE, PBE+D3 and RPBE+D3 show a noticeable improvement of the ORR catalytic activities (lowering of the overpotential) after the OH co-adsorbs on the Fe site by $\sim$ 0.14, $\sim$ 0.25, and $\sim$ 0.17 V, respectively. On the other hand, on the Co site, no significant change can be found by using each functionals.

The slight improvement of the ORR activity in the Fe-N₄-C is mainly due to the formation of the $\sigma$-bonding between $2p_{z}$ and $3d_{z^{2}}$ orbitals upon termination of the $\mathrm{{{}^{\ast}}OH}$ on the Fe site (see Supplemental Materials [40]). Accordingly, the proportion of $\sigma$-bonding state on the Fe(OH)-$\mathrm{N_{4}}$-C site decreased, and the interaction between the ORR intermediates is dominated by the $\pi$-bonding states (i.e., the filling degree of $3d_{zx}$ and $3d_{zy}$). Therefore, the adsorption of ORR intermediates on the Fe(OH)-$\mathrm{N_{4}}$-C becomes weaker (see Supplemental Materials [40]), and the ORR performance is improved. The same behavior is observed for the Co-$\mathrm{N_{4}}$-C with the OH-termination. However, since all the orbitals of $3d_{z^{2}}$, $3d_{zx}$, and $3d_{zy}$ are almost filled after the $\mathrm{{{}^{\ast}}OH}$ co-adsorption at the Co site, and also due to the half-metallic nature of the Co-$\mathrm{N_{4}}$-C, no significant change is observed in the adsorption strength of the ORR intermediates as well as the ORR catalytic activity.