Effects of Doped N, B, P, and S Atoms on Graphene toward Oxygen Evolution Reactions.

Molecular oxygen and hydrogen can be obtained from the water-splitting process through the electrolysis technique. However, harnessing energy is very challenging in this way due to the involvement of the 4e- reaction pathway, which is associated with a substantial amount of reaction barrier. After the report of the first N-doped graphene acting as an oxygen reduction reaction catalyst, the scientific community set out on exploring more reliable doping materials, better material engineering techniques, and developing computational models to explain the interfacial reactions. In this study, we modeled the graphene surface with four different nonmetal doping atoms N, B, P, and S individually by replacing a carbon atom from one of the graphitic positions. We report the mechanism of the complete catalytic cycle for each of the doped surfaces by the doping atom. The energy barriers for individual steps were explored using the biased first-principles molecular dynamics simulations to overcome the high reaction barrier. We explain the active sites and provide a comparison between the activation energy obtained by the application of two computational methods. Observing the rate-determining step, that is, oxo-oxo bond formation, S-doped graphene is the most effective. In contrast, N-doped graphene seems to be the least useful for oxygen evolution catalysis compared to the undoped graphene surface. B-doped graphene and P-doped graphene have an equivalent impact on the catalytic cycle.

Introduction

With the dawn of technological evolution class="Chemical">n, certainpan> electrochemical reactionpan>s such as pan> class="Chemical">oxygen reduction reaction (ORR),[1] oxygen evolution reaction (OER),[2] and hydrogen evolution reaction[3] can be achieved by electrolysis[4] and photocatalysis[1] to address the energy-related issues. The obstacles in these multielectron transfer pathways are the high overpotential and energy barrier of the rate-determining step (RDS). Though precious metal complexes provide a good catalytic size-performance ratio, they are bound with instability and high cost and prone to poisoning. Later in the 21st century, these metal-based electrocatalysts were modified, assisted, or replaced by newly introduced graphene, a carbon-based 2D material.[5] The sp2 carbon atoms arranged in a honeycomb structure with π–π conjugation, aided by high mechanical strength and electrical conductivity, make graphene the most effective alternate catalyst.[6] Since then, other 2D materials such as graphdiyne,[7−9] graphyne,[10,11] graphitic nitride,[12−14] black phosphorous,[15−17] and boron nitride[18,19] play an essential role in improving metal-free heterogeneous catalyst. Graphene in its bare form is not significantly responsive toward electrocatalysis, though the artificially created distortion leads to improvement in performance. Another instance of enhanced catalytic activity is observed by doping metals and nonmetals that create an electron and spin density distortion on the graphene surface. Nonmetallic atoms such as N,[20,21] B,[22,23] P,[24] and S[25,26] act as doping material. The experimental aspects of these materials deal with engineering methodology,[27] industrial-scale production, abundant starting materials, and excellent electrochemical performance in ambient conditions.[28] However, the theoretical[29−31] investigations with atomistic details focus on active sites, the catalytic mechanism with energetically favorable pathways, and modeling new catalysts using available information. Previous studies shed light on the remarkable effect of singly doping N, B, S, and P atoms on the ORR,[23,32−35] without much detail about the detailed mechanism of OER.[25] The 4e– transfer pathway of OER[36] is a complex and energy-consuming process, and it is observed that even metal-free doped material became unstable after a few cycles.[37−39] Theoretical investigations revealed the charge distortion on the graphene sheet leading to ORR and OER catalytic activities.[30] Thorough inspection indicated that the carbon atom adjacent to the dopant was the active catalytic site in N-doped graphene.[30,35] However, the dopant itself acted as the active site for ORR and OER in the case of boron-,[40] phosphorous-,[41] and sulfur-[25] doped materials. From the previous studies, it was evident that the OER at 0 V is an energetically uphill process[42] and the last step of the catalytic cycle went downhill only after a specific potential application.[43] However, the reaction barrier for the individual steps has not been explored extensively, which can provide an opportunity to explore the role of various dopants on the mechanism and feasibility of reactions. The climbing image nudged elastic band method was used to explore the reaction barrier for the RDS of OER, that is, the formation of an oxo–oxo bond[44] using density functional theory (DFT). Lack of detailed mechanistic study of the catalytic cycle in finding out the activation energy for the individual steps of OER has motivated us to perform the current systematic and explicit study on various single-atom-doped graphene surfaces.

