Importance of the gas-phase error correction for O 2 when using DFT to model the oxygen reduction and evolution reactions

DFT modelling of the oxygen reduction and evolution reactions (ORR and OER) habitually makes use of semiempirical corrections to oxygen in the gas phase. Although such corrections are tacit in the model, they should not be overlooked. In this article, we calculate the errors in the total energy of oxygen for commonly used exchange-correlation functionals, PW91, RPBE, PBE, and BEEF-vdW, to show that, for all functionals tested, the error is at least 0.3 eV. We discuss the impact this sizeable error in oxygen has on the modelling of the ORR and the OER. The error due to oxygen affects not only the overall equilibrium potential of the reaction, but also the energies of individual mechanistic steps. This illustrates that understanding the reasoning behind the semiempirical corrections for oxygen is important for researching new catalysts which may have different potential limiting steps.


Introduction
The oxygen reduction reaction, ORR, is one of the fundamental reactions in electrochemistry [1][2][3].The ORR takes place at the cathode of polymer electrolyte membrane fuel cells (PEMFCs) and so it is essential in the development of efficient new energy technologies.However, the reaction is sluggish, and is by far the largest contribution to voltage losses for these devices [4].Platinum is the most active catalyst for the reaction but since it is expensive and scarce, research is being carried out to try and decrease catalyst loading whilst increasing activity [3,5,6].
Based on e.g.cyclic voltammetry, X-ray photoelectron spectroscopy, Raman spectroscopy, and computational models [7][8][9][10], adsorbed hydroxyl (*OH) is a key ORR intermediate.For different catalytic materials, the overpotential of the reaction can be described as a function of the *OH free energy of adsorption, ΔG OH , giving rise to a Sabatier-type volcano plot of reactivity, where surfaces with intermediate binding energies have the lowest overpotential [3,10,11].Indeed, for Pt-based catalysts, the optimal binding of *OH is approximately 0.10-0.15eV weaker than for Pt(1 1 1) terrace sites [12,13].To modify the *OH adsorption energies of Pt-based catalysts, one can employ compressive strain and/or ligand effects introduced by alloying and dealloying treatments [14][15][16][17].
Assessing the active sites for the reaction on the surface is an important step in catalyst development.Since the original work from Feliu's group [18], various experimental studies have shown that introducing step defects to single-crystal surfaces increases the catalytic activity in acid [19][20][21].Conversely, other studies show that decreasing the size of nanoparticles, and therefore increasing the proportion of step edge sites, decreases the activity [22,23].In an attempt to understand these opposing results, a number of DFT studies have been carried out [13,16,24].It is observed that the ORR intermediates bind increasingly strong to Pt surfaces as the numbers of firstand second-nearest neighbours of the active sites decrease [13,16,24].Therefore, if ideal Pt sites for the ORR bind *OH more weakly than Pt(1 1 1), undercoordinated step sites are not likely to be responsible for the enhancement.In fact, Pt(1 1 1) sites with increased number of second-nearest neighbours are predicted to bind *OH in the range of interest (0.10-0.15 eV weaker than Pt(1 1 1)), which is found at the bottom of step edges and similar over-coordinated sites [24][25][26] present also in hollow nanoparticles, nanoframes and akin nanostructures [6,27].
For computational modelling to be an effective method of screening catalysts, the results must be in line with experiments and obtained within reasonable timeframes.Huge efforts have gone into increasing the accuracy of DFT by improving the functionals used.For instance, the exchange-correlation energy can be swiftly estimated by means of the generalised gradient approximation (GGA), which uses the local density and its gradient at a given point [28][29][30].They can be extended further to meta-GGAs which include an approximate of the second derivative [31,32], and long-range interactions can also be accounted for in a variety of ways [33][34][35].Hybrid functionals are computationally demanding because they include exact, non-local, Fock exchange, the proportion of which is often determined in an empirical way.Such exchange is combined with a GGA functional for the evaluation of correlation contributions to the total energy [36][37][38].In general, the more accurate the functional used, the more expensive the calculation and so a balance between accuracy and cost needs to be found.To illustrate this, for a group of 148 molecules, the mean absolute errors with respect to experiments for a GGA functional, a meta-GGA and a hybrid are, respectively, 0.76, 0.26 and 0.22 eV [32].Thus, for molecules, hybrid functionals are commonly used.Conversely, for metals, which have delocalised electron density, GGAs are often used, as they provide a good comparison to experiment for bulk and surface properties using reasonable computational resources [39][40][41].
When modelling heterogeneous reactions, where metal surfaces and gas-phase molecules are involved, there is a question of which type of functional to use.In general, GGAs are used in order to get an accurate description of the metal surface, however this can introduce unacceptable levels of error for the molecules involved [31,32].These errors can be mitigated by the addition of semiempirical corrections to the total energy of the gases, which is regularly done for instance, in the modelling of reactions belonging to the carbon cycle [42][43][44].
There is enormous interest in electrochemical reactions involving molecular oxygen, such as the ORR and the OER [2][3][4]10,45], but this molecule is particularly poorly described by GGAs and other types of exchange-correlation functionals [31] and so a semiempirical correction is obtained from the equilibrium redox potential of water [46,47].In this paper we use four common GGA functionals, Perdew-Wang-1991 (PW91) [28], Perdew-Burke-Ernzerhof (PBE) [29], revised Perdew-Burke-Ernzerhof (RPBE) [30], and Bayesian error estimation functional with van der Waals correlation (BEEF-vdW) [33].According to Web of Science, to date, the most commonly used functionals to model the ORR or the OER are PBE and RPBE [48].We use the functionals to evaluate the error in the calculated energy of O 2 and discuss how it affects the ORR and OER predictions.

