A DFT study of the oxygen reduction reaction mechanism on be doped graphene

Graphene despite its high surface area has very limited activity towards the oxygen reduction reaction (ORR), demonstrating selectivity towards the unfavorable two-electron mechanism. We have employed the spin polarized density functional theory method to investigate the oxygen reduction reaction activity of the heteroatom p-type beryllium (Be) doped graphene in gas and aqueous media. The preferred doping sites, active sites and reaction mechanisms available on the doped graphene surfaces were investigated with increasing Be concentrations of 0.03 ML, 0.06 ML and 0.09 ML. Our results reveal that oxygen is physisorbed on bare graphene, however, Be at the lattice sites provides site for oxygen adsorption and ORR. Generally, Be concentration increase in a single honey-comb of graphene from 1 to 2, increasing oxygen activation and dissociation on graphene, compared to their isolated defective counterparts. Considering the conjugated Be defective surfaces, oxygen binds dissociatively on the doped surfaces preferentially in the order 0.06 ML > 0.09 ML > 0.03 ML. Whereby strong binding at the 0.06 ML doped surface corresponds to least charge loss to oxygen. Surfaces with the least affinity for O*, OH* and H2O*, and can transfer the most charges to O2 leading to dissociation (as seen on the 0.03 ML surface), shows the highest reaction activity towards the ORR. From these results surfaces that bind O2 weakly and dissociatively could hold promise for the ORR reaction. Single atom catalyst of Be, based on graphene is promising for the ORR.


Introduction
The global interest in electrochemical energy technologies, such as fuel cell and battery technologies, is rising. This is towards attaining energy sustainability and to reduce the global carbon economy. Traditional solid-oxide fuel cells (SOFCs) have the highest performance and durability but operate at temperatures above 600 °C. Lowering the operating temperature is critical. However, the energy-conversion efficiency of the low-temperature SOFC is mainly limited by the slow kinetics of the oxygen reduction reaction (ORR) at the cathode (Borup et al. 2007;Marques 2013).
Pt and its alloys are effective cathodes (Lim et al. 2009), however, due to the high cost, low abundance and poor durability of the Pt-based catalysts, alternative low-cost efficient metal and carbon-based cathodes have been of interest worldwide in recent years. For instance, the interest in carbon sheets like graphene is due to their low cost, environmental friendliness, outstanding activity and stability (Ma et al. 2016).
Be has an excellent ability to fine-tune the physicochemical properties and catalytic activity of graphene. Lattice Be heteroatom-doped graphene has received attention for several applications including band gap tuning for solar cells (López-Urías et al. 2015;Olaniyan et al. 2019), magnetic induction (Song et al. 2017), as a single metal catalyst support (Groves et al. 2012), an anodic material for lithium ion batteries (Ullah et al. 2017;Fatheema and Rizwan 2020), mainly owing to Be belonging to the family of light atoms. However, other heteroatom-doped graphene cathodic fuel cell materials have attracted considerable interest in very recent years. For example, materials like p-type doped graphene (e.g., boron), n-typed doped graphene (e.g., nitrogen and phosphorus) and halogen doped graphene have received significant attention for the ORR. This includes their ternary and quaternary doped compounds (Wu et al. 2010;Ozaki et al. 2010;Qu et al. 2010;Choi et al. 2012;Yang et al. 2012a, b;Yang et al. 2012a, b;Choi et al. 2013;Kaukonen et al. 2013;Li et al. 2014;Rao et al. 2014;Kakaei and Balavandi 2016) which have shown high catalytic activity and selectivity towards the ORR.
The ORR is reported to occur along two possible pathways; via the 2-electron pathway which could lead to hydrogen peroxide formation or via 4-electron pathway which will directly transform into water. Understanding the detailed mechanisms have focused on other heteroatom doped graphene and not beryllium. Using both experiment and density functional theory calculations, the detailed mechanisms on doped graphene has been studied, i.e., nitrogen (Ni et al. 2012), phosphorus (Zhang et al. 2015;Bai et al. 2016), boron (Fazio et al. 2014), boron and nitrogen (Zhao et al. 2013), nitrogen, sulfur and phosphorus (Choi et al. 2013), MnN 4 (Lu et al. 2015) and FeN 4 (Kattel and Wang 2014).
To the best of our knowledge, although p-type heteroatom doped graphene are efficient for the ORR in fuel cells, work on p-type dopants have focused only on boron (B). However, no work systematic studies have been conducted on the effect of beryllium dopant on graphene for the ORR. Thus, we have conducted for the first time, a comprehensive study of the Be-doped graphene, to understand the impact of Be doping on surface active sites, adsorption modes of adsorbates and the detailed reaction mechanism of the oxygen reduction reaction. In this study, graphene has been doped with beryllium at varying monolayer concentrations of 0.03 ML, 0.06 ML and 0.09 ML. This work provides a basis for the study of binary, ternary and quaternary Be-doped systems and their efficiencies as a cathodic fuel cell material.