Different precursors, conditions, and preparation techniques lead to different doping sites and consistency. Effectively n class="Chemical">graphene undergoes dopinpan>g at bulk, edge, anpan>d Stonpan>e–Wales defect positionpan>s. Sinpan>ce the very first reported pan> class="Chemical">N-doped graphene, X-ray photoelectron spectroscopy (XPS) data confirm the presence of three different N-doped positions, namely pyridinic, pyrrolic, and graphitic.[45−48] Considering the active site for adsorption of the reactive species, some studies revealed pyridinic nitrogen and others reported graphitic nitrogen to be responsible for effective catalysis for N-doped graphene(NGr).[49−53] Microwave plasma-oriented synthesis of free-standing N-doped graphene shows maximum at. % of graphitic nitrogen.[20] B-doped graphene prepared for Li-ion battery application through green synthesis method shows XPS peak for BC3, otherwise known as graphitic B.[54] XPS analysis of S-doped graphene exfoliated through electrochemical method shows S 2p peak for graphitic S.[55] P is also known to occupy edge as well as bulk position.[56] Thus, it is evident that along with N, B/S/P atoms can also occupy the graphitic position by replacing the carbon atom in bulk graphene. For an appropriate comparative study, we chose all doping sites to be graphitic positions, and the doping was performed on a graphene sheet containing 72 carbon atoms, which leaves us at 71 carbons and 1 dopant. So, for NGr, BGr, SGr, and PGr, the dopant concentration is 1.62, 1.23, 3.62, and 3.50% atomic weight, respectively. The concentration of the dopants may vary depending on the preparative process. Various synthesis methods provide N-doped graphene in the concentration range of 0.4–15.0 wt % in monolayer or multilayer graphene.[20] The successful synthesis of BGr through chemical vapor deposition and later through growth technique on a solid surface leads to 2.5 and 5.0 at. wt %, respectively. In a recent study,[57] P is reported to occupy 6.40 at. wt % on the graphene surface and is significantly higher than the previously reported value of 1.94. S-doped graphene prepared for ORR[58] or bifunctional catalysis[59] was examined with S atomic weight concentration of 2.2 and 0.8%, respectively. So, the method applied to synthesize the doped material plays a vital role in the substituent loading.

In this study, we have compared the effect of single doping of n class="Chemical">N, B, S, and P atoms onpan> a pan> class="Chemical">graphene sheet on OER by replacing graphitic carbon as edge selectivity behaves differently than the bulk selectivity. Previous DFT calculations[60] have shed light on the effect of molecular oxygen and oxygen-containing reactive species on surface adsorption and binding energy. The binding strength of OH– to the surface seems to have an inverse relation with surface coverage.[61] However, in our study, we created an ideal experimental condition, following which we added 1 OH– to make the system 1 M basic medium content. Thus, we ignored the surface oxidation, which might occur from an abundance of reactive species. We explain the active site and will provide a comparative study of the effect of various heteroatoms on the OER in the alkaline medium and activation energy when singly doped on the graphene surface using two different computational methods. The energy barriers were calculated using the metadynamics-based first-principle molecular dynamics (FPMD) method using one of the DFT functionals. For the interfacial reactive systems, the molecular dynamics simulations within DFT functionals are relatively computationally expensive. We also plan to present a comparative finding of the study using a cost-effective density functional tight-binding (DFTB) method, which seems to be an appropriate method for similar systems for future investigations.

Results and Discussion

