Theoretical Investigations of the Electrochemical Reduction of CO on Single Metal Atoms Embedded in Graphene.

Single transition metal atoms embedded at single vacancies of graphene provide a unique paradigm for catalytic reactions. We present a density functional theory study of such systems for the electrochemical reduction of CO. Theoretical investigations of CO electrochemical reduction are particularly challenging in that electrochemical activation energies are a necessary descriptor of activity. We determined the electrochemical barriers for key proton-electron transfer steps using a state-of-the-art, fully explicit solvent model of the electrochemical interface. The accuracy of GGA-level functionals in describing these systems was also benchmarked against hybrid methods. We find the first proton transfer to form CHO from CO to be a critical step in C1 product formation. On these single atom sites, the corresponding barrier scales more favorably with the CO binding energy than for 211 and 111 transition metal surfaces, in the direction of improved activity. Intermediates and transition states for the hydrogen evolution reaction were found to be less stable than those on transition metals, suggesting a higher selectivity for CO reduction. We present a rate volcano for the production of methane from CO. We identify promising candidates with high activity, stability, and selectivity for the reduction of CO. This work highlights the potential of these systems as improved electrocatalysts over pure transition metals for CO reduction.

Introduction

Electrochemical conversion and reduction of CO2 (CO2RR) coupled with a renewable energy source (such as wind or solar) has attracted particular attention as an approach for generating carbon-neutral fuels and industrial chemicals.[1,2] From the perspective of catalyst development, improvements in activity, selectivity over hydrogen evolution, and stability under reducing conditions remain major challenges.[3] Recent candidates such as transition metal dichalcogenides,[4,5] metal-functionalized porphyrin-like structures,[6−8] and metal–organic frameworks[9,10] have demonstrated promise for reducing CO2 to CO. Of the catalysts explored for CO2RR to further reduced hydrocarbons and alcohols, Cu metal and Cu-based catalysts[11−14] have shown the best activity and stability, but these catalysts still require nearly 1 V in overpotential to produce 1 mA/cm2 current density (when normalized to the electrochemically active surface area).

Carbon-based materials are particularly interesting potential electrocatalysts for CO2RR. Due to their low cost and ability to form a wide range of nanostructures, these materials have recently been studied for applications in electrochemical energy conversion and storage, and they represent one promising class of alternatives to commercial precious-metal catalysts.[15−19] Single metal atoms embedded in a conductive graphene network present unique active sites and distinct properties. These materials have been the subject of numerous theoretical[20−23] and experimental studies[21,24] and have been shown to be good catalysts for electrochemical reactions such as oxygen reduction.[25−28]

In this work, we use density functional theory (DFT) calculations to study a wide range of transition metals embedded in single vacancies (M@SV) of graphene for their activity and selectivity toward the electrochemical reduction of CO (CORR). Computational screening of new catalysts for CO2 reduction is particularly challenging in that activation energies for key reaction steps and the corresponding kinetic analyses are necessary to predict CO reduction activity, which presents many open theoretical challenges.[29] Previous work has generally focused on the thermodynamics of the coupled proton–electron transfer steps involved in the reaction which result in necessary but not sufficient criteria for high activity.[29−32] Here, we present the first theoretical screening study of novel materials for COR that uses a fully explicit description of the solvent to determine the energies of critical electrochemical transition states. We present the full reaction pathway  on Pt@SV. We then present scaling relations for the transition state of the rate determining CO to CHO reduction step as well as for the elementary steps for hydrogen evolution. Using a microkinetic analysis, we predict systems with improved activity and selectivity over transition metals.

Computational Details

