Density Functional Theory Study of the Oxygen Reduction Reaction Mechanism on Graphene Doped with Nitrogen and a Transition Metal

The active centers of carbon nonplatinum catalysts doped with cobalt, iron, nickel, and copper have been simulated by quantum-chemical density functional theory methods. The thermodynamics of the electrochemical oxygen reduction reaction (ORR) on model catalysts has been determined. It was found that among the studied catalysts, graphene doped with cobalt and iron showed the best properties. A two-state reactivity effect has been found on a cobalt-containing catalyst, and a more detailed reaction mechanism has been proposed, including the stages of charging by an extra electron and association with water. The proposed mechanism explains several effects that have arisen during the modeling in relation to the classical mechanism.


INTRODUCTION
The modern world every day becomes more and more dependent on electricity, which can be generated in various ways. Fossil fuels are a convenient but nonrenewable resource that pollute the atmosphere. Alternative energy seems to be a solution to this problem; however, electricity must be not only generated but also stored until the moment of use; the mismatch in the time of consumption and generation of electricity is one of the main problems of alternative energy. 1 Both these problems can be resolved by hydrogen energy, namely, by the use of fuel cells. The main obstacle to their widespread use is the presence of platinum as a catalyst, due to which the price of fuel cells increases significantly. Thus, the aim of this work was to search for alternative catalysts for one of the main reactions occurring in fuel cellsthe electrochemical oxygen reduction reaction (ORR). Many works have been devoted to the study of ORR catalysts. Platinum and materials based on it are often used as an electrode. 2−4 However, in recent years, there has been an increasing interest in nonplatinum catalystscarbon materials including various metals and their compositions (mixtures, alloys, and metals embedded in the support structure). 5,6 Of particular interest are materials with active sites included in the structure of a carbon material, formed from a metal atom and surrounding nonmetal atoms (nitrogen, oxygen, sulfur, phosphorus, etc.). 7 These nonmetallic dopants with different electronegativities can cause polarization of the carbon skeleton and create charged active centers, which can significantly improve the slow kinetics of the ORR. 8−12 Recent experimental studies and DFT calculations have shown that monoatomic catalysts including a transition metal (TM) and nonmetals such as TM−NC (TM = Fe, Co, Mn, Ni, Zn, etc.) 8,11−18 exhibit excellent electrocatalytic activity toward the ORR. Graphene, especially defective graphene with vacancies, can be an excellent support for monoatomic catalysts. 19 The delectrons of the central transition metal can be regulated by nitrogen coordination, and hence, their electrocatalytic performance for the ORR is improved. 20,21 Researchers have reported 22 a high-performance iron-based electrocatalyst coated with Fe cations coordinated by the N (Fe−N 4 ) pyridine type. Through DFT calculations, they confirmed that the density of states of the central Fe cations was regulated by coordinated N atoms at the active sites of Fe−N 4 , resulting in an efficient four-electron process and reduced overpotential for the ORR. There are also DFT data 23 reporting that the ORR occurs via the 4e pathway on the FeN 4 catalysts with an activation energy comparable to that for modern Pt-based catalysts.
To accelerate the process of searching for the optimal material, theoretical methods are employed, such as quantumchemical modeling using density functional theory (DFT) methods. 24,25 This approach allows us to determine the adsorption properties of the catalyst and the effect of the composition and structure of the material on the kinetics of the catalytic process and to compare various catalytic materials without resorting to experimental analysis. However, modeling can be complicated by the possibility of the reaction proceeding through states with different multiplicities in the low-and highspin states, 17 accompanied by transitions from one state to another. The possibility of such processes necessitates additional calculations in the simulation.
In this work, the properties of carbon materials doped with nitrogen and metals: cobalt, iron, nickel, and copper are studied using quantum-chemical modeling.