Each step of OER is involved with OH– as active species, and the first step is the adsorn class="Chemical">ptionpan> of the OH–. As a nucleophile, the OH– will be adsorbed onpan> the positive charge dense regionpan>. Thus, by observinpan>g the radial distributionpan> functionpan>s (RDF) of pan> class="Chemical">oxygen of water molecules to the dopant and the adjacent carbon atoms from classical molecular dynamics (CMD) simulations, we determine the active site for further investigations through FPMD simulations. The doped graphene sheet, along with the water molecules placed along the z-direction, is figuratively depicted in Figure a. The energy-minimized configuration of the doping site and neighboring atoms are depicted in Figure b–e for NGr, BGr, SGr, and PGr, respectively. The RDFs, calculated from CMD simulations, are depicted in Figure f–i for NGr, BGr, SGr, and PGr, respectively. We observe a sharp peak at 3.2 Å, followed by a minimum and a broad peak ranging from 4 to 8 Å for the RDF of N–Ow. The RDFs of Ow with three carbon atoms adjacent to N are broad peaks ranging from 3 to 8 Å. This suggests that all carbon atoms next to the dopant (N) are similar. Since N is more electronegative, hence interacts strongly with the hydrogen atoms of water molecules and effectively shows a simultaneous interaction with oxygen. This observation is along with the phenomenon of charge distortion and excess positive charge accumulation on the carbon atoms next to the nitrogen. Following the facts mentioned above, we choose a random carbon (among C1, C2, and C3) to be the active site. All these carbon atoms show similar behavior toward oxygen atoms of water molecules. For BGr, a p-type dopant being less electronegative than the carbon primarily creates a negatively dense charge environment around the neighboring carbon atoms. Also, from RDF data, we observe a stronger correlation of Ow with the doped material as compared to the adjacent carbon atoms. Bader charge analysis of B-doped graphene has confirmed the increase of charge density in the surrounding carbon atoms and B acquiring extreme positive charge density;[62] also, because of its electron-donating ability (p-type), it can transfer 0.47 electrons to the C on the graphene surface.[63] Thus, B dopant can act as an active site, confirming our finding. Similarly, plotting the Ow RDF concerning the P and adjacent carbon atoms, a peak is observed at 3.5 Å for P–Ow, which seems missing for the adjacent carbon atoms. P atom, owing to its oxygen affinity, has been employed as the active site for reactive species adsorption.[41] S atom having similar electronegativity as carbon is known to distort the spin density and creates a defect. The RDF of S–Ow produces a sharp peak at 3.4 Å. For reference, the Ow–C in pure graphene is provided in Figure S6. The atomic orbital mismatch can induce a positive charge on the S-atom, which becomes capable of OH– adsorption.[64] In N, S co-doped graphene, increasing S atomic percentage increases the OER catalytic efficiency.[65] Thus, S becomes the most probable choice as the active site for all catalytic steps with OH*, O*, and OOH* as the adsorbate. We observed a peculiar trait for the Ow RDF concerning the C atoms adjacent to the dopant when compared with the undoped graphene. For dopants, that overly alter charge density (NGr and BGr), the Ow–C RDF seems to deviate more from the usual behavior in comparison to PGr and SGr; the later ones alter spin density due to orbital mismatch. The RDFs confirm the presence of a defect in the bulk graphene site. The following results are calculated from the FPMD simulations.

Figure 1

(a) Represents the simulation box with a singly nonmetal-doped graphene sheet and 100 water molecules. (b–e) depict the system consisting of NGr, BGr, SGr, and PGr graphene sheets. (f–i) represent the RDF of oxygen of water molecule (Ow) to the N, B, S, and P atom of the graphene sheet. For comparison and detection of the active site, the RDFs of Ow are also calculated to the three adjacent carbon atoms and presented in the same panel as the doped material.

(a) Represents the simulation box with a singly nonn class="Chemical">metal-pan> class="Chemical">doped graphene sheet and 100 water molecules. (b–e) depict the system consisting of NGr, BGr, SGr, and PGr graphene sheets. (f–i) represent the RDF of oxygen of water molecule (Ow) to the N, B, S, and P atom of the graphene sheet. For comparison and detection of the active site, the RDFs of Ow are also calculated to the three adjacent carbon atoms and presented in the same panel as the doped material.

