Oxygen Evolution Reaction Kinetic Barriers on Nitrogen-Doped Carbon Nanotubes

We investigate kinetic barriers for the oxygen evolution reaction (OER) on singly and doubly nitrogen-doped single-walled carbon nanotubes (NCNTs) using the climbing image nudged elastic band method with solvent effects represented by a 45-water-molecule droplet. The studied sites were chosen based on a previous study of the same systems utilizing a thermodynamic model which ignored both solvent effects and kinetic barriers. According to that model, the two studied sites, one on a singly nitrogen-doped CNT and the other on a doubly doped CNT, were approximately equally suitable for OER. For the four-step OER process, however, our reaction barrier calculations showed a clear difference in the rate-determining *OOH formation step between the two systems, with barrier heights differing by more than 0.4 eV. Thus, the simple thermodynamic model may alone be insufficient for identifying optimal OER sites. Of the remaining three reaction steps, the two H2O forming ones were found to be barrierless in all cases. We also performed solvent-free barrier calculations on NCNTs and undoped CNTs. Substantial differences were observed in the energies of the intermediates when the solvent was present. In general, the observed low activation energy barriers for these reactions corroborate both experimental and theoretical findings of the utility of NCNTs for OER catalysis.


INTRODUCTION
With the looming threat of climate change, impending rise in global population, and increasing energy demands, the development of alternative and sustainable energy sources becomes essential. 1−3 Because of the intermittent availability of sources such as solar and wind energy, an efficient storage system is required to realize their potential fully. 1,2 Long-term storage capability, the absence of carbon emissions during conversion, high energy density, and flexibility of use both in mobile and stationary applications make hydrogen an appealing energy vector. 4 Perhaps the most promising candidate for largescale hydrogen production is water splitting through electrolysis, as it enables the clean conversion of water into oxygen and hydrogen gases. 3,5 Unfortunately, the oxygen evolution reaction (OER) occurring at the anode requires high overpotentials and is rather slow because of the relatively large molecular rearrangement involved in the formation of the O 2 bond. 5,6 Consequently, the development of durable highactivity OER catalysts is of paramount importance for the conversion and storage of renewable energy.
Using both experimental and computational approaches, dozens of potential OER catalysts have been investigated to date. Many studies have focused on noble metals and their oxides such as Pd, 7,8 Pt, 9, 10 Au, 11 Ru, 10 RuO 2 , 12 Ir, 10 and IrO 2 , 13 but the scarcity and high cost of these compounds hinder their application. In principle, a cheaper alternative would be to employ transition-metal oxides such as Co 2 O 3 , 14,15 MnO 2 , 16,17 Ni−Fe oxide, 18 and Fe 100−y−z Co y Ni z O x . 19 Although they can show activities that are comparable to their RuO 2 and IrO 2 counterparts, 20 the low conductivity of these oxides presents a major obstacle in practical applications. 21 In an attempt to circumvent the issues commonly associated with noble and transition-metal oxides, recent investigations have focused on carbon-based structures such as the carbon nanotubes (CNTs) discovered by Iijima. 22 Their superb mechanical and electrical properties and large specific surface area make CNTs excellent catalyst supports, 15,21,23−26 even though the pristine tubes themselves are considered to be inactive toward OER. Nitrogen doping of CNTs produces some distinct defects and nitrogen functionalities that can alter the structure of the tube, with the proportion of different defects determined by the fabrication technique and the amount of nitrogen present in doping. 27 For example, graphitic nitrogens, which merely replace carbons in the hexagonal lattice, do not significantly distort the structure of the nanotube. In contrast, pyridinic defects introduce both fiveand nine-ring defects into the nanotube. With larger amounts of nitrogen, pyrrolic five-ring defects start to appear which facilitate the closing of the nanotube through cap formation. 27 Importantly, many of these nitrogen-doping modifications can significantly improve the OER activity of the CNT in alkaline solutions. 28,29 A number of alternative reaction schemes have been proposed for OER depending on pH and the catalytic material. 30 The two main mechanisms for water oxidation catalysis are known as the radical oxo-coupling mechanism (ROC) and the water nucleophilic attack mechanism (WNA). 31 The defining feature of the ROC mechanism is a bimolecular coupling of two oxo-species that have been formed at adjacent surface sites. In contrast, the WNA mechanism is characterized by the presence of three distinct surface intermediates: *OH, *O, and *OOH, where "*" indicates the single attachment surface site. Of these two mechanisms, WNA is thought to be predominant on most metal oxide surfaces and molecular catalysts such as CNTs and graphene. 32 Its exact mechanism depends on pH, but in alkaline medium, it can be written as 33 The net reaction is the formation of two water molecules and one oxygen from four OH − and the transfer of four electrons. Crucially, the intermediates are identical in both alkaline and acidic media.
In recent decades, density functional theory (DFT) has become a promising tool for the design of new electrochemical catalysts by helping to identify and prescreen promising materials and by shedding light on the atomistic mechanisms underlying the experimentally observed chemical processes. 34−36 For example, utilizing a simple thermodynamical model, Li et al. 37 were able to demonstrate a linear relationship between the binding energies of the surface-adsorbed *OOH and *OH species for OER and its reverse, the oxygen reduction reaction (ORR), on N-doped graphene. On the basis of DFT calculations, they predicted theoretical overpotentials comparable to those measured with Pt catalysts, suggesting that Ndoped carbon nanoribbons could be highly functional OER/ ORR catalysts.
In addition to reaction energetics, DFT can also be used to predict kinetic barriers by coupling it to climbing image nudged elastic band (CI-NEB) calculations. In a CI-NEB calculation, the reaction path is explored by a set of spring-connected beads corresponding to adjacent molecular configurations along the reaction coordinate. Contrary to the typical nudged elastic band (NEB) calculation, 38,39 in CI-NEB, these beads can also move in the direction of the reaction path, progressively improving the estimate of the barrier heights by increasing bead density in these regions. 40 For example, in the case of the cathodic watersplitting reaction, that is the hydrogen evolution reaction (HER), Holmberg et al. recently utilized such calculations to show that on nitrogen-doped CNTs, the Volmer−Heyrovsky mechanism dominates over the Volmer−Tafel one with the Heyrovsky step as rate-determining. 41 In this work, we examined the energetics and barrier heights of WNA-type OER on nitrogen-doped CNTs in an alkaline solution through the four-step reaction mechanism (1)−(4). We modeled the solvent by a droplet consisting of 45 water molecules positioned directly above the active nanotube site. By combining DFT and CI-NEB calculations, we studied the reaction energetics and pathways on a semi-quantitative level to gauge the usefulness of different nitrogen-doped CNTs for practical applications. The systems and active sites employed in this study were based on our previous publication 42 where we carried out a thermodynamic analysis of the catalytic activity of various N-doped CNTs.
The remainder of this article is organized as follows: an overview of the systems studied and the computational methods employed is given in section 2. In section 3, we first present and discuss the major findings of our simulations for solvent-free reaction barrier calculations and then for the ones employing our water droplet model. The results are summarized in section 4 with considerations for possible future research directions.

