Sensitivity of Density Functional Theory Methodology for Oxygen

Article pubs.acs.org/JPCC

Sensitivity of Density Functional Theory Methodology for Oxygen Reduction Reaction Predictions on Fe−N4‑Containing Graphitic Clusters Joshua P. McClure,*,† Oleg Borodin,† Marco Olguin,† Deryn Chu,† and Peter S. Fedkiw‡ †

Electrochemistry Branch, U.S. Army Research Laboratory, 2800 Powder Mill Road, Adelphi, Maryland 20783, United States, Department of Chemical and Biomolecular Engineering, North Carolina State University, Raleigh, North Carolina 27695-7905, United States

‡

S Supporting Information *

ABSTRACT: Density functional theory (DFT) was used to examine the O2 reduction reaction on Fe−N4-containing graphitic carbon clusters (Fe−N4−G) modeled after recent experimentally identified active sites, Mössbauer spin-state predictions and electrochemical reaction behavior in alkaline media. A detailed analysis of the O2, O, H2O, OOH, and OH adsorbate interactions on the Fe−N4−G cluster with solvation and/or dispersion corrections are considered. The total and partial density of states for the α- and β-spin orbitals are compared for the adsorbate of interest, Fe atom and surrounding graphitic cluster. Relative free-energy diagrams are constructed, which allow us to compare DFT predictions to experimental results for O2 reduction on systems containing embedded Fe−N4 clusters. For all reaction steps, different DFT functionals are explored and the respective geometries, energetics, and spin-states for each adsorbate interaction are reported for six commonly used functionals including B3LYP, M062X, M06-L, PBE, TPSSh, and ωB97X-D. Functionals with high fractions of exact exchange were found to favor higher spin-states, as well as stronger binding of the adsorbates, making these methodologies less feasible for Fe−N4-containing electrocatalysts when compared to experimental data. Pure functionals with and without empirical correlation exhibit different ground spin-states and geometries, however the free energy diagrams yield similar conclusions at relevant overpotentials. The activation energy for the O−OH bond scission step, as well as OH desorption from the Fe−N4−G cluster are discussed since the barrier could prohibit a pure 4e− ORR. Finally, we discuss the energetically unfavorable steps for select overpotentials, which provides the experimentalist with a tuning knob for electrocatalytic design.

■

INTRODUCTION Non-Pt-group metal (non-PGM) electrocatalysts have gained considerable attention for the electro-reduction of O2 due to the potential to decrease the costs of power and energy devices. For example, non-PGM electrocatalysts have gained interest for use in anion-exchange membrane (AEMFC), hybrid acidalkaline membrane (HMFC) and bicell membrane fuel cells (BCMFC).1−3 There are many opportunities for the use of non-PGM catalysts in applications other than fuel cell devices, however non-PGM electrocatalysts are commonly used as oxygen electro-reduction catalysts. Various non-PGM electrocatalysts have been considered for the oxygen reduction reaction (ORR) and include metal-containing macrocycles (i.e., metals such as Fe, Co, Mn, Ni, etc.),4 heat-treated metalcontaining macrocycles,3,5,6 and metal-free nitrogen-doped carbons,7−9 among others. The heat-treated metal-containing materials are particularly interesting since it has been observed experimentally that these materials maintain their electrocatalytic activity for long durations necessary for device operation.10 However, the drawback to heat-treated metal© 2016 American Chemical Society

containing electrocatalysts is a lack of understanding of the active site(s) formed during the aggressive synthesis. Moreover, progress in developing and optimizing new non-PGM electrocatalysts is hindered by the lack of fundamental understanding of the mechanism of the electro-reduction of O2 on these materials. In a number of studies, Fe-containing macrocycles are decomposed on or throughout, for example, carbonaceous nanomaterials.5,11 Alternatively, various transition metal precursors may be added to coordination polymers for complexation which is then followed by heat-treatment.12,13 Despite the synthesis method or approach, the resulting formation products and the active site(s) are only potentially related to metal species.14−18 For example, previous reports suggest metal-free nitrogen-doped electrocatalysts have been shown to promote the 3−4e − reduction of O 2 in alkaline media both Received: August 23, 2016 Revised: November 28, 2016 Published: November 29, 2016 28545

DOI: 10.1021/acs.jpcc.6b08498 J. Phys. Chem. C 2016, 120, 28545−28562

Article

The Journal of Physical Chemistry C experimentally and theoretically.7,8,19,20 Recent, and previous, experimental efforts support the generation of an ORR-active site related to Fe-coordinated moieties. For example, Yeager et al. reported enhanced electrocatalytic activities for catalysts containing heat-treated metal salts in the presence of a coordination polymer.21 Also, Ramaswamy et al.22 used a Δμ spectral subtraction technique and suggested that the active site(s) for heat-treated FeTPP supported on carbon black were likely to involve Fe−N4-like coordination with the Fe atom in a centrosymmetric carbonaceous environment. Recently, an experimental method was reported for preparing Fe−N4containing electrocatalysts with specificity for four-coordinate iron.23 The authors found that after acid-washing there were three dominant spin-states for the Fe site, which were likely due to local environmental alterations resulting from the method of preparation. A complementary 57Fe Mossbauer spectroscopy study revealed up to five different Fe-species after pyrolyzing iron acetate adsorbed on carbon black.15 The authors attributed three of the five Fe species to electrocatalytic activity and assigned these to Fe−N4-like sites with the Fe atoms in low, intermediate and high spin states. In addition to fourcoordinate Fe species, Lefevre et al.24 reported the presence of two active sites after heat-treatment (i.e., Fe−N2 and Fe−N4 sites), and proposed the Fe−N2 sites were more ORR active than Fe−N4 sites. However, in more recent experimental15,17 and theoretical25,26 studies, it was reported that Fe−N2 sites are less chemically stable than Fe−N4 sites, which suggest that Fe− N4 species embedded in carbonaceous materials are reasonable active moieties for consideration under prolonged operation. The identification of active sites related to heat-treated nonPGM electrocatalysts and their mechanistic activity related to the ORR is a flourishing area of research, which requires complementary theoretical methods to develop a complete picture. Density functional theory (DFT) is a relatively inexpensive method for probing the ORR on catalytic systems (i.e., > 50 atoms) compared to, for example, the complete active space multiconfiguration self-consistent field methods that are computationally prohibitive simulations for our system because they scale ∼ N7 (i.e., where N is the number of electrons or basis functions).27,28 Systems involving open-shell electronic structures with multireference character such as transitionmetal complexes, clusters and extended systems are challenging for DFT studies, and the construction of a nonempirical/ empirical DFT methodology that accurately describes openshell systems is considered a highly important but challenging undertaking in theoretical chemistry.27,29,30 With that in mind, DFT has been widely used in studying one or more steps believed to take place in the ORR. For example, the ORR on Co−Nx species embedded in a graphene sheet,31 Co- and Fecontaining polyaniline composites,32 nitrogen-doped graphene,19,20 Fe−Nx species embedded within a carbon framework,25,26 other transition metal-containing species embedded within a graphene sheet, 25 and Fe- or Co-containing phthalocyanines have been considered.34−36 Most of the previous studies limited their investigations of the ORR to one particular DFT methodology leaving the question of the reliability of these predictions and sensitivity of the results to DFT functional largely unanswered.25,26,33,37 It is important to note that multireference and multiconfiguration ab initio methods on heme-related systems with tetra-, penta- and hexa-coordinate species have been previously explored in terms of the appropriate employed active space.28 Specifically, it was

found that the double-shell effect (i.e., 3d orbitals of Fe reinforced with a shell of 4d orbitals) for heme-related systems is important in describing the static and dynamic correlation of d orbitals with respect to spin-state energetics of Fe-containing systems. Considering previously reported experiments and ab initio investigations of heme-related systems, we explored various DFT methodologies for graphitic Fe−N4 containing clusters with and without implicit solvation for various thermodynamically feasible electron transfer steps, as well as the energy barrier of the O−O bond scission step. The investigation of different DFT methodologies for the ORR on Fe−N4 graphitic clusters will start filling this knowledge gap by considering the sensitivity of predicted spin-state ordering, adsorbate/substrate geometries, thermodynamic energies and barriers that are important for understanding the ORR mechanism with respect to the choice of an empirical or nonempirical DFT functional. As a contribution, we explored the ORR on non-PGM electrocatalysts with Fe−N4-containing carbon clusters (denoted as Fe−N4−G). Different empirical and nonempirical DFT methods are used to predict adsorbate interactions on Fe−N4-containing aromatic carbon clusters. For each methodology, we determined the ground spin-states before and after the adsorption of H2O, O2, O, OH, and OOH species on the Fe active site. Free energy (ΔG) reaction diagrams are developed from the adsorbate binding energies at different applied potentials (U) and are compared to reported experimental data. We further discuss the barrier for O−OH bond scission for various DFT methods with and without solvation. We discuss the predicted ground spin-states and degeneracies for each method and compare the total and partial density of states for the α- and β-spin orbitals. In addition, we compare our findings to literature relevant graphitic periodic systems containing embedded Fe−N4 species and experimentally relevant systems with electrocatalysts that possess Fe− N4−G moieties.

■

EXPERIMENTAL SECTION Density Functional Theory Calculation Methodology. Density functional theory calculations were performed using the Gaussian 09 package.38 Six different functionals as standalone or in combination form were considered and are as follows: two different functionals from Truhlar’s M06 family including the meta-GGA M06-L39 (0% HF exchange) and the global-hybrid meta-GGA M06-2X40 (54% HF exchange), the combination form GGA PBE for exchange and correlation41 (0% HF exchange), the hybrid GGA B3LYP42,43 (20% HF exchange), the long-range-corrected GGA + MM ωB97X-D44 functional (22% HF exchange short-range and 100% at longrange) and the global-hybrid meta-GGA TPSSh45 (10% HF exchange). (Note: GGA, HF and MM refer to generalized gradient approximation, Hartree−Fock and molecular mechanics, respectively) Geometry optimizations and vibrational analyses were performed for the Fe−N4-containing graphitic carbon clusters with the following chemical formula C36H16FeN4 (i.e., Fe−N4−G) that are shown in Figure 1. In all studies, the edge carbon atoms of the graphitic cluster were terminated with hydrogens.46 The overall system charge was set to +2 for simulations before and after placing an adsorbate molecule near the metal site following previous work.47 The Fe−N4-containing graphitic carbon cluster with and without adsorbate (ads) molecules (i.e., O2, O, OH, OOH, and H2O) were optimized at the 6-31G(d) basis set48 level using a 28546


Article

