TM atoms on B/N doped defective graphene as a catalyst for oxygen reduction reaction: a theoretical study

Abstract

According to its high specific surface area and unique electronic properties, graphene with single transition metal atoms attached to defects in the graphene sheets is attractive for use in hydrogen fuel cells for oxygen reduction reaction on the cathode. It is motivated by the experimental observations for oxygen reduction reaction, and we use density function theory to systematically study the single transition metal atom-(B/N doped) vacancy complexes in graphene. The binding energies between single transition metal atoms and vacancies are calculated, along with adsorption energies of O₂, OOH, HOOH, O and OH. Our results indicate that N-doping can effectively improve the binding strength of metal atoms with divacancies. According to the adsorption energies of oxygen reduction reaction intermediates, it is found that Fe–N doped divacancy, Co–N doped divacancy and Zn–N doped divacancy complexes are promising candidates for use in hydrogen fuel cell cathodes for oxygen reduction reaction.

---

Introduction

The need for clean, sustainable sources of energy for transportation and portable application has inspired research in fuel cells over recent years. The oxygen reduction reaction (ORR), which attracts tremendous attention due to the slow kinetics that take place at the cathode electrode,^1,2 is a key electrocatalytic reaction for fuel cell applications. Up to now, the most widely used catalysts for ORR are still Pt and Pt-based alloys,^3–6 yet several problems such as high cost, high overpotential and low durability issues of Pt catalysts hinder its large scale application. Therefore, great efforts have been devoted to search for low-Pt and even metal free catalysts^7–13 aiming to overcome these problems.

Graphene based single metal atom–vacancy complexes have been successfully prepared in the laboratory.^14–19 Gan et al. had observed the metal atoms (Au and Pt) replaced carbon atoms and were located in single or multiple vacancies. The activation energies for the in-plane migration of both Au and Pt atoms in graphene are around 2.5 eV, indicating covalent bonding between metal and carbon atoms.¹⁴ The vacancies had been successfully created by high-energy atom/ion bombardment, which could be filled with different single atoms (Pt, Co, and In).¹⁸ The structures of a monovacancy (V1), a divacancy (V2) and a Pt atom trapped in a divacancy were shown in the HRTEM images. Robertson et al. used focused electron beam irradiation to create monovacancy and divacancy in graphene within a defined area, which then acted as trap sites for mobile Fe atoms initially resident on the graphene surface.¹⁹ They had confirmed that these Fe containing defect structures are stable in comparison to Fe atoms incorporated into graphene edge and the structures of a Fe@V1 and a Fe@V2 were shown in false collar AC-TEM images. The chemical doping has been proved to be an effective approach to tailor the property of the graphene. The substitutional dopants, for example, B^20–22 and N^21–25 atoms, have been used commonly to substitute carbon atoms in graphene and carbon nanotubes and to tune their properties. In recent years, lots of work has been done about the TM-(B/N doped) vacancy complexes.^26–40 Federico et al. discussed the oxygen reduction reaction (ORR) and oxygen evolution reaction (OER) activity of transition metal atom (belong to groups seven to nine in the periodic table)-N doped divacancy complex. They found that the transition metals do not have intrinsic catalytic activities since their adsorption behavior can be severely altered by changes in the local geometry of the active site, the chemical nature of the nearest neighbours.²⁷ The Fe atom-N doped divacancy complex had been investigated, and its profile suggested that O₂ molecules were inherently favorable for reduction into H₂O following a four-electron process.³⁵ Kattel et al. calculated the chemisorptions of ORR species O₂, O, OH, OOH, and H₂O and O–O bond breaking in peroxide occur on both Fe atom-N doped divacancy complex and Fe atom-two N atoms doped divacancy complex.^39,40 A complete ORR was predicted to occur via a single site 4e⁻ mechanism on both complexes. The Fe atom-two N atoms doped divacancy complexes were predicted to be prime candidate because it's higher stability. The Mn atom-N doped divacancy complex was predicted to also have novel catalytic activity for ORR using dispersion-corrected density function theory (DFT), which will be comparable to that of the Pt catalyst.³⁴ Kaukonen et al. studied the possibility of using various impurity atoms embedded in graphene with monovacancy and divacancy for oxygen reduction reaction (ORR) with DFT theory.²⁹ According to the calculation results, single metal atom (Ni, Pd, Pt, Sn)-divacancy complexes can be good candidates for the use in fuel cell cathodes for ORR. However, the binding mechanisms of metal atom with vacancies or doped vacancies have not been investigated in detail, and the chemisorptions of ORR species have not been studied systematically.

First-principles methods have been proven to be extremely successful in a rational design of solid catalysts.⁴¹ In the present work, we employ DFT to study the feasibility of defective graphene with single metal atom impurities for operation in cathodes of FCs. We stress that the aim of our study is not to explain a specific reaction or catalytic mechanism. Our focus is on pointing that the high stability is an important factor in determining a material's suitability as catalyst. In the following, we report how the single metal atoms are anchored on monovacancy or divacancy graphene, whose interaction should be as strong as possible to warrant high stability of the system. The monovacancy or divacancy doped with B or N atoms are further explored. In the end, the various ORR intermediates adsorption on the complexes have been investigated and compared so that the appropriate adsorption energy is expected to can be obtained in some complexes.