Computational details
All calculations were carried out using the spin-polarized density functional theory (DFT) method where periodic boundary conditions were applied as implemented within the Quantum-Espresso suite of programs (Giannozzi et al. 2009). This performs fully self-consistent DFT calculations to solve the Kohn-Sham equations (Kohn and Sham 1965). The generalized gradient approximated exchange correlation functional of Perdew, Burke, and Ernzerhof (Perdew et al. 1996) as well as the ultra-soft Vanderbilt pseudopotentials were used. The Grimme's D2 correction was implemented for dispersion corrections. The plane wave basis set was employed with a kinetic energy cutoff of 40 Ry, as obtained from convergence test. The Fermi-surface effects were treated by the smearing technique of Fermi-Dirac, using a smearing parameter of 0.003 Ry. An energy convergence threshold defining self-consistency of the electron density was set to 10 -6 eV and a beta defining mixing factor for selfconsistency of 0.2.
A p(3 × 3) supercell of graphene was employed in all adsorption studies which is made up of 32 carbon atoms, whereby adsorbates were more than 6 Å apart. Also, a vacuum region of 20 Å perpendicular to each surface was introduced to avoid interactions between periodic units. The Monkhorst-Pack k-point Brillouin Zone sampling method was implemented (Pack and Monkhorst 1977), with k-point sampling grids of (7 × 7 × 1) for the unit cell and (3 × 3 × 1) for the p(3 × 3) supercell. The graphics of the atomic structures have been prepared with the XCrysDen software (Kokalj 1999). Transition states were searched using the CI-NEB method, where structures at the saddle point are characterized by a single imaginary frequency.

Results and discussions
The resulting C-C bond of pure graphene has been reported to be 1.43 Å using density functional theory (Wu et al. 2010). Our calculated C-C bond is seen to be 1.42 Å, in good agreement with previous studies. In Fig. 1, Be impurities were introduced at three (3) unique doping sites, i.e. (D) lattice point, where carbon at the lattice site is replaced with impurity, site (E) where Be bridges two carbon atoms and site (F) where there is bridge doping and carbon vacancy at the neighboring lattice site.
In Fig. 1, three (3) surface oxygen adsorption sites (A, B and C) were considered for the adsorption of the O 2 molecule. The molecules were adsorbed at three unique sites on the surface, i.e., hollow (A), bridge (B) and top (C) sites.

Oxygen adsorption on bare graphene
Molecular oxygen was first adsorbed onto clean graphene at the surface adsorption sites i.e., (A), (B) and (C) as indicated in Fig. 1. The binding energies are calculated as shown below: where E ads is the adsorption energy indicating the binding strength and stability of products relative to reactants.E G+O 2 is the energy of adsorbed system, E G is the energy of isolated graphene and E O 2 is the energy of isolated oxygen molecule.
Oxygen can be bound to a transition metal in several ways, namely, dissociatively (typical O-O distances are around 2.5-2.8 Å), or associatively in either a peroxo (typical O-O distances are around 1.4-1.55 Å) or a superoxo (typical O-O distances are around 1.3 Å) form (Gutsev et al. 2000). Upon geometry optimization, oxygen adsorption on graphene is thermodynamically unfavoured, with a distance of 3.64 Å from graphene and a binding energy of 0.96 eV. Super-oxo (O 2− ) species are formed on pure graphene. This being a diatomic oxygen molecular species accounts for pure graphene promoting oxygen transformation through the twoelectron pathway (H 2 O 2 intermediate) and not the 4-electron pathway. We observed that the O-O bond is polarized resulting in bond activation from 1.23 Å in isolated oxygen molecule to about 1.25 Å in the adsorbed state. (See Table 1).