The Journal of Physical Chemistry C

the multiplicities m = 2, 4, 6 were considered. The stability of ground-state wave function was verified for each method and ads system with respect to spin-state. Activation energies for O−OH dissociation were performed using the intrinsic reaction coordinate (IRC) calculations for select DFT methods. A threemolecule specification (i.e., reactants, products, and guess for transition state (TS)) was specified, and a search for the TS was performed using the synchronous transit-guided quasi-newton method. The initial guess for the TS was found by systematically scanning over the O−O bond followed by reoptimizing the geometry. Charge-transfer and delocalization between the ads and Fe active site in the Fe−N4−G cluster was evaluated using an orbital-by-orbital analysis with compositions of molecular orbitals, overlap populations between molecular fragments, bond orders and density-of-states spectra calculated using the AOMix program after performing a full population analysis with Gaussian 09.55,56 The total and partial density of states (i.e., TDOS and PDOS, respectively) of α- and β-spin orbitals were further investigated by partitioning the Fe−N4−G system and ads into fragments with Gaussian peak shapes. The PDOS and overlap population density of states analysis (OPDOS) were determined for the respective orbitals summed into each atom or fragment of interest and using the Mulliken population analysis (MPA) to determine the orbital character. Since the Mulliken population analysis is basis set dependent, we included additional TDOS and PDOS with the TZVP basis set as shown in Figures S3 and S4 for the isolated Fe−N4−G cluster, which for the energy levels considered did not produce any significant difference.

Figure 1. DFT images of representative optimized geometries for the (a) isolated Fe−N4−G and (b) H2O−Fe−N4−G in step 1, (c) SO− O2−Fe−N4−G in step 2a, (d) EO−O2−Fe−N4−G in step 2b, (e) HOO−Fe−N4−G in step 3, (f) O−Fe−N4−G in step 4, and (g) HO− Fe−N4−G in step 5. Reaction step 6 is the regeneration of the isolated Fe−N4−G system after desorption of OH-bound species. Note: EO and SO refer to end-on and side-on, respectively.

quadratically convergent self-consistent field procedure49 with subsequent vibrational analysis in order to confirm that no imaginary frequencies were present in the optimized structure. Additional optimizations and vibrational analysis were performed on the isolated ads molecules at the same level of theory. Single-point (SP) energy calculations with the larger triple-ζ quality basis set 6-311+G(2df,2p) were performed using geometries optimized with the smaller 6-31G(d) basis set. A select set of SP calculations with the TZVP basis set was performed (Figure S1) to further probe sensitivity of DFT predictions to the choice of basis set.50 The influence of surrounding solvent on the ORR was investigated by comparing SP energies in the gas-phase, i.e., no implicit solvent, with the SP energies for the cluster surrounded by the solvent cavity reaction fields (SCRF) using the polarized continuum model (PCM) and a dielectric constant (ε) of 78.3553 to approximate H2O solvation.51 Grimme’s empirical dispersion (GD3) correction was used in conjunction with the B3LYP, M06-L, M06-2X and PBE functionals to improve the description of van-der Waals interactions.52 (Note: TPSSh and ωB97X-D functionals have built-in empirical dispersion) Basis set superposition errors (BSSE) were determined using the counterpoise (CP) method.53 For BSSE simulations, the optimized bound ads-Fe−N4−G bound clusters were separated into two fragments; fragment one, the isolated Fe−N4−G with the lowest energy spin state; and fragment two, the ads molecule. The adsorption energies at the 6-311+G(2df,2p) level showed little difference before and after the BSSE correction as shown in Figures S2a,b and the representative O2bound system. All systems were open-shell, thus we studied different spin states (S) for the isolated Fe-containing system, as well as different S for the clusters after placing the ads ∼2 Å from the Fe active site with subsequent geometry optimization. Following previous calculations for O2 bonded to nonheme groups, we examined energies of the different possible states.28,54 Herein, multiplicity (m) is adopted instead of S, where m = 2S + 1, where m = 1, 2, 3, 4, 5, 6, and 7 are interchangeably referred to as singlet, doublet, triplet, quartet, quintet, sextet, and septet, respectively. For the isolated Fe− N4−G cluster, the multiplicities considered were m = 1, 3, 5. After placing the O2, O or H2O molecule near the Fe site within the Fe−N4−G cluster, the multiplicities considered were m = 1, 3, 5, 7. After adsorption of OH and OOH on Fe−N4−G,

■

RESULTS AND DISCUSSION Proposed Mechanism of the O2 Reduction Reaction in Alkaline Media. The reduction of O2 occurs at the cathode compartment during AEMFC, HMFC, or BCMFC operation on the Fe−N4−G electrocatalyst. The isolated Fe−N4−G graphitic cluster is modeled as shown in Figure 1a (isolated). The chosen chemical formula with a Fe:N:C atomic ratio of 1:4:36 is similar to previously reported experimental data,57 and a recent study that used low-temperature scanning tunneling microscopy to image a 2 nm × 2 nm area with a Fe−N4−G containing graphene sheet.58 The ORR may undergo a number of different electron-transfer and chemical reaction steps.59,60 The broad range of non-PGM electrocatalysts suitable for the ORR in high pH has been explained by outer-sphere and innersphere reaction mechanisms, whereby the latter mechanism involves the enhancement and activation of the peroxide (HO2−) intermediate which may lead to completely breaking the O2 molecule (i.e., 4e− transfer).1,60 The inner-sphere ORR mechanism suggests direct adsorption of the O2 molecule to the active site(s) followed by subsequent reaction steps that may or may not involve O2 bond breakage with a decrease in overpotential for an initial outer-sphere O2/O2− step by ∼0.7 V. The exact mechanism of the ORR on non-PGM electrocatalysts is not known, however general descriptors and species involved for the ORR provide insight into the common reaction limitations and barriers.61−63 A potential reaction scheme for the ORR in alkaline media has been previously discussed31 and is shown below according to the following reaction equations/ steps with modifications: 28547


Article


end-on and side-on configurations throughout the study as EO and SO, respectively. Next, we systematically examine the ORR for the six DFT functionals and discuss differences for all reaction steps 1−6 as illustrated in Figure 1. Isolated Fe−N4−G Aromatic Carbon Cluster. Figure 1a shows a representative optimized isolated Fe−N4−G model clusters in a top- and side-view. The Fe−N4−G clusters were determined to be ground state triplets for all methods except the M06-2X functional with a quintet ground state. In the triplet ground-state, all structures adopt a flat morphology, as opposed to that seen for the quintet state which shows distortion of the Fe atom out-of-plane (Note: not shown). Therefore, the isolated Fe−N4−G cluster is in-plane before adsorption of an ads molecule except in the M06-2X case, whereby the Fe−N bond distance is increased due to the higher ground-state spin prediction. Figure 2 shows that the Fe−N

H 2O (or OH−) + *Fe−N4−G → G−N4−Fe−H 2O (or OH−)

(1)

O2 + 2H 2O + G−N4−Fe−H 2O → G−N4−Fe−O2 + 3H 2O

(2)

G−N4−Fe−O2 + 3H 2O + 4e− → G−N4−Fe−OOH + OH− + 2H 2O + 3e−

(3)

G−N4−Fe − OOH + OH− + 2H 2O + 3e− → G−N4−Fe−O + 2(OH−) + 2H 2O + 2e− or G−N4−Fe−OOH + H 2O + e− → G−N4−Fe−H 2O + OOH−

(4)

G−N4−Fe−O + 2(OH−) + 2H 2O + 2e− → G−N4−Fe−OH + 3(OH−) + H 2O + e−

(5)

G−N4−Fe−OH + 3(OH−) + H 2O + e− → 4(OH−) + H 2O + *Fe−N4−G

(6)

where the adsorption takes place on the Fe active site within the Fe−N4−G. The first reaction, step 1, shows adsorption of H2O (or OH) (Figure 1b) on a bare Fe−N4−G followed by displacement of the preadsorbed H2O (or OH) molecule by an O2 adsorbate in reaction steps 2a or 2b (i.e., Figure 1, parts c and d, respectively).64,65 (Note: step 1 may involve a preadsorbed OH-bound group instead of H2O since a pH = 13 is considered.) Protonation of the adsorbed O2 creates an intermediate OOH species in reaction step 3 (Figure 1e) which may or may not lead to O−O breakage. If the O−O bond ruptures this results in an adsorbed O radical in step 4 (Figure 1f) followed by protonation to a bound OH step 5 (Figure 1g). Alternatively, reaction step 3 may preferentially undergo a 2e− transfer reaction and desorb from the active site to form a soluble HO2− species (i.e., peroxide species for pH > 12).60 The successful removal of the OH-bound species in step 6 completes the overall reduction of O 2 to OH − with regeneration of the Fe−N4−G active site, and is referred to as the complete 4e− transfer reaction. Many experimental efforts for non-PGM electrocatalysts in alkaline media have reported the number of electrons transferred during the electroreduction reaction to be between 2 and 4 with Fe−Nxcontaining systems mostly possessing between 3 and 4.17,66,67 For the present study, we consider the ORR mechanism shown for steps 1−6 and evaluate the electron transfer behavior. Other potential electron transfer steps have been previously discussed, however the experimental determination of 3−4 electrons implies activation of the HOO adsorbate and subsequent protonation and electron-transfer steps.26,36 The steps here-in are general descriptors and provide a foundation for evaluating the different DFT methods and describing ground spin-states in terms of α- and β-spin orbitals. It is worth noting that in addition to the above reaction steps, the O2 and OOH adsorbate reactions steps were optimized for two different starting configurations; end-on (Pauling model) or side-on (Griffiths model) configurations to assess initial bonding orientation effects on the reaction energies.21 We refer to the

Figure 2. Longest Fe−N bond distances extracted from the optimized geometries for each reaction step adsorbate (HOH, O2−SO, O2−EO, OOH, O, OH) for all DFT functionals studied. Note: SO and EO refer to side-on and end-on, respectively.