Energetics of Catalytic Cycle

Step 1: The first step of the OER involves adsorn class="Chemical">ptionpan> of OH– as our reactionpan> takes place inpan> the alkalinpan>e medium. OH– is formed by completely removinpan>g onpan>e pan> class="Chemical">hydrogen atom away from the water molecule closest to the chosen active site in the final configuration obtained from the FPMD simulations of the doped graphene and water system. The Na+ was added to maintain system neutrality. For the adsorption of the reactive species, the active site has to come out of the surface. Hence, a change in the dihedral angle is expected. Considering these, we have chosen the coordination number of the active site (Cα, B, S, and P for NGr, BGr, SGr, and PGr, respectively) to the oxygen of OH– as CV1; see Figure b,c. The dihedral angle was designated as CV2 defined by considering the C–C–N–Cα angle for NGr and the C–C–C–X angle for XGr (X = B/S/P) as depicted in Figure a,f, respectively. As the distance between the reactive species and surface varies from system to system, we have to consider different fixed distance cutoff values (d0) while defining CV1. The p and q values for all systems are provided in Table S2. The wall potential to trace the product state properly, along with all simulation parameter details are provided in Table S3. The CV1 ranges from 0 to 1, where the lowest and the highest range corresponds to the nonbonded and chemisorbed state of the OH– on the active site, respectively. The activation barrier is found to be 5.11 ± 0.33, 3.71 ± 0.27, 4.61 ± 0.46, and 3.92 ± 0.43 kcal mol–1 for NGr, BGr, SGr, and PGr, respectively.

Figure 2

(a,f) panels depict the dihedral angle chosen as the 2nd collective variable (CV) of step 1 metadynamics for NGr and BGr, respectively. (b–e) represent the CVs chosen for step (1–4) of NGr-aided OER, respectively. (g–j) pictorially depict the CVs defined for step (1–4) of BGr-aided OER, respectively. For PGr- and SGr-mediated OER, the CVs are the same as BGr by replacing B with S and P.

The energy required for adsorption of OH– on the active site is found to be least for B-doped graphene. Phosphorous, because of its larger size and better electron dispersion ability, shows better performance as compared to nitrogen as a dopant.

The activation barrier of OER step 1 with unn class="Chemical">doped pan> class="Chemical">graphene will be relatively higher than N-doped graphene. Because of electron and spin density defects created by doping foreign material, a reduction in the energy barrier was observed. The reactant, transition state, and product snapshots of step 1 for NGr, BGr, SGr, and PGr are depicted in Figure a–c,e–g,i–k,m–o, respectively. The averaged free-energy surfaces obtained from three independent metadynamics simulations are depicted in Figure d,h,l,p for NGr, BGr, SGr, and PGr, respectively. At the end of the simulations, sodium cation stays well separated from the OH– and does not participate in the reaction.

Figure 3

Represents the snapshots of reactant (R), transition state (TS), product (P), and surface contour plot of the free energy obtained from the metadynamics simulations of step 1, that is, the adsorption of OH–. Panels (a–d), (e–h), (i–l), and (m–p) depict the reactant, transition state, product, and the free-energy surface for N, B, S, and P, respectively. The free-energy values are presented in kcal mol–1.

Represents the snapshots of reactant (R), transition state (TS), product (P), and surface contour plot of the free energy obtained from the metadynamics simulations of step 1, that is, the adsorn class="Chemical">ptionpan> of OH–. Panpan>els (a–d), (e–h), (i–l), anpan>d (m–p) depict the reactanpan>t, tranpan>sitionpan> state, product, anpan>d the free-energy surface for pan> class="Chemical">N, B, S, and P, respectively. The free-energy values are presented in kcal mol–1.