Be-doped graphene
On doping at the three investigated sites i.e., D, E and F (see Fig. 1), and upon optimization, the ground state structures were obtained (see Fig. 2). The preferred dopant site was investigated, where the defect formation energies were calculated with the formula below; where E D is the energy of the Be-doped graphene, E C is the energy of carbon atom in gas phase, E G is the energy of isolated graphene, E B is the energy of beryllium atom in gas phase and n is the number of dopants i.e., 1, 2 and 3 for 0.03 ML, 0.06 ML and 0.09 ML doping, respectively.
At 0.03 ML Be doping concentration, it is observed that the dopant prefers the point defective site (D) over the interstitial site (E) (Fig. 2). The bridge site doping with vacancy, i.e., structure (F) optimizes to (D). For  Beryllium binds at the point lattice site more favorable with a defect formation energy of 7.68 eV for D and F, and binds least at the interstitial site with a defect formation energy of 8.66 eV at the bridge (E) site. In structure E, around the defective point, graphene C-C bond shrinks from an inter-nuclear distance of 1.42 Å to about 1.34 Å and C-Be interaction bond is at 1.33 Å. In structure D and F, the C-Be bond is weaker at 1.56 Å and C-C shrinks less to 1.38 Å. The structural distortions are due to the higher radius of the impurity causing strain in the material. Beryllium doping of graphene is an endothermic process requiring energy at both the lattice and interstitial sites.

Oxygen adsorption on Be-doped graphene
Beryllium absorbs at the point defective site, more Be is introduced into the point lattice sites of graphene, to produce higher concentrations of 0.06 ML and 0.09 ML. Subsequently, oxygen adsorption studies were investigated on the 0.03 ML, 0.06 ML and 0.09 ML doped graphene sheets (see Table 2).
At 0.03 ML Be doping, the preferred binding of O 2 was investigated at the hollow, bridge and top sites using Eq. 1. It was observed that the oxygen molecule binds most strongly when absorbed at the top site of Be (see Fig. 3). At the hollow site, the oxygen molecule moved from the hollow position to the top site of Be, where oxygen is absorbed at the edge with a binding energy of 1.27 eV. On the bridge site, oxygen aligns to surface carbon atoms, with a binding energy of 2.16 eV. However, at the top site, oxygen binds most strongly leading to oxygen decomposition onto the   Table 2). The most active site and preferred binding site for oxygen on doped graphene at 0.03 ML concentration is seen to be the top site. Dissociation of the oxygen species is indicative of favorable transformation via the 4-electron pathway. At the dopant concentration of 0.06 ML, the binding energy of O 2 were investigated on two surfaces, an isolated defective surface (a) and a conjugated defective surface (b) as depicted in Fig. 4. This is to investigate the preferred binding site of O 2 when introduced to the surface, either to a populated defective site or an isolated defective site at 0.06 ML concentration. Increasing Be concentration generally increases the binding strength of oxygen on graphene. The binding of oxygen on the conjugated defective surface is preferred as oxygen is decomposed, whiles an activated per-oxo-like species is formed on the isolated defective surface. At the conjugated defective site, the binding energy is − 7.62 eV with a binding coordination number of 4. In the isolated case, a binding energy of − 1.63 eV is obtained.
Be dopant concentration in graphene lattice was increased to 0.09 ML (see Fig. 5), in the isolated (a) and conjugated (b) defective cases. This is to investigate the preferred binding site of O 2 when introduced to the surface, either to a populated defective site or an isolated defective site at 0.09 ML concentration. From Table 2, oxygen prefers to adsorb at the conjugated defective surface with a binding energy of − 1.63 eV over the isolated defective surface with energy of 0.98 eV. Activated peroxo-like species is obtained on the isolated defective surface whiles on the conjugated defective surface (where binding is preferred) a decomposed oxygen molecule is seen.
These results indicate that on 0.03 ML, 0.06 ML and 0.09 ML Be-doped graphene, oxygen is activated to dissociation into atomic oxygen which will favor transformation via the 4-electron pathway and not the 2-electron pathway as seen on pure graphene. Be impurities oriented in an isolated fashion are generally less efficient for molecular oxygen binding as compared to the conjugated defects. The oxygen adsorbed at the conjugated defective site favors decomposition, which in in turn favors the 4-electron pathway. This also reduces the possibility of reduction via the 2-electron pathway which is undesirable and competitively leads to hydrogen peroxide formation. From this study introduction of Be impurities in a single honeycomb ring of graphene has more impact on the binding efficiency of doped graphene.