Methods and models

Spin-polarized DFT calculations are performed using the Dmol³ code⁴² embedded in the Materials Studio software. The generalized gradient approximation (GGA) with the Perdew and Wang (PW91)^43–45 is used to describe the exchange-correlation functional. A double numerical plus polarization (DNP) is used as the basis set, while the ion–electron interaction is described using DFT semicore pseudopotentials (DSPPs).⁴⁶ The convergence tolerance energy of 1.0 × 10⁵ Ha, maximum force of 0.002 Ha Å⁻¹, and maximum displacement of 0.005 Å are adopted in the geometry optimization. The smearing techniques are used to achieve self-consistent field convergence with a smearing value of 0.005 Ha. The graphene sheet is represented using a hexagonal supercell containing 72 atoms, with a p (6 × 6) structure in the x–y plane and a vacuum layer of 15 Å along the z direction placed between the sheets, which lead to negligible interaction between the periodic image.^47–49 For geometric optimization, the Brillouin zone integration is performed with a (5 × 5 × 1) k-point sampling.

The single transition metal atom (TM = Sc, Ti, V, Cr, Mn, Fe, Co, Ni, Cu, Zn, Pd, Pt) binding at the vacancies were created through three steps:^50,51 (i) the creation of a monovacancy (V1) and divacancy (V2), as shown in Fig. 1(a) and (b); (ii) the substitution of under-coordinated carbons by B or N atoms (VB1/VN1/VB2/VN2) and (iii) TM atom incorporation in the center of vacancy (TM–V1/TM–VB1/TM–VN1/TM–V2/TM–VB2/TM–VN2), as depicted in Fig. 1(c) and (d). Following the method as previous reports,^52–54 the formation energy E^f_V of vacancy and the binding energy (E_b) of a TM–vacancy complex, is calculated as:

E^f_V = E_V + xμ_C − E_Gra − yμ_N(B) (1)

E_b = E_TM–V − E_V − E_TM (2)

where E_V and E_Gra is the total energy of vacancy graphene (doped or not) and primitive graphene, respectively. μ_C is the chemical potential of C defined as the energy per C atom obtained from infinite graphene sheet,^52,55 μ_N is the chemical potential of N defined as half of the total energy of N₂ molecule,^56–58 μ_B is the chemical potential of B defined as energy of one atom in a metallic α–B phase.⁵⁹ E_TM–V, E_V and E_TM are the energies of TM–vacancy complex, pristine vacancy and the free TM atom,^52,60,61 respectively. The x and y is the number of removed C atoms and added N(B) atoms in both monovacancy and divacancy.

Fig. 1 Optimized geometry of V1 (a), V2 (b), and binding geometries of TM atom on the V1 (c), and V2 (d). Yellow, carbon atom; purple, TM atom.

After the geometry optimization, there appear to be slight deformation of the vacancy substrate. Thus, the binding energy is classified into two parts: one is the interaction energy and the other is the deformation energy, which are defined by eqn (3) and (4) respectively

E_inter = E_TM–V − E_V–a − E_TM (3)

E_defor = E_b − E_inter = (E_TM–V − E_V − E_TM) − (E_TM–V − E_V–a − E_TM) = E_V–a − E_V (4)

where E_V–a is the single point energy of the vacancy substrate, the structure of which comes from the optimized geometry of the TM–vacancy system.

The adsorption energy of X molecule, E_ads(X), on TM–vacancy complex is calculated as

E_ads(X) = E_X/TM–V − E_TM–V − E_X (5)

where E_X/TM–V, and E_X are the energies of the TM–vacancy complex adsorbed with X and the free X molecule respectively. All three types of energy will be calculated by using the same calculation setting. With this definition, a negative value indicates an exothermic adsorption. It is well known that the standard DFT could not properly treat the long-term dispersion interactions in which electronic overlap is very small, wherein the local density approximation has a tendency to overestimate the binding energy, compared with experimental results, while GGA tends to underestimate the binding energy. However, in our study, as the interaction between the absorbed and surface atoms was in short range, the DFT calculation was adopted, which could predict the physical properties of many atom systems to a reasonable accuracy. In order to evaluate the reliability of the absolute value, we also compared our results of O₂ adsorption energy on Fe–VN2 (−1.02 eV) with the previous results^{33,35,39,40,62} which have been confirmed with high activity. According to different calculated methods, the adsorption energy varies from −0.88 eV to −0.98 eV, agreeing well with our results. Therefore, it is rational that the molecular adsorption energies on TM–vacancy complex are calculated by DFT method. In particular, in our study, the adsorption energies are used to compare for which we need only relative values in the same method.

Results and discussions

Defective graphene sheet as substrate