Method
The DFT calculations were carried out using the VASP code [49].For each of the functionals (PW91 [28], PBE [29], RPBE [30], and BEEF-vdW [33]) each gas phase molecule was relaxed using the conjugate gradient method until the maximum force on any atom was below 0.01 eV Å −1 .The effect of the core electrons on the valence electron density was described using the Projector Augmented Wave (PAW) method [50] and a finite temperature Gaussian smearing was used to facilitate convergence of the self-consistent field process.The Fermi level was smeared with a width of 0.001 eV, and all energies were extrapolated to 0 K.The valence electron density was described with a plane-wave basis set with a cut-off of 450 eV.A convergence test for the free energy of O 2ðgÞ þ 2H 2ðgÞ ⇌ 2H 2 O ðlÞ with cutoffs between 300 and 1000 eV with PBE is included in the Supporting Information to show that this value provides converged reaction energies.Spin unrestricted calculations were performed for the oxygen molecule and the final magnetisation for all functionals was 2 μ B , as expected for a triplet state.The total Gibbs energies of gas-phase species were calculated using ΔG°¼ ΔE DFT þ ΔZPE À TΔS°.Zero-point energies (E ZPE ) were determined using the harmonic oscillator approximation.For molecular hydrogen and oxygen, entropy contributions at 298.15 K (S°) were taken from available thermodynamic data [51].A suitable liquid-phase correction was applied to the entropic correc-tion for water [10,46,47,51].Heat capacity effects were not included, in view of their small contributions to formation energies in the range from 0 to 298.15 K [52].