Step 2: This step involves the formation of an oxo complex by proton abstraction from the OH* by another OH–. The OH– was prepared by removing one of the n class="Chemical">hydrogens from the pan> class="Chemical">water molecule forming a hydrogen bond with H of the adsorbed OH*. The CVs CV1 and CV2 represent the coordination number between Oa and Ha and Ob and Ha. CV1 and CV2 are presented in Figure c,h for NGr and BGr, respectively. The CVs defined in Figure h are used for PGr and SGr by replacing B with P and S. The required parameters p and q and the Gaussian parameters are tabulated in Tables S2 and S3. Pictorially, the free-energy contour plots are depicted in Figure a,c,e,g for NGr, BGr, SGr, and PGr, respectively, along with the snapshots of reactant, transition state, and product in the inset figures. The energy barriers thus calculated are 2.66 ± 0.24, 1.79 ± 0.10, 2.56 ± 0.57, and 1.69 ± 0.05 kcal mol–1 for NGr, BGr, SGr, and PGr, respectively. The d0 values chosen and the energy barriers obtained from independent simulations are tabulated in Table S4. All values are less than the energy barrier value for O* formation calculated for undoped graphene, that is, 2.68 ± 0.15 kcal mol–1. The cation maintains neutrality and attains a stable state by forming a hydrogen bond with the water molecules in its vicinity.

Figure 5

(a,b), (c,d), (e,f), and (g,h) represent the contour plots of the free-energy surface obtained from metadynamics simulations of step 2, that is, oxo complex formation on the surface and step 4, that is, molecular oxygen desorption for NGr-, BGr-, SGr-, and PGr-assisted OER catalysis. The snapshots of the reactant (R), transition state (TS), and product (P) are provided in the inset figure, and the arrow marks their respective positions on the contour surface.

Step 3: Both homogeneous and heterogeneous catalysis involved in OER claim the oxo–oxo bond formation reaction as the RDS.[66−68] Here, we have approached the direct mechanism, where the reactive species OH– directly forms a bond with the O*. The OH– was consumed in the previous step; another one was produced by removing the n class="Chemical">hydrogen atom from the pan> class="Chemical">water molecule closest to the O*. CV1 is simply defined as the Oa coordination number concerning Oc. For CV2, the coordination number between Oc and Od is considered, where Od belongs to an external water molecule. Physically CVs for NGr and BGr are defined in Figures d,i, respectively. For PGr and SGr, the CVs can be redefined by replacing B with P and S in Figure i. The position of the external water molecule and O* is fixed for the constant d0 defined for CV1 and CV2. The physical parameters p/q and the K values for wall potential along with metadynamics simulation parameters are given in Tables S2 and S3 of the Supporting Information.

The activation barriers for oxo–oxo bond formation obtained from three metadynamics simulations were calculated, and the values are given in Table S4. CV1 ranges from 0 to 1 as the distance between Oa and Oc, whereas CV2 ranges between 1 and 0 due to an increase in distance between Oc and Od for all systems. The activation barriers are 18.23 ± 0.48, 15.97 ± 0.51, 14.73 ± 0.92, and 17.88 ± 0.76 kcal mol–1 for NGr, BGr, SGr, anpan>d PGr, respectively. As observed, SGr exhibits the least value of energy barrier for the oxo–oxo bonpan>d formatiopan> class="Chemical">n, whereas NGr has the highest. There are very few experimental and theoretical studies concerning the effect of graphene surface solely doped with P or S on the OER catalysis. Sulfur,[25] phosphorous,[69] and boron[62,70] are proven to be excellent dopants on graphene-based materials for bifunctional electrocatalysis. The comparison of these observed values with the energy barrier of the 3rd step of the undoped graphene catalytic cycle, that is, 21.19 ± 0.51 kcal mol–1, unveils the fact that nitrogen is least appropriate for increasing the efficiency of graphene as an OER catalyst. Boron and sulfur atoms provide equivalent data for the energy barrier of hydroperoxo complex formation and seem to be the most useful. The snapshots of reactant, transition state, product, and the surface contour plot obtained from the metadynamics simulations are presented in Figure , panel (a–d), (e–h), (i–l), and (m–p) for NGr-, BGr-, SGr-, and PGr-assisted OER catalytic cycle, respectively.

Figure 4

