Structure of Fe–Nx–C Defects in Oxygen Reduction Reaction

Article pubs.acs.org/JPCC

Structure of Fe−Nx−C Defects in Oxygen Reduction Reaction Catalysts from First-Principles Modeling Edward F. Holby,*,† Gang Wu,‡ Piotr Zelenay,‡ and Christopher D. Taylor†,§ †

Materials Science and Technology Division and ‡Materials Physics and Applications Division, Los Alamos National Laboratory, Los Alamos, New Mexico 87545, United States S Supporting Information *

ABSTRACT: The structure of active sites in Fe-based nonprecious metal oxygen reduction reaction catalysts remains unknown, limiting the ability to follow a rational design paradigm for catalyst improvement. Previous studies indicate that N-coordinated Fe defects at graphene edges are the most stable such sites. Density functional theory is used for determination of stable potential oxygen reduction reaction active sites. Clusters of Fe−Nx defects are found to have N-coordination-dependent stability. Previously reported interedge structures are found to be significantly less stable than in-edge defect structures under relevant synthesis conditions. Clusters that include Fe−N3 defects are found to spontaneously cleave the O−O bond.

S

metal atoms in close association and that this facilitates the multielectron reduction steps in ORR.13−15 In this paper, we use density functional theory (DFT) and ab initio molecular dynamics (MD) to characterize the properties of the active sites in NPMCs. We focus on their molecular configuration, surface accessibility, sensitivitiy to N and Fe chemical potentials, and response to an aqueous environment. Particular attention is paid to the clustering tendencies of different N-coordinated structures to illuminate the structure of the MxNy centers. We then explore the role that chemical potentials of Fe and N play in stabilizing/destabilizing such sites. Finally, we consider the adsorption of O2 to various MxNy centers embedded within graphene edge defect sites and the implications for ORR.

uccessful development and possible implementation of nonprecious metal catalysts (NPMCs) for the oxygen reduction reaction (ORR) promise to significantly decrease the materials cost associated with proton exchange fuel cells1 but require detailed characterization of the chemical and structural composition of the ORR active site. A scientific understanding of active site structure will provide synthesis targets that optimize activity, stability, and selectivity, with the ultimate goal of increasing active site density, durability, and efficiency. Graphene nanoribbon (GNR)a two-dimensional carbon sheet with edgeshas previously been employed as a model system for NPMCs2 since it can host a variety of potential ORR active sites and is consistent with the observation of high-levels of sp2 (graphitic) carbon in synthesized NPMCs with high ORR activity. It has been shown that nitrogen, metal, and vacancy defects (which, in some combination, must constitute the active site or sites3) are more stable at edges than in bulk graphene.2,4 The edge nature of these defects is an important requirement for product and reactant accessibility, facilitating mass-transport efficiency. Nanoribbon edge sites also have unique electronic and magnetic properties4,5 that differentiate them from their counterparts in bulk graphene.6−8 To maintain an optimal configuration under operating conditions, it is likely that the active sites are covalently embedded within the planar structure of the graphitic matrix. Nitrogen has been shown to coordinate the nonprecious metal atoms directly,3 but the nature of the active MxNy complexes is actively debated in the literature.6−12 In fact, the stability of any particular MxNy complex will be dependent upon the M- and N-chemical potentials realized during synthesis conditions, a fact which has not been addressed previously in the literature. It has been suggested that the active site is composed of multiple © 2014 American Chemical Society

■

METHODOLOGY To study the formation of FexNy active sites, three geometries are considered using DFT: 2N-coordinated interedge, 3Ncoordinated intraedge, and 4N-coordinated intraedge. These geometries are composed of Fe atoms between the edges of Nterminated zigzag edges, Fe atoms situated above a monovacancy coordinated by 3 N atoms (two of which are at the zigzag edge termination), and Fe atoms in divacancy positions coordinated by 4 N atoms (two of which are at the zigzag edge termination) (see Figure 1). A graphene nanoribbon eight-C-pairs-long and five-C-pairs-wide, with two FeNx defects of a given geometry, is used for the 3N and 4N cases, allowing for possible defect clustering distances of 1, 2, 3, and 4 Received: April 2, 2014 Revised: June 5, 2014 Published: June 6, 2014 14388