Oxygen reduction to water via the associative pathway
At the cathode of PEMFCs, O 2 is reduced by four proton-electron transfer steps to water The standard free energy of the overall reaction, ΔG 0 ORR , is −4.92 eV.This can be converted to a standard equilibrium potential for the reaction, E°= 1.23 V vs RHE.It is generally assumed that all catalysts follow the same 4-step associative mechanism [3,10] described in Eqs. ( 1)-( 4), since studies have shown that there is a high kinetic barrier towards direct O 2 dissociation on Pt [10], see section 3.2.However, it has been proposed that the dissociative pathway might hold the key to enhance O 2 electrocatalysis [53].
Hereon, * represents a free surface site and *X are adsorbed species (X = O, OH, OOH).The calculated free energy of each step in the catalytic pathway can be written as in Eqs. ( 5)-( 8), where ΔG y H2O is zero, as water and proton-electron pairs are used as free-molecule references for the adsorption energies and the formation of O 2 .
Every Gibbs free energy calculated using DFT, ΔG y i , has some positive or negative error, ε i , associated with it.The error is calculated with respect to the ideal (experimental) energy, ΔG i , by means of Eq. ( 9) [43]: The total reaction energy can be determined from the overall reaction and is also the sum of each individual step, Eq (10).
In principle, DFT-GGAs model adsorbed states accurately (or at least with an error comparable to that of clean surfaces), and H 2 O and H 2 are also well described [31].Hence, it is generally assumed that ɛ i ¼ 0 for H 2 and H 2 O and that ɛ ÃOH ≈ɛ ÃOOH ≈ɛ ÃO ≈0.Therefore, using Eqs.( 9) and (10) it can be shown that errors in the total reaction energy are due only to oxygen (Eq.( 11)).
We calculated ΔG y ORR for four commonly used functionals and, using ΔG ORR = −4.92eV (that is, the experimental value), determined ɛ O2 , the results of which are shown in Table 1.Values of 1.08 and 1.01 V for the ORR equilibrium potential have previously been recorded for PBE and RPBE, respectively [54].The systematic shift of 0.03 V likely stems from the fact that the DFT calculations were carried out using different codes and basis sets.For all functionals, ɛ O2 is significantly larger than 0.1 eV, which is generally accepted as the minimum error in DFT-GGAs.In electrochemical experiments, when the pH is varied by one unit , the potential is shifted by 0.059 V.When ΔG y ORR is converted to E°vs RHE, and the error in eV is converted to V by dividing by the number of transferred electrons (ɛ O2 to ɛ 0 O2 ), all values are greater than 0.059 V, illustrating that the error is significant not only in terms of energies, but also in terms of potentials.For all the functionals tested, the error is conspicuous, but the largest error is for BEEF-vdW and it has been suggested that ɛ H2 = 0.1 eV should also be used [55].Since the exchange-correlation energies calculated with PW91, PBE, RPBE and BEEF-vdW are based on the same general equation using different coefficients for the various terms in the equation (see further details in [33]), the differences among these functionals stem from the choice of those coefficients.BEEF-vdW also includes non-local interactions, which might also affect its predictions.In Table 1 we observe that the difference between the experimental result and the result predicted by DFT without the O 2 correction is so large that the common practice of O 2 correction in the gas phase is essential.Indeed, 1.23 V vs RHE and −4.92 eV are ubiquitously written instead of the actual DFT-calculated free energies and equilibrium potentials.
In nature, oxygen reduction catalysed by enzymes happens at very low overpotentials (~0.1 V) and each of the steps in the enzymatic pathway has a very similar energy [56].An ideal electrocatalyst, as shown with the black line in Fig. 1, mimics enzymes and every step has an equal energy, ΔG 0 ORR =4 ≡ −1.23 eV.The same criterion holds for an ideal OER electrocatalyst, so that every step takes 1.23 eV [57][58][59].For real electrocatalysts, the overpotential required for each step is calculated from the adsorption energies of the intermediates (*O, *OH and *OOH) and the equilibrium potential and it is essential to consider adsorbate-solvent interactions, especially hydrogen bonding [54,[60][61][62].In these systems the assumption that ɛ ÃOH ≈ɛ ÃOOH ≈ɛ ÃO ≈0 may not always be applicable.The method of calculating the equilibrium potential for each of the reaction steps affects the error in *O and *OH [54] and there seems to be a systematic error in how the OAO bond in *OOH is described [63].Using BEEF-vdW, some approaches allow for estimation of these errors [64].However, the ideal catalyst is affected only by ɛ O2 , since the energetics of its isometric electrochemical steps are obtained from the quotient of the equilibrium potential and the number of electrons transferred.
Fig. 1 illustrates that at each step of the mechanism, ɛ O2 causes the energy released to be underestimated.Furthermore, ɛ O2 does not just affect the overall energy of the reaction, the modelling of any individual mechanistic step involving O 2 will also be affected.For ORR this affects the first step in the reaction, Eq. ( 1).From a thermodynamic perspective, the potential-limiting step of any reaction is the one of greatest interest.Accurate modelling of this step influences the design of catalysts with controlled shapes and design of surface atom arrangements to maximise the number of active sites for the reaction [25,65].For a large number of electrocatalytic materials active for the ORR, including porphyrins, perovskites, monoxides, doped TiO 2 and functionalised graphitic materials, the potential limiting step was calculated [59].The percentage of these catalysts limited by each step is reported in Table 2.The results show that for 61% of the catalysts analysed, the first step, Eq. ( 1), is the potential limiting step.If a greater number of metal surfaces [10,13,16,24,66], which tend to bind *OH more strongly than oxidized materials, were included in the analysis, the percentage limited by the final step, Eq. ( 4), would likely increase until it is approximately 50%.However, the number of catalysts limited by the first step, Eq. ( 1), would still be significant and so it is particularly important to account for ɛ O2 for modelling the ORR since O 2 is involved in the potential-limiting step of the reaction for a majority of current catalysts, see section S3 for the effect of ɛ O2 on the ORR overpotentials for real catalysts.

Oxygen reduction to water via the dissociative pathway
If the O 2 molecule dissociates directly, instead of forming *OOH, (i.e. 2 Ã þ O 2 ⇌ 2 Ã O), it is said that the ORR proceeds through the Table 1 For each functional, the calculated reaction energy, ΔG † ORR , and the equilibrium potential, E o vs RHE, are reported.The errors in these values due to oxygen (where ɛ ORR ¼ Àɛ O2 ) are also presented.2, where the same errors observed in Fig. 1 are present for the different exchange-correlation functionals analysed.Since the ideal catalyst is obtained from thermodynamic considerations, it is straightforward to conclude that the free energy for the dissociative adsorption of O 2 should be ΔG diss = 0. Furthermore, in terms of Gibbs energy, the Brønsted-Evans-Polanyi (BEP) relation reported by Nørskov et al. [67] for O 2 dissociation on stepped transition metals is the following (in eV): Table 2 The percentage of catalysts [59] analysed which are limited by each step in the associative ORR and OER reaction pathways.The ORR proceeds from step (1) to (4), whereas the OER proceeds from step (4) to (1).

% of catalysts limited by step
Pathway step for ORR for OER Fig. 2. Free diagram for the 4-electron ORR reaction on an ideal catalyst via the dissociative pathway.The free energy for an ideal experimental electrocatalyst (black line) is shown.The calculated free energy profiles for ideal catalysts calculated with PW91 (blue), PBE (green), RPBE (lilac) and BEEF-vdW (red) are also shown.
Fig. 3. Free energy diagram for the 4-electron OER reaction on an ideal catalyst.The free energy for an ideal experimental electrocatalyst (black line) is shown.
The calculated free energy profiles for ideal catalysts calculated with PW91 (blue), PBE (green), RPBE (lilac) and BEEF-vdW (red) are also shown.
If the ideal catalyst obeys Eq. ( 12), its kinetic barrier would be G TS ≈ 1.11 ± 0.26 eV, which is large for a room-temperature process.In fact, the limit for surmountable barriers at room temperature has been set at 0.75 eV, corresponding to a turnover frequency of 1 s −1 per site [68].This justifies the fact that most ORR analyses focus on the associative pathway.However, we note that the ideal thermodynamic catalyst does not necessarily need to follow Eq. ( 12) and might follow another with a smaller offset or may follow no BEP relations at all.

Oxygen evolution reaction
The arguments in the previous section can be extended to other reactions involving O 2 , especially the oxygen evolution reaction, OER, which is the opposite overall reaction to the ORR.Fig. 3 shows that, for the OER mechanism (which is assumed to be opposite of the associative ORR pathway [2,69]), ɛ O2 affects the overall energy needed for the reaction and the energetics of ideal catalysts in the same way.However, Table 2 shows that 87% of catalysts are limited by the second and third steps, the reverse of Eqs. ( 2) and (3).Neither of these steps involve O 2 and so for most catalysts, modelling of the OER potential limiting step is unaffected by ɛ O2 .Nevertheless, electrocatalysis of the OER is an area of intense research with new catalysts being developed regularly [45,65].Since the ultimate aim is to equalise the energy needed for each step, modelling of all four steps is needed and so ɛ O2 needs to be considered in the search for highly active materials.Finally, we note that the simple considerations made here for the total energy of O 2 , the ORR/OER reaction energies and redox potentials, and ideal ORR/OER catalysts, can be coupled with more advanced methods for the quantification of uncertainty in activity and stability plots [64,70,71].

Conclusion
It is important to understand that tacit semiempirical corrections are nearly always used for molecular oxygen in DFT-based electrocatalysis models.This paper shows that for common functionals, PW91, PBE, RPBE and BEEF-vdW, using the computationally calculated value for O 2 would introduce significant errors (greater than 0.3 eV) into the overall energy for the ORR and OER.The equilibrium potentials are also affected by these errors, and the DFT-calculated values are far from the experimental value of 1.23 V vs RHE.Until an exchange-correlation functional is developed which can describe catalyst surfaces and molecules at the same level of accuracy, a semiempirical correction for oxygen cannot be discarded.Understanding the reasoning for this correction is important in developing new ORR and OER catalysts.In general, matching the experimental equilibrium potential should be a necessary, initial step in computational electrocatalysis studies.

Declaration of Competing Interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.
pathway.After this initial chemical step, *O + *OH, *O + H 2 O, *OH, and * + H 2 O are formed sequentially.The free energy diagram for the ideal catalyst following the dissociative pathway is shown in Fig.
Fig. 1.Free energy diagram for the 4-electron ORR reaction on an ideal catalyst via the associative mechanism.The free energy for an ideal experimental electrocatalyst (black line) is shown.The calculated free energy profiles for ideal catalysts calculated with PW91 (blue), PBE (green), RPBE (lilac) and BEEF-vdW (red) are also shown.dissociative