Our geometry-optimized graphene sheet has an in-plane C–C distance of 1.42 Å, which agrees well with previous theoretical observations.^52,63,64 The formation energies and other properties of these six vacancies graphene structures were calculated as listed in Table S1.† The formation of V1 in Fig. 1(a) and V2 in Fig. 1(b) is energetically unfavorable (Table S1†), consistent with previous calculations.^16,52 The formation energies of B/N doped vacancies are compared, and it is found the formation energy in the order VN < VB < V in Table S1.† B/N dopant can reduce the number of the half-filled dangling electronic orbits of vacancies which are electronically unfavourable, thus leading to more favorable energetic configuration. In the cases of primitive monovacancy, the dangling electronic orbits are only half-filled and unsymmetrical, which destabilizes the defective carbon structure. In contrast, B or N doped monovacancies, the dangling electronic orbits of N are filled and the dangling electronic orbits of B atoms are empty, both of which are less destabilizing and symmetrical, consistent with DFT-predicted order of the vacancies formation energies in Table S1.† In the cases of divacancy, the dangling electronic orbits are also half-filled but symmetrical. The relaxed V1 and V2 geometries show that a Jahn–Teller distortion occurs and a weak bond is formed between carbon atoms bordering defects labelled by dotted line in Fig. 1(a) and (b).⁶⁵ The weak bonds are also confirmed by the deformation electron density as shown in Fig. S1.† The picture shows the electron accumulation in the weak bond site, meaning the formation of weak covalent bond.

TM adsorbed on the defective graphene complexes

We now turn to metal atoms adsorbed on vacancies in a defective graphene sheet. Such a structure can be associated with a substitutional impurity in a graphene. The typical atomic configuration of a TM–vacancy complex is shown in Fig. 1(c) and (d). Because the metal atomic radii is larger than that of the carbon atom, the metal atoms displace outwards from the graphene surface.^60,66 The binding energy E_b of the metal atom with vacancy is shown in Fig. 2. For the TM–V1 and TM–V2 complexes, the values of E_b are about −7 eV, except for the atoms with almost full d-shells, such as Cu and Zn. The binding energies of TM–V1 are stronger than the corresponding TM–V2, except for Cu–V1 and Zn–V1, which have a slight weaker binding strength than corresponding TM–V2. However, when the under-coordinated C atoms are replaced by B or N atoms, the binding strength of TM–V(N/B)2 is stronger than the corresponding TM–V(N/B)1. For the monovacancy, the introductions of B or N atoms into the vacancy decrease the binding energy of metal atom with monovacancy. However, the introductions of N atoms into the divacancy enhance the binding energy of metal atoms. The introductions of B or N into the vacancies have an important influence on the binding strength of TM. The trends of the binding strength E_b are similar for all TM–vacancy complexes, except for Sc and Ti. As shown in Fig. 2, it is observed that the most stable systems are TM–VN2 complexes, followed by TM–V1 complexes.

Fig. 2 The binding energy, E_b of the defective graphene sheet with TM atoms adsorbed on the vacancies.

The binding properties of TM adatom on vacancy are analyzed by calculating the charge density difference in real space. The charge density difference of Sc(Ti)–V1(V2/VN2/VB2) is shown in Fig. 3. The charge density difference of TM–V1 is defined as following equation:

Δρ = ρ(TM–V) − ρ(V) − ρ(TM) (6)

where the charge density is calculated in the same supercell for all quantities, ρ(TM–V) is the total charge density of the TM–vacancy system, and ρ(V) and ρ(TM) are the charge densities of noninteracting pristine vacancy and TM atom, respectively. The charge density difference quantifies the redistribution of electron charge due to the interaction between adatom and monovacancy. The accumulated charge density around the TM–C bonds indicates that TM is strongly bound on the site vacancy of graphene, in accordance with the high exothermic adsorption energy in Fig. 2, implying that the covalent bonding is formed between TM and C atoms. In addition, Fig. 3 shows a charge transfer from TM to C atom, which is consistent with the electronegativity difference between C and TM atoms. The Hirshfeld charge of the TM–vacancy systems is calculated in Table S2,† which also indicates that a great amount of electrons are transferred from TM to substrate.

Fig. 3 Top view and side view of the charge density difference for Sc/Ti atom adsorbed on vacancy complexes: Sc–V1 (a), Sc–V2 (b), Sc–VN2 (c), Sc–VB2 (d), Ti–V1 (e), Ti–V2 (f), Ti–VN2 (g); Ti–VB2 (h). The isosurface value is 0.004 e Å⁻³. The red colour represents electron accumulation, while the blue colour represents electron depletion. The intensity of colour depends on the amount of electron change: the darkest red marks the most accumulation; the darkest blue marks the most depletion. Yellow, carbon atom; blue, nitrogen atom; pink, boron atom; purple, Sc atom; light red, Ti atom.