The simulations were carried out using the Atomic Simulation Environment (ASE).[33] The electronic structure calculations were performed using the Quantum ESPRESSO program package.[34] The electronic wave functions were expanded in series of plane waves with a cutoff energy of at least 500 eV and an electron density cutoff corresponding to 5000 eV for all M@SV systems (for Fe@SV, an energy cutoff of 800 eV and an electron density cutoff of 8000 eV were used). Core electrons were approximated with ultrasoft pseudopotentials.[35,36] We used the BEEF-vdW exchange-correlation functional, which has been shown to accurately describe chemisorption as well as physisorption properties on graphene and on transition metals.[37]

For adsorption energy calculations, a supercell of lateral size (4 × 4) was used, and a vacuum region of at least 18 Å was used to decouple the periodic images. The first Brillouin zone was sampled with a (4 × 4 × 1) Monkhorst–Pack[38]k-point grid. All the atoms were allowed to relax until the maximum force on each atom was below 0.02 eV/Å for each calculation. All calculations were performed with spin polarization. We applied the model developed by Krasheninnikov et al.[20] as an initial guess for the spin state of the metal atom in a graphene single vacancy. We did not find that the energetics we calculated were very sensitive to spin state, as shown in Supplementary Figure 1.

As GGA-level DFT is known to poorly describe the magnetic states, we benchmarked these results against hybrid calculations employing the HSE06 exchange-correlation functional[39] for the two cases of Cu@SV and Pt@SV to ensure that the calculated adsorption energies with GGA level functionals were consistent. Single-point HSE06 calculations were carried out in VASP[40−42] on the geometries optimized in GGA-level DFT. Ultrasoft pseudopotentials[43] were used with a planewave kinetic energy cutoff of 400 eV. The free energy of adsorption was corrected from the electronic energy with zero point energy and entropy contributions estimated using the harmonic approximation.[44]

The CHE model was used to calculate the free energy levels of all adsorbates.[45] In this model, the chemical potential of the proton–electron pair is equal to gas-phase H2,and the electrode potential is taken into account by shifting the electron energy by –eU,where e and U are the elementary charge and the electrode potential relative to RHE, respectively.

The effect of solvation was determined using one monolayer of an icelike water structure[46,47] determined through minima hopping (Supplementary Figure 2).[48,49] To model coupled proton–electron transfers, an excess hydrogen was added to the water monolayer, resulting in a spontaneous separation of charge at the interface where the excess hydrogen atoms transferred electrons to the surface.[50] The entire system is therefore electroneutral, and no compensating homogeneous background charge is added. Transition states were converged using the climbing image nudged elastic band (CI-NEB) method.[51] Forces on the climbing image were converged to <0.02 eV/Å. The plane wave cutoff, electron density cutoff, functional, and calculator parameters were the same as those used in adsorbate geometry optimizations.

The potential dependence of the activation energies was determined using the same assumptions as in Liu et al. for transition metals.[29] Both molecular dynamics simulations of water on Pt(111)[52] and the correlations between surface work functions and potentials of zero charge of metal surfaces suggest the water dipole at the interface to be small close to the potential of zero charge.[53,54] The net dipole for the water structures used in this work was found to be ∼0.8 eV, which is similar to that found for transition metals.[29] This value was subtracted from the calculated work function to correct for the use of an oriented water structure[54,55] The charge extrapolation scheme[56] developed by Chan and Nørskov was used to extrapolate the potential-dependent electrochemical barriers to 4.0 eV, which corresponds to 0 VRHE at pH 7. Transition states were referenced to the aqueous initial state as determined using the CHE.

To evaluate the feasibility of coupled proton-electron transfers, we evaluated the width of the adsorbate-induced states for the H–CO transition state on Pt@SV using a smearing of 0.01 eV and k-point sampling of (16 × 16 × 1) (Supplementary Figure 3). Since the width of the adsorbate-induced states is on the order of 0.1 eV, the Heisenberg uncertainty principle gives an electron lifetime of 10–14 s. Electron transfer is therefore instantaneous on the time scale of proton transfer which occurs on the ns to ps time scale.[57,58]