SYSTEMS AND METHODS
We employed spin-polarized DFT in our calculations as incorporated in the Quickstep module 43 of the publicly available CP2K simulation suite. 44 The Quickstep package utilizes the hybrid Gaussian/plane wave method where the Gaussian-type basis is used to represent the wave functions, and the electronic density is described by an auxiliary basis set. As in our previous studies on nitrogen-doped CNTs, 41,45 we combined the Perdew−Burke−Ernzerhof generalized-gradient approximation functional 46 with the Grimme D3-dispersion correction. 47 The Kohn−Sham orbitals were expanded in the molecularly-optimized double-valence polarized basis set (MOLOPT-SR-DZVP) with core electrons represented by the Goedecker−Teter−Hutter pseudopotentials. 48−50 In all calculations, the plane-wave kinetic energy cutoff was set to 600 Ry. The force convergence criterion of the geometry optimizations was left to its default value of 0.023 eV/Å. To alleviate spin contamination issues, 51,52 the diagonalization algorithm 43 with electronic smearing at 300 K was utilized for the NEB calculations.
Our system consisted of a six-unit (14,0)-zigzag CNT initially with 336 carbon atoms and a water droplet of 45 molecules. The simulations were carried out in a 40.0 × 40.0 × 25.6 Å 3 cell with periodic boundary conditions imposed in all three directions. The droplet was carved out from a previously equilibrated 348.15 K NVT pristine CNT calculation with full water coverage 41 using a cutoff distance criterion from the adsorbing atomic site on the nanotube. After this, a 4 ps NVT simulation was performed on the water droplet−CNT system at 348.15 K. To enable the use of a 1 fs timestep, all the water molecules were deuterated for the NVT calculation.
In the nitrogen-doped CNTs, either one or two carbons in the tube were replaced with nitrogens, as depicted in panels C and D of Figure 1. These graphitically N-doped CNTs model a system under low nitrogen conditions with approximately 0− 1% of the dopant present. Sites 2 and 3 in Figure 1 were chosen for investigation based on our previous study of OER thermodynamics on both nitrogen-doped CNTs and graphene, 42 employing the methodology developed by the Nørskov group. 20,33,53 For the singly and doubly N-doped CNTs (NCNT and N 2 CNT), sites 2 and 3, respectively, were predicted to be most active. As pristine single-walled CNTs are known to be poor catalysts for OER, 21 our focus here is on the nitrogen-doped tubes, even though the pristine CNT has been included in our calculations for completeness.
The CI-NEB method 40 was used to find the minimum energy path and to calculate the energy profile of the four-step The Journal of Physical Chemistry C Article OER process for the chosen reaction sites. A set of eight beads was used in each of the NEBs reported here, while some of the calculations in the Supporting Information utilized up to 12. The convergence of the NEB calculations was investigated by testing three different values for the maximum force criterion: 0.023, 0.1, and 0.21 eV/Å. As explained in the Supporting Information, in all tested cases, the changes observed in the activation energies were below 0.05 eV. Following the approach of Tuomi et al., 54 we settled on an intermediate convergence criterion of 0.10 eV/Å for the water droplet calculations, whereas the looser convergence criterion of 0.21 eV/Å was used for the solvent-free NEB calculations to conserve computational resources.
In each OER step of the solvated NEB calculations, a hydrogen atom was removed from a water molecule at least 2.7 Å distant from the active site to yield the reacting OH − in the initial structure and the charge-depleted CNT. The water droplet was shifted so that it was always directly above the active site. For the initial image of the NEB calculations, a core of five to six water molecules surrounding the OH − was first relaxed together with surface atoms up to the two nearest neighbors of the active CNT site. Next, the outer water shell was optimized while keeping the inner shell fixed. In step 1 of the OER, one water molecule neighboring the OH − was frozen throughout the calculation to stop the topological hole defining the OH − group from migrating to the surface. For the final image of the NEB calculations, a shell of 18 water molecules was allowed to relax around the reacting species, while keeping the rest of the droplet fixed to the geometry obtained during the construction of the initial image. This larger amount of free water molecules was required to stabilize the *OOH-species in step (3). The initial guesses for the intermediate images of the NEB calculation were obtained by linear interpolation of atomic positions of the initial image to those of the final image using the open-access ASE package (https://wiki.fysik.dtu.dk/ ase/). Only the possible surface adsorbates, such as *O in step (3), the reacting OH − , and a portion of the CNT corresponding to the two nearest neighbors of the active CNT site, were allowed to move to curtail the computational cost of the NEB calculations. Consequently, the atoms in our systems were effectively categorized into three groups based on their permitted motion during the NEB calculation: atoms whose positions were allowed to vary freely during the simulation; atoms whose positions changed from image to image because of the linear interpolation; and atoms whose positions remained fixed. The vast majority of the atoms in the system were included in the last category.
In addition to the water-containing NEBs discussed above, we also ran simulations on the same sites on solvent-free Ndoped CNTs and on a pristine tube. We also performed a series of preliminary calculations to investigate the effect of different surface sites, various NEB settings, such as the convergence criteria, and different types of approaches for how to fix the water molecules during the NEB. As the suitability of the utilized DFT method for these kinds of systems has been verified by several previous publications, 41,42,45,54 we did not explore other functionals or basis sets. The bulk of these tests are summarized in the Supporting Information (see, e.g., Figures S1 and S2).

RESULTS AND DISCUSSION
The most important initial discovery of this study was that the processes where water is created, namely, reactions (2) and (4), are barrierless. This observation is in qualitative agreement with the findings by Fortunelli et al. that regardless of whether solvent effects are included, the water-forming steps in OER and ORR on a Pt surface possess low barriers. 55 Indeed, this lack of barrier was observed both in solvent-free and solvated NEBs, on several different sites, and in test calculations where the OH − was initially one water molecule away from the reactive site. In this case, the NEB calculation contained a proton-hopping event and the actual OH − reaction with the surface site or some surface-adsorbed species. The barrierlessness manifested itself in that the system relaxed to the final geometry of the NEB process directly during the geometry optimization of the initial image, even in cases where the distance between the approaching OH − and the reacting surface species was large. Because of this lack of barrier in steps (2) and (4), the next subsections will focus solely on steps (1) and (3), namely, the production of *OH and *OOH.
3.1. Solvent-Free CNTs. Because the thermodynamically focused approach of the Nørskov model often employs solventfree surface calculations for estimating catalytic activity, we wanted to investigate barrier heights for these kinds of CNT systems as well. In this case, NEB calculations were performed on pristine and singly and doubly N-doped CNTs. On the pristine tube, all sites are equivalent, whereas on NCNT and N 2 CNT, sites 2 and 3 (see Figure 1) were used, respectively.
The NEB results for the NCNT and N 2 CNT are shown in Figure 2 for step (1). It can be seen that without the solvent, the adsorption of the OH − is a barrierless process. Analogous to the observations made during our HER studies, 41,54 when the OH − adsorbs on the nanotube, one sees a distortion of the tube structure, with the carbon atoms near the reaction site moving toward the adsorbate, reflecting a change from sp 2 to sp 3 hybridization as the coordination of the C-atom active site becomes tetrahedral. This rise of the C-atom active site along the tube normal was 0.44 Å for the pristine CNT, for example.
The greatest adsorption energy was observed for site 2 in NCNT because of its proximity to the nitrogen atom. The adsorption energy decreases from N 2 CNT to CNT, which is in The Journal of Physical Chemistry C Article agreement with the adsorption energy study by Murdachaew and Laasonen, 42 and we see that the energy difference between the neighboring curves is about 0.5 eV in the end. As expected from Figure 2, the C−OH bond has its smallest value of about 1.45 Å in the NCNT case and its largest value of 1.48 Å for the CNT case. Additional bond lengths for this and other NEB calculations are reported in Table S2 of the Supporting Information.
In contrast to step (1), for step (3), a small barrier of about 0.1 eV is observed in all studied cases, as depicted in Figure 3.
This barrier is partly due to the presence of the water molecule formed in step (2), as one sees a twist in the approaching OH during the reaction. Note that in this case, the *O in the initial state had attached to a bridge site instead of on top of a carbon atom (see snapshots of the initial and transitions states in Figure 3). The O−O bond length in *OOH is similar in all three cases and has a value of 1.47−1.48 Å, but there are significant differences in the C−O bond lengths. In N 2 CNT and CNT, this bond is about 1.54 Å, whereas in NCNT it is only 1.47 Å. This difference is also reflected in Figure 3 in the stability of the *OOH-species. For N 2 CNT, image 8 is about −0.9 eV lower in energy than image 1. In contrast, for NCNT, this number is about −1.6 eV with the corresponding value for CNT located between these two. The obtained barrier heights can also be compared with the OER results obtained by Fortunelli et al. through DFT NEB calculations with a continuum solvent model on Pt(111) surfaces. 55 While platinum is not a particularly suitable catalyst for OER, 3,56,57 it is one of the best catalysts for both HER and ORR. 3 In the most favorable reaction mechanism studied for the platinum surface in a vacuum, the calculated kinetic reaction barrier for the rate-determining step (RDS) for *OOH formation from surface-adsorbed *OH and *O was 1.04 eV. 55 As expected, this is substantially higher than the barrier for the RDS of reaction step (3) observed here.
3.2. Solvated CNTs. The solvated NEB results for the first OH − attachment (step (1)) on NCNT and N 2 CNT are shown in Figure 4. In contrast to the solvent-free nanotube results displayed in Figure 2, we see that for both NCNT and N 2 CNT, a small barrier of less than 0.1 eV is observed. As in the solventfree CNT case, site 2 in NCNT seems to bind *OH more strongly than site 3 in N 2 CNT, with the C−O bond length at 1.43 Å compared to 1.49 Å for N 2 CNT. Regardless, the overall shapes of the two curves are very similar. In both cases, the formed *OH intermediate accepts two hydrogen bonds and donates one, resulting in a substantially weaker adsorption than that observed for the solvent-free calculations, as can be seen by comparing Figures 2 and 4.
Unlike in step (1) of Figure 4, there are clear differences in the reaction barriers for the third OH − attachment producing *OOH [step (3)], as seen from Figure 5. The barrier for N 2 CNT is around 0.5 eV, while the same reaction in NCNT is virtually barrierless, at around 0.1 eV. This finding is in line with our thermodynamic calculations on pristine and graphitically N-doped CNT systems, where step (3) was observed to be the RDS. 42 Both of these RDS barriers are, however, relatively low when compared, for example, to the 0.92 eV reaction barrier calculated by Luckling et al. for OER on TiO 2 . 58 In both cases, the product *OOH species donates one and accepts two hydrogen bonds, one for each of its oxygens, and is thus again more weakly bound to the surface relative to the solvent-free CNT. The C−O bond lengths are 1.60 Å for N 2 CNT and 1.51  (1)] calculated using CI-NEB. The green, red, and blue curves refer to the pristine, singly N-doped, and doubly N-doped CNTs, respectively. The two snapshots display the molecular configurations for the initial and final states of the N 2 CNT reaction (i.e., images 1 and 8, respectively).   (1) and (3).
Comparing the solvated calculations with each other, again, *OOH is more strongly bound to the surface in the NCNT case compared to the N 2 CNT one, with a difference of about 0.6 eV between the relative energies of the final images. As mentioned, the C−O bond for the *OOH on NCNT is around 0.1 Å shorter than the N 2 CNT one. This is reasonable because site 2 is the nearest neighbor of N on NCNT, whereas site 3 on N 2 CNT is the second and third nearest neighbor of the two Natoms as shown in Figure 1.
Looking at Figures 3 and 5, it is evident that the solvent can have a deciding effect on the reaction barrier. This finding agrees with previous research by Sha et al. 59 and Fortunelli et al. 55,60 who observed substantial solvent effects in their continuum solvent studies of OER and ORR on platinum surfaces. Similar effects were also reported by Fester et al. in their study of OER and water-assisted water dissociation on cobalt oxide nanoislands. 61 The reaction paths for the solvated NEB calculations of this study have been combined with our energy calculations in Figure 6 to form a summary of the reaction path. The details of how the individual energies have been computed can be found in the Supporting Information together with a summary of the relative initial, transition, and final-state energies in Table S1. As in our OER thermodynamic study of nitrogen-doped CNTs, 42 the reference potential has been set to that of the standard hydrogen electrode. The slightly arbitrary nature of incorporating the bystander OH − and H 2 O species in our solvated calculations makes direct comparison difficult with the typical thermodynamically obtained step plots where the free energy of the whole process sums up to the experimental value of 4.92 eV (see, e.g., Man et al. 20 or Murdachaew and Laasonen 42 ). As in the DFT calculations of OER on SrCoO 3 by Tahini et al., 62 the RDS for the nitrogen-doped CNTs is the formation of *OOH at step (3). As mentioned, this result is also in agreement with our thermodynamic calculations 42 and, furthermore, with the Brønsted−Evans−Polanyi principle of surface chemistry, 63,64 which predicts that in a multistep reaction, the step with the largest thermodynamic energy change will also possess the largest kinetic barrier. On the basis of the differences in the calculated energy paths, our findings indicate that while the thermodynamic approach is useful for initial surveys of many catalysts and sites, subsequent NEB or other reaction barrier calculations are needed to ascertain their functionality further. Our findings here and in our adsorption energy study point to the same direction as the recent experimental results for multiwalled nitrogen-doped CNTs by Davodi et al. 65 that these kinds of systems can exhibit good performance for OER on par with the latest metal-based catalysts. 66−68 As these systems show excellent catalytic activity for HER as well, 65 they serve as promising candidates for bifunctional metal-free catalysts for full water splitting.