TM–monovacancy complexes. We found that all metal atoms considered form chemical bonds with the under-coordinated C, N or B atoms at the vacancy by breaking the weak C–C or B–B bond at the pentagon in the reconstructed vacancy. Therefore, it is instructive that we chose the Sc–V and Ti–V complexes as examples to investigate binding mechanisms in detail. First, TM–V1 complexes are considered. The Sc atom (3d4s²) has three valence electrons. One can assume that two of them go into the Sc–C covalent σ bonds, while the third replaces the π electron of the missing carbon atom. Thus, as shown in Fig. 4(a), the bonding states will appear below the Femi level. It is shown that the bonding states are concentrated around the Sc–C bond while the antibonding states are concentrated around the TM atom site in Fig. 4(a). For the Ti atom (3d³4s²), it has four valence electrons and three of them go into the Ti–C covalent σ bonds, while the fourth replaces the π electron of the missing carbon atom. As shown in Fig. 4(b), the bonding states present below the Femi level, which are deeper than those in Sc–V1 complex. It is shown that the bonding states are concentrated around the Ti–C bond. The orbits between TM atom and adjacent carbon atoms can be seen clearly. The tetravalent Ti atom is a perfect substitution for a tetravalent C atom and four valence electrons are saturated. The three bonding electrons rather than two and the deeper location of bonding states of Ti–C confirm the stronger binding of Ti–V1.

Fig. 4 PDOS of Sc–V1 (a), Ti–V1 (b), Sc–VN1 (c), Ti–VN1 (d), Sc–VB1 (e), Ti–VB1 (f). The corresponding orbits at specific energy levels are shown in the insert pictures.

When the C atoms are substituted by N atoms, the dangling N atom will have lone pair electrons and a dangling electron. After removing and substituting of four C atoms, the number of remained C atoms in graphene sheet is even as well as the π electrons. The Sc atom (3d4s²) has three valence electrons and all of them will go into the Sc–N σ bond. However, the Ti atom (3d³4s²) has four valence electrons and three of them go into the Ti–N σ bond and the fourth one will act as an π electron. Four valence electrons are not saturated and the DOS is antisymmetric. Noticeably, it is observed that the filled antibonding bonding state (−0.37 eV) occurs below the Fermi level in Fig. 4(d), which will weaken strength of the bond between TM and graphene. In addition, the trivalent Sc atom is a more perfect substitution for a trivalent N atom. Therefore, Sc atom binds stronger than Ti atom.

When the C atoms are substituted by B atoms, the dangling B atom will have only one dangling electron. After removing and substituting of four C atoms, the number of remained C atoms is even. The situation is similar to N-doped case, and we also observe the antibonding states in Ti–VB1, locates at −0.37 eV. For both B and N doping, the number of electrons participating in bonding states is small compared to primitive monovacancy. This may explain the weaker binding of TM–VN(B)1 than TM–V1.

TM-divacancy complexes. The similar ideas can be applied to TM–V2 complexes. In a naked V2 complex, every dangling bond atom has one potential σ electrons and one π electron. After removing two C atoms, the number of remained C atoms is even as well as the π electrons. The Sc atom has three valence electrons, and all of them go into the Sc–C covalent σ bonds while Ti atom has four valence electrons which enter the Ti–C covalent σ bonds. The four bonding electrons in Ti–V2 complex rather than three in Sc–V2 complex, and the deeper location of bonding states of Ti–C, as shown in Fig. 5(a) and (b), confirm the stronger binding of Ti–V2.

Fig. 5 PDOS of Sc–V2 (a), Ti–V2 (b), Sc–VN2 (c), Ti–VN2 (d), Sc–VB2 (e), Ti–VB2 (f). The corresponding orbits at specific energy levels are shown in the insert pictures.

When the carbon atoms are substituted by N atoms, the dangling N atom will have lone pair electrons and a dangling electron. After removing and substituting of six C atoms, the number of remained C atoms is even as well as the π electrons. The binding mechanism is similar to the monovacancy complex, except that there are antibonding states in the deeper energy level for both Sc–VN2 and Ti–VN2, as shown in Fig. 5(c) and (d). We observed the similar conditions in Sc–VB2 and Ti–VB2, as shown in Fig. 5(e) and (f).

The enhancement of binding between TM and N-doped divacancy. Above has analyzed the interaction between TM and defective graphene sheet in detail, which indicates that the bonds between TM and VN2 are weak due to the deeper antibonding states. On the other hand, our calculations gives the contrary result that the binding strength of TM–VN2 is the strongest among the complexes in Fig. 2. The binding energy of monovacancy increases in the order E_b(VB1) < E_b(VN1) < E_b(V1) while the binding energy of divacancy increases in the different order E_b(VB2) < E_b(V2) < E_b(VN2). In order to explore the enhanced binding strength between TM and VN2, we have calculated the interaction energy between TM atoms and vacancy as well as the deformation energy of substrate after binding with TM atoms. The energies are shown in Table S3.† For interaction energy, it increases in the order E_inter(VB1) < E_interb(VN1) < E_inter(V1) and E_inter(VB2) < E_interb(VN2) < E_inter(V2). It indicates that TM atoms interact most strongly with C atom in primitive vacancies and the interaction energy between TM atoms and divacancy is stronger than with monovacancy.