Represents the snapshots of reactant (R), transition state (TS), and product (P), as well as the surface contour plot of the free energy obtained from the metadynamics simulations of step 3, that is, desorption of molecular oxygen from the graphene surface. Panel (a–d), (e–h), (i–l), and (m–p) depict the reactant, transition state, product, and the free-energy surface for N, B, S, and P, respectively. The unit of free energy depicted is kcal mol–1.

Represents the snapshots of reactant (R), transition state (TS), and product (P), as well as the surface contour plot of the free energy obtained from the metadynamics simulations of step 3, that is, desorn class="Chemical">ptionpan> of molecular pan> class="Chemical">oxygen from the graphene surface. Panel (a–d), (e–h), (i–l), and (m–p) depict the reactant, transition state, product, and the free-energy surface for N, B, S, and P, respectively. The unit of free energy depicted is kcal mol–1.

Step 4: The final step of the 4e– transfer pathway is the abstraction of a proton from the n class="Chemical">hydroperoxo complex (pan> class="Chemical">OOH*) adsorbed on the graphene surface along with the desorption of molecular oxygen. The final structure obtained from step 3 is adopted to produce the initial geometry for step 4. The closest and properly oriented water molecule toward the hydrogen atom of OOH* was traced, and one hydrogen atom was removed to produce reactive OH–. The three different nuclear coordinates for independent metadynamics simulations are obtained by altering the distance between the hydrogen of OOH* and OH–. CV1 and CV2 were defined by considering the coordination number of the active sites with Oa and Hc with Oe, respectively, depicted in Figures e,j for NGr and BGr, respectively. The same CVs as BGr can be used for SGr and PGr. The necessary parameters defined concerning CV1 and CV2, along with the electronic parameters for running metadynamics simulations, are tabulated in Tables S2 and S3, respectively. The CV1 values range from 1 to 0, representing the covalently bonded and completely desorbed state, respectively, between the active site and Oa. CV2 ranges between 0 and 1, depicting the progress of proton transfer reaction. We plotted the free-energy surfaces and calculated the activation barrier for oxygen departure. The values hence obtained for NGr, BGr, SGr, and PGr are 2.51 ± 0.21, 2.42 ± 0.16, 3.40 ± 0.93, and 2.74 ± 0.41 kcal mol–1, respectively, which are lower than 4.89 ± 0.63 for undoped graphene. In Table S4, we assimilated the do values and corresponding energy barriers for step 4. The barriers for the 2nd and 4th steps are significantly close, hinting that the proton transfer event contributes mostly toward the activation barrier of the 4th step. The desorption of molecular oxygen has a very low barrier or proceeds barrier less, as explained by previous studies.[71,72] The reactant, transition state, and product snapshots along with the free-energy contour plot are provided in Figure b,d,f,h for NGr-, BGr-, SGr-, and PGr-assisted OER catalytic cycle, respectively.

(a,b), (c,d), (e,f), and (g,h) represent the contour plots of the free-energy surface obtained from metadynamics simulations of step 2, that is, oxo complex formation on the surface and step 4, that is, molecular n class="Chemical">oxygen desorpan> class="Chemical">ption for NGr-, BGr-, SGr-, and PGr-assisted OER catalysis. The snapshots of the reactant (R), transition state (TS), and product (P) are provided in the inset figure, and the arrow marks their respective positions on the contour surface.

Relative Energy Profile and Charge Density of Reactive Species

The 2D energy profiles for n class="Chemical">NGr, BGr, SGr, and PGr were calculated and are plotted in Figure a–d, respectively. For plotting all free-energy values simultaneously, we report relative free energy. The free energy (ΔG) of the reactant of step 1 is considered 0 kcal mol–1, and the rest of the species in the intermediate steps were scaled up accordingly. For the same species, the values of ΔG for reactants for different steps vary because of the energy biasing. While plotting the 2D profile by scaling the ΔG values, we calculated the energy barriers. Mulliken and Lowdin charge analyses were performed regarding the locality of the extra electron produced in each step. The trends obtained for the change in charge of all n class="Chemical">dopants are provided in Figure S5.