ORR reaction mechanisms on Be-doped graphene
On the most active Be-doped surfaces for 0.03 ML, 0.06 ML and 0.09 ML doped surfaces the ORR was studied and compared. In the oxygen reduction reaction, two elementary steps were considered for the complete reduction of dissociated oxygen in to water i.e., Eqs. 3 and 4.
(3) 1. O * + H * → OH * For electrochemical processes, the barriers for proton transfer to adsorbates from solution are expected to be quite small. Hence, the reaction pathway is analyzed taking into account the hydrogen electrode potential as utilized by Peterson et al. (Peterson et al. 2010). Each reaction step involves an electrochemical proton electron transfer, and the free energy change of each step is determined as a function of the applied electrical potential known as the computational hydrogen electrode (CHE) model. In the CHE, the chemical potential of a proton electron pairs, µ(H + ) + µ(e − ) is equal to half of the chemical potential of gaseous hydrogen (1/2 µ(H 2 )) at a potential of 0 V. The chemical potential of each adsorbed species is calculated as free energies. (See Table 3).
∆Free energy for each compound (intermediates and products) is calculated as; However, entropy contributions were considered at 298 K using standard table. ∆ZPE is the change in ground state zero-point vibrational energy. The reaction energies (ΔE) of adsorbed systems are calculated relative to the reactant molecules (i.e., energy of isolated O 2 , H 2 and doped graphene) where; where E (products) is the energy of the intermediate and E (reactants) is the energy of reactants i.e., energy isolated O 2 , H 2 and doped graphene.
The least negative potential serves as the limiting potential. The change in free energy (ΔG n ) of each step will change with a simple linear function of applied potential (U) to eliminate the limiting potential.
where n is the number of proton-electron pairs transferred relative to CO 2 and e is the number of elementary positive charges transferred (number of electrons) when an external potential U was applied.
From Fig. 6, it is seen that the reaction opens up on 0.03 ML doped graphene at 0 V. On the 0.06 ML and 0.09 ML doped surfaces, the limiting step is the hydroxyl hydrogenation step at 0 V. Using Eq. 8, at 298 K, the reaction is exergonic at an onset potential of − 0.5 V and − 1.5 V on 0.09 ML and 0.06 ML respectively (see Fig. 7). The reaction is more feasible on both 0.03 ML and 0.09 ML surfaces at 0 V than on the 0.06 ML surface.
Using the Langmuir-Hinshelwood mechanism, we also considered hydrogenation of surface oxygen with where E TS is the absolute energy of the transition state and ∆E IS is the absolute energy of the intermediate state leading to product formation. From the reaction energy profile diagram in Fig. 8, kinetically, OH* formation is most feasible in the order 0.06 ML (1.62 eV) > 0.09 ML (1.69 eV) > 0.03 ML (1.84 eV). The subsequent step, i.e., water formation is most favored on 0.03 ML (0.59 eV) > 0.09 ML (1.95 eV) > 0.06 ML (2.53 eV) surfaces. Strong binding of intermediate species e.g., O* and OH* at active sites have been reported to control the reaction kinetics on transition metals where strong binding increases the barrier for ORR reactions by reducing surface activity.
As seen in Table 3, our results show that stronger binding of O* is consistent with lower barrier for OH formation. The O* binding strength shows the order 0.06 ML (− 7.42 eV) > 0.09 ML (-1.42 eV) > 0.03 ML (0.89 eV) and the reverse trend is seen for the O* hydrogenation barriers i.e., 0.03 ML > 0.09 ML > 0.06 ML with energies of 1.84 eV > 1.69 eV > 1.62 eV, respectively. This indicates that strong binding of atomic oxygen on Be-doped graphene does not impede its hydrogenation as expected.
Oxygen is rather hydrogenated easily as the hydrogenation product formed, i.e., OH* species, is not a leaving group. However, OH* hydrogenation leads to the formation of water, which is a leaving group, hence the strong OH* binding energies are consistent with higher energy barriers. On both 0.06 ML and 0.09 ML, OH* interactions are greater at 4.40 eV compared to the − 2.41 eV on 0.03 ML surface, hence hydrogenation and leaving of OH* is challenging on those surfaces. The hydrogenation of OH* is challenging and the rate determining step in the ORR (9) E a = E TS − E IS process on both 0.06 ML and 0.09 ML surfaces. On the 0.03 ML surface however, the hydrogenation of the O* species is the rate limiting step.
Overall, the 0.03 ML surface favors the ORR reaction kinetically due to the rate limiting barrier of 1.84 eV for the hydrogenation of O* compared to 2.5 eV and 1.9 eV on 0.06 ML and 0.09 ML surfaces respectively for the hydrogenation of OH*. Although substantial charge transfer to the O 2 molecule is desirable to favor the O 2 dissociation and the 4-electron pathway, additionally, surfaces that can preferentially bind O* strongly and OH* weakly will be kinetically suitable for the ORR. The Be point doped graphene at dilute concentrations like 0.03 ML is the most promising candidate that can be explored further at much lower concentrations for fuel cell applications.
Comparing our results to experimental and DFT results on P and N co-doped graphene, (Han and Chen 2020) surfaces with the highest net charge transfer to the O 2 molecule (as seen for 0.03 ML and 0.09 ML doped systems in this work), has the least onset potential of 0.419 V reported on PN 3 co-doped ratio. PN 3 co-doped graphene is reported to provide the lowest onset potential with the highest electrochemical ORR activity using the CHE model. The most active surface G-PN 3 adsorbs O 2 dissociatively, as opposed to associative adsorptions on the other surfaces. These results suggest that the Be doped surfaces (provide dissociate adsorptions) and are promising candidates for the ORR. The onset potential recorded on the 0.03 ML doped graphene i.e., 0.0 V is lower than 0.4 V predicted on the PN 3 doped surface, this implies Be dilute of graphene holds a great promise for the ORR reaction.

Conclusion
Using spin-polarized DFT calculations and the CHE method, the effect of p-type beryllium dopant on the ORR reaction was studied. Beryllium binds at the point lattice site of graphene more favorably than at the interstitial site. The presence of Be at the lattice site of graphene tuned the Fig. 8 Reaction Profile Diagram for oxygen reduction into water on Be-doped Graphene at 0 K in gas media selectivity of graphene towards the 4-electron ORR pathway by causing oxygen dissociation. Through the 4-electron pathway, two steps are encountered for O* and OH* hydrogenation into water. Kinetically, the reaction is seen to be controlled by the binding strength of O* and O*/H*. The stronger the binding of O*, O*/H* and H 2 O* the more kinetically challenging the process is as seen on the 0.06 ML doped surface. The ORR is favored most on the surface that transfers most charges to oxygen and binds oxygen species the weakest as observed for the single atom catalyst at 0.03 ML isolated doping. These results suggest that doping at lower Be concentrations below 0.03 ML, could be desirable to improve the rate of reaction. This work provides a basis for exploring more dilute Be doped materials and Be based ternary compounds as cathodic material for the ORR.