After the adsorption of TM on the vacancies complex, the graphene wills go under deformation. The deformation energy of TM–V2 is nearly 2.5 eV, which is much higher than the TM–VB2 and TM–VN2 (less than 1 eV). The large deformation energy of TM–V2 complex makes this complex less stable, which also explains the smaller binding energy of TM–V2 than TM–V1. However, the comparable high interaction energy as well as low deformation energy contributes to the large binding energy of TM–VN2 complex. Therefore, it is concluded that the enhanced binding between TM and VN2 mainly comes from the decreased deformation energy. For the monovacancy complex, their deformation energies are comparable and effect is equal.

The oxygen reduction reaction intermediates adsorbed on the TM–vacancy complexes

To catalyze the ORR, the catalyst must be able to adsorb an oxygen molecule. However, the adsorption energy of O₂ on the catalyst's surface should be in an appropriate range. If the adsorption strength between the catalyst and O₂ molecule is very weak, the coverage of O₂ may be too low to have a high reaction rate and the O–O bond would be too strong to be broken. On the other hand, if O₂ binds very strong to the catalyst, it would poison the surface and reduce its activity. Previous works have shown that platinum demonstrates a high ORR activity. Theoretical studies of O₂ adsorbed on the Pt surface suggest that the adsorption energy is between −0.53 and −1.02 eV,^51,67,68 based on different models and methods. Therefore, the adsorption energy on the Pt surface can be used as an important reference to evaluate ORR activity of other catalysts. In our systems, TM–VN2 complexes have the strongest binding strength in comparison with the other five complexes. Furthermore, it is well known that the TM–N4 fragments have been widely investigated.^{35,39,40,69,70} Therefore, O₂ adsorbed on TM–VN2 will be investigated systematically. According to previous study, the oxygen reduction reaction (ORR) on no-metal catalyst at the acidic conditions, the number of electron transfers is generally less than 4, indicating that both 4e and 2e reduction occur simultaneously in ORR. It is believed that the 4e reduction involves sequentially the*O₂, *OOH, *OH and *O intermediates (* indicates the adsorbed species) as every elementary step of the reaction is exothermic. For the 2e reduction, it could involve *OOH and *HOOH intermediates (without O–O bond breaking) and H₂O₂ is as the final product. The proposed mechanisms can be summarized as eqn (7)–(15). Eqn (7)–(13) are the 4e reduction pathway and eqn (14) and (15) are the 2e pathway. The *OOH intermediate appears in both pathways and thus is the key species determining the selectivity. The main purpose of the current work is to evaluate the potential catalyst that can promote the ORR. The equations are shown as:

* + O₂ → ^*O₂ (7)

^*O₂ + H_(aq)⁺ + e⁻ → *OOH (8)

*OOH + H_(aq)⁺ + e⁻ → 2*OH (9)

*OOH + H_(aq)⁺ + e⁻ → *O + H₂O (10)

*O + H_(aq)⁺ + e⁻ → *OH (11)

*OH + H_(aq)⁺ + e⁻ → *H₂O (12)

*H₂O → H₂O_(aq) + * (13)

*OOH + H_(aq)⁺ + e⁻ → *HOOH (14)

*HOOH → H₂O_2(aq) + * (15)

where * indicates an empty adsorption site and *X indicates an adsorbed molecule. These reactions consist of three main ORR pathways:

I. (7) → (8) → (9) → (12) → (13)

II. (7) → (8) → (10) → (11) → (12) → (13)

III. (7) → (8) → (14) → (15)

The reaction pathways I and II mainly differ by the formation of 2*OH (9) and *O + H₂O (10), respectively. Both of the pathways have four electrons transportation. The reaction pathway III is totally different from reaction pathway I and II because only two electrons transportation occurs in reaction. For reaction pathway I and II, the rate-determining step may be O–O break or the remove of OH.

The adsorption of O₂, OOH, HOOH, O and OH are shown in Fig. 6. Upon binding of O₂ the bond length TM–N increase slightly by ∼0.00–0.09 Å, as shown in Table S4.† The TM–N bonds gain length either equally or two bonds become longer compared with the two other ones according to different adsorption geometry. For most cases, the dioxygen prefers to bind horizontally along the graphene plane, Fig. 6(a), with the exception of Fe, Co, Pd and Pt, as shown in Fig. 6(b), via one oxygen atom mode. Upon binding of OOH the bond length TM–N increase by ∼0.01–0.08 Å, similar to O₂ adsorption. The OOH prefers to bind via one oxygen atom along the graphene plane for most cases, as shown in Fig. 6(d), instead of a two oxygen mode (Sc, Ti and V, in Fig. 6(c)). Upon binding of HOOH, the bond length TM–N increase by ∼0.01–0.33 Å. The HOOH will be dissociated into two OH molecular fragments for TM–VN2 (TM = Sc, Ti, V, Cr, Mn, Fe, Zn), as shown in Fig. 6(e), or dissociated into one O atom and one H₂O molecular for Co–VN2, as shown in Fig. 6(f) or keep molecular structures for TM (TM = Ni, Cu, Pd, Pt)–VN2, as shown in Fig. 6(g). We also considered the adsorption of O atom (Fig. 6(h) and (i)) and OH (Fig 6(j) and (k)). The O atom was adsorbed on top of TM when TM belongs to Sc, Ti, V, Cr, Mn, Fe and Co, as shown in Fig. 6(h). For the other TM–VN2 complexes, the most stable O atom adsorption site was on top of the next nearest neighbor C atom, as shown in Fig. 6(i). Upon binding of O atom, the bond length TM–N increase by ∼0.00–0.06 Å.