bond distance is ∼1.9−1.93 Å for M06-L, B3LYP, PBE, TPSSh, and wB97X-D, whereas the Fe−N bond distance is ∼2.0 Å for the M06-2X method. Moreover, for DFT functionals with HF ≤ 10% the Fe−N bond distance is ≤∼1.91 Å, whereas the Fe− N bond distance increases with HF > 10%. The predicted ground-state triplet for the Fe−N4−G system agrees with previous DFT calculations for tetra-coordinated heme systems, as well as experimental reports that suggest either a 3A2g or 3Eg triplet state.28 However, previous experimental results were found to disagree with CASPT2 predictions, which provides an example where the more computationally expensive approach fails to reproduce the experimentally determined spin-state. In addition, previous studies found that the Fe−N bond distance played a key role in the spin-state quintet-triplet splitting predictions.28 For example, increased Fe−N bond distance as shown in Figure 2 results in high ground spin-state prediction. The calculated triplet ground spin state for Fe−N4−G is like that observed for other four-coordinate Fe-macrocycles.17 A recent theoretical investigation found a triplet spin-state for a FeIIN2+2/C system in agreement with our study using DFT with periodic boundary conditions and the PBE exchangecorrelation within GGA formalism using double-ζ basis sets.37 Therefore, the higher spin-state prediction for M06-2X is possibly related to the increased amount of HF % exchange, which forces the Fe atom out-of-plane due to the unpaired electrons. It was previously observed that pseudo-octahedral 28548


Article


difference of the ⟨S2⟩ for all ORR steps per functional). The % difference of the ⟨S2⟩ was less than ∼2.1% after averaging all functionals per ORR step. (Note: The exact ⟨S2⟩ values are in Table S1) The small % difference for each reaction step and methodology indicates that the predicted lowest energy state for the penta-coordinate systems in this study do not suffer from a high level of spin-contamination. A more rigorous approach would be to correct the ground-state energy of the broken symmetry solution with respect to the high-spin state. However, it was previously found that penta-coordinate systems were less affected by spin-contamination when compared to the spin-contamination of hexa-coordinate systems.70 Moreover, since our percent difference between the calculated and expected ⟨S2⟩ are low for the DFT methods in this study, we believe the trends and evaluation to be sufficient. Experimental results showed that for isolated Fe−N4−G systems all three spin-states were observed, which was reported to be due to slight differences in the local environment (i.e., distal ligands) and potential five-coordinate Fe species.23 Specifically, a FeN2+2-like micropore-hosted site was found that exhibited a ground-state triplet multiplicity, which closely resembles the predicted ground-state multiplicity in our model system, whereas a five-coordinate N−FeN2+2-like site was found experimentally to be a quintet. Therefore, slight perturbations in the environment may result in detection of multiple spinstates for the Fe−N4-containing systems due to the disordered surroundings. For the Fe−N4−G system, the M06-L, B3LYP, PBE, and TPSSh methods appear to be more appropriate since these methods converge on a similar ground spin-state and do not significantly favor high or low spin-states compared to the M06-2X and ωB97X-D functionals. The α- and β-spin orbital TDOS and PDOS for the isolated Fe−N4−G system were compared for each method as shown in Figures 4 and S5. The populations are shown for the total density of states, as well as PDOS based on partitioning the Fe−N4−G into three fragments; Fe, N4, and the remaining CxHy atoms. The highest occupied population on the Fe atom is shown for the β-spin orbitals for both the M06-L and PBE functionals compared to their lower energy α-spin orbitals. Interestingly, the HOMO β-spin for PBE and M06-L are almost entirely attributed to the Fe atom. The TPSSh method shows a mixed PDOS for the Fe atom and the remaining N4−G layer with the HOMO β-spin energy value of ∼0.89 and ∼0.77 eV lower than that observed for M06-L and PBE methods, respectively. The relative contribution of the Fe atom for the TPSSh method is lower in comparison to the M06-L and PBE functionals yielding approximately equal PDOS contributions (i.e., when compared at ∼−6.0 eV) from the remaining N4−G layer and the Fe atom. In contrast, B3LYP and TPSSh show a small PDOS for the Fe atom for the α-spin and corresponding HOMO β-spin that include a small Fe atom contribution with additional mixed populations resulting from the remaining N4− G layer. The M06-2X and ωB97X-D functionals show a significantly less Fe atom PDOS in comparison to the M06-L, PBE and TPSSh functionals. The HOMO position of the Fe dorbitals correlates with the strength of interaction with the O2 adsorbate with the shortest Fe−O bond observed for the PBE functional. In contrast, the longest Fe−O bond was observed for the TPSSh, B3LYP and ωB97X-D functionals, which show significantly less Fe orbital contribution. Next, we considered the approximate HOMO−LUMO gaps for α-spin, which are 6.16 eV > 4.64 eV > 2.67 eV > 2.17 eV >

iron complexes systematically favor low-spin states for pure functionals and high-spin states for hybrid functionals that contain larger fractions of HF exchange.68 A thorough study of isolated Fe-containing macrocycles such as FeP, FePc, FePz, FeP(Cl), and FeP(THF)2 was performed for a broad range of DFT methods.69 It was found that DFT was the only adequate method suitable for description of isolated iron porphyrins and other Fe-containing compounds since the ground-state was accurately represented as opposed to an inaccurate ground-state representation for the expensive CI and CASPT2 methods.28,69 However, the choice of DFT method was found to be of utmost importance since it was determined that GGA and meta-GGA functionals failed to correctly predict high spin-state systems; meta-GGA methods outperformed GGA methods for high spin-state predictions.69 For this study, the spin preference for each method is also seen when comparing the high spin (HS) and low spin (LS) state energies for the Fe−N4−G model clusters. (Note: the LS and HS states for the isolated systems are m = 1 and 5.) The energy difference for the HS and LS states was calculated using eq 7 as follows: ΔEHS − LS = EHS − ELS

(7)

Figure 3 shows the ΔEHS‑LS for all DFT methods. For the isolated Fe−N4−G system, the higher spin state is favored for

Figure 3. Energy difference between high and low spin states for all functionals with the 6-311+G(2df,2p) basis set and SCRF polarized continuum model (PCM) with ε = H2O to include implicit solvation effects.