A steady-state, mean-field microkinetic model was used to determine theoretical turnover frequencies. The Catalysis Microkinetic Analysis Package (CATMAP)[59] was used to solve the set of coupled differential equations that describe the reaction rates. We assume a CO pressure of 1 mbar based on the experimental Faradaic efficiency under CO2 reaction conditions.[29,60]

Results and Discussion

Model System

Transition metals can be trapped in the single and double vacancies of the graphene and carbon nanotube (CNT) structures at various levels of nitrogen doping.[20,21,24,25] Depending on the size of the vacancy, single metal atoms or clusters of several atoms can be localized in nanotubes or graphene layers.[21] For simplicity, we will not consider the possibility of the transition metal in the double vacancy. Stable M@SV systems were found for Cr, Mo, W, Mn, Fe, Ru, Os, Co, Rh, Ir, Ni, Pd, Pt, and Cu. Figure a shows the formation energies of these transition metals doped into a graphene single vacancy site, referenced to the bulk metal and bare graphene with a single vacancy.[27]

Figure 1

(a) Formation energies with respect to bulk chemical potentials for single impurity atoms at a single vacancy on graphene (M@SV). The inset shows the atomic structure of M@SV. (b) Comparison between the projected density of states for the Pt atom in Pt@SV and a Pt surface atom on Pt (111).

Almost all the transition metal dopants form strong bonds with the single vacancy of graphene. This phenomenon has been observed in experiment for Au-doped graphene[21] and has been explored theoretically[20,25,31] for metal-doped graphene and metal-functionalized porphyrin-like graphene. As an example, Figure b shows the projected density of states (PDOS) for the Pt atom in both Pt(111) and Pt@SV. It can be seen that the isolated Pt@SV atom has a very different electronic structure than a surface atom in Pt(111): the Pt states are much more localized, and there is overlap with the carbon states. The unoccupied states corresponding to the narrow peak at ∼1 eV above the Fermi energy (EF) originate from the hybridization of the in-plane d-states of Pt with the sp2 states of its neighboring C atoms. It is also important to note that not only d-states but also s- and p-states of Pt contribute to the hybridization and form sharp peaks close to the Fermi level. Clearly, such a significant amount of unoccupied antibonding states reveals a strong Pt–C bonding. The narrow peak also indicates dangling bonds, which introduces a degree of covalency into the binding with adsorbates and rationalizes the differences in scaling of the associated binding energies, as detailed below.[61]

Benchmarking against HSE06

HSE06 calculations were carried out in VASP[40−42] to evaluate the validity of applying GGA-DFT to these metal-doped graphene systems, where the electronic structure of the active catalytic center may be localized. As illustrated in Figure , truncated systems were used in order to reduce computational time, and we considered Cu and Pt transition metal centers. Overall, we obtained less than 0.2 eV variation between HSE and GGA level functionals, as illustrated in Figure .

Figure 2

Comparison of CHO/CO binding on Cu@SV and Pt@SV with different functionals. Hybrid (HSE06) calculations gave similar energies to GGA (RPBE and BEEF) calculations.

Free Energy Diagram for CO Reduction to CH4

Figure shows a free energy diagram for the reduction of CO to CH4 on Pt@SV. Here, we focus on further reduced products from CO, so we use CO as the reactant species. We considered both pathways of CO going to CHO and COH, and the COH intermediate was found to be much higher in energy. Additionally, we consider both pathways of CHO going to CHOH and CH2O. Although the CHOH and CH intermediates were found to be stable on transition metals, we found that these intermediates were much less stable on the M@SV surfaces than the CH2O and CH2OH intermediates.[29] We therefore expect the reaction to proceed through the CHO and CH2O intermediates as shown in the solid line in Figure .

Figure 3

(a) Free energy diagram for the reduction of CO to CH4 on Pt@SV. All barriers and intermediate binding energies were calculated with an explicit solvent. (b) Configurations of the intermediates along the reaction path on the active site in the presence of the solvent.