dx.doi.org/10.1021/jp503266h | J. Phys. Chem. C 2014, 118, 14388−14393

The Journal of Physical Chemistry C

Article

implemented in the Vienna ab initio simulation package (VASP) code.16−19 The generalized gradient approximation, as parametrized by Perdew, Burke, and Erzerhoff (PBE),20,21 is used for the exchange and correlation component of the energy calculations. Monkhorst−Pack k-point meshes of 5 × 1 × 1 (five along the length of the GNR) are used as in previous work on graphene nanoribbon systems.2 Large vacuum spaces (∼20 Å) are introduced between ribbons and their periodic images to limit undesired self-interactions. Formation energies are calculated as follows Gdefect = Edefect + ΔnCμC + ΔnHμH − ΔnNμ N − ΔnFeμFe − EGNR

(1)

where μx represents the chemical potential of species x; Edefect is the internal energy calculated for the specific Fex−Ny structure; and EGNR is the internal energy calculated for the Fe and N free nanoribbon. The reference state energies employed for the chemical potentials are the energy of a single C atom in bulk graphene for C, 1/2 energy of the N2 molecule for N, the energy of a single Fe atom in bulk BCC iron for Fe, and 1/2 the energy of the H2 molecule for H. Δnx represents the change in the number of species x in the defected vs the defect-free nanoribbon. In eq 1, E is the ground state internal energy value output from VASP. Due to the exploratory nature of this study, entropy and zero-point energy (ZPE) are not included (and are unlikely to affect the obtained trends). Other factors that will contribute to accurate energetics include solvation and dispersion effects. Ab initio molecular dynamics simulations are performed with temperature scaling and the Verlet algorithm, as implemented in VASP. Time steps of 1 fs are used. The temperature is held constant at 300 K. MD simulations are used to investigate the response of the in vacuo active sites to additional opportunities for coordination with O2 and water under aqueous and thermalized conditions. To accelerate the molecular dynamics simulations, coarser energy cut-offs and k-point meshes are used. Geometry optimizations of final snapshots from the MD simulations use the same parameters as the optimization calculations mentioned above.

■