all functionals except for the ωB97X-D functional. The most negative ΔEHS‑LS is observed for the M06-2X functional with the largest HF exchange fraction. The TPSSh and B3LYP functionals favor higher spin-states to a lesser degree than M06L, PBE and M06-2X. Therefore, the M06-2X and ωB97X-D methods appear to be outliers with certain spin-states significantly favored over others. Our results are consistent with the known trends for DFT functionals in predicting iron ground spin-states with the exception of results observed for the ωB97X-D functional with 22% HF exchange at short-range and 100% at long-range, which significantly favors lower spinstates.68,69 In addition, it is common for open shell systems to exhibit spin contamination from higher spin states.70 Therefore, we verified the calculated versus expected ⟨S2⟩ values to be less than ∼2.5% for all DFT methodologies (i.e., averaging the % 28549


Article


Figure 4. Total density of states and partial density states of the α- and β-spin orbitals for the isolated Fe−N4−G systems for the M06-L, PBE, and TPSSh DFT methodologies. Note: Inset images are generated from the canonical orbitals for the highest occupied molecular orbitals.

1.74 eV > 1.55 eV and for the β-spin are 5.96 eV > 4.63 eV > 2.43 eV > 1.77 eV > 0.69 eV > 0.63 eV for the ωB97X-D, M062X, B3LYP, TPSSh, M06-L and PBE methods, respectively. The difference in HOMO−LUMO gap is interesting and is important for the interaction of the respective frontier orbitals with the adsorbate of interest. For example, ωB97X-D, M062X, B3LYP and TPSSh methods have a lower Fe PDOS at the HOMO−LUMO levels compared to the pure PBE and M06-L methods. On the basis of these considerations, we explored in detail all the ads steps to better understand the reaction parameters and the method space and how the difference in the PDOS of the central Fe site and spin-state preferences affected the overall reaction energetics. O2 Binding, Orientation and Spin-State Preference. The O2 binding energies (ΔEBE‑O2) for the O2−Fe−N4−G system are reported in Table 1 and were calculated using eq 8 as follows: ΔEBE − O2 = EO2 − Fe − N4 − G − (E Fe − N4 − G + EO2)

where EO2−Fe−N4−G, EFe−N4−G, and EO2 are the total system energies of the O2−Fe−N4−G cluster, the energy of the isolated Fe−N4−G cluster, and the energy of the isolated O2 molecule, respectively. Figure 1 steps 2a and 2b show representative images of the optimized O2−Fe−N4−G systems for either the side-on (SO) or end-on (EO) configurations, respectively. The isolated O2 molecule (i.e., zero charge with a ground triplet spin state), and H2O/OH− molecules were optimized separately at the same level of theory. The Fe−O bond distances for the optimized O2−Fe−N4−G clusters are plotted for each DFT method in Figure 5. The final optimized Fe−O bond distances vary with the method, multiplicity and starting O2 configuration. The PBE, M06-L, M06-2X, and ωB97X-D show that despite the starting configuration (i.e., O2 molecule placed SO or EO above Fe site), the final optimized Fe−O distance for the ground-state is approximately equal with the Fe−O distance for PBE < M06-2X < M06-L < ωB97X-D. In contrast, B3LYP and TPSSh predict different Fe−O bond distances for the different starting configurations. Figure 6 shows that the PDOS of the Fe atom calculated using the

(8) 28550


Article


Table 1. O2 Binding Energies (eV) and Bond Distances at 0 K without ZPE Correction for Lowest Energy Multiplicity for Fe− N4−G Systems from Single Point DFT Calculations Using 6-311+G(2df,2p) Basis Set (Unless Otherwise Noted) with and without the Polarized Continuum Model (PCM) and ε = H2O and/or Grimme’s Dispersion Correction (GD3) at Geometries Optimized Using the 6-31G(d) Basis Seta sample

method

PCM

Fe−N4−G Fe−N4−G Fe−N4−G Fe−N4−Gc Fe−N4−G Fe−N4−G Fe−N4−G Fe−N4−G Fe−N4−G Fe−N4−Gc Fe−N4−G

B3LYP M06-2X M06-L PBE TPSSh ωB97X-D B3LYP M06-2X M06-L PBE TPSSh TPSShd ωB97X-D ωB97X-Dd B3LYP B3LYPd M06-2X M06-2Xd M06-L M06-Ld PBE PBEd

x x x x x x

Fe−N4−G Fe−N4−G Fe−N4−G Fe−N4−G Fe−N4−Gc

GD3

x x x x x x x x x x

M 5, 7, 7, 3, 5, 5, 5, 7, 7, 3, 5, 5, 5, 5, 5, 5, 7, 7, 7, 7, 3, 3,

4 6 6 2 2 2 4 6 6 2 2 2 4 4 4 4 6 6 6 6 2 2

O2 (eV)

O−O (Å)b

O−OH (Å)b

0.34 −0.23 −0.29 −0.38 0.067 0.10 0.0035 −0.53 −0.49 −0.55 −0.055 −0.16 −0.86 −0.29 0.23 −0.10 −0.51 −0.43 −0.47 −0.72 −0.44 −0.68

1.316 1.301 1.320 1.268 1.287 1.266 1.316 1.301 1.320 1.268 1.287 1.287 1.266 1.266 1.316 1.316 1.301 1.301 1.320 1.320 1.268 1.268

1.426 1.441 1.442 1.449 1.439 1.411 1.426 1.441 1.442 1.449 1.439 1.439 1.415 1.415 1.426 1.426 1.441 1.441 1.442 1.442 1.449 1.449

% inc 8.8 8.4 8.2 3.3 5.5 4.6 8.8 8.4 8.2 3.3 5.5 5.5 4.6 4.6 8.8 8.8 8.4 8.4 8.2 8.2 3.3 3.3

(17.9) (20.1) (18.2) (17.8) (18.0) (16.6) (17.9) (20.1) (18.2) (17.8) (18.0) (18.0) (16.6) (16.6) (17.9) (17.9) (20.1) (20.1) (18.2) (18.2) (17.8) (17.8)

a The M (multiplicity) and % inc (% Increase) columns refer to the ground spin-states for the O2−Fe−N4−G and HOO−Fe−N4−G clusters, respectively. bBond distance determined with 6-31G(d) basis set. cDenotes O2 end-on configuration. dBinding energies predicted with the lower 631G(d) basis set.

However, ΔEBE‑O2 is binding for B3LYP at the 6-31G(d) basis set as shown for each functional in Table 1. The M06-2X and M06-L methods predict similar ground spin-states and O2 SO binding modes with a binding energy difference less than 0.1 eV for all systems including gas, solvation-corrected, and dispersion-corrected methods. The ωB97X-D method shows a significant difference for dispersion- and solvation-corrected methods with ∼1 eV positive change, which implies that adding both dispersion and solvation-corrections for this method essentially has a canceling effect in terms of binding energy, thus resulting in more positive binding energies. The negation is also observed for TPSSh but to a lesser degree, whereby the added solvation-correction resulted in nonbinding. Interestingly, the M06-L and PBE methods have different ground spin-states and binding orientations with the M06-L method ground-state binding O2 in the SO mode and the PBE method ground-state binding O2 in the EO mode despite the initial starting configuration. Moreover, multiple O2 binding modes were found for M06-L and PBE as shown in the inset in Figure 5. That is, the binding mode was preferred to be SO after increasing the number of unpaired electrons for the M06L method. In contrast, an EO binding mode for the PBE method was observed despite the spin-state with no stable configuration found for the septet system. This difference is interesting since the choice of methodology leads to different predicted results in binding orientation and spin-state orderings. Despite the subtle differences in spin-state and binding mode configurations, Figure 2 shows a relative increase in Fe− N bond distance for all methods after O2 binding, which suggests all methods geometrically account for the O2 addition. However, the degree to which the O2 adsorption pulls the Fe

Figure 5. Shortest Fe−O bond distances extracted from the optimized geometries for each reaction step adsorbate (HOH, O2−SO, O2−EO, OOH, O, OH) for all DFT functionals studied. Note: SO and EO refer to side- on and end-on, respectively. Inset images shown for alternate O2-binding configurations for the PBE and M06-L methodologies.

B3LYP and TPSSh methods is lower than that seen for the M06-L and PBE methods. (See Figure S6 for PDOS and TDOS for B3LYP, M06-2X, and ωB97X-D methods.) Therefore, placing an O2 molecule near the Fe center does not result in a strong binding affinity and instead predicts a preferable nonbinding or weak binding state, which is observed and reported in Table 1 at 0 K before zero-point correction. 28551


Article


Figure 6. Total density of states and partial density states of the α- and β-spin orbitals for the ground-state O2−Fe−N4−G systems for each DFT methodology. Note: Inset images are generated from the canonical orbitals for the highest occupied molecular orbitals.

out-of-plane varies with each functional, whereby the B3LYP and M06 family methodologies exhibit the largest increase in Fe−N distance. The Fe−N distance in this study is slightly lower than that reported for macrocyclic compounds as shown in Table 2 except for the bond distances found in this study with the M06 functional family. For comparison, Table 2 lists the absolute O2 binding energies from several theoretical reports for macrocyclic molecules and embedded Fe−Nx sites in periodic systems. For the Fe−N4-containing molecules/slabs the O2 binding energy energies ranged from −1.61 to −0.4 eV using different DFT methods, which are close to the gas phase binding energies in this study for the PBE, M06-L, M06-2X, and ωB97X-D methods. We note that the results reported in the literature are dependent on methodology and tend toward lower binding energies when using nonempirical GGA formalism methodologies with 0% HF exchange. This effect is also observed for macrocyclic systems shown in Table 2, which have, for example, a ∼0.7 eV lower binding energy for the FePc system when using the GGA-PW91 method compared to the BLYP method. In addition to binding energy,

the preferred binding mode and spin-state affects the % elongation of the O−O bond. For example, B3LYP, M06-2X, and M06-L methods have O−O elongations of ∼8−9%, whereas TPSSh and B97X-D elongations are increased to a lesser degree by ∼5−6%. Interestingly, the PBE functional shows the lowest percentage increase at ∼3%, which is expected for the predicted increased electron pairing for the lower spin system. Table 2 includes results for the heme−O2 (i.e., FeP(O2)) systems, where it was determined that the binding energy of O2 (i.e., penta-coordinated system) increased in the following order: B3LYP < PBE0 < CASPT2 < BP86.70 This finding highlights the variation in predictions per DFT methodology, which obviously is important for O2 binding and spin-state predictions. In addition, the authors found that the increase in binding energy corresponded with an increase in the O−O bond distance relative to the unbound O2 adsorbate that correspond to the trends found in this study. Despite the range of O−O bond distances found, the O−O bond increase observed for the hypothetical O2−Fe−N4−PAN clusters in this study is consistent with elongation of the O−O bond for initial 28552


Article


Table 2. O2 Binding Energies (eV) and Bond Distances at 0 K without Zero-Point Energy (ZPE) Corrections for Different FeContaining Macrocycles, Graphene, and Carbon Cluster Systems from Previous Worka sample

method

state

multiplicity

O2 (eV)

O−O (Å)

M−N (Å)

M−O (Å)

FePc71 FePcb,71 FePc72 FeP73 FeTPP73 FeTPPb,73 heme-O2c,70 heme-O2c,70 heme-O2c,70 heme-O2c,70 Fe−N4-C25 Fe−N4-Cb,25 Fe−N4−Gr26 Fe−N474 Fe−N4,b,74 FeN2+2/C37 FeN2+2/Cb,37

GGA-PW91 GGA-PW91 BLYP VWN-BP VWN-BP VWN-BP PBE0 B3LYP BP86 CASPT2 GGA-PBE GGA-PBE GGA-PBE BLYP BLYP PBE PBE

gas gas gas gas gas gas gas gas gas gas gas gas gas gas gas gas gas

LE LE LE 1 1 NR LE LE LE LE NR NR NR NR NR NR NR

−1.16 −0.40 −0.45 −0.58 −0.78 −0.15 −0.03 0.07 −0.70 −0.43 −1.45 −1.61 −0.98 −0.56 −0.28 −1.18 −1.90

1.274 1.356 1.300 1.267 1.27 1.349 ∼1.23 ∼1.23 >∼1.23 >∼1.23 1.29 1.37 NR NR NR 1.29 1.39

1.949 NR 1.958 1.981 1.972 NR NR NR NR NR NR NR NR NR NR NR NR

1.744 1.893 NR 1.733 1.732 NR ∼2.5 ∼2.5 ∼1.9 ∼1.9 1.69 1.83 NR NR NR NR NR

Note: LE, NR, and BSSE are lowest energy, not reported, and basis set superposition error, respectively. bDenotes final side-on configuration. Corrected for ZPE, BSSE, and broken symmetry spin contamination.

a c

displacing the preadsorbed molecules for each functional using the following equation:

adsorption along the ORR pathway. The O−O bond increase for the O2−Fe−N4−PAN cluster indicates bond rupture is likely following the protonation step. Wang et al.71 reported an increase in O−O bond distance for iron phthalocyanines (FePc) compared to cobalt phthalocyanines (CoPc) in an effort to describe the difference in electrocatalytic activity of the two Pc molecules with different transition metal centers. All methods account for O−O increases, however there are discrepancies in the binding modes and spin-state configurations. Therefore, we suggest that complementary methods be utilized when investigating the first O2 binding step. Despite the differences, we were interested if the choice of methodology would affect subsequent reaction steps in a similar manner. Therefore, in addition to the initial O2 adsorption, we investigated the adsorbate interactions for reaction steps 3−6. Preadsorbed H2O or OH Displacement with O2. In an AEMFC, HMFC or BCMFC with a high pH cathode, the Fe site within Fe−N4−G may be preadsorbed with species such as H2O or OH prior to O2 adsorption. Preadsorbed H2O has been discussed for nonprecious metal catalysts for proton-exchange membrane fuel cells.64,65 Therefore, to bind O2 (i.e., initiate the O2 electro-reduction), bound H2O or OH must be displaced as would be encountered in an alkaline environment. The energy of such a displacement reaction is calculated as the difference between O2 binding energy and adsorbate binding energy: ΔENet − ads = EBE − O2 − EBE − ads

ΔGSystem = ΔE + ΔZPE − T ΔS + ΔGU + ΔGpH

(10)

where ΔE is the binding energy calculated according to reaction steps 1−6, and ZPE and TΔSare the zero-point energy and entropy terms obtained from the vibrational analysis, respectively. ΔGU = −eU where U is the potential at the electrode and e is the charge transferred. ΔGpH = kBT × ln(10) × pH, where kB is the Boltzmann constant, T = 298.15 K, and pH = 13 for alkaline media. After correcting for zero-point, thermal, entropic, and pH effects, the displacement of H2O and/or OH by O2 is preferred for all methods except M06-2X, which was unable to displace H2O (i.e., for the implicit solvation model step 2). Even after increasing ΔGU the M062X method did not result in favorable displacement of H2O. This may be a result of overbinding the H2O to the Fe−N4−G in the quintet state as opposed to other methods that predict triplet spin states before and after preadsorbing H2O. The displacement of H2O or OH by the O2 molecule initiates the ORR, where it is worth noting that the displacement of OH is necessary in step 6 to regenerate the catalyst site. Figure 7 shows the PDOS of the Fe atom with the preadsorbed H2O molecule for α- and β-spin orbitals. Except for M06-2X, all methodologies show the Fe atom has approximately the same population of states before and after preadsorption of H2O, which suggests that the Fe states are available for bonding and displacing the prebound molecule(s) by an incoming O2 molecule. However, the M06-2X functional shows a shift in the Fe atom orbital character to a lower energy state due to a strong ads interaction. This is further shown in Figure 2 (step 2), where the Fe−N distance is similar for the isolated Fe−N4− G system and H2O−Fe−N4−G system prior to O2 displacement except for the M06-2X method. For all methods, the incorporation of solvation and/or dispersion lowers the overall displacement energy with the solvation-adjusted energy less than the dispersion-adjusted energy. This difference in energy is discussed in more detail in the subsequent ORR Free Energy Diagram section. With ΔGU = 0 V, the solvent-corrected

(9)

where ΔENet‑ads, EBE‑O2, and EBE‑ads refer to the net, O2, and H2O or OH binding energies, respectively. A negative ΔENet‑ads indicates that O2 is likely to displace the prebound adsorbate. The ΔENet‑H2O was negative for the PBE functional without any corrections, however the TPSSh, ωB97X-D, M06-2X, and B3LYP functionals show positive energies before zero-point and free energy corrections. Alternatively, if the Fe site is preadsorbed with an OH group, then ΔENet−OH for all DFT methods predicts that the O2 molecule will readily displace the OH group. Next, we compared the free energy (ΔG) values for 28553


Article


Figure 7. Total density of states and partial density states of the α- and β-spin orbitals for the preadsorbed HOH−Fe−N4−G systems for the M06L, PBE and TPSSh DFT methodologies. Note: Inset images are generated from the canonical orbitals for the highest occupied molecular orbitals.

energies for reaction Step [2] are lowest and increase in energy in the following order PBE < M06-L < TPSSh < B3LYP < ωB97X-D < M06-2X. This ordering corresponds to the TDOS HOMO−LUMO gap for the α- and β-spin previously discussed for the isolated Fe−N4−G system with the displacement reaction being more favorable as the HOMO−LUMO gap decreases. Similar conclusions are made for the preadsorbed OH species. Since the HO−Fe−N4−G system is encountered again in reaction step 5, the TDOS and PDOS are discussed in more detail in subsequent sections in parallel to the displacement in step 2. We note that the preadsorption and coverage of Fe−N4 sites with OH was previously described in terms of water activation and thus requires a Fe2+ state to allow efficient adsorption of the O2 molecule.17 Moreover, it was shown a + 3 to +2 reduction of the Fe center is needed to favorably displace the prebound species.47 Therefore, the isolated Fe−N4−G system and subsequent displacement step are described most aptly for the B3LYP, TPSSh, M06-L, and PBE methods (i.e., with respect to the implicit solvation model), whereas the ωB97X-D and M06-2X methods are not

suitable for the conditions considered. In all cases, the displacement of a preadsorbed H2O or OH molecules by O2 initiates the ORR in a series of proton−electron transfer steps. ORR Adsorbates (O2, O, OOH and OH) Bound to Fe− N4−G. After O2 binding to the Fe site, subsequent reactions take place as shown in reaction steps 3−6 and Figures 1e−g. These include protonation of O2−Fe−N4−G to generate a bound HOO−Fe−N4−G. This may be followed by breaking of the O−O bond, thus leading to radical-bound O−Fe−N4−G, which then leads to an additional protonation step that creates the HO−Fe−N4−G system. The binding energies (ΔEBE‑ads) of each step were determined as follows: ΔE BE − ads = Eads − Fe − N4 − G − (E Fe − N4 − G + Eads)

(11)

where Eads−Fe−N4−G, EFe−N4−G, and Eads are the total system energies of the ads−Fe−N4−G cluster, the energy of the isolated Fe−N4−G cluster, and the energy of the isolated ads molecule, respectively. (Note: ads refers to the O2, OOH, O, or OH adsorbates) After protonation of the bound O2 molecule to generate the HOO−Fe−N4−G, the orientation for all systems 28554


Article


Figure 8. Total density of states and partial density states of the α- and β-spin orbitals for the lowest energy HOO−Fe−N4−G systems for each DFT method. Note: Inset images are generated from the canonical orbitals for the highest occupied molecular orbitals.

O−O bond length elongations greater than ∼17%. This suggests that the bound HOO leads to a break in the O−O bond to desorb an OH− molecule and create a O−Fe−N4−G system as shown in reaction step 4. An O−O bond distance of 1.45 Å was previously reported for a bound HOO−Fe−Nxcontaining system, which is close to the M06-2X, M06-L, and PBE methods as shown in Table 1.37 This is further illustrated in Figures 8 and S7 that show the PDOS and TDOS for the OOH−Fe−N4−G systems. The α- and β-spin show larger Fe PDOS at the HOMO level for the PBE and TPSSh methods compared to the M06-L method. In contrast, the M06-L method shows higher Fe PDOS at the LUMO level for the βspin compared to the PBE and TPSSh methods. The lower Fe population available at the HOMO−LUMO gap indicates that the site will not favorably compete against other ads at this step and will instead proceed through either desorption or elongation to break. After the O−O bond breaks in reaction step 4, a radical oxygen atom is attached to Fe to generate the O−Fe−N4−G system. This results in changes to the Fe−N with a slight

adopted an EO binding mode with similar Fe−O bond distances regardless of the DFT method. However, the Fe−N distance for M06-L and M06-2X were ∼2.1 Å compared to ∼1.94 Å for all other methods. This is correlated with the spinstate preference with M06-L and M06-2X methods predicting higher sextet spin-states. The PBE, TPSSh and ωB97X-D methods predict doublets, whereas B3LYP predicts a quartet ground spin-state. This is shown in Figure 3 where, for example, PBE and M06-2X possess preferable low and high spin-states, respectively. Figure 3 also shows that the ΔEHS‑LS for the PBE and M06-2X methods have nearly a 4−5 eV difference. This is significant and highlights the low- and highspin predictions for pure and hybrid DFT methods, respectively. Moreover, the predicted ground spin-state per DFT functional creates systems that are not geometrically equivalent. For example, Figure 2 shows that the M06 methods have increased Fe−N distance, whereas the remaining methods have smaller Fe−N distance due to increased spin-pairing. Despite the spin-state difference, all methods have similar Fe− O bond distance as shown in Figure 5, and all methods result in 28555


Article


Figure 9. Comparison of relative Free energy change (ΔG) per reaction step for the entire reaction profile with (a) B3LYP, (b) M06-L, (c) M06-2X, (d) ωB97X-D, (e) PBEPBE, and (f) TPSSh methods and the 6-311+G(2df,2p) basis set, SCRF polarized continuum model (ε = H2O), pH = 13, T = 298 K, and a range of applied potentials from U = 0 to 1.23 V (Note: U = 0 V is designated as the open-circuit potential, i.e., no overpotential).

the overall free energy diagrams to understand if these differences alter the reaction profiles. ORR Free Energy Diagrams. The ΔGSystem for reaction steps 1−6 were calculated using eq 10 for a range of ΔGU, where the overpotential is defined as η = U0 − U, and U0 refers to 1.23 V vs SHE at standard conditions.63 The free energy diagrams are discussed with respect to the applied potential and are only discussed with respect to overpotential (i.e., pHcorrected to SHE) when making comparisons to experimental reports. The zero-reference state in this system is set by reaction step 1, which essentially preadsorbs a H2O or OH molecule followed by displacement in step 2 after pH correction for O2/O2−, which is considered a surfaceindependent outer-sphere redox process as previously described with a decrease in η of ∼0.76 V.31,60 Steps 3−6 undergo single electron-transfer steps assuming an inner-sphere reaction.60 Figure 9 shows the difference per reaction steps 1−6 for all six methods in this study with solvation-corrected energies. The Grimme dispersion and gas-phase calculations are shown in Figures S8 and S9, respectively. When ΔGU = 0 V, no method

upshift for the B3LYP and PBE functionals, whereas the M06-L and M06-2X methods show a significant decrease and increase, respectively. Moreover, the high-spin state is significantly favored for M06-2X and M06-L methods compared to low spin-state preference for the PBE, ωB97X-D, and TPSSh methods as shown in Figure 3. Figure 3 shows that M06-L is unique in that it straddles the zero-line indicating the method only slightly favors higher spin-states, and may suggest the method will be useful for examining reactions where degeneracy is known to exist, and where multiple spin-states are found for similar Fe-coordinations.15 Moreover, the Fe−O bond distance for the O−Fe−N4−G system is decreased compared to the HOO−Fe−N4−G bound system. Protonation of the oxygen in the O−Fe−N4−G cluster leads to a HO−Fe− N4−G system, which results in an increased Fe−O bond distance and a substantial Fe−N bond distance change for M06-L, M06-2X and ωB97X-D methods in contrast to the TPSSh, PBE, and B3LYP methods that show little change. The use of different pure and hybrid DFT methods has subtle changes in individual ads systems, which led us to investigate 28556


Article


Figure 10. Total density of states and partial density states of the α- and β-spin orbitals for the lowest energy O−Fe−N4−G systems for each DFT method. Note: Inset images are generated from the canonical orbitals for the highest occupied molecular orbitals.

as less Fe states are available near the HOMO−LUMO gap the reactivity decreases as shown in Figure 10. Parts a and b of Figure 11 overlay reaction steps 1−6 when the ΔGU is 0.2 and 0.5 V, respectively. It is possible that a high activation energy for the O−O bond scission could prohibit the 4e− ORR. Therefore, barriers were investigated for the PBE, M06-L and B3LYP methods for the OOH dissociation step. The influence of solvent on the O−OH reaction was also examined by contrasting the reaction energies for the gas-phase with those obtained using implicit solvent represented via the PCM (i.e., ε = H2O) model. The PBE, M06-L, and B3LYP methods were chosen per Figure 9. Figure 9 shows qualitatively similar trends for B3LYP when compared to M06-L for steps 4−6. Similarly, the PBE and TPSSh methods reveal similar trends for the same ORR steps. Figure 11 inset shows that during O−O bond scission, the OH breaks away from the bound OOH molecule and readsorbs on a nearby carbon atom. The bond distances and geometries are shown in Figure S11 for the PBE/6-31+G(d,p) calculations with PCM (i.e., ε = H2O) for implicit solvent. The IRC

predicts downhill energies for all reaction steps with the M06-L, PBE, B3LYP, and TPSSh methods showing uphill energies for the O−Fe−N4−G and HO−Fe−N4−G reaction steps. The number of electrons transferred is ≥3 per O2 molecule for B3LYP, M06-L, PBE and TPSSh when the ΔGU > 0.2, 0.5, 0.5, and 0.75 V, respectively. When all reaction steps are downhill then the number of electrons transferred per O2 molecule is referred to as the complete reduction (4-electron transfer), which approximately occurs for B3LYP, M06-L, PBE and TPSSh when the ΔGU > 0.5, 0.75, 0.5, and 0.75 V, respectively. The uphill energies are most readily seen when transitioning from step 4 to step 5 for all methodologies in Figure 9 except M06-2X and ωB97X-D, which have an uphill transition from step 3 to step 4. To better understand the protonation of the O−Fe−N4−G system, Figure 10 compares the PDOS and TDOS for the M06-L, PBE, and TPSSh methods. (The PDOS and TDOS for M06-2X, ωB97X-D and B3LYP functionals are shown in Figure S10.) The increasing uphill energy from steps 4 to 5 correlates with the increase in the HOMO−LUMO difference with the gaps for PBE < M06-L < TPSSh. Moreover, 28557


Article


shown in Figure S11. The B3LYP functional yields less exothermic reaction energy and higher barrier than the PBE functional, while the M06-L functional yields a slightly more exothermic reaction energy and lower barrier. These results further emphasize the importance of methodology on overall ORR predictions69,70 Finally, the OH-desorption from the O− OH bound group was examined. The OH adsorption energy and free energies were found to be 1.26 and 0.73 eV from PBE/ 6-31+G(d,p) PCM (i.e., ε = H2O) calculations, however, after the next electron is transferred to the complex, the adsorption free energy is reduced to 0.07 eV indicating that the OHdesorption from the nearby carbon atom is not likely to be a limiting step. All methods show downhill energies for most steps at a certain overpotential, except desorption of the subsequent hydroxide species into solution. Therefore, at low overpotentials (Figure 11a) the ORR proceeds through reaction steps 3−4, and requires an additional overpotential (Figure 11b) to proceed through a 3−4 electron-transfer reaction. The ORR proceeds more readily per B3LYP, M06-L, TPSSh and PBE methodologies. However, only B3LYP predicts that the ORR will be a near 4-electron transfer when ΔGU = 0.5 V. The M06-L and PBE methods predict that when ΔGU = 0.5 V there will only be a 3−4 electron-transfer reaction. These methods are in good agreement with experimental results that show a potential-dependent ORR with the number of electrons ∼4 for heat-treated non-PGM Fe−Nx-containing materials.11,57,60 Figure 12 shows the PDOS and TDOS for the HO−Fe− N4−G before desorption, which is still shown to be slightly uphill even with a ΔGU = 0.5 V. (Figure S12 shows the PDOS and TDOS for the M06-2X, B3LYP, and ωB97X-D for comparison) This agrees with a previous study that suggests that the HO-bound metal may be used as an activity descriptor for a wide-range of non-PGM electrocatalysts.34 All methods show similar Fe−O bond distances, however the expansion of the Fe−N ring plays an important role in removing the OHbound species. Interestingly, the M06-L method allows for the OH ads to come out-of-plane with an increase in Fe−N distance. Also, there is an increase in the HOMO−LUMO gap compared to the PBE and TPSSh methods that show smaller Fe−N bond distances and lower Fe states near the HOMO− LUMO gap. However, despite the slight uphill energy difference, the free energy profiles predict that the transfer will be 3 to 4 electrons. Recently, it was shown that heat-treated Fe−N4-containing electrospun fibers exhibited ORR selectivities with calculated 3.9 and 3.6 electrons (i.e., different heattreatment conditions) transferred per O2 molecule in a pH = 13 media at 0.6 V, respectively.57 The number of electrons transferred per O2 molecule was potential-dependent and only approached a complete reduction of O2 at potentials ≥0.6 V. The increase in electron transfer was potentially related to the increase in graphitic carbon content and porosity as witnessed by a decrease in the full-width at half-maximum of the C 1s peak in the X-ray photoelectron spectra and an increase in the surface area, respectively. It has also been reported that a slight change in the carbon atom configuration surrounding the Fe− N4 site will alter the overall free energy diagram; for example, increasing the graphitic sites has been shown to decrease the overall allowed distortion of the Fe site yielding lower binding energies for the O2 molecule.26,37 The DFT methodologies in this study reveal subtle changes for each ads step. Moreover, the choice of DFT functional has a significant effect on the groundstate geometries, which allow the adoption of various spin-state

Figure 11. Comparison of relative Free energy change (ΔG) per reaction step for entire profile for all functionals with the 6311+G(2df,2p) basis set, SCRF polarized continuum model (ε = H2O), pH = 13, T = 298 K and an applied potential of (a) U = 0.2 V (η ∼ 1.0 V vs SHE) and (b) U = 0.5 V (η ∼ 0.7 V vs SHE). Inset (top) shows the dissociation of the OOH adsorbate on the Fe site. TS and * refer to transition state and Fe−N4−G site, respectively.

calculations confirmed the connection of the reaction path between the barrier and reactant geometry. The ΔE barrier for dissociation for PBE is 0.41 eV for the 6-311+G(2df,2p) basis set with implicit solvation, which is very similar to the barrier 0.38 eV obtained from PBE calculations using smaller basis set 6-31+G(d,p) level, which suggests that the smaller 6-31+G(d,p) basis set is adequate for investigating this reaction. (See Figure S11) After free-energy corrections, the barrier is further reduced from 0.38 to 0.27 eV, which would increase the overpotential necessary to drive the 4e− reduction of O2. Inclusion of the implicit solvent reduced the barrier from 0.71 to 0.38 eV at the PBE/6-31G(d,p) level indicating a beneficial effect of the polarizing solvent on the reaction kinetics and a need to include solvent effects in DFT calculations for the ORR. The PBE/6-31+G(d,p) barrier in the gas-phase of 0.71 eV is slightly higher than the barrier of 0.56 eV for the O−OH scission reported from periodic plane-wave DFT study for a similar substrate.26 A dependence of the O−OH scission reaction energy and barriers on the density functional is also 28558


Article


Figure 12. Total density of states and partial density states of the α- and β-spin orbitals for the lowest energy HO−Fe−N4−G systems for each DFT method. Note: Inset images are generated from the canonical orbitals for the highest occupied molecular orbitals.

functionals, respectively, then the individual adsorbate systems favor higher spin-states and overbind the adsorbate making these functionals less feasible for non-PGM Fe−N4-containing systems. This is most evident by the removal of OH− in the final HO−Fe−N4−G cluster desorption step, which positively increases up to 2 eV in energy after protonation of the O−Fe− N4−G cluster. However, the pure functionals with no HF exchange such as M06-L and PBE showed the smallest HOMO−LUMO gaps for the isolated Fe−N4−G system, which correlated with the ORR reactivity and ability to proceed through subsequent transfer steps including removal of the bound HO-group with increasing overpotentials. Interestingly, pure DFT methodologies showed the highest PDOS for the Fe atom located in the frontier orbitals, whereas the hybrid and meta-hybrid DFT methods, with increased HF quantities, showed smaller or little PDOS for the Fe atom with states delocalized onto the surrounding N4−G cluster. Qualitatively,

preferences. If one disregards the spin-state energy preference and geometric differences, and only considers the overall ORR profile then many DFT functionals (B3LYP, M06-L, TPSSh, and PBE) appear to arrive at similar conclusions after considering subtle binding issues for the initial O2 adsorption step. Moreover, the DFT methodologies in this study reach qualitatively similar conclusions compared to previous theoretical work and support the findings reported for several experimental studies.

■

CONCLUSIONS The ORR was studied on Fe−N4−G clusters with different DFT methodologies using all-electron basis sets. Different ground spin-state preferences were found for different DFT functionals at each adsorbate step for the proposed innersphere ORR. When the HF quantity is increased locally and/or increased nonlocally, as for the M06-2X and ωB97X-D 28559


Article

The Journal of Physical Chemistry C the free energy profiles obtained were similar for pure empirical, pure nonempirical and hybrid methodologies with low fractions of HF exchange. However, the pure DFT methods differed in their prediction for the lowest energy spin-states for the different adsorbates. The O2−Fe−N4−G system with the PBE methodology was found to favor a low spin-state in an end-on mode and the M06-L method favored a higher spin state in a side-on mode. In contrast, the binding modes for O, OH and OOH did not vary for the DFT methods, however the predicted ground spin-states and binding energies were different. Ultimately, the difference in spin-state preference and binding energy mode did not affect the overall energy profile for the O2, OOH, O, and OH steps, which showed that the reactions are approximately downhill in ground-state energies for select methodologies when the potential is ∼0.7−0.8 V vs SHE for thermodynamic considerations. In addition, the activation energy for the O− OH scission was investigated since it could prohibit a true 4e− transfer. It was found that the free-energy barrier for O−OH scission is 0.27 eV (i.e., with solvation correction) obtained from PBE calculations, whereas B3LYP and M06-L functionals yielded slightly higher and lower barriers, respectively. Moreover, desorption of the OH-bound on a nearby carbon atom following O−OH scission is removed after electron addition and is not thought to be a limiting step. All methods reveal uphill energies in the final protonation and OH desorption steps, which varied from 0.5 to 2 eV depending on the DFT method. The inclusion of dispersion or solvation corrections resulted in lower overall free energy profiles for each ORR step, however the inclusion of both corrections had a near canceling effect. Our results suggest that when studying the ORR on nonPGM materials and electrocatalytic sites embedded in disorganized environments that multiple DFT methodologies be utilized with low fractions of HF exchange. This includes hybrid methods with low HF quantities, as well as pure DFT methodologies. Moreover, when considering pure DFT methodologies it may be more informative to include a comparison between empirical and nonempirical methods. The qualitative comparisons between the six methodologies in this study allow one to predict potentially useful DFT methods. The exact role that solvation, dispersion, and HF exchange quantity play in influencing the overall ORR profile is not currently known and will be the subject of a future study. Lastly, the findings here support previous theoretical and experimental work reported for Fe−N4-containing extended systems with an added emphasis on spin-state dependencies and geometric considerations within each DFT formalism and the inclusion of implicit solvation. Our work serves as a first step toward understanding the sensitivity and reliability of DFT predictions of the ORR on non-PGM electrocatalysts. The insights gained through the comparison of the DFT methods in this study also provide experimentalists with additional knowledge for developing new non-PGM electrocatalysts.

■

■

spin contamination tables, transition-state geometries, and activation energies (PDF)

AUTHOR INFORMATION

Corresponding Author

*(J.P.M.) E-mail: [email protected]. Telephone: (301)394-4375. ORCID

Joshua P. McClure: 0000-0002-2442-1139 Oleg Borodin: 0000-0002-9428-5291 Author Contributions

The manuscript was written through contributions of all authors. All authors have given approval to the final version of the manuscript. Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS Thanks are expressed to the U.S. Department of the Army and U.S. Army Materiel Command for supporting this work. We would like to thank the DoD High Performance Computing Modernization Program (HPCMP) for computing facilities and support.

■

REFERENCES

(1) Unlu, M.; Abbott, D.; Ramaswamy, N.; Ren, X. M.; Mukerjee, S.; Kohl, P. A. Analysis of Double Layer and Adsorption Effects at the Alkaline Polymer Electrolyte-Electrode Interface. J. Electrochem. Soc. 2011, 158 (11), B1423−B1431. (2) Unlu, M.; Zhou, J.; Kohl, P. A. Self-Humidifying Hybrid AnionCation Membrane Fuel Cell Operated Under Dry Conditions. J. Electrochem. Soc. 2010, 157 (10), B1391−B1396. (3) Jiang, R. Z.; Chu, D. Comparative Study of CoFeNx/C Catalyst Obtained by Pyrolysis of Hemin and Cobalt Porphyrin for Catalytic Oxygen Reduction in Alkaline and Acidic Electrolytes. J. Power Sources 2014, 245, 352−361. (4) Liao, M. S.; Scheiner, S. Comparative Study of Metal-Porphyrins, -Porphyrazines, and -Phthalocyanines. J. Comput. Chem. 2002, 23 (15), 1391−1403. (5) Jiang, R. Z.; Tran, D. T.; McClure, J. P.; Chu, D. Increasing the Electrochemically Available Active Sites for Heat-treated Hemin Catalysts Supported on Carbon Black. Electrochim. Acta 2012, 75, 185−190. (6) McClure, J. P.; Devine, C. K.; Jiang, R.; Chu, D.; Cuomo, J. J.; Parsons, G. N.; Fedkiw, P. S. Oxygen Electroreduction on Ti- and FeContaining Carbon Fibers. J. Electrochem. Soc. 2013, 160 (8), F769− F778. (7) McClure, J. P.; Thornton, J. D.; Jiang, R. Z.; Chu, D.; Cuomo, J. J.; Fedkiw, P. S. Oxygen Reduction on Metal-Free Nitrogen-Doped Carbon Nanowall Electrodes. J. Electrochem. Soc. 2012, 159 (11), F733−F742. (8) Ratso, S.; Kruusenberg, I.; Vikkisk, M.; Joost, U.; Shulga, E.; Kink, I.; Kallio, T.; Tammeveski, K. Highly Active Nitrogen-Doped FewLayer Graphene/Carbon Nanotube Composite Electrocatalyst for Oxygen Reduction Reaction in Alkaline Media. Carbon 2014, 73, 361−370. (9) Wang, Z. J.; Jia, R. R.; Zheng, J. F.; Zhao, J. G.; Li, L.; Song, J. L.; Zhu, Z. P. Nitrogen-Promoted Self-Assembly of N-Doped Carbon Nanotubes and their Intrinsic Catalysis for Oxygen Reduction in Fuel Cells. ACS Nano 2011, 5 (3), 1677−1684. (10) Faubert, G.; Lalande, G.; Cote, R.; Guay, D.; Dodelet, J. P.; Weng, L. T.; Bertrand, P.; Denes, G. Heat-treated Iron and Cobalt Tetraphenylporphyrins Adsorbed on Carbon Black: Physical Characterization and Catalytic Properties of these Materials for the Reduction

ASSOCIATED CONTENT

S Supporting Information *

The Supporting Information is available free of charge on the ACS Publications website at DOI: 10.1021/acs.jpcc.6b08498. Additional total and partial density of states plots for the different density functionals and methodologies, freeenergy profiles for gas-phase and dispersion calculations, 28560


Article

The Journal of Physical Chemistry C of Oxygen in Polymer Electrolyte Fuel Cells. Electrochim. Acta 1996, 41 (10), 1689−1701. (11) Jiang, R. Z.; Tran, D. T.; McClure, J.; Chu, D. Heat-treated Hemin Supported on Graphene Nanoplatelets for the Oxygen Reduction Reaction. Electrochem. Commun. 2012, 19, 73−76. (12) Wu, J.; Li, W. M.; Higgins, D.; Chen, Z. W. Heat-Treated NonPrecious Catalyst Using Fe and Nitrogen-Rich 2,37,8-Tetra(pyridin-2yl)pyrazino[2,3-g]quinoxaline Coordinated Complex for Oxygen Reduction Reaction in PEM Fuel Cells. J. Phys. Chem. C 2011, 115 (38), 18856−18862. (13) Wu, G.; More, K. L.; Johnston, C. M.; Zelenay, P. HighPerformance Electrocatalysts for Oxygen Reduction Derived from Polyaniline, Iron, and Cobalt. Science 2011, 332 (6028), 443−447. (14) Birry, L.; Zagal, J. H.; Dodelet, J. P. Does CO Poison Fe-Based Catalysts for ORR? Electrochem. Commun. 2010, 12 (5), 628−631. (15) Kramm, U. I.; Herranz, J.; Larouche, N.; Arruda, T. M.; Lefèvre, M.; Jaouen, F.; Bogdanoff, P.; Fiechter, S.; Abs-Wurmbach, I.; Mukerjee, S.; et al. Structure of the Catalytic Sites in Fe/N/CCatalysts for O2-Reduction in PEM Fuel Cells. Phys. Chem. Chem. Phys. 2012, 14 (33), 11673−88. (16) Rao, C. V.; Cabrera, C. R.; Ishikawa, Y. In Search of the Active Site in Nitrogen-Doped Carbon Nanotube Electrodes for the Oxygen Reduction Reaction. J. Phys. Chem. Lett. 2010, 1 (18), 2622−2627. (17) Tylus, U.; Jia, Q.; Strickland, K.; Ramaswamy, N.; Serov, A.; Atanassov, P.; Mukerjee, S. Elucidating Oxygen Reduction Active Sites in Pyrolyzed Metal-Nitrogen Coordinated Non-Precious-Metal Electrocatalyst Systems. J. Phys. Chem. C 2014, 118 (17), 8999−9008. (18) Yang, J. B.; Liu, D. J.; Kariuki, N. N.; Chen, L. X. Aligned Carbon Nanotubes with Built-in FeN4 Active Sites for Electrocatalytic Reduction of Oxygen. Chem. Commun. 2008, 3, 329−331. (19) Yu, L.; Pan, X. L.; Cao, X. M.; Hu, P.; Bao, X. H. Oxygen Reduction Reaction Mechanism on Nitrogen-Doped Graphene: A Density Functional Theory Study. J. Catal. 2011, 282 (1), 183−190. (20) Zhang, L. P.; Xia, Z. H. Mechanisms of Oxygen Reduction Reaction on Nitrogen-Doped Graphene for Fuel Cells. J. Phys. Chem. C 2011, 115 (22), 11170−11176. (21) Yeager, E. Dioxygen Electrocatalysis: Mechanisms in Relation to Catalyst Structure. J. Mol. Catal. 1986, 38 (1−2), 5−25. (22) Ramaswamy, N.; Mukerjee, S. Fundamental Mechanistic Understanding of Electrocatalysis of Oxygen Reduction on Pt and Non-Pt Surfaces: Acid versus Alkaline Media. Adv. Phys. Chem. 2012, 2012, 491604. (23) Kramm, U. I.; Herrmann-Geppert, I.; Behrends, J.; Lips, K.; Fiechter, S.; Bogdanoff, P. On an Easy Way to Prepare Metal-Nitrogen Doped Carbon with Exclusive Presence of MeN4-type Sites Active for the ORR. J. Am. Chem. Soc. 2016, 138 (2), 635−40. (24) Lefevre, M.; Dodelet, J. P.; Bertrand, P. Molecular Oxygen Reduction PEM Fuel Cells: Evidence for the Simultaneous Presence of Two Active Sites in Fe-Based Catalysts. J. Phys. Chem. B 2002, 106 (34), 8705−8713. (25) Kattel, S.; Wang, G. F. A Density Functional Theory Study of the Oxygen Reduction Reaction on Me-N4 (Me = Fe, Co, or Ni) Clusters Between Graphitic Pores. J. Mater. Chem. A 2013, 1 (36), 10790−10797. (26) Kattel, S.; Wang, G. F. Reaction Pathway for Oxygen Reduction on FeN4 Embedded Graphene. J. Phys. Chem. Lett. 2014, 5 (3), 452− 456. (27) Jacob, C. R.; Reiher, M. Spin in Density-Functional Theory. Int. J. Quantum Chem. 2012, 112 (23), 3661−3684. (28) Chen, H.; Lai, W. Z.; Shaik, S. Multreference and Multiconfiguration Ab Initio Methods in Heme-Related Systems: What Have We Learned So Far? J. Phys. Chem. B 2011, 115 (8), 1727−1742. (29) Cramer, C. J.; Truhlar, D. G. Density Functional Theory for Transition Metals and Transition Metal Chemistry. Phys. Chem. Chem. Phys. 2009, 11 (46), 10757−10816. (30) Cohen, A. J.; Mori-Sanchez, P.; Yang, W. T. Challenges for Density Functional Theory. Chem. Rev. 2012, 112 (1), 289−320. (31) Kattel, S.; Atanassov, P.; Kiefer, B. Catalytic Activity of CoN(x)/C Electrocatalysts for Oxygen Reduction Reaction: A Density

Functional Theory Study. Phys. Chem. Chem. Phys. 2013, 15 (1), 148− 153. (32) Chen, X.; Sun, S. R.; Wang, X. Y.; Li, F.; Xia, D. G. Density Functional Theory Study of the Oxygen Reduction on a CobaltPolypyrrole Composite Catalyst. J. Phys. Chem. C 2012, 116 (43), 22737−22742. (33) Orellana, W. Catalytic Properties of Transition Metal-N4 Moieties in Graphene for the Oxygen Reduction Reaction: Evidence of Spin-Dependent Mechanisms. J. Phys. Chem. C 2013, 117 (19), 9812−9818. (34) Guo, J. S.; He, H.; Chu, D.; Chen, R. R. OH−-Binding Effects on Metallophthalocyanine Catalysts for O2 Reduction Reaction in Anion Exchange Membrane Fuel Cells. Electrocatalysis 2012, 3 (3−4), 252− 264. (35) He, H.; Lei, Y.; Xiao, C.; Chu, D.; Chen, R.; Wang, G. Molecular and Electronic Structures of Transition-Metal Macrocyclic Complexes as Related to Catalyzing Oxygen Reduction Reactions: A Density Functional Theory Study. J. Phys. Chem. C 2012, 116 (30), 16038− 16046. (36) Chen, R. R.; Li, H. X.; Chu, D.; Wang, G. F. Unraveling Oxygen Reduction Reaction Mechanisms on Carbon-Supported Fe-Phthalocyanine and Co-Phthalocyanine Catalysts in Alkaline Media. J. Phys. Chem. C 2009, 113 (48), 20689−20697. (37) Szakacs, C. E.; Lefèvre, M.; Kramm, U. I.; Dodelet, J. P.; Vidal, F. A Density Functional Theory Study of Catalytic Sites for Oxygen Reduction in Fe/N/C Catalysts Used in H2/O2 Fuel Cells. Phys. Chem. Chem. Phys. 2014, 16 (27), 13654−61. (38) Frisch, M. J.; Trucks, G. W.; Schlegel, H. B.; Scuseria, G. E.; Robb, M. A.; Cheeseman, J. R.; Scalmani, G.; Barone, V.; Mennucci, B.; Petersson, G. A., et al. Gaussian 09, Revision A.02; Gaussian, Inc.: Wallingford CT, 2009. (39) Zhao, Y.; Truhlar, D. G. A New Local Density Functional for Main-Group Thermochemistry, Transition-Metal Bonding, Thermochemical Kinetics, and Noncovalent Interactions. J. Chem. Phys. 2006, 125 (19), 194101. (40) Zhao, Y.; Truhlar, D. G. The M06 Suite of Density Functionals for Main Group Thermochemistry, Thermochemical Kinetics, Noncovalent Interactions, Excited States, and Transition Elements: Two New Functionals and Systematic Testing of Four M06-Class Functionals and 12 Other Functionals. Theor. Chem. Acc. 2008, 120 (1−3), 215−241. (41) Perdew, J. P.; Burke, K.; Ernzerhof, M. Generalized Gradient Approximations Made Simple. Phys. Rev. Lett. 1996, 77 (18), 3865− 3868. (42) Becke, A. D. Density-Functional Thermochemistry. III. The Role of Exact Exchange. J. Chem. Phys. 1993, 98 (7), 5648−5652. (43) Lee, C. T.; Yang, W. T.; Parr, R. G. Development of the ColleSalvetti Correlation-Energy Formula into a Functional of the Electron Density. Phys. Rev. B: Condens. Matter Mater. Phys. 1988, 37 (2), 785− 789. (44) Chai, J. D.; Head-Gordon, M. Long-Range Corrected Hybrid Density Functionals with Damped Atom-Atom Dispersion Corrections. Phys. Chem. Chem. Phys. 2008, 10 (44), 6615−6620. (45) Tao, J. M.; Perdew, J. P.; Staroverov, V. N.; Scuseria, G. E. Climbing the Density Functional Ladder: Nonempirical MetaGeneralized Gradient Approximation Designed for Molecules and Solids. Phys. Rev. Lett. 2003, 91 (14), 146401. (46) Kostov, M. K.; Santiso, E. E.; George, A. M.; Gubbins, K. E.; Nardelli, M. B. Dissociation of Water on Defective Carbon Substrates. Phys. Rev. Lett. 2005, 95 (13), 136105. (47) Anderson, A. B.; Sidik, R. A. Oxygen Electroreduction on FeII and FeIII Coordinated to N4 Chelates. Reversible Potentials for the Intermediate Steps from Quantum Theory. J. Phys. Chem. B 2004, 108 (16), 5031−5035. (48) Rassolov, V. A.; Ratner, M. A.; Pople, J. A.; Redfern, P. C.; Curtiss, L. A. 6-31G* Basis Set for Third-Row Atoms. J. Comput. Chem. 2001, 22 (9), 976−984. 28561


Article

The Journal of Physical Chemistry C (49) Bacskay, G. B. A Quadratically Convergent Hartree-Fock (QCSCF) Method. Application to Closed Shell Systems. Chem. Phys. 1981, 61 (3), 385−404. (50) Curtiss, L. A.; McGrath, M. P.; Blaudeau, J. P.; Davis, N. E.; Binning, R. C.; Radom, L. Extension of Gaussian-2 Theory to Molecules Containing Third-Row Atoms Ga-Kr. J. Chem. Phys. 1995, 103 (14), 6104−6113. (51) Scalmani, G.; Frisch, M. J. Continuous Surface Charge Polarizable Continuum Models of Solvation. I. General Formalism. J. Chem. Phys. 2010, 132 (11), 114110. (52) Grimme, S.; Antony, J.; Ehrlich, S.; Krieg, H. A Consistent and Accurate Ab Initio Parameterization of Density Functional Dispersion Correction (DFT-D) for the 94 Elements H-Pu. J. Chem. Phys. 2010, 132 (15), 154104. (53) Simon, S.; Duran, M.; Dannenberg, J. J. How Does Basis Set Superposition Error Change the Potential Surfaces for Hydrogen Bonded Dimers? J. Chem. Phys. 1996, 105 (24), 11024−11031. (54) Chen, H.; Cho, K. B.; Lai, W.; Nam, W.; Shaik, S. Dioxygen Activation by a Non-Heme Iron(II) Complex: Theoretical Study Toward Understanding Ferric-Superoxo Complexes. J. Chem. Theory Comput. 2012, 8 (3), 915−26. (55) Gorelsky, S. I. AOMix, version 6.82; Program for Molecular Orbital Analysis. (56) Gorelsky, S. I.; Lever, A. B. P. Electronic Structure and Spectra of Ruthenium Diimine Complexes by Density Functional Theory and INDO/S. Comparison of the Two Methods. J. Organomet. Chem. 2001, 635, 187−196. (57) McClure, J. P.; Jiang, R. J.; Chu, D.; Fedkiw, P. S. Oxygen Electroreduction on Fe- or Co-Containing Carbon Fibers. Carbon 2014, 79, 457−469. (58) Deng, D.; Chen, X.; Yu, L.; Wu, X.; Liu, Q.; Liu, Y.; Yang, H.; Tian, H.; Hu, Y.; Du, P.; et al. A Single Iron Site Confined in a Graphene Matrix for the Catalytic Oxidation of Benzene at Room Temperature. Sci. Adv. 2015, 1 (11), e1500462. (59) Watanabe, M.; Tryk, D. A.; Wakisaka, M.; Yano, H.; Uchida, H. Overview of Recent Developments in Oxygen Reduction Electrocatalysis. Electrochim. Acta 2012, 84, 187−201. (60) Ramaswamy, N.; Mukerjee, S. Influence of Inner- and OuterSphere Electron Transfer Mechanisms During Electrocatalysis of Oxygen Reduction in Alkaline Media. J. Phys. Chem. C 2011, 115 (36), 18015−18026. (61) Sha, Y.; Yu, T. H.; Merinov, B. V.; Goddard, W. A. Prediction of the Dependence of the Fuel Cell Oxygen Reduction Reactions on Operating Voltage from DFT Calculations. J. Phys. Chem. C 2012, 116 (10), 6166−6173. (62) Guo, J. S.; Li, H. X.; He, H.; Chu, D.; Chen, R. R. CoPc- and CoPcF16-Modified Ag Nanoparticles as Novel Catalysts with Tunable Oxygen Reduction Activity in Alkaline Media. J. Phys. Chem. C 2011, 115 (17), 8494−8502. (63) Norskov, J. K.; Rossmeisl, J.; Logadottir, A.; Lindqvist, L.; Kitchin, J. R.; Bligaard, T.; Jonsson, H. Origin of the Overpotential for Oxygen Reduction at a Fuel-Cell Cathode. J. Phys. Chem. B 2004, 108 (46), 17886−17892. (64) Zhu, H.; Paddison, S. J.; Zawodzinski, T., Jr. Effect of Ligand, Support and Solvent on the O2 Binding of Non-Precious Metal Catalysts: An Ab Initio Study. ECS Trans. 2012, 45 (2), 97−107. (65) Sidik, R. A.; Anderson, A. B.; Subramanian, N. P.; Kumaraguru, S. P.; Popov, B. N. O2 Reduction on Graphite and Nitrogen-Doped Graphite: Experiment and Theory. J. Phys. Chem. B 2006, 110 (4), 1787−1793. (66) Bonakdarpour, A.; Lefevre, M.; Yang, R. Z.; Jaouen, F.; Dahn, T.; Dodelet, J. P.; Dahn, J. R. Impact of Loading in RRDE Experiments on Fe-N-C Catalysts: Two-or-Four-Electron Oxygen Reduction? Electrochem. Solid-State Lett. 2008, 11 (6), B105−B108. (67) Jaouen, F.; Charreteur, F.; Dodelet, J. P. Fe-Based Catalysts for Oxygen Reduction in PEMFCs. J. Electrochem. Soc. 2006, 153 (4), A689−A698.

(68) Bowman, D. N.; Jakubikova, E. Low-Spin versus High-Spin Ground State in Pseudo-Octahedral Iron Complexes. Inorg. Chem. 2012, 51 (11), 6011−6019. (69) Liao, M. S.; Watts, J. D.; Huang, M. J. Assessment of the Performance of Density-Functional Methods for Calculations on Iron Porphyrins and Related Compounds. J. Comput. Chem. 2006, 27, 1577−1592. (70) Radon, M.; Pierloot, K. J. Phys. Chem. A 2008, 112, 11824− 11832. (71) Wang, G. F.; Ramesh, N.; Hsu, A.; Chu, D.; Chen, R. R. Density Functional Theory Study of the Adsorption of Oxygen Molecule on Iron Phthalocyanine and Cobalt Phthalocyanine. Mol. Simul. 2008, 34 (10−15), 1051−1056. (72) Sun, S. R.; Jiang, N.; Xia, D. G. Density Functional Theory Study of the Oxygen Reduction Reaction on Metalloporphyrins and Metallophthalocyanines. J. Phys. Chem. C 2011, 115 (19), 9511−9517. (73) Shi, Z.; Zhang, J. Density Functional Theory Study of Transition Metal Macrocyclic Complexes’ Dioxygen-Binding Abilities and Their Catalytic Activities Toward Oxygen Reduction Reaction. J. Phys. Chem. C 2007, 111 (19), 7084−7090. (74) Chen, X.; Li, F.; Zhang, N.; An, L.; Xia, D. Mechanism of Oxygen Reduction Reaction Catalyzed by Fe(Co)-Nx/C. Phys. Chem. Chem. Phys. 2013, 15, 19330.

28562


Sensitivity of Density Functional Theory Methodology for Oxygen

Recommend Documents