Overall, on Pt@SV, high overpotential for the CO2RR originates from the barrier for proton–electron transfer to CO to form CHO:On transition metals, it is the linear scaling between the energies of the H–*CO transition state and *CO that poses a fundamental limit on the activity toward C1 products.[29] Assuming that this step is limiting in all other M@SV surfaces, we focus on this barrier as an indicator of CO reduction activity across our metal-doped graphene surfaces.

Transition State Scaling

Improvements in CO2RR catalytic activity therefore require materials that follow a different scaling of the CO* binding energy and the H–*CO transition state energy that that of transition metals. Herein, we show that single metal active sites in metal-doped graphene behave significantly differently from pure transition metals. Figure a displays the energy of the H–*CO transition state complex vs the binding energy of *CO. All energies are given at 0 VRHE and pH 7. CO does not generally bind without transition metal substitution on N-doped graphene materials, but, as a result of incorporating transition metal substitution, the CO adsorption energy becomes negative in all of the M@SVs considered.[61] It can be seen from Figure a that the free energies of *CO and H–*CO for M@SVs deviate significantly from the scaling relation for transition metals, and a few candidates have much lower transition state energies, approaching an activation energy near 0.7 eV, which corresponds to a turnover frequency of 1 site–1 s–1 at room temperature (298 K) and atmospheric pressure (1 atm). Moreover, for most of the M@SVs, the CO* adsorption energies fall near that of copper in the intermediate binding strength region in the transition metal scaling, which is located closest to the peak of the 211 activity volcano.[29]

Figure 4

(a) Activation energy for eq vs *CO binding energy on M@SV (black), M(111) (blue), and M(211) (green) surfaces. (b) Energy of *CO rotation at the active site for Cu(111), Cu (211), and Co@SV. M(111) and M(211) data is taken from Liu et al.[29]

In order to understand this deviation from transition metal scaling, we have explored the behavior of the H–*CO transition state at the active site. It has been previously demonstrated that the improved scaling of the H–*CO transition state on (211) transition metal surfaces over (111) transition metal surfaces is due to a lowered energy of rotation of the *CO intermediate. Similarly this rotation energy was shown to be lowered on top over 3-fold sites on Ni(111).[62] A rotation of the adsorbed *CO initial state to transition state configuration is less energetically costly on the M@SV surfaces than for Cu(111) or Cu(211) transition metal surfaces, contributing to the improved transition state scaling (Figure b).

Microkinetic Modeling

Microkinetic models are important in the prediction of trends in catalytic activity. Although DFT calculations may not quantitatively predict the absolute rate for a given surface, the trends in activity for a set of surfaces can be well described.[63] As demonstrated by this study, transition state energies and the corresponding reaction rate constants can be important descriptors for predicting activity over thermodynamic intermediate binding energies. Here we develop a steady-state, mean-field microkinetic model to find the potential dependent rate of the reduction of CO (Figure a). We use the scaling relations of binding energies of intermediates with the descriptors of the *CO binding energy and H–*CO transition state energy. In our microkinetic model, only one reaction intermediate is allowed to occupy the single metal atom active site. We do not consider any intermediates binding to any sites other than the metal atom.[61] At any point in our descriptor space, the coverages of the considered intermediates (including the empty sites) always sum to 1.

Figure 5

(a) 2D activity map of CO reduction to CH4 with descriptors of CO binding energy and H–*CO transition state energy at −0.5 VRHE. (b) Activation free energies for highest barriers in CO reduction and hydrogen evolution, as a function of the binding energy of CO. In all cases, the largest barrier for CO reduction is smaller than the largest barrier for hydrogen evolution. M@SV surfaces that include a nitrogen atom near the active site are denoted with a filled circle.