MODEL AND CALCULATION METHODS
The calculations were carried out using the Gaussian 09 software. 26 The simulation results were visualized using the GaussView 6 software. The modeling was carried out by the density functional theory (DFT) method using the B3LYP functional 27,28 and the 6-31G* basis set, similar to the previous studies. 29,30 The effect of the solvent (water) was taken into account in the model of self-consistent reaction field (SCRF). 31,32  To analyze the catalytic parameters of materials, the surface regions of carbon materials were modeled using a fragment of graphene (10 conjugated aromatic rings) and fragments of graphene with nitrogen atoms and the studied metal incorporated in the structure, namely, a metal atom surrounded by four N atoms inside 10 conjugated aromatic rings ( Figure 1). Experimentally, such structures can be obtained using phthalocyanines of the corresponding metals as precursors. 33 The choice of structures was based on the rationale for the available experimental data, which will be published later.
The reaction of electrochemical oxygen reduction can proceed mainly in two ways 34 the two-and four-electron reduction schemes according to the associative and dissociative mechanisms. The reactions 1a, 1b, and 2 represent the ORR scheme in an alkaline solution, where 1 and 2 correspond to the associative and dissociative mechanisms, respectively.
The four-electron reduction associative mechanism 1b was chosen as the predominant one for this type of active center. 35 The following intermediates were chosen as the key points for modeling 36 ( To compare the energies of the optimized structures, the Gibbs free energy was used as the electron energy corrected with the thermal free energy, calculated as follows 26 where E ele is the total electronic and nuclear repulsion energy at 0 K, E ZPE is the zero-point vibrational energy, E therm is the total thermal internal energy, k B is Boltzmann's constant, S tot is the system entropy, and T is the temperature (T = 298.15 K).
where E t is the thermal internal energy due to translation, E r is the internal energy due to rotational motion, E v is the internal energy due to vibrational motion, and E e is the internal energy due to electronic motion. The sum of the energies of the final optimized substances and the catalyst was taken as a zero level.
The adsorption energy was calculated using the formula 37 To determine the energy effects of the ongoing processes, the following elementary reactions occurring on the catalyst surface were identified 38 ( OH (ads) + e − → OH − which correspond to the transitions of intermediates 2 → 3, 3 → 4, 4 → 5, and 5 → 6, respectively (see Table 1). The free energy of elementary reactions (used to construct all graphs and tables) was calculated as follows taking into account the pH of the solution and the change in the electrode potential 39 where ΔG s is the free energy change of the system, eU is the contribution of free energy due to a change in the values of the electrode potential U, and k B T·ln 10·pH is the free energy contribution due to changes in the pH values. The value pH = 14 was considered for further calculations. According to the calculated results, the total change in free energy at the above stages, that is, for the overall O 2 + 2H 2 O + 4e − → 4OH − reaction, is 8.10 eV at zero electrode potential and pH, which is close to the value computed in a recent work using DFT calculations with the PBE functional; 40 this indicates that the accuracy of the present calculations is adequate.
To analyze the efficiency of the process of electrochemical oxygen reduction, the overpotential ηORR was calculated as follows 41 where 1.20 is 1/4 of the total free energy change and ΔG i is the free energy change for each stage.

RESULTS AND DISCUSSION
Considering the optimized structures of the intermediates simulating the catalytic process (Table 1), it can be concluded that the adsorption of oxygen is localized on the surface, and therefore, the active center of the catalyst is a metal atom embedded in the carbon structure. The energy profiles of the investigated electrochemical reaction of oxygen reduction both in the free form and on the studied carbon materials were obtained as a result of modeling (see Figures 2, 3, and 4).
A carbonaceous material doped with one nitrogen atom was modeled. Although the nitrogen atom was an adsorption center, it did not significantly affect the activity of the catalyst ( Figure  2). Possibly, such a discrepancy with the published works 16,18 demonstrating successful catalytic applications of N-doped graphene can be explained by the location of nitrogen in the model in the graphene layer rather than on the edge surface. A low catalytic activity of such an active site (basal-type graphitic nitrogen) has been confirmed by other researchers. 32 Analyzing the thermodynamics of the elementary reactions on the catalyst containing cobalt, we observed the effect of two-state reactivity, 42 in which the minimum was observed alternately in the low-and high-spin states (in doublet and quartet states, respectively) ( Figure 3).
Thus, during the adsorption of OOH and OH radicals, a lower energy of the system was observed in the quartet state, while for other stages, the doublet state was lower in energy. Therefore, this effect makes the reaction easier to proceed due to less overpotential on the electrodes, which can improve the efficiency of the catalyst. The total profile of the states with the minimum energy was taken to compare with the other catalysts, namely, a quartet for adsorbed OOH and OH and a doublet for the rest. Certainly, the mechanism and kinetics of the two-state reactivity would depend on the structure and spin− orbit coupling of the pertinent minimal energy points on the doublet-quartet seam of crossing, 43,44 but its detailed analysis is beyond the scope of the present work.
Analyzing the profiles of the reaction on the catalysts, a gradual decrease in energy could be seen, which indicates the possibility of the reaction proceeding according to the selected mechanism ( Figure 4).
Analyzing the thermodynamics of oxygen sorption on the active sites of the catalyst (Table 2), we observed an increase in the Gibbs free energy due to the entropy factor according to the calculations. The energy of oxygen adsorption on the active center is in the range of 0.17−0.29 eV (16−28 kJ/mol) when summing only the electronic energy and the energy of zero-point vibrations, that is, considering only the energy of adsorption at 0 K. Cobalt and iron show the greatest affinity for oxygen, which manifests itself in the largest decrease in the energy of the system. However, the contribution of the thermal enthalpy and entropy factors, which amounts to ∼0.4 eV, makes the adsorption process slightly unfavorable in terms of the free energy. A slight increase in the free energy of the system at the stage of O 2 adsorption was noted in the literature too. 31 In    The energies of the intermediate compound adsorption on the active sites of the catalysts were calculated ( Table 3). The adsorption of OOH was characterized by a relatively low adsorption energy (<0.75 eV), while a high affinity for the catalysts was found for O and OH (>1 eV). Weak adsorption of the peroxide radical could lead to an increase in the fraction of the byproduct of the reaction, hydrogen peroxide, that is, an increase in the proportion of the process proceeding according to the two-electron scheme.
The constancy of the total thermal, enthalpy, and entropy contributions (∼0.30−0.46 eV) was observed when calculating the adsorption energy based only on the electronic energy and the energy of zero-point vibrations.
On the studied catalysts containing cobalt, the following feature of OOH adsorption was observed: adsorbate displacement and a smaller distance between the metal and more distant oxygen, which can prevent the breaking of the O−O bond and the desorption of OH − ( Figure 5).
However, one or more strong adsorptions of intermediates can be detrimental to the ORR process. An overly strong adsorption energy means that the release of free energy at this stage is huge, which can lead to a small release of free energy at other stages, since the overall change in the free energy of the system is fixed. The energies of each of the stages and overpotential are shown in Table 4.
The efficiency of the process is determined by the stage with the least energy release, determining the overpotential on the catalyst. The lower the overpotential, the closer the catalyst properties are to ideal. The catalysts containing cobalt and iron are characterized by the lowest overpotential. It should be noted that it is these catalysts that have become the most frequent subjects for experimental studies. 18,33,35 According to theoretical calculations, 41 cobalt is one of the best metals in terms of catalytic properties in the ORR, surpassed only by iridium and rhodium, which, however, are precious metals, making them less promising for commercial use.
It can be noted that for different catalysts, different stages are rate-limiting; for cobalt, copper, and nickel, this corresponds to the reduction of adsorbed molecular oxygen to the peroxide radical, whereas for iron, this is the process of charging and desorption of the hydroxide ion. This is probably due to the fact that the iron atom in the models is characterized by the lowest positive charge as compared to the other metals (Table 5).
Analyzing the energy profile of the studied reaction and taking into account the external voltage (U = 1.20 V) (Figure 6), it can be seen that the path of the reaction proceeds and greatly changes when the doping metal is varied. The closest to the ideal "zero" line are the catalysts that include cobalt and iron.
An attempt was made to expand the reaction mechanism ( Figure 7) by including the stages of charging (2 → 3, 5 → 6, 7 → 8, and 10 → 11) and association with water (3 → 4 and 8 → 9). The original mechanism corresponds to the transitions 1 → 2 → 5 → 7 → 10 → 12. Here, the stage of association was modeled by adding one water molecule to the model. To simulate the charging stage, the charge of the system was changed to −1. The energy of an individual electron was taken as zero since it was assumed that electrons originate from the cathode and not from the solution and that graphene is a conductor and can transfer charge without significant energy consumption. Figure 7 shows the reaction profile for a Co-containing catalyst. As it was revealed above, for cobalt, a change in the spin state is possible; therefore, both the low-(doublet) and high-spin (quartet) states for the intermediates were calculated. The plot shows the states with the lowest energy: high spin for 5 and 10 and low spin for all others. Modeling the reaction according to the extended mechanism, a strong influence of the system charge on its energy is observed. A negative charge of the system leads to a decrease in its energy, which, in particular, can facilitate the process of oxygen adsorption on the catalyst. The presence of such a charge solves the problem of negative energy of oxygen adsorption on the catalyst. However, in this mechanism, energy barriers appear in the elementary reactions O 2 − ads + H 2 O → OOH ads + OH − (4→5), O − ads + H 2 O → OH ads + OH − (9 → 10), and OH − ads → OH − (11 → 12). These barriers are quite understandable and correspond to the cleavage of the H−O bond and the desorption of the hydroxide ion, respectively. On the other hand, it can be assumed that the barrier value can be reduced if the charging processes take place with partial charges or in parallel with other chemical processes.

CONCLUSIONS
As a result of the study, it was revealed that graphene doped with nitrogen and metals can exhibit high catalytic characteristics in the reaction of electrochemical oxygen reduction and the use of