RESULTS A. Active Site Clustering from DFT. The relaxed structures for each of the different geometries obtained at different Fe−Fe spacing are shown in Figure 1. The spacing between Fe−Nx clusters is increased, and formation energies are calculated following eq 1. Calculated Fe−Fe distances and formation energies are given in Table 1. The formation energies with respect to the nearest neighbor cluster (clustering energies) as a function of Fe−Fe spacings are plotted in Figure 2 for the three considered geometries. As can be seen from Figure 2, the different active site geometries have distinct formation energies as a function of Fe−Nx defect spacing. Both the 2N interedge and 3N geometries show a drastic thermodynamic preference for Fe−Fe clustering, while the 4N geometry formation energy is relatively insensitive to Fe−Fe spacing. At higher temperatures, entropic effects will alter these formation energies. The N2 reference molecule, with higher entropy due to increased degrees of freedom as a gas, will be more stable at higher temperatures, leading to an effective stabilization of lower-N structures (such as the NN

Figure 1. Energy minimized structures for (a) 2N interedge geometry, (b) 3N geometry, and (c) 4N geometry. Each geometry shows nearest neighbor (NN), 2nd nearest neighbor (2NN), 3rd nearest neighbor (3NN), and 4th nearest neighbor (4NN) clusters, respectively. Note that (a), (b), (c), and (d) are periodic horizontally and (a) is also periodic vertically.

carbon-pair units situated along the zigzag edge. In the case of the 2N-coordinated structure, lower formation energies are found by using a six-C-pair-wide ribbon. Such an arrangement allows the Fe to sit between two edge N atoms (Figure 1a), which is found to be more stable than the staggered arrangement realized with a five-C-pair-wide ribbon. By calculating formation energies of the defect clusters at all possible Fe−Fe spacings, the clustering tendency of defects of a given geometry can be evaluated. Formation energies are computed using density functional theory (DFT), as 14389

dx.doi.org/10.1021/jp503266h | J. Phys. Chem. C 2014, 118, 14388−14393

The Journal of Physical Chemistry C

Article

poor conditions. Figure 3 shows the free energy of each lowestenergy structure for each of the three considered geometries as

Table 1. Structure Fe−Fe Bond Distances and Formation Energies at Reference State Fe and N Conditions structure FeN2 FeN2 FeN2 FeN2 FeN3 FeN3 FeN3 FeN3 FeN4 FeN4 FeN4 FeN4

NN 2NN 3NN 4NN NN 2NN 3NN 4NN NN 2NN 3NN 4NN

Fe−Fe distance (Å)

formation energy (eV)

2.17 4.89 7.34 9.83 2.14 4.83 7.88 9.88 2.25 5.05 7.46 9.85

9.52 11.3 11.1 11.1 2.37 3.00 3.11 3.86 0.75 0.69 0.82 0.76

Figure 3. Stability plot of different N and Fe/N structures as a function of N chemical potential.

a function of chemical potential of N. Additionally, an isolated N-defect (N replacing carbon at the graphene edge) and N-free nanoribbon are also shown for comparison. Based in Figure 3, a purely thermodynamic interpretation suggests that there are three different stability regions. At low chemical potential of N, the defect-free nanoribbon is most stable. At relative chemical potential of N between −0.72 and 0.21 eV, a single edge N atom is the most stable structure. Above 0.21 eV, the separated FeN4 geometry is found to be most stable. In addition to the three stable structures found for the bulk Fe reference state, a fourth structure (the FeN3 clustered structure) is stable under certain synthesis conditions with high chemical potential of Fe and low chemical potential of N. It is likely that different synthesis conditions also lead to variations in the Fe chemical potential. As such, it is possible to plot a predominance diagram of phases as a function of chemical potentials of N and Fe. Since all structures have a fixed number of Fe and N, each structure has a formation energy plane in the Fe and N chemical potential space. The slope of this plane in a given dimension is dictated by the number of Fe or N atoms in the structure. The structure with the lowest-energy plane at a given chemical potential of Fe and N is the most thermodynamically stable of the considered structures at those conditions. The calculated difference in formation energy between a free Fe atom and one bound in a BCC lattice is ∼5 eV, and so an Fe chemical potential range from 0−5 eV is considered. The stability regions of different structures are plotted in Figure 4. One important finding that comes from considering formation energies as a function of chemical potential of N and Fe is just how unstable the 2N interedge structures are compared to the other considered geometries (∼8.77 eV less stable than the lowest-energy FeN4 structure). While every possible interedge structure cannot be accounted for, care was taken to minimize the energy of the 2N interedge structures, including the use of a wider nanoribbon to ensure proper registry of the zigzag edges. It is found that having the edge N

Figure 2. Clustering energies as a function of Fe−Fe spacing following the structures depicted in Figure 1.

clustered 3N and 4N structures) versus higher-N, non-nearestneighbor structures. These results imply that the tendency of a given active site geometry to form multimetallic (i.e., clustered) configurations depends minimally on whether it is 2N, 3N, or 4N. Multimetallic active sites are more likely to catalyze ORR through unique pathways (as we demonstrate in Section C). Furthermore, a potential pathway for tunable catalysis is opened by specifically designing heterometallic clusters. This phenomenon may explain the increased activity of synthesized Fe−Co catalysts.22 B. Thermodynamic Phase Stability from DFT. By fixing the C, Fe, and H chemical potential and varying that of N, the formation energies of the different geometries under varied synthesis conditions (varied N-source such as N2 or NH3 and varied activities/partial pressures of the N-source reactant) are compared. As the clustered Fe structures have less N than their separated Fe counterparts (for the 3N and 4N cases) and different geometries have varied amounts of N, the chemical potential of N is shown to have an important role in determining which structures are most stable at different synthesis conditions. By plotting formation energies as a function of N chemical potential (following eq 1 with the N2 reference state set as μ = 0 eV) it is possible to compare relative active site structure thermodynamic stability under N-rich or N14390

dx.doi.org/10.1021/jp503266h | J. Phys. Chem. C 2014, 118, 14388−14393

The Journal of Physical Chemistry C

Article

adsorption process. O2 adsorption structures are presented in Figure 5.

Figure 4. Stability plot of different N and Fe/N structures as a function of chemical potential of N and Fe. Fe chemical potentials greater than zero are possible for Fe reference states less stable than bulk Fe (such as iron salts or lone Fe atoms). N chemical potentials less than zero are possible when entropic and partial pressure considerations stabilize a N2 molecule vs the 0 K molecule reference state utilized.