Fig. 6 The adsorption geometries of O₂ (a) and (b), OOH (c) and (d), HOOH (e)–(g), O (h) and (i) and OH (j) and (k) on TM–VN2 complexes. Yellow, carbon atom; blue, nitrogen atom; purple, TM atom; red, oxygen atom; white, hydrogen atom.

For the adsorption of OH, the top of TM was the most stable site for most of TM, as shown in Fig. 6(j), except for Pd and Pt, which prefer to top of the next nearest neighbor C atom, as shown in Fig. 6(k). Upon binding of OH, the bond length TM–N increase by ∼0.00–0.17 Å. For all the adsorption fragments, the OH or HOOH (two OH) have the most significant influence on the structure of the substrate.

The adsorption energies of O₂, OOH, HOOH, O and OH on TM–VN2 were calculated and shown in Fig. 7(a) and Table S5.† The adsorption energy increases with the increasing d electrons of TM atoms. The adsorption energies of O₂, OOH, HOOH, O and OH have same trends and we give the linear relationship in Fig. 7(b). A good defect-candidate for ORR should first have a high binding energy to insure stability, and second it should bind O₂ appropriate to insure the add H and the remove of OH. We use the adsorption energies of O₂ in the scale from −0.50 to −1.0 eV as the criterion to estimate the activity. The adsorption energies of O₂ on TM–VN2 are compared with the situation on previous results.^{35,39,40,69,70} We find that the adsorption energies of various ORR intermediates of Fe–VN2 and Co–VN2 are comparable with previous results. They have been confirmed to be active for ORR. According to our calculations, HOOH will be dissociated into two OH spontaneously for Fe–VN2 structure and the ORR will be occur under 4e reduction pathways (I) and the rate-determining step will be the remove of OH. For Co–VN2, however, HOOH will break into O and H₂O spontaneously and the ORR will be under 4e reduction pathways (II) and the rate-determining step will be the remove of OH too. So the ORR will be under 4e reduction pathways for both Fe–VN2 and Co–VN2 and the rate-determining step will be the remove of OH. In the TM–VN2 (TM = Ni, Cu, Pd, Pt) systems, the HOOH will keep molecular state and the ORR will be under 2e reduction pathways (III) and the rate-determining step will be the break of O–O bond. However, our results of Mn–VN2 are different from others.³⁴ The PDOS of O₂ adsorption on TM–VN2 (TM = Mn, Fe, Co, Ni) are shown in Fig. 8. It is observed that the bonding states appear below Fermi level for four TM–VN2 complexes. The deeper the bonding states, the stronger the binding strength. However, the bonding states of the Co–VN2 are deeper than the Fe–VN2, but antibonding states also appear near the Fermi level, which is the reason why the O₂ adsorption on Co–VN2 is weaker than on Fe–VN2.

Fig. 7 The adsorption energies of O₂, OOH, HOOH, O and OH on the TM–VN2 complexes (a) and their scale relationship (b).

Fig. 8 PDOS of O₂ adsorb on Mn–VN2 (a), Fe–VN2 (b), Co–VN2 (c) and Ni–VN2 (d) complex. The corresponding orbits at specific energy levels are shown in the insert pictures.

Conclusion

Using first principles calculations, we have studied the possibility of using various impurity metal atoms embedded in graphene with vacancies for oxygen reduction reaction (ORR). A good single metal atom–vacancies complexes candidate for ORR should first have a high binding energy between metal atom and vacancy to insure stability, and second it should bind O₂ appropriate to ensure the break of O–O bond and remove of H₂O molecular. TM–VN2 is the most stable one among the six kinds of complexes. We considered the various ORR intermediates adsorption on TM–VN2 complexes, and it was found that Fe–VN2, Co–VN2 and Zn–VN2 should be potential ORR catalyst.

Acknowledgements

The authors would like to thank the funding from National Natural Science Foundation of China (No. 51372095) and the “211” and “985” projects of Jilin University. The DFT calculations utilized resources at the High Performance Computing Center, Jilin University.

Notes and references