In addition to the rate of CO reduction on a given catalyst, the selectivity over the hydrogen evolution reaction (HER) is an important factor in determining the efficiency of the catalyst. To determine HER activity on the M@SV graphene surfaces, we calculated the transition state energies of the two-step HER.Across all the M@SV systems studied, the first proton adsorption step (eq ) presented a higher activation energy than the formation of H2(g). All HER barriers and binding energies are included in the microkinetic model. The highest barrier for CO reduction is lower than the highest barrier for HER (Figure b). Therefore, we expect these surfaces to be selective for CO reduction. Selectivities for HER and COR are given in Supplementary Figure 5. Energetics of *OH adsorption was also considered and included in the microkinetic model: coverages are shown in Supplementary Figure 6. The M@SV systems that we have explored in this work fall into the region of high *CO coverage due to the higher ΔGH–*CO on these surfaces. Surfaces with the same ΔGCO but much lower ΔGH–*CO could be covered by *OH, but we do not expect *OH coverage for these surfaces. We do not expect *OH poisoning for the surfaces and applied potentials considered here.

As nitrogen doping is considered as a successful strategy for creating vacancies in graphene, we expect that nitrogen may be present near the single atom transition metal site under experimental conditions. There is an interplay between the extent of nitrogen doping and the number of created defects, the effects of which have been explored theoretically[64] and experimentally.[65,66] To explore nitrogen doping in these systems, we calculated the transition state energies of Rh@SV and Co@SV with a single nitrogen atom replacing one of the four carbon atoms coordinated to the metal atom. We find that although the *CO binding energy and the H–*CO transition state energy were altered in the presence of nitrogen compared to the pure metal, the scaling found for the pure M@SV systems was still followed (Figure a).

It is worth mentioning that several experimental efforts have been made to synthesize single metal atoms doped in graphene and nanotube.[21,24,67,68] One typical strategy is to first create vacancies in graphene by energetic Ar ions or nitrogen doping and predeposit metal clusters or in situ evaporate metal impurities.[21,24] Most recently, ball milling[67] and atomic layer deposition[68,69] have been successfully used to experimentally synthesize single Pt- and Pd-doped graphene.

We note that there are many other catalytic sites that have been found to exist experimentally on metal-doped graphene materials. Metal atoms and small clusters can be used to decorate edge sites on graphene,[67] and nanoparticles can be inserted at vacancies in graphene.[21,70,71] At high nitrogen doping, porphyrin-like structures can be formed.[6−8] While the exploration of the activation energies and kinetics on all such sites is out of the scope of this work, we recognize that the presence of different sites could convolute the experimentally measured activity and selectivity. This work proposes that single vacancies in graphene doped with single transition metal atoms can selectively reduce CO. Future experimental work could explore the feasibility of synthesizing such well-defined systems.

Conclusion

In this work, we performed density functional theory calculations to study transition metals embedded in graphene as dopants (M@SV) for the electrochemical reduction of CO. We evaluated the stability of M@SV systems, and we found that all the considered M@SVs form stable structures. The complete reduction pathways were further considered on examples of these surfaces by constructing the lowest energy pathway for CORR to gain insight into catalytic activities, reaction mechanisms, and reaction products. We include activation energies for electrochemical steps determined using a fully explicit description of the solvent. Our results indicate that these materials present markedly different scaling relations between key intermediates in the electrochemical reduction of CO (i.e., *CO and H–*CO transition state) from pure transition metal scaling. A complete microkinetic analysis including competition with the HER was performed. We identified Pd- and Ni-doped graphene as the best candidates with high activity and selectivity for the reduction of CO. This study is, to date, the first computational screening study for CO reduction based on ab initio calculations of electrochemical transition states and associated mean field kinetics. The results point to the potential of these metal-doped graphene systems as highly active and selective CORR electrocatalysts. We note that many different sites can exist in experimentally synthesized metal-doped graphene, which may pose a challenge in the experimental evaluation of the theoretically predicted activity. Future theoretical work will also consider alternative coordinations and substrates for single metal sites and metal clusters.