directly across from another edge N, as shown in Figure 1(a) with a 6-C wide ribbon, lowers the system energy versus the edge N not being directly across from another edge N as is the case with a 5-C wide ribbon. Removing the edge-passivating H atoms is found to result in recombination of the zigzag edge, thus effectively turning the edge defect into a bulk defect (at the center of an infinite graphene sheet). In previous work, such defects were shown to be considerably less stable than edge defects.2 This finding is in contrast to the previously proposed NPMC active site model with Fe between graphene edges. Such structures are considerably less stable than other defects occurring at graphene edges and, as such, are unlikely candidates for NPMC active sites. C. Relation between O2 Dissociation and Active Site Geometry from DFT. The initial step in oxygen reduction is the adsorption of oxygen to the active site. It is therefore instructive to consider this adsorption step and compare the behavior of different active sites. DFT and geometry optimization indicate that placing an O2 molecule above a clustered FeN3 site (i.e., Fe2N5) leads to cleavage of the O2 bond. This spontaneous behavior is confirmed using the nudged elastic band method,23 which shows that the bond breaking reaction is barrierless. The possibility for the two undercoordinated Fe to each bind an adsorbed O atom appears to provide the strong driving force for O2 dissociation. Since oxygen dissociation is spontaneous, it is likely that a dissociative ORR pathway will proceed on clustered bimetallic FeN3 structures. In the case of the lone FeN3 and FeN4 sites, as well as the clustered FeN4 site, the O2 binds to the Fe atom(s), and addition of an H atom leads to the formation of OOH without dissociating the O−O bond. This suggests that an associative mechanism is likely at such sites. This mechanism is generally less preferable as it can result in the undesired formation of H2O2 instead of H2O. A mixed bimetallic cluster FeN3/FeN4 site also led to O2 dissociation suggesting that the presence of the FeN3 structure aids in the dissociative

Figure 5. Optimized geometries for (a) FeN3 structure with adsorbed O2; (b) FeN4 structure with adsorbed O2; (c) bimetallic FeN3 cluster structure with adsorbed O2; (d) bimetallic FeN4 cluster structure with adsorbed O2; (e) bimetallic FeN3/FeN4 cluster structure with adsorbed O2. Each structure is shown from above and from the side.

These findings show that active-site geometry plays a key role in ORR pathway and hence the four- vs two-electron selectivity. Multimetal defects at nanoribbon edges are postulated to be ideal active sites for reactions where initial bond breaking is an important step (such as ORR and, potentially, electrochemical ammonia synthesis). D. Ab Initio Molecular Dynamics Simulations of the Response of the Active Site to Solvation. The metal− nitrogen clusters determined from first-principles thermodynamics may or may not be stable upon exposure to an aqueous environment. To address this question directly, ab initio molecular dynamics of this system is performed in the presence of solvent molecules and at thermal temperatures (300 K). Water molecules are placed in a grid up to a density of 1.0 g/ cm3 in the space exclusive of the van der Waals volume around the graphene sheet. The system is then gradually heated to 300 K over a series of 0.5 ps, and then this temperature is maintained for a total of 2.5 ps. A time step of 1 fs is used for these simulations. While this is not sufficient time to reach complete equilibration of the catalyst/water interface (roughly 14391

dx.doi.org/10.1021/jp503266h | J. Phys. Chem. C 2014, 118, 14388−14393

The Journal of Physical Chemistry C

■

Article

CONCLUSIONS Using DFT, a range of potential NPMC ORR active site/ graphene nanoribbon edge defect geometries are, for the first time, compared in a self-consistent fashion. On the basis of the formation energies of these FexNy sites, it is found that they are thermodynamically driven to cluster. Furthermore, depending upon synthesis conditions, either the 4N or 3N structures are the most stable Fe-containing defects. It is postulated that both these structures may coexist. The 2N interedge structures are so much higher in energy than the other defects studied that their existence is postulated to be either thermodynamically prohibited or metastable due to the kinetics of their decomposition. Clustered FeN3 structures (Fe2N5) seem to excel at cleaving the O2 bond with zero barrier and thus are likely to follow a dissociative pathway for ORR. This pathway is expected to be more selective, avoiding formation of H2O2 at the cost of potentially overbinding ORR intermediates. Ab initio molecular dynamics indicates that this spontaneous reaction is likely to be unaffected by solvation. Solvent does not appear to affect the stability of these edge defects. Future work will consider the full ORR pathway.

100 ps is required for this), it is sufficient time to observe reconstruction of water molecules about the active site. The results of the ab initio molecular dynamics simulations are shown in Figure 6. The water molecules are equilibrated