Figure 6

2D energy profile of the energy barriers of the 4e– reaction pathway obtained from metadynamics simulations for OER assisted by (a) NGr, (b) BGr, (c) SGr, and (d) PGr.

2D energy profile of the energy barriers of the 4e– reaction pathway obtained from metadynamics simulations for OER assisted by (a) n class="Chemical">NGr, (b) BGr, (c) SGr, and (d) PGr.

We followed the same simulation protocol for SCC-n class="Chemical">DFTB as Perdew–Burke–Ernzerhof (PBE). All the simulationpan> details such as p, q specified for CVs, K for wall potential, anpan>d electronpan>ic parameters are tabulated inpan> Tables S2 anpan>d S3. Different d0 values alonpan>g with the activationpan> barriers obtainpan>ed from inpan>dependent metadynamics simulationpan>s for different steps are tabulated inpan> Table S4. The results obtainpan>ed from SCC-pan> class="Chemical">DFTB-based metadynamics simulation of OER steps for all dopants are provided in the Supporting Information, Section 4. The ΔG values hence obtained are comparable to the PBE-D3 method. For all dopants, the 1st and 3rd steps are energy-consuming, with the 3rd being the RDS. We observe the free-energy barrier for SGr to be 14.94 ± 0.41 kcal mol–1, which is in close agreement with that from PBE, that is, 14.73 ± 0.92 kcal mol–1. S-doped graphene seems to be the most effective as a heterogeneous catalyst for both cases, followed by BGr. The activation barrier for oxo–oxo bond formation obtained from PBE and SCC-DFTB are 15.22 ± 0.87 and 15.97 ± 0.51 kcal mol–1, respectively. The RDS barrier for PGr from both levels of theory also mildly deviates with values 17.34 ± 1.06 and 17.88 ± 0.76 kcal mol–1. NGr is least efficient in enhancing the catalytic ability of undoped graphene as the RDS has an activation barrier of 17.91 ± 0.56 kcal mol–1 from SCC-DFTB as compared to 18.23 ± 0.48 kcal mol–1 obtained from PBE simulations. From our observations of the energy barrier of oxo–oxo bond formation, we conclude that boron doping, which increases the negative charge density, and sulfur doping, which alters spin density on the adjacent carbons, are relatively useful as doping materials, followed by P and N. The 2nd and 4th steps, which involve proton transfer, are low energy demanding steps, with activation barriers for 2nd step ranging from 1.80 to 4.00 kcal mol–1. The 4th step has values within 2.00–4.00 kcal mol–1 range for all dopants irrespective of the DFT-MD methods. The energy barrier for the 4th step is the combined effect of oxygen desorption and proton transfer from OOH* to OH–. As depicted in the earlier thermodynamic barrier studies, molecular oxygen desorption is a barrier-less process as found from our kinetic studies. The activation barrier for the 4th step is primarily due to the proton transfer. The results obtained for DFTB are in line with PBE except for the 2nd step of NGr and PGr. We observe that for NGr and PGr, the 4th step requires even less energy than the 2nd step. The pictorial representation of reactant (R), transition state (TS), product (P), and surface contour plots for steps 1–4 and 2D energy profile are depicted in Figures S1–S4 for NGr, BGr, SGr, and PGr, respectively. The complete 2D energy profile with scaled energy barrier values is depicted in Figures S1q, S2q, S3q, and S4q for NGr, BGr, SGr, and PGr, respectively.