SUMMARY AND CONCLUSIONS
In this study, we have looked at the kinetic barriers for the fourstep OER on pristine and nitrogen-doped CNTs, both with and without the solvent. The studied sites were chosen based on results from our previous thermodynamic investigation of these systems. 42 We employed a droplet model containing 45 water molecules to account for solvent effects. We observed that the presence of the solvent had a substantial impact on the reaction barrier. In all of our calculations, the third reaction step, which forms *OOH, was rate-determining, whereas the H 2 O-forming second and fourth steps were found to be barrierless.
Judging from our thermodynamic calculations, 42 the two studied sites on our singly and doubly N-doped CNT systems were approximately equally preferable for OER. However, we observed a clear difference between the reaction barriers at the NCNT and N 2 CNT sites. For the solvated calculations, the barriers of the *OH-producing first step are below 0.1 eV, but for the rate-determining *OOH formation step, the N 2 CNT NEB barrier is 0.54 eV, which is substantially higher than the corresponding 0.12 eV obtained in the NCNT case. On the basis of these results, while the solvent-free thermodynamic model is very useful for the prescreening of potential catalysts, it might not be sufficient for identifying the best systems for OER.  . Calculated reaction energy path for the OER process. The snapshots display the reaction process on solvated NCNT. Note that the number of water molecules in the CNT-water droplet system indicated by * in the figure changes at each step so that the total number of atoms in the system at each step remains the same (see Section S6 of the Supporting Information for more details). The kinetic barriers have been included in the step energies in all cases.
The Journal of Physical Chemistry C Article It should be noted that the conclusions of this article rely solely on energetic considerations and have not taken into account any entropic or zero-point energy effects because of the excessive computational resources required for such calculations. 69,70 We did, however, perform some rudimentary tests to study the impact of applied electrical potential on the reaction barriers by manipulating the total charge of the system. As shown in Table S3 of the Supporting Information, we found that the relative energies of the final NEB states and the approximate barriers decreased with the increasing system charge, as is expected for OER. In addition to a more careful investigation of the impact of overpotential, entropic, and zeropoint energy effects for the systems studied here, it would be interesting to extend the barrier calculations to other promising CNT systems such as tubes containing Stone−Wales-defects. 42 In general, the observed low activation energy barriers further solidify the emerging scientific consensus that nitrogen-doped CNTs can serve as excellent metal-free catalysts for the OER.

Notes
The authors declare no competing financial interest.

■ ACKNOWLEDGMENTS
We would like to thank Finland's Center for Scientific Computing (CSC) for the computational resources provided for this study. We are also grateful to Nico Holmberg for insightful advice and useful discussions throughout this project. Finally, we thank the Academy of Finland project 13292520 for funding.