■

ASSOCIATED CONTENT

S Supporting Information *

Relaxed atomic positions for the structures in Figures 1 and 5. This material is available free of charge via the Internet at http://pubs.acs.org.

■

AUTHOR INFORMATION

Corresponding Author

*Phone: (505) 665-0034. E-mail: [email protected]. Present Address §

Strategic Research & Innovation, DNV GL, 5777 Frantz Road, Dublin, OH 43017 (C.D.T.) 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.

■

Figure 6. Final, geometry-optimized snapshots taken from the ab initio molecular dynamics simulations of the Fe2N5 cluster in the presence of water showing (a) coordination to a constrained O2 molecule held over the active site and (b) the optimized 2O* state, indicating spontaneous dissociation starting from state (a).

ACKNOWLEDGMENTS The authors wish to thank the Los Alamos National Laboratory for funding under the Laboratory Directed Research and Development (LDRD) program and for institutional computing resources. Los Alamos National Laboratory is operated by Los Alamos National Security LLC for the National Nuclear Security Administration of the U.S. Department of Energy under contract DE-AC528-06NA25396. This work used the Extreme Science and Engineering Discovery Environment (XSEDE), which is supported by National Science Foundation grant number OCI-1053575. This work was facilitated by CNM User Proposal 28858. Use of the Center for Nanoscale Materials was supported by the U.S. Department of Energy, Office of Science, Office of Basic Energy Sciences, under Contract No. DE-AC02-06CH11357.

over the Fe2N5 active site, and an O2 molecule is introduced to replace the nearest water molecule over the active site. The geometry of the active site and the water molecule is held fixed while the other degrees of freedom are relaxed, and then the active site itself is relaxed while the oxygen molecule is held fixed, leading to the geometry shown in (a) of Figure 6. When the constraint on O2 is relaxed, O2 spontaneously dissociates in the presence of water, just as is observed for the in vacuo case. The energy released by dissociation of the O2 bond is 3.96 eV (1.98 eV/O atom). Further investigation of the potential energy surface for ORR using a combination of ab initio molecular dynamics and geometry optimization in these solvent environments is deferred to a separate paper. In all cases, however, the active site geometry is maintained and does not appear to be immediately unstable in an aqueous environment.

■

REFERENCES

(1) Chen, Z.; Higgins, D.; Yu, A.; Zhang, L.; Zhang. A Review On Non-Precious Metal Electrocatalysts For PEM Fuel Cells. J. Energy Environ. Sci. 2011, 4, 3167−3192.

14392

dx.doi.org/10.1021/jp503266h | J. Phys. Chem. C 2014, 118, 14388−14393

The Journal of Physical Chemistry C

Article