We extracted the charge for the structures of n class="Chemical">dopant anpan>d the reactive species obtainpan>ed from metadynamics simulationpan>s. We have presented the chanpan>ge inpan> the Mulliken anpan>d Löwdinpan> charge distributionpan>s for pan> class="Chemical">NGr, BGr, SGr, and PGr in Figure S5(a1–a4),(c1–c4),(e1–e4),(g1–g4) and S5(b1–b4),(d1–d4),(f1–f4),(h1–h4), respectively. For each step, we observe that the oxygen of the reactive species, that is, OH–, eventually attains positive charge density when it approaches the product state as it acts as a nucleophile and shares its electron to form the chemisorbed species, the oxo complex, or the hydroperoxide complex or for molecular oxygen formation. Another significant change is observed for the doping site for which the change in the charge shows a linear trend. This proves the fact of homogeneous charge distribution on the graphene sheet. There is no accumulation of electron density on the reactive species. The continuous rise and fall of the charge values in the case of the reactive species oxygen can be attributed to the nature of metadynamics simulations as the reaction occurs multiple times in a course of the total number of steps. The same trend is observed for all steps irrespective of the dopant.

Conclusions

Because of a lack of research studies of singly n class="Chemical">doped pan> class="Chemical">graphene as a heterogeneous catalyst, its effect on catalysis is not well explored. We found very few instances of solely SGr and PGr used as OER catalyst. So, here we have provided an energetic and dynamic study using metadynamics simulation technique, employing PBE and SCC-DFTB level of theories, for the complete 4e– transfer pathway of OER using singly doped graphene with 4 most studied nonmetal dopant atoms. The active site for adsorption is the C atoms adjacent to nitrogen for NGr and the dopant itself for the rest and was confirmed from the Ow and sheet elements RDF. The free-energy change for NGr, BGr, PGr, and SGr are close in value. The oxo–oxo bond formation leading to the hydroperoxo complex in the 3rd step is the RDS. In terms of catalytic efficiency, S-doped graphene is followed by BGr and PGr. PGr despite being an n-type dopant shows better performance as OER catalysis than NGr due to the size effect. Sulfur, boron, and phosphorous outperform nitrogen as a doping material. We observed enhanced catalytic performance of graphene due to doping foreign material. Also, the charge analysis of the species involved in the catalysis hints at the distribution of the extra electron produced in each step rather than accumulation on the doping site. The step-wise activation barriers obtained utilizing SCC-DFTB align with the PBE generated ones, with a deviation of less than 1 kcal mol–1 in general. It turns out that DFTB may be one of the reliable levels of theory for reactive interaction of water molecules on the surface and can be used as an alternate for computationally expensive DFT functionals.

Computational Methods

The initial structure for FPMD was obtained from n class="Disease">CMD simulationpan>s usinpan>g the large-scale atomic molecular massively parallel simulator.[73] To perceive various molecular inpan>teractionpan>s, adapan> class="Chemical">ptive intermolecular reactive empirical bond order,[74] SPC/E,[75] and general amber force field[76,77] parameters were used. The details of the CMD simulations are provided in the Supporting Information along with nonbonded interaction parameters between sheet elements and water molecules tabulated in Table S1. After obtaining the energy-minimized stable doped-graphene water system, we performed FPMD simulation using the Quick step[78] module available within the CP2K software[79] suite. Each system was simulated within the NVT ensemble for 2 ps with a timestep of 0.5 fs using the Nosé–Hoover thermostat.[80,81] We used PBE[82] functional with Grimme’s 3rd order dispersion correction[83,84] and 600 Ry plane wave cutoff. For explaining core electrons, Goedecker–Teter–Hutter pseudopotentials were used.[85] Hence, the obtained equilibrated final structure from FPMD simulations was considered for the calculations of free energy by the metadynamics method.[86−88] The details of metadynamics simulations are depicted in the Supporting Information. For various steps of the catalytic cycle, we considered the 4e– reaction pathway for OER, which is defined as follows

Metadynamics simulations were performed for each of the above steps individually by defining two CVs. The calculation for finding the energy barrier was also performed using n class="Chemical">DFTB,[89] a cost-effective method as compared to the DFT method. The parameters associated with pan> class="Chemical">DFTB simulations are detailed in the Supporting Information. To check the convergence of our adopted metadynamics method, we performed three independent simulations for each step, where we varied the initial coordinates of the reactive species, hence altering the fixed cutoff distance (d0) to the reproducibility of the product state. The d0 values along with the obtained activation barriers are tabulated in Table S4. The snapshots of the simulated systems were prepared using VMD.[90]