1. J. K. Nørskov, J. Rossmeisl, A. Logadottir, L. Lindqvist, J. R. Kitchin, T. Bligaard and H. Jonsson, J. Phys. Chem. B, 2004, 108, 17886–17892.
2. Z. Feng, H. Falk, S. Uwe, S. Fritz, B. Peter and H. Iris, Environ. Sci. Technol., 2006, 40, 5193–5199.
3. U. A. Paulus, A. Wokaun, G. G. Scherer, T. J. Schmidt, V. Stamenkovic, N. M. Markovic and P. N. Ross, Electrochim. Acta, 2002, 47, 3787–3798.
4. H. R. Colón-Mercado and B. N. Popov, J. Power Sources, 2006, 155, 253–263.
5. J. Greeley, I. E. L. Stephens, A. S. Bondarenko, T. P. Johansson, H. A. Hansen, T. F. Jaramillo, J. Rossmeisl, I. Chorkendorff and J. K. Nørskov, Nat. Chem., 2009, 1, 552–556.
6. C. H. Cui, L. Gan, M. Heggen, S. Rudi and P. Strasser, Nat. Mater., 2013, 12, 765–771.
7. S. Kattel, P. Atanassov and B. Kiefer, J. Mater. Chem. A, 2014, 2, 10273–10279.
8. J. Wu, Z. Yang, X. Li, Q. Sun, C. Jin, P. Strasser and R. Yang, J. Mater. Chem. A, 2013, 1, 9889–9896.
9. Z. Chen, D. Higgins, A. Yu, L. Zhang and J. Zhang, Energy Environ. Sci., 2011, 4, 3167–3192.
10. R. Bashyam and P. Zelenay, Nature, 2006, 443, 63–66.
11. R. Othman, A. L. Dicks and Z. H. Zhu, Int. J. Hydrogen Energy, 2012, 37, 357–372.
12. F. Jaouen, E. Proietti, M. Lefèvre, R. Chenitz, J.-P. Dodelet, G. Wu, H. T. Chung, C. M. Johnston and P. Zelenay, Energy Environ. Sci., 2011, 4, 114–130.
13. M. Lefevre, E. Proietti, F. Jaouen and J. P. Dodelet, Science, 2009, 324, 71–74.
14. Y. J. Gan, L. T. Sun and F. Banhart, Small, 2008, 4, 587–591.
15. J. A. Rodriguez-Manzo, O. Cretu and F. Banhart, ACS Nano, 2010, 4, 3422–3428.
16. J. K. Florian Banhart and A. V. Krasheninnikov, ACS Nano, 2011, 5, 26–41.
17. A. W. Robertson and J. H. Warner, Nanoscale, 2013, 5, 4079–4093.
18. H. Wang, et al., Nano Lett., 2012, 12, 141–144.
19. A. W. Robertson, et al., Nano Lett., 2013, 13, 1468–1475.
20. Y. G. Zhou, X. T. Zu, F. Gao, J. L. Nie and H. Y. Xiao, J. Appl. Phys., 2009, 105, 014309.
21. S. S. Yu, W. T. Zheng, C. Wang and Q. Jiang, ACS Nano, 2010, 4, 7619–7629.
22. S. S. Yu and W. T. Zheng, Nanoscale, 2010, 2, 1069–1082.
23. Y. Fujimoto and S. Saito, Phys. Rev. B: Condens. Matter Mater. Phys., 2011, 84, 245446.
24. H. M. Jeong, J. W. Lee, W. H. Shin, Y. J. Choi, H. J. Shin, J. K. Kang and J. W. Choi, Nano Lett., 2011, 11, 2472–2477.
25. L. P. Zhang, J. B. Niu, L. Dai and Z. H. Xia, Langmuir, 2012, 28, 7542–7550.
26. Y. Tang, Z. Yang and X. Dai, J. Chem. Phys., 2011, 135, 224704.
27. F. Calle-Vallejo, J. I. Martinez and J. Rossmeisl, Phys. Chem. Chem. Phys., 2011, 13, 15639–15643.
28. M. J. López, I. Cabria and J. A. Alonso, J. Phys. Chem. C, 2014, 118, 5081–5090.
29. M. Kaukonen, A. V. Krasheninnikov, E. Kauppinen and R. M. Nieminen, ACS Catal., 2013, 3, 159–165.
30. K. Li, Y. Li, H. Tang, M. Jiao, Y. Wang and Z. Wu, RSC Adv., 2015, 5, 16394–16399.
31. Y. Tang, Z. Yang and X. Dai, Phys. Chem. Chem. Phys., 2012, 14, 16566–16572.
32. M. Z. Yun-Hao Lu, C. Zhang and Y.-P. Feng, J. Phys. Chem. C, 2009, 113, 20156–20160.
33. A. G. Saputro and H. Kasai, Phys. Chem. Chem. Phys., 2015, 17, 3059–3071.
34. Z. Lu, G. Xu, C. He, T. Wang, L. Yang, Z. Yang and D. Ma, Carbon, 2015, 84, 500–508.
35. J. Zhang, Z. Wang and Z. Zhu, J. Power Sources, 2014, 255, 65–69.
36. C. E. Szakacs, M. Lefevre, U. I. Kramm, J. P. Dodelet and F. Vidal, Phys. Chem. Chem. Phys., 2014, 16, 13654–13661.
37. J. Sun, Y. H. Fang and Z. P. Liu, Phys. Chem. Chem. Phys., 2014, 16, 13733–13740.
38. W. Liang, J. Chen, Y. Liu and S. Chen, ACS Catal., 2014, 4, 4170–4177.
39. S. Kattel and G. Wang, J. Phys. Chem. Lett., 2014, 5, 452–456.
40. S. Kattel, P. Atanassov and B. Kiefer, Phys. Chem. Chem. Phys., 2014, 16, 13800–13806.
41. J. K. Norskov, T. Bligaard, J. Rossmeisl and C. H. Christensen, Nat. Chem., 2009, 1, 37–46.
42. B. Delley, J. Chem. Phys., 2000, 113, 7756.
43. J. P. Perdew, K. A. Jackson, M. R. Pederson, D. J. Singh and C. Fiolhais, Phys. Rev. B: Condens. Matter Mater. Phys., 1992, 46, 6671–6687.
44. J. P. Perdew and Y. Wang, Phys. Rev. B: Condens. Matter Mater. Phys., 1992, 45, 13244–13249.
45. J. A. White, D. M. Bird, M. C. Payne and I. Stich, Phys. Rev. Lett., 1994, 73, 1404–1407.
46. B. Delley, Phys. Rev. B: Condens. Matter Mater. Phys., 2002, 66, 155125.
47. Y. Gao, N. Zhao, J. Li, E. Liu, C. He and C. Shi, Int. J. Hydrogen Energy, 2012, 37, 11835–11841.
48. L.-J. Zhou, Z. F. Hou and L.-M. Wu, J. Phys. Chem. C, 2012, 116, 21780–21787.
49. L. Wu, T. Hou, Y. Li, K. S. Chan and S.-T. Lee, J. Phys. Chem. C, 2013, 117, 17066–17072.
50. W. I. Choi, S.-H. Jhi, K. Kim and Y.-H. Kim, Phys. Rev. B: Condens. Matter Mater. Phys., 2010, 81, 085441.
51. Y. Feng, F. Li, Z. Hu, X. Luo, L. Zhang, X.-F. Zhou, H.-T. Wang, J.-J. Xu and E. G. Wang, Phys. Rev. B: Condens. Matter Mater. Phys., 2012, 85, 155454.
52. S. Kattel, P. Atanassov and B. Kiefer, J. Phys. Chem. C, 2012, 116, 8161–8166.
53. A. T. Lee, J. Kang, S.-H. Wei, K. J. Chang and Y.-H. Kim, Phys. Rev. B: Condens. Matter Mater. Phys., 2012, 86, 165403.
54. S. Kattel, RSC Adv., 2013, 3, 21110.
55. Y. Li, Z. Zhou, P. Shen and Z. Chen, ACS Nano, 2009, 3, 1952–1958.
56. F. Emanuele, D. V. Cristiana, S. Annabella and P. Gianfranco, J. Phys. Chem. C, 2007, 111, 9275–9282.
57. Y. Kesong, D. Ying, H. Baibiao and H. Shenghao, J. Phys. Chem. B, 2006, 110, 24011–24014.
58. A. G. Garcia, S. E. Baltazar, A. H. R. Castro, J. F. P. Robles and A. Rubio, J. Comput. Theor. Nanosci., 2008, 5, 2221–2229.
59. S. Azevedo, Phys. Lett. A, 2006, 351, 109–112.
60. A. V. Krasheninnikov, P. O. Lehtinen, A. S. Foster, P. Pyykkö and R. M. Nieminen, Phys. Rev. Lett., 2009, 102, 126807.
61. Y. Shang, J.-X. Zhao, H. Wu, Q.-H. Cai, X.-G. Wang and X.-Z. Wang, Theor. Chem. Acc., 2010, 127, 727–733.
62. W. Orellana, J. Phys. Chem. C, 2013, 117, 9812–9818.
63. D. E. Jiang, B. G. Sumpter and S. Dai, J. Chem. Phys., 2007, 126, 134701.
64. C. Rajesh, C. Majumder, H. Mizuseki and Y. Kawazoe, J. Chem. Phys., 2009, 130, 124911.
65. O. V. Yazyev, Rep. Prog. Phys., 2010, 73, 056501.
66. H. Sevinçli, M. Topsakal, E. Durgun and S. Ciraci, Phys. Rev. B: Condens. Matter Mater. Phys., 2008, 77, 195434.
67. L. Qi, X. Qian and J. Li, Phys. Rev. Lett., 2008, 101, 146101.
68. M. P. Hyman and J. W. Medlin, J. Phys. Chem. C, 2007, 111, 17052–17060.
69. S. Kattel, P. Atanassov and B. Kiefer, Phys. Chem. Chem. Phys., 2013, 15, 148–153.
70. S. Kattel and G. Wang, J. Mater. Chem. A, 2013, 1, 10790–10797.

---

Footnote

† Electronic supplementary information (ESI) available. See DOI: 10.1039/c5ra15315h

---