(2) Holby, E. F.; Taylor, C. D. Control Of Graphene Nanoribbon Vacancies By Fe And N Dopants: Implications For Catalysis. Appl. Phys. Lett. 2012, 101, 064102−1−064102−4. (3) Wu, G.; Chen, Z.; Artyushkov, K.; Garzon, F. H.; Zelenay, P. Polyaniline-derived Non-Precious Catalyst for the Polymer Electrolyte Fuel Cell Cathode. ECS Trans. 2008, 16, 159−170. (4) Longo, R. C.; Carrete, J.; Ferrer, J.; Gallego, L. J. Structural, Magnetic, And Electronic Properites Of Nin And Fen Nanostructures (N = 1−4) Adsorbed On Zigzag Graphene Nanoribbons. Phys. Rev. B 2010, 81, 115418−1−115418−10. (5) Ritter, K. A.; Lyding, J. W. The Influence Of Edge Structure On The Electronic Properties Of Graphene Quantum Dots And Nanoribbons. Nat. Mater. 2009, 8, 235−242. (6) Kattel, S.; Atannassov, P.; Kiefer, B. Catalytic Activity Of Co-Nx/ C Electrocatalysts For Oxygen Reduction Reaction: A Density Functional Theory Study. Phys. Chem. Chem. Phys. 2013, 15, 148−153. (7) Calle-Vallejo, F.; Martinez, J. I.; Rossmeisl, J. Density Functional Studies Of Functionalized Graphitic Materials With Late Transition Metals For Oxygen Reduction Reactions. Phys. Chem. Chem. Phys. 2011, 13, 15639−15643. (8) Studt, F. The Oxygen Reduction Reaction on Nitrogen-Doped Graphene. Catal. Lett. 2013, 143, 58−60. (9) Titov, A.; Zapol, P.; Kral, P.; Liu, D.-J.; Iddir, H.; Baishya, K.; Curtiss, L. A. Catalytic Fe-xN Sites in Carbon Nanotubes. J. Phys. Chem. C 2009, 113, 21629−21634. (10) Lee, D. H.; Lee, W. J.; Lee, W. J.; Kim, S. O. Theory, Synthesis, and Oxygen Reduction Catalysis of Fe-Porphyrin-Like Carbon Nanotube. Phys. Rev. Lett. 2011, 106, 175502-1−175502-4. (11) Chen, R.; Li, H.; Chu, D.; Wang, D. Unraveling Oxygen Reduction Reaction Mechanism on Carbon-Supported Fe-Phthalocyanine and Co-Phthalocyanine Catalysts in Alkaline Solutions. J. Phys. Chem. C 2009, 113, 20689−20697. (12) Kattel, S.; Atanassov, P.; Kiefer, B. Stability, Electronic and Magnetic Properties of In-Plane Defects in Graphene: A FirstPrinciples Study. J. Phys. Chem. C 2012, 116, 8161−8166. (13) Holby, E. F.; Wu, G.; Zelenay, P.; Taylor, C. D. Metropolis Monte Carlo Search for Non-Precious Metal Catalyst Active Site Candidates. ECS Trans. 2012, 50, 1839−1845. (14) Tributsch, H.; Koslowski, U. I.; Dorbrandt, I. Experimental And Theoretical Modeling Of Fe-, Co-, Cu-, Mn-Based Electrocatalysts For Oxygen Reduction. Electrochim. Acta 2008, 53, 2198−2209. (15) Yeager, E. Dioxygen Electrocatalysis: Mechanisms In Relation To Catalyst Structure. J. Mol. Catal. 1986, 38, 5−25. (16) Kresse, G.; Hafner, J. Ab Initio Molecular-Dynamics Simulation Of The Liquid-Metal-Amorphous-Semiconductor Transition In Germanium. Phys. Rev. B 1994, 49, 14251−14269. (17) Kresse, G. Furthmuller, Efficiency Of Ab-Initio Total Energy Calculations For Metals And Semiconductors Using Plane-Wave Basis Set. J. Comp. Mater. Sci. 1996, 6, 15−50. (18) Kresse, G.; Furthmuller, J. Efficient Iterative Schemes For Ab Initio Total-Energy Calculations Using Plane-Wave Bsis Set. Phys. Rev. B 1996, 54, 11169−11186. (19) Kresse, G.; Hafner, J. Ab Initio Molecular Dynamics For Liquid Metals. Phys. Rev. B 1993, 47, 558−561. (20) Perdew, J. P.; Burke, K.; Ernzerhof, M. Generalized Gradient Approximation Made Simple. Phys. Rev. Lett. 1996, 77, 3865−3868. (21) Perdew, J. P.; Burke, K.; Ernzerhof, M. Errata: Generalized Gradient Approximation Made Simple. Phys. Rev. Lett. 1997, 78, 1396−1396. (22) Wu, G.; More, K. L.; Johnston, C. M.; Zelenay, P. HighPerformance Electrocatalysts For Oxygen Reduction Derived From Polyaniline, Iron, And Cobalt. Science 2011, 332, 443−447. (23) Jonsson, H.; Mills, G.; Jacobsen, K. W. Classical and quantum dynamics in condensed phase simulations; Berne, B. J., Cicotti, G., Coker, D. F., Ed.; World Scientific: Singapore, 1998; pp 385−404.

14393

dx.doi.org/10.1021/jp503266h | J. Phys. Chem. C 2014, 118, 14388−14393