Open Access Article
This Open Access Article is licensed under a
Creative Commons Attribution 3.0 Unported Licence

Mapping free energy regimes in electrocatalytic reductions to screen transition metal-based catalysts

Srinivasan Ramakrishnan , Ross A. Moretti and Christopher E. D. Chidsey *
Department of Chemistry, Stanford University, Stanford, CA 94305, USA. E-mail:

Received 11th April 2019 , Accepted 26th June 2019

First published on 27th June 2019


The free energy landscape of catalytic intermediates in the two-electron reduction of proton donors and/or CO2 to H2, CO and HCO2 is mapped with density functional theory to screen catalyst candidates from a library of different transition metals and ligands. The goal is to minimize the free energy corrugations between reactants, catalytic intermediates and each desired product, simultaneously screening against intermediates with low free energy that would be traps, and against necessary intermediates with high free energy. Catalysts are initially screened for those with: (a) standard state free energy of the metal hydride intermediate ergoneutral with HCO2, which is the lowest energy product with weak proton donors, and (b) standard free energy of the metal carbonyl intermediate sufficiently high to avoid trapping. The design method is tested on a diverse range of ligands including cyclopentadienyl, polypyridyl, amino, phosphino and carbonyl ligands, around three earth-abundant d6 transition metal ions, Mn(I), Fe(II) and Co(III), using the BP86 density functional, the double-zeta 6-31+G* basis, LANL2DZ effective core potential on the metals and SMD polarizable continuum model for acetonitrile as solvent, which have previously provided chemically accurate values of several redox potentials, pKa's and ligand exchange equilibria for transition metal complexes. Among the 36 complexes screened, an Fe(II) center ligated to two bipyridines and a pyridine with a solvent-bound sixth coordination site for hydride formation from phenol as the proton donor is identified as a promising candidate for ergoneutral hydride formation without trapping by CO. The redox-active bipyridine ligands are predicted to provide near ergoneutral sites for accumulating the two electrons needed to form the hydride. To test the predictions, an Fe(II) complex was prepared with the desired ligand environment using a pentadentate ligand to prevent ligand exchange. The synthesized complex was indeed found to be active towards electrocatalytic proton reduction as well as CO2 reduction at the predicted redox potentials with no trapping by CO. However, contrary to the in silico predictions, we found electrochemical evidence of CO2 binding after the first reduction leading to CO production. Mapping the free energies of key catalytic intermediates such as the metal hydride and metal carbonyl species by using density functional theory (DFT) serves as a first step in catalyst screening spanning large libraries of metals and ligands. In order to screen against all the intermediates in the catalytic pathway, such as reduced metal-bound CO2 intermediates, further refinement and validation of the DFT methods are needed.


The multi-electron reduction of abundant feedstocks such as protons and CO2 to fuels and value-added products1–3 requires catalysts that can rapidly mediate these reductions with high selectivity and energy efficiency. Several earth-abundant heterogeneous electrocatalysts such as metallic copper can catalyze CO2 reduction, but often suffer from poor selectivity and energy efficiency.4–7 Furthermore, the challenge of knowing the exact chemical nature of the catalytic active site on metal surfaces makes the rational optimization of such catalysts difficult. On the other hand, with transition metal complexes, one can in principle tune the relative energies of the catalytic intermediates by simply changing the metal and ligand combination, choosing from a vast library to favor one stoichiometric product over another. While heterogeneous catalysts may continue to be more practical, molecular complexes serve as excellent model systems to advance single-site catalyst design.

The multi-electron reduction of CO2 to methanol or C2 products is the ultimate goal, but currently no strategies exist to precisely design catalysts that can achieve these transformations rapidly with optimal selectivity. Therefore, as a first step towards this process, we choose the two-electron reduction of protons to H2 and of CO2 to HCO2 and CO (eqn (1)–(3)) as the target reactions for catalyst design, given the availability of mechanistic information for these individual reactions. Furthermore, these transformations involve common intermediates, which provide useful handles to reduce the complexity of the design process. Interestingly, there are many more known single-site molecular electrocatalysts for H2 and CO production8–11 than for HCO2 production,12–17 and even the few examples of the latter involve precious metals and require high overpotentials.

HA + CO2 + 2e → HCO2 + A(1)
2HA + 2e → H2 + 2A(2)
HA + 2CO2 + 2e → CO + HCO3 + A(3)

In order to optimize catalytic activity, our two main strategies involve (a) establishing the relative free energies of the stoichiometric products for a given proton donor and (b) minimizing the free energy corrugation of the catalytic intermediates relative to zero driving force as illustrated in Scheme 1. The latter strategy is central to achieving the desired catalytic activity as intermediates that are high energy or deep trap states in the catalytic pathway will deactivate the catalysts. While the relative rates of different processes will be ultimately governed by the relative heights of transition states, the first screen of catalysts described in this work aims to understand how free energy corrugations of key catalytic intermediates are determined by the transition metal and the ligand environment.

image file: c9sc01766f-s1.tif
Scheme 1 Catalyst design strategies.

A significant recent development in transition metal-mediated catalysis has been the application of density functional theory-based computations to model the reactivity of transition metal complexes.18–20 DFT-based methods have proven to predict with reasonable accuracy redox potentials, pKa's and ligand-exchange equilibria of transition metal complexes in aprotic polar solvents that are well described by simple polarizable continuum models.18,21–24 These methods have been employed to model post-facto the catalytic landscape in a single family of catalysts with relatively small perturbations to the ligand field.25–27 While the resulting changes in reactivity are extremely useful for the system under study, there is little insight offered into the reactivity of a larger set of complexes. Accelerating the screening of metal and ligand choices towards predicting the optimal catalyst candidate with such methods requires the identification of suitable thermodynamic and kinetic descriptors.

The free energy of the transition metal hydride intermediate relative to reactants, products and other intermediates, an important thermodynamic quantity that determines the driving force towards hydride transfer to CO2 and protons to produce HCO2 and H2 respectively, has been well studied as an important descriptor of catalytic activity.28–31 For instance, we previously showed, using experimentally calibrated DFT, how specific ligand environments change the free energy of CO2 insertion into a ruthenium hydride complex.24 The hydride complex we were studying despite being an active transfer hydrogenation catalyst for ketones32 was not found to be a catalyst for CO2 reduction due to product inhibition. In another case, a cyclopentadienyl Ru complex33 was predicted by using DFT to have near ergoneutral CO2 insertion into the corresponding Ru-H intermediate. However, CO2 coordinated to the singly reduced Ru-center prior to any hydride formation resulting in a Ru-bound CO intermediate that was a trap even after a third electrochemical reduction preventing turnover at reasonable potentials. This example highlighted the need to model the energetics of off-path intermediates during catalyst design. While in certain cases, the metal-bound CO can be labilized with a third reduction,34 the effect of the metal or the nature of reduction (ligand vs. metal-centered) on the dissociation is not well understood. Perhaps the most well-studied metal hydride complex for CO2 reduction is the [Ru(tpy)(bpy)H]+ complex (tpy = 2,2,6,6-terpyridine, bpy = 2,2-bipyridine) by Creutz and others.28,35–37 While experimental as well as computational studies have unequivocally established that the hydride donating ability of this complex is well suited for ergoneutral CO2 insertion to produce HCO2, under electrocatalytic conditions, CO is the major product.38 Therefore, the intermediates leading to CO production thermodynamically and kinetically compete with metal hydride chemistry, motivating the need for computational catalyst design approaches to go beyond hydricity as a single descriptor.

In this work, we explore a diverse library of transition metals and ligands and show the broader utilization of DFT in the catalyst design process. We compute the standard state free energies of the metal hydride and metal-bound CO intermediates relative to the stoichiometric reactants and products for metal complexes with diverse ligand environments around three earth-abundant d6 transition metal ions (Mn(I), Fe(II) and Co(III)), to screen for optimal metal–ligand combinations. This coarse level of screening based on just two standard state thermodynamic descriptors, using the BP86 density functional,33,39 weeds out bad candidates effectively and identifies two bipyridines (bpy) and a pyridine (py) ligand around Fe(II) as promising candidates among the 36 complexes studied. Fortuitously the two redox-active bipyridyl units are predicted to provide optimal reduction potentials for accumulating the two electrons needed to form the hydride. Because the simple [FeII(bpy)2(py)(CH3CN)]2+ complex would not be experimentally stable due to ligand scrambling at ambient temperatures, a pentadentate iron complex, [Fe(bpy2PYMe)S]2+ (bpy2PYMe = 1-(2-pyridyl)-1-bis-(6-2,2′-bipyridyl)ethane, S = CH3CN), reported previously by Long et al.40 was synthesized. [Fe(bpy2PYMe)S]2+ was indeed found to electrocatalytically reduce protons as well as CO2 at modest overpotentials with no catalyst trapping by CO. While the calculated reduction potentials and the free energies of the metal hydride and carbonyl intermediates are validated by the experimental results, the binding energy of CO2 to the singly reduced Fe complex was found to be energetically favorable in contrast to the predictions, resulting in the uphill formation of CO instead of the thermodynamically favored HCO2. This work, therefore, provides an efficient first step of catalyst design, by computationally screening against bad candidates via the mapping of free energies of key intermediates in the catalytic pathway. More work is needed to further narrow the candidate space through better modelling of transition metal–CO2 interactions as a function of the redox state of the transition metal complex. To the best of our knowledge, this is the first reported application of DFT for the a priori design of molecular electrocatalysts.

Results and discussion

Thermodynamics of the stoichiometric reactions

For the two-electron reduction of CO2 to HCO2 and of protons to H2 in the standard states of reactants and products, the relative driving force in acetonitrile as a function of the pKa of the stoichiometric proton donor, HA, is shown in eqn (4). They are based on the thermochemical analyses shown in Table S1.
image file: c9sc01766f-t1.tif(4)

Given the difference in proton stoichiometry between eqn (1) and (2), plotting image file: c9sc01766f-t2.tifversus the pKa of the proton donor HA in acetonitrile yields a straight line that crosses zero at a pKa of ca. 24 (Fig. S1). With weaker acids in acetonitrile (pKa > 24) there will therefore be a thermodynamic bias for CO2 reduction to HCO2versus proton reduction.41 For example, phenol (pKa = 29.14)42 should disfavor H2 while acetic acid (pKa = 23)42 should favor H2.

Under weakly acidic conditions, the two-electron, one-proton reduction of CO2 to CO produces bicarbonate (eqn (3)). Due to the lack of experimental free energies for this reaction in acetonitrile, we estimate an upper bound based on the two-proton reduction of CO2 to CO and water (eqn S1). We estimate this latter process (eqn S1) to be roughly uphill by ca. 4.8 kcal mol−1 relative to H2 production under water levels of ∼1 mM in acetonitrile (eqn S2), in agreement with reported values.43

Catalytic intermediates in the two-electron reductions

The reactants, products and catalytic intermediates involved in the two-electron reduction of CO2 and protons as expected from chemical precedent are shown in Scheme 2, mediated by a generic transition metal complex [M–S]n+, where M is a transition metal in a six-coordinate ligand environment in which S is a labile solvent ligand (CH3CN in this case). The chemical potential of electrons is set to a value of −1.6 V vs. Fc+/Fc with phenol as the proton donor that makes CO2 and HCO2 ergoneutral in the standard state, H2 at a relative free energy of 6.7 kcal mol−1 (eqn (4)), and CO with an upper bound of ca. 11.5 kcal mol−1 (see the ESI).44
image file: c9sc01766f-s2.tif
Scheme 2 Two-electron reduction pathways for a generic transition metal electrocatalyst [M–S]n+, HA = proton donor. The free energies of catalytic intermediates are chosen arbitrarily.

Two-electron reduction of the metal complex [M–S]n+ with the dissociation of S is followed by protonation at the metal center to form the metal hydride intermediate. CO2 insertion into the metal hydride bond leads to HCO2 while direct protonation of the metal hydride yields H2. Additionally, the singly or doubly reduced metal complex, [M](n−1)+ and [M](n−2)+ can directly bind CO2 at the metal center33 subsequently leading to the potential release of CO via an [M–CO]n+ intermediate.

In silico screening of catalyst candidates

Within this class of reactivity, we first calculate two key thermodynamic quantities, viz. the free energy of the metal hydride image file: c9sc01766f-t3.tif and the metal carbonyl intermediates image file: c9sc01766f-t4.tif relative to the thermodynamically favorably product with weak proton donors, HCO2. Given the same proton stoichiometry in eqn (1) and (3), the comparison of the two thermodynamic quantities image file: c9sc01766f-t5.tif and image file: c9sc01766f-t6.tif relative to HCO2 is independent of the choice of the proton donor, which will merely shift the relative energies of the stoichiometric products (Table S3). Scheme 3 lists the families of complexes we chose for a comparison of their reactivity within our design framework. They encompass d6 earth-abundant metal ions (Mn(I), Fe(II) and Co(III)) and a diverse set of strong field ligands including polypyridines, phosphines, amines and carbonyls to enhance metal–ligand binding.
image file: c9sc01766f-s3.tif
Scheme 3 Families of d6 transition metal complexes screened (S = CH3CN, bpy = 2,2′-bipyridine, dppe = bis(diphenylphosphino) ethane, and en = 1,2-ethylene diamine).

The two thermodynamic quantities image file: c9sc01766f-t7.tif and image file: c9sc01766f-t8.tif relative to HCO2 are plotted versus each other in Fig. 1 for all the complexes studied. These values are independent of the choice of the proton donor due to the equal proton stoichiometry in eqn (1) and (3). For selective and ergoneutral HCO2 production, the optimal value of image file: c9sc01766f-t9.tif is zero and image file: c9sc01766f-t10.tif is greater than zero, and vice versa for CO production. image file: c9sc01766f-t11.tif values greater than 10 kcal mol−1, for example, would inhibit the formation of [M–CO]n+ at room temperature.45

image file: c9sc01766f-f1.tif
Fig. 1 2D plot of the calculated values of image file: c9sc01766f-t17.tif and image file: c9sc01766f-t18.tif (in kcal mol−1 units) for the d6 transition metal complexes listed in Scheme 3. Level of theory: BP86/LANL2DZ(M)/6-31+G*(C,H,N,O,P,Cl)/SMD(MeCN).

First, independent of the choice of metal, it is evident that all the data points fall roughly on a straight line for the cyclopentadienyl, pincer and bis-bipyridine frameworks. The slopes of the linear fits are all about −0.5, meaning that as the free energy of [M–CO]n+ decreases and the free energy of [M–H](n−1)+ goes up. This is reasonable as an increase in electron density at the metal center is expected to increase the [M–H](n−1)+ free energy but lower the [M–CO]n+ free energy due to increased back-bonding into the π* orbital of the carbonyl ligand. For the same metal, however, changing the ligands causes minor deviations from the best fit line, presumably due to different extents of sigma and pi interactions in these complexes. Second, all the Co complexes are clustered in the upper-left quadrant of the plot while all the Mn complexes are clustered in the lower-right quadrant across all the families of ligands. This suggests that while CO does not thermodynamically trap the Co(III) complexes, the corresponding hydrides are all very stable relative to HCO2 and therefore unreactive towards CO2. Conversely, in the Mn(I) systems, the hydrides lie higher than HCO2 in relative energy, and therefore will react with CO2, but the corresponding carbonyl intermediates are very stable potentially leading to catalytic trap states. It is evident that of the three metals, the complexes of Fe are the best potential catalyst candidates as their image file: c9sc01766f-t12.tif and image file: c9sc01766f-t13.tif values lie closest to the origin.

The slopes of the fits in all three plots are roughly the same but the y-intercepts vary. Notice that changing the overall charge only moves you along the best fit line. In the cyclopentadienyl and pincer complexes (Fig. 1A and B), for example, the y-intercept is negative (ca. −5 kcal mol−1) suggesting that [M–CO]n+ will be dominant irrespective of the other ligands. The bis-bipyridine family of complexes, on the other hand, have a positive y-intercept of ca. 5 kcal mol−1 (Fig. 1C), resulting in a ca. 10 kcal mol−1 spread, which is greater than the ca. 5 kcal mol−1 error associated with such calculations. The bis-bipyridine iron complex with L = pyridine (image file: c9sc01766f-t14.tif, image file: c9sc01766f-t15.tif) is therefore an interesting candidate to test the predictions for ergoneutral metal hydride formation and the lack of catalyst trapping by CO, within an error of ±5 kcal mol−1 in the calculations (vide infra).

While the exact origin of the different intercept in the case of the bis-bipyridine family compared to the cyclopentadienyl and pincer families is unclear, we hypothesize that having a high degree of pyridyl ligation could play a role, presumably from a good balance between the σ-donating and π-accepting ability of pyridyl ligands. The former is responsible for conferring nucleophilicity on the metal center leading to a higher image file: c9sc01766f-t16.tif, and the latter for destabilizing [M–CO]n+ through competition for π-back donation from the metal center.

To further probe this hypothesis, we modeled the effect of increasing the π-accepting ability of the pyridyl ring on image file: c9sc01766f-t19.tif and image file: c9sc01766f-t20.tif for the [Fe(bpy)2(py)S]2+ system (Fig. 2). Firstly, conjugating the pyridine with one of the bipyridine ligands to make a terpyridine (tpy) ligand reduces the pyridine to Fe-CO dihedral angle (C-N-Fe-C) from 30° to zero, and aligns the π* orbital of the pyridine with the Fe-CO axis, which destabilizes the carbonyl ligand through competition for back-bonding. Consistent with this prediction, image file: c9sc01766f-t21.tif for the [Fe(bpy)(tpy)]2+ system increases by ca. 3 kcal mol−1 with almost no effect on image file: c9sc01766f-t22.tif. Secondly, substitution at the para-position of the pyridine ring with a nitro group has a similar effect on image file: c9sc01766f-t23.tif and image file: c9sc01766f-t24.tif. Consistent with the known trapping of the [Ru(bpy)(tpy)]2+ system by CO,38 the Ru analogs have much lower values of image file: c9sc01766f-t25.tif owing to the greater back-bonding ability of Ru compared to Fe. While the π-accepting nature of the ligand trans to CO in a metal complex is regularly invoked as a determinant for thermodynamic stability (the trans influence), our results highlight a substantial cis influence on the energy of the carbonyl (but not of the hydride), similar to effects we have seen previously.24

image file: c9sc01766f-f2.tif
Fig. 2 Effect of conjugation and substitution of the pyridine ring on image file: c9sc01766f-t26.tif and image file: c9sc01766f-t27.tif in the [M(bpy)2(py)]2+ system (M = Fe, Ru).

Electrocatalytic activity of [Fe(bpy2PYMe)(CH3CN)]2+

In order to experimentally test these predictions, we require a pentadentate ligand to prevent ligand scrambling. Long et al. recently reported the synthesis and characterization of the pentadentate ligand ‘bpy2PYMe’ (Scheme 4).40 This ligand closely matches the desired ligand environment in Fig. 1C for iron (L = pyridine), and therefore we synthesized and studied the electrocatalytic activity of the complex [Fe(bpy2PYMe)(CH3CN)]2+.
image file: c9sc01766f-s4.tif
Scheme 4 Structure of [Fe(bpy2PYMe)S](CF3SO3)2, where bpy2PYMe = 1-(2-pyridyl)-1-(6-2,2′- bipyridyl)ethane and S = CH3CN.

The calculated free energies of the catalytic intermediates for the different two-electron reduction processes (eqn (1)–(3)) mediated by the complex [Fe(bpy2PYMe)(CH3CN)]2+ are shown in Fig. 3, with phenol as the stoichiometric proton donor. In addition to possessing low corrugations with respect to image file: c9sc01766f-t28.tif and image file: c9sc01766f-t29.tif, viz. 1.7 kcal mol−1 and 7.4 kcal mol−1 respectively, the most endergonic on-path intermediates for this system are within 5 kcal mol−1. This complex is therefore likely to avoid high thermodynamic barriers and operate close to the standard potential for the reduction reactions. We note that in Fig. 3, the metal formate intermediate is predicted to be a product inhibitor. Previous work suggested that formate can be labilized with the modification of the solvent composition to provide H-bond stabilization.12,24

image file: c9sc01766f-f3.tif
Fig. 3 DFT-calculated energy level diagram (standard state) for CO2 reduction to HCO2 (in black) and CO (in grey), and proton reduction to H2 (in grey), mediated by [Fe(bpy2PYMe)(CH3CN)]2+, with phenol as the proton donor.

The cyclic voltammogram of [Fe(bpy2PYMe)(CH3CN)]2+ in acetonitrile is shown in Fig. 4. The complex undergoes two one-electron reductions at −1.55 V and −1.66 V vs. Cp2Fe+/0 (Cp2Fe = ferrocene). The experimental redox potentials match up perfectly with those reported by Long and coworkers.40,46 Moreover, the calculated redox potentials (Table 1) are within about 100 mV of these measured values. For further validation, CO gas was passed through a CD3CN solution of [Fe(bpy2PYMe)(CH3CN)]2+ in a J-Young NMR tube for ca. 30 minutes. No changes in the 1H NMR spectrum of [Fe(bpy2PYMe)(CH3CN)]2+ were observed, suggesting that the corresponding [Fe(bpy2PYMe)(CO)]2+ complex is in negligible concentration, which is in agreement with our calculations (Fig. 3).

image file: c9sc01766f-f4.tif
Fig. 4 Cyclic voltammograms of 2 mM [Fe(bpy2PYMe)(CH3CN)]2+ at 100 mV s. (A) [Fe(bpy2PYMe)(CH3CN)]2+ only (black solid), with the addition of 0.3 M phenol (red dashed) followed by satd. CO2 (ca. 0.28 M, blue dotted). (B) [Fe(bpy2PYMe)(CH3CN)]2+ only (black solid), with the addition of 0.3 M acetic acid (red dashed) or 0.5 M acetic acid (blue dotted). Insets show the linear fit of the ratio of the current with the substrate (icat) to the peak current in the absence of any substrate (ipc) at ca. −1.75 V as a function of the inverse square root of the scan rate (eqn (5)).
Table 1 Calculated and measured redox potentials of [Fe(bpy2PYMe)(CH3CN)]2+
E°/V vs. Cp2Fe+/0
Calculated Experimental
1st reduction −1.46 −1.55
2nd reduction −1.65 −1.66

Upon the addition of 0.3 M phenol (pKa = 29.14 in acetonitrile42) to the CV solution containing 2 mM [Fe(bpy2PYMe)(CH3CN)]2+, there is a slight enhancement in the current at the second reduction potential suggesting that H2 evolution is taking place through the formation of a metal hydride and subsequent protonation. Saturating this solution with CO2 (ca. 0.28 M)47 results in a higher (ca. 2.5-fold) current enhancement near the second reduction potential (ca. −1.66 V vs. Cp2Fe+/0), indicating that CO2 is a substrate for an electrocatalytic process mediated by [Fe(bpy2PYMe)(CH3CN)]2+. Assuming a 2 mM concentration of phenoxide and a 300 mM concentration of phenol48 we calculate that this electrocatalytic process is occurring at an overpotential of ca. 200 mV, from eqn (4) with pKa = 29.14 for phenol. Scan rate dependent ratios of the current in the presence of CO2 to the peak current in its absence (eqn (5), Fig. 4A (inset) and S2) yield a pseudo-first order rate constant of ca. 2.4 s−1 for the two-electron electrocatalytic process involving the reduction of CO2. Unexpectedly, there is also a positive shift of the onset potential of the first peak in the CV with no corresponding current increase at this potential. Upon the addition of CO2 in the absence of phenol (Fig. S3) a similar positive shift in the onset potential is observed in the CV. Both these results suggest that there is energetically favorable CO2 binding accompanying the first electrochemical reduction of [Fe(bpy2PYMe)(CH3CN)]2+, giving rise to the positive shift in the onset in the presence of CO2. On this basis, we assign the oxidative feature around ca. −1.35 V vs. Cp2Fe+/0 to the [Fe(bpy2PYMe)(CO2)]2+/+ couple. Fig. S4 shows that the experimental data can be approximately simulated by CO2 binding and a small amount of reductive deoxygenation at extremely negative potentials, presumably due to carbonate formation.33 This unexpected result of favorable CO2 binding to the singly reduced intermediate [Fe(bpy2PYMe)]+, contrary to the DFT-based predictions in Fig. 3, is discussed below.

The effect of proton strength on the electrocatalytic behavior of [Fe(bpy2PYMe)(CH3CN)]2+ was explored with acetic acid (pKa = 23.51 in acetonitrile)42 as the proton donor instead of phenol. In the cyclic voltammograms (Fig. 4B), the electrocatalytic currents near the second reduction potential of [Fe(bpy2PYMe)(CH3CN)]2+ are over three times greater than with phenol. However, the overpotential, as compared to the phenol case, is about 200 mV higher, as image file: c9sc01766f-t30.tif in the standard state (Table S1). With 0.5 M acetic acid, a pseudo-first order rate constant greater than 150 s−1 (the actual rate is limited by the cell resistance) was estimated from scan-rate dependent CVs (Fig. S2). Upon the addition of CO2 (ca. 0.28 M) there is a greater current enhancement and positive shift of the onset potential (Fig. S5), qualitatively similar to what was observed with phenol. Bulk electrolysis of a 2 mM solution of [Fe(bpy2PYMe)(CH3CN)]2+ was first performed with 0.3 M phenol and 0.1 M KPF6 as the supporting electrolyte in CO2-saturated CD3CN (3 mL solution) at a potential of −1.65 V vs. Cp2Fe+/0 near the onset of the second reduction wave (Fig. 4A, see the Experimental details). Decomposition of the Fe species in solution preceded any appreciable turnover, and therefore no products were detected. The electrolysis experiment was then repeated with 0.3 M acetic acid as the proton donor instead of phenol. A 1 mL aliquot of the headspace was injected into a gas chromatograph to detect and quantify gaseous products (Fig. S6). Based on the calibration of the GC peak areas with 1% gas standards, the amounts of gaseous products were estimated to be 8.1 × 10−6 moles of CO and 2.6 × 10−5 moles of H2. They correspond to 9% and 30% of the total charge passed in the electrolysis experiment respectively for two-electron stoichiometry, suggesting that other reduced products are present. In the 1H NMR spectrum of the electrolyte solution, the peak corresponding to HCO2 at ca. 8.2 ppm49 was not observed. Decomposition of the Fe species in solution prevented comprehensive analyses of the liquid phase products. With acetic acid as the proton donor, H2 and HCO2 are predicted to be roughly ergoneutral in the standard state (Fig. S1). However, under the electrolysis conditions, there is a ca. 200 mV extra driving force towards H2 from the higher proton activity due to the absence of a stoichiometric conjugate base, as well as from homoconjugation effects.42,48 These effects along with the second reduction potential of [Fe(bpy2PYMe)(CH3CN)]2+, −1.66 V vs. Cp2Fe+/0 (and the electrolysis potential), being ca. 400 mV more negative than the potential for H2 production under the electrolysis conditions (vide supra), together explain the observed H2 in the electrolysis.

We hypothesize that the lack of the thermodynamically favored product, HCO2, under the electrolysis conditions, is due to the rate of CO2 insertion into the metal hydride being significantly slower than the rate of protonation of the metal hydride. CO2 insertion into hydrides is typically a slow process (ca. 1.82 × 10−2 M−1 s−1 for [Ru(tpy)(bpy)H]+ at 298 K in acetonitrile50). The incorporation of proton-directing groups into the ligand backbone to accelerate the rate of CO2 insertion into the metal hydride, as we have explored in previous work,24 would favor the HCO2 pathway.

The production of CO as a minor product indicates direct binding of CO2 during the electrochemical reduction of [Fe(bpy2PYMe)(CH3CN)]2+. The direct CO2 binding pathway, evidenced by the positive shift of the first reduction wave in the CV in Fig. S3, competes with hydride formation (Scheme 2). However, DFT calculations at the current level of theory suggest that the binding of CO2 after the first reduction is uphill by ca. 7 kcal mol−1 (Fig. 3). The same calculation using the hybrid functional B3LYP instead of BP86 suggests that the CO2 binding step after the first reduction is downhill by about 4 kcal mol−1. It is, therefore, possible that while BP86 captures redox potentials and the free energy of key intermediates such as the metal hydrides and metal carbonyls with reasonable accuracy, it underestimates the free energy of CO2 binding to the reduced Fe-complexes. Therefore, a catalyst search algorithm that employs the free energies of the metal hydride and carbonyl intermediates alone (Fig. 1) will effectively weed out bad candidates based on these two descriptors. However, within this narrowed space of candidates, favorable CO2 binding to the metal center could lead to CO production, as evidenced in the case of [Fe(bpy2PYMe)(CH3CN)]2+. This motivates further refinement and experimental validation of DFT methods to adequately model the thermodynamics of CO2 binding to reduced metal complexes, in order to map the free energies of the whole range of intermediates in a catalytic pathway.


In this work we have shown how the calculation of two thermodynamic parameters, viz. the relative free energies of the metal hydride and the CO-bound intermediates by using DFT streamlines the search for an appropriate transition metal and ligand environment for catalyzing the multi-electron electrochemical reduction of CO2 or protons at low driving forces. We applied this catalyst screening approach on a library of common ligands around three earth-abundant d6 metal ions and tested the in silico predictions via the synthesis and electrochemical studies of an Fe-based complex, predicted to avoid trapping by CO, form a metal hydride of appropriate free energy, and possess optimal redox potentials. These predictions were validated, and we found the iron-based electrocatalyst to be active towards CO2 reduction as well as H2 production at the predicted potentials. The model, however, underestimated the free energy of CO2 binding to the reduced metal complex. We highlight the need for refinements to the DFT-methods to adequately capture the free energy of this latter step that leads to CO production, in order to further narrow the candidate space. This work is, to the best of our knowledge, the first application of DFT as a screen for molecular electrocatalysts across diverse molecular environments and paves the way forward for computations to take the lead in identifying promising catalysts for electrocatalytic transformations.

Experimental details

The ligand bpy2PYMe and its Fe(II) complex were prepared as per reported protocols.401H NMR of bpy2PYMe (CDCl3, 400 MHz, ppm): δ 8.62 (3H, d), 8.29 (2H, d), 8.22 (2H, d), 7.73 (2H, t), 7.68 (2H, d), 7.57 (1H, t), 7.24–7.20 (4H, m), 7.13 (2H, m), 2.53 (3H, s). bpy2PYMe was then stirred with Fe(CF3SO3)2 for 24 hours at room temperature in acetonitrile, followed by filtration through a Celite plug. The solvent was removed, and the crude residue was recrystallized with CH3CN–diethyl ether. The resulting dark-red/black crystals were filtered and dried in a vacuum oven at 50 °C for three days. 1H NMR (CD3CN, 500 MHz, ppm): δ 9.24 (1H), 8.53 (3H), 8.32 (2H, d), 8.13 (2H, d), 8.03 (2H), 7.97 (2H, d), 7.87 (2H, t), 7.38 (1H), 7.30 (1H, t), 7.21 (1H, d), 6.89 (1H, d), 2.80 (3H, s).

All electrochemical experiments were performed inside a N2 glovebox fitted with a CO2 feedthrough using a SP-200 potentiostat from Bio-Logic Co. For the cyclic voltammetry experiments, a glassy carbon electrode (3 mm diameter) from BASi Inc. was used in combination with a Ag/AgNO3 (10 mM)/TBAPF6 (100 mM) reference electrode and a Pt wire counter electrode in a sealed cell. The pseudo-first-order rate constant kcat for the catalytic waves in Fig. 4 was determined using eqn (5)12 (ν is the scan rate).

image file: c9sc01766f-t31.tif(5)

Bulk-electrolysis was performed in an ‘H-cell’ with two compartments, one for the working and reference electrodes and the other for the counter electrode, separated by a porous glass frit. The working electrode was a glassy carbon rod (BASi Inc.), and the reference and counter electrodes were the same as used in the cyclic voltammetry experiments. Both compartments were sealed with rubber septa and stirred vigorously during electrolysis. The entire process was carried out in a N2-filled glove box. Headspace product analyses were performed using a Shimadzu gas chromatograph fitted with a 100 μL sample injection loop, an FID, a TCD and a CarbonPLOT column with N2 carrier gas. Percent H2 and CO were determined with one-point calibration with a 1% gas standard purchased from Sigma-Aldrich Co.

Computational details

All calculations based on Kohn–Sham density functional theory were performed using the Gaussian 0951 (Rev D. 01) software package. The BP86 (ref. 52 and 53) functional was used for all the calculations in combination with a double-zeta 6-31+G*54–56 basis set on all the p-block elements, and the LANL2DZ57 effective core potential on the transition metals. The relevant intermediates were first optimized in the gas phase followed by a harmonic analysis on the stationary point to obtain enthalpic, entropic and zero-point energy corrections to the electronic energy in the standard state, as implemented in Gaussian 09. The solvation energy was then determined with a single point calculation with a polarizable continuum model (SMD) for acetonitrile.58 Standard state corrections were made to the free energies to account for the change in going from 1 mol per 24.46 L (gas phase) to 1 M (solution phase). This level of theory provided excellent agreement with experimental values of the reduction potentials versus the ferrocene/ferrocenium couple in acetonitrile (eqn (6)–(8)), as well as hydricities for representative complexes in each family of complexes listed in Scheme 3.37,59,60
Oxn+ + Cp2Fe ⇌ Red(n−1)+ + Cp2Fe+(6)
image file: c9sc01766f-t32.tif(7)
image file: c9sc01766f-t33.tif(8)

For all the complexes reported in this work, the lowest spin states (singlets for the even-electron intermediates and doublets for the odd-electron intermediates) were found to be the thermodynamically favorable states. This is expected for the local density functional BP86 given the choice of strong-field ligands.20 The B3LYP52,61–63 functional was employed in one instance for comparison of the CO2 binding energies to the singly reduced Fe complex. Simulations of the electrochemical data were performed with Digielch64 using a planar semi-infinite 1D diffusion model (see the ESI).

Conflicts of interest

The authors declare no competing financial interests.


This work was partially supported by the Global Climate and Energy Project (GCEP) at Stanford University, through a grant titled “Electrohydrogenation: Enabling Science for Renewable Fuels”. The authors would like to thank Professors Todd J. Martinez, Robert M. Waymouth and T. Daniel P. Stack for useful discussions. S. R. would also like to thank the Center for Molecular Analysis and Design at Stanford University for a graduate fellowship.


  1. M. Mikkelsen, M. Jørgensen and F. C. Krebs, Energy Environ. Sci., 2010, 3, 43–81 RSC.
  2. D. L. Dubois, Inorg. Chem., 2014, 53, 3935–3960 CrossRef CAS PubMed.
  3. A. M. Appel, J. E. Bercaw, A. B. Bocarsly, H. Dobbek, D. L. Dubois, M. Dupuis, J. G. Ferry, E. Fujita, R. Hille, P. J. A. Kenis, C. A. Kerfeld, R. H. Morris, C. H. F. Peden, A. R. Portis, S. W. Ragsdale, T. B. Rauchfuss, J. N. H. Reek, L. C. Seefeldt, R. K. Thauer and G. L. Waldrop, Chem. Rev., 2013, 113, 6621–6658 CrossRef CAS PubMed.
  4. C. W. Li and M. W. Kanan, J. Am. Chem. Soc., 2012, 134, 7231–7234 CrossRef CAS PubMed.
  5. Y. Chen, C. W. Li and M. W. Kanan, J. Am. Chem. Soc., 2012, 134, 19969–19972 CrossRef CAS PubMed.
  6. Y. Hori, in Modern Aspects of Electrochemistry, Springer New York, New York, NY, 2008, pp. 89–189 Search PubMed.
  7. K. P. Kuhl, E. R. Cave, D. N. Abram and T. F. Jaramillo, Energy Environ. Sci., 2012, 5, 7050–7059 RSC.
  8. K. A. Grice and C. P. Kubiak, Adv. Inorg. Chem., 2014, 66, 163–188 CrossRef CAS.
  9. C. Costentin, M. Robert and J.-M. Savéant, Chem. Soc. Rev., 2013, 42, 2423–2436 RSC.
  10. M. L. Helm, M. P. Stewart, R. M. Bullock, M. R. DuBois and D. L. DuBois, Science, 2011, 333, 863–866 CrossRef CAS PubMed.
  11. J. R. McKone, S. C. Marinescu, B. S. Brunschwig, J. R. Winkler and H. B. Gray, Chem. Sci., 2014, 5, 865–878 RSC.
  12. P. Kang, C. Cheng, Z. Chen, C. K. Schauer, T. J. Meyer and M. Brookhart, J. Am. Chem. Soc., 2012, 134, 5500–5503 CrossRef CAS PubMed.
  13. J. R. Pugh, M. R. M. Bruce, B. P. Sullivan and T. J. Meyer, Inorg. Chem., 1991, 30, 86–91 CrossRef CAS.
  14. C. Caix, S. Chardon-Noblat and A. Deronzier, J. Electroanal. Chem., 1997, 434, 163–170 CrossRef CAS.
  15. A. Taheri and L. A. Berben, Chem. Commun., 2016, 52, 1768–1777 RSC.
  16. S. Roy, B. Sharma, J. Pécaut, P. Simon, M. Fontecave, P. D. Tran, E. Derat and V. Artero, J. Am. Chem. Soc., 2017, 139, 3685–3696 CrossRef CAS PubMed.
  17. B. M. Ceballos and J. Y. Yang, Proc. Natl. Acad. Sci. U. S. A., 2018, 115, 12686–12691 CrossRef CAS PubMed.
  18. S. Raugei, D. L. Dubois, R. Rousseau, S. Chen, M. H. Ho, R. M. Bullock and M. Dupuis, Acc. Chem. Res., 2015, 48, 248–255 CrossRef CAS PubMed.
  19. C. J. Cramer and D. G. Truhlar, Phys. Chem. Chem. Phys., 2009, 11, 10757–10816 RSC.
  20. F. Neese, Coord. Chem. Rev., 2009, 253, 526–563 CrossRef CAS.
  21. L. E. Roy, E. Jakubikova, M. Graham Guthrie and E. R. Batista, J. Phys. Chem. A, 2009, 113, 6745–6750 CrossRef CAS PubMed.
  22. B. H. Solis and S. Hammes-Schiffer, Inorg. Chem., 2011, 50, 11252–11262 CrossRef CAS PubMed.
  23. Y. Zhu, Y. Fan and K. Burgess, J. Am. Chem. Soc., 2010, 132, 6249–6253 CrossRef CAS PubMed.
  24. S. Ramakrishnan, K. M. Waldie, I. Warnke, A. G. De Crisci, V. S. Batista, R. M. Waymouth and C. E. D. Chidsey, Inorg. Chem., 2016, 55, 1623–1632 CrossRef CAS PubMed.
  25. M. L. Pegis, B. A. McKeown, N. Kumar, K. Lang, D. J. Wasylenko, X. P. Zhang, S. Raugei and J. M. Mayer, ACS Cent. Sci., 2016, 2, 850–856 CrossRef CAS PubMed.
  26. A. J. Göttle and M. T. M. Koper, J. Am. Chem. Soc., 2018, 140, 4826–4834 CrossRef PubMed.
  27. B. H. Solis and S. Hammes-Schiffer, Inorg. Chem., 2014, 53, 6427–6443 CrossRef CAS PubMed.
  28. J. T. Muckerman, P. Achord, C. Creutz, D. E. Polyansky and E. Fujita, Proc. Natl. Acad. Sci. U. S. A., 2012, 109, 15657–15662 CrossRef CAS PubMed.
  29. S. I. Johnson, R. J. Nielsen and W. A. Goddard, ACS Catal., 2016, 6, 6362–6371 CrossRef CAS.
  30. C. Tsay, B. N. Livesay, S. Ruelas and J. Y. Yang, J. Am. Chem. Soc., 2015, 137, 14114–14121 CrossRef CAS PubMed.
  31. M. S. Jeletic, M. T. Mock, A. M. Appel and J. C. Linehan, J. Am. Chem. Soc., 2013, 135, 11533–11536 CrossRef CAS PubMed.
  32. W. Baratta, G. Chelucci, S. Gladiali, K. Siega, M. Toniutti, M. Zanette, E. Zangrando and P. Rigo, Angew. Chem., Int. Ed., 2005, 44, 6214–6219 CrossRef CAS PubMed.
  33. S. Ramakrishnan and C. E. D. Chidsey, Inorg. Chem., 2017, 56, 8326–8333 CrossRef CAS PubMed.
  34. Y. Huang, C. C. Neto, K. A. Pevear, M. M. Banaszak Holl, D. A. Sweigart and Y. K. Chung, Inorg. Chim. Acta, 1994, 226, 53–60 CrossRef CAS.
  35. C. Creutz and M. H. Chou, J. Am. Chem. Soc., 2009, 131, 2794–2795 CrossRef CAS PubMed.
  36. C. Creutz, M. H. Chou, H. Hou and J. T. Muckerman, Inorg. Chem., 2010, 49, 9809–9822 CrossRef CAS PubMed.
  37. Y. Matsubara, E. Fujita, M. D. Doherty, J. T. Muckerman and C. Creutz, J. Am. Chem. Soc., 2012, 134, 15743–15757 CrossRef CAS PubMed.
  38. Z. Chen, C. Chen, D. R. Weinberg, P. Kang, J. J. Concepcion, D. P. Harrison, M. S. Brookhart and T. J. Meyer, Chem. Commun., 2011, 47, 12607–12609 RSC.
  39. K. M. Waldie, S. Ramakrishnan, S. K. Kim, J. K. Maclaren, C. E. D. Chidsey and R. M. Waymouth, J. Am. Chem. Soc., 2017, 139, 4540–4550 CrossRef CAS PubMed.
  40. M. Nippe, R. S. Khnayzer, J. A. Panetier, D. Z. Zee, B. S. Olaiya, M. Head-Gordon, C. J. Chang, F. N. Castellano and J. R. Long, Chem. Sci., 2013, 4, 3934–3945 RSC.
  41. E. S. Wiedner, M. B. Chambers, C. L. Pitman, R. M. Bullock, A. J. M. Miller and A. M. Appel, Chem. Rev., 2016, 116, 8655–8692 CrossRef CAS PubMed.
  42. B. D. McCarthy, D. J. Martin, E. S. Rountree, A. C. Ullman and J. L. Dempsey, Inorg. Chem., 2014, 53, 8350–8361 CrossRef CAS PubMed.
  43. M. L. Pegis, J. A. S. Roberts, D. J. Wasylenko, E. A. Mader, A. M. Appel and J. M. Mayer, Inorg. Chem., 2015, 54, 11883–11888 CrossRef CAS PubMed.
  44. These values neglect the homoconjugation of phenol with phenoxide which is discussed in greater detail later.
  45. The addition of 43 kcal mol−1 (hydricity of formate in acetonitrile) to the values in Fig. 1 will yield hydricity values of the metal hydrides in Scheme 3.
  46. Ref. 40 does not report electrocatalytic H2 evolution with the Fe(II) complex.
  47. E. Fujita, C. Creutz, N. Sutin and D. J. Szalda, J. Am. Chem. Soc., 1991, 113, 343–353 CrossRef CAS.
  48. Assuming a 1 mM concentration of the conjugate base, the proton activity is ca. 300 Ka, which corresponds to a 2.5 unit lower pH. This corresponds to a 150 mV extra driving force relative to the standard state. Homoconjugation of acetic acid at these concentrations (ref. 42) leads to a further decrease of 1 pH unit, which corresponds to another 60 mV extra driving force.
  49. S. T. Ahn, E. A. Bielinski, E. M. Lane, Y. Chen, W. H. Bernskoetter, N. Hazari and G. T. R. Palmore, Chem. Commun., 2015, 51, 5947–5950 RSC.
  50. H. Konno, A. Kobayashi, K. Sakamoto, F. Fagalde, N. E. Katz, H. Saitoh and O. Ishitani, Inorg. Chim. Acta, 2000, 299, 155–163 CrossRef CAS.
  51. M. J. Frisch, G. W. Trucks, H. B. Schlegel, G. E. Scuseria, M. A. Robb, J. R. Cheeseman, G. Scalmani, V. Barone, B. Mennucci, G. A. Petersson, H. Nakatsuji, M. Caricato, X. Li, H. P. Hratchian, A. F. Izmaylov, J. Bloino, G. Zheng, J. L. Sonnenberg, M. Hada, M. Ehara, K. Toyota, R. Fukuda, J. Hasegawa, M. Ishida, T. Nakajima, Y. Honda, O. Kitao, H. Nakai, T. Vreven, J. A. Montgomery, J. E. Peralta, F. Ogliaro, M. Bearpark, J. J. Heyd, E. Brothers, K. N. Kudin, V. N. Staroverov, R. Kobayashi, J. Normand, K. Raghavachari, A. Rendell, J. C. Burant, S. S. Iyengar, J. Tomasi, M. Cossi, N. Rega, J. M. Millam, M. Klene, J. E. Knox, J. B. Cross, V. Bakken, C. Adamo, J. Jaramillo, R. Gomperts, R. E. Stratmann, O. Yazyev, A. J. Austin, R. Cammi, C. Pomelli, J. W. Ochterski, R. L. Martin, K. Morokuma, V. G. Zakrzewski, G. A. Voth, P. Salvador, J. J. Dannenberg, S. Dapprich, A. D. Daniels, O. Farkas, J. B. Foresman, J. V Ortiz, J. Cioslowski and D. J. Fox, Gaussian 09, Revision D.01, Gaussian, Inc., Wallingford CT, 2009 Search PubMed.
  52. A. D. Becke, Phys. Rev. A, 1988, 38, 3098–3100 CrossRef CAS.
  53. J. P. Perdew, Phys. Rev. B, 1986, 33, 8822–8824 CrossRef PubMed.
  54. R. Ditchfield, W. J. Hehre and J. A. Pople, J. Chem. Phys., 1971, 54, 724–728 CrossRef CAS.
  55. P. C. Hariharan and J. A. Pople, Theor. Chim. Acta, 1973, 28, 213–222 CrossRef CAS.
  56. V. A. Rassolov, M. A. Ratner, J. A. Pople, P. C. Redfern and L. A. Curtiss, J. Comput. Chem., 2001, 22, 976–984 CrossRef CAS.
  57. P. J. Hay and W. R. Wadt, J. Chem. Phys., 1985, 82, 270–283 CrossRef CAS.
  58. A. V. Marenich, C. J. Cramer and D. G. Truhlar, J. Phys. Chem. B, 2009, 113, 6378–6396 CrossRef CAS PubMed.
  59. D. P. Estes, A. K. Vannucci, A. R. Hall, D. L. Lichtenberger and J. R. Norton, Organometallics, 2011, 30, 3444–3447 CrossRef CAS.
  60. K. M. Waldie, A. L. Ostericher, M. H. Reineke, A. F. Sasayama and C. P. Kubiak, ACS Catal., 2018, 8, 1313–1324 CrossRef CAS.
  61. A. D. Becke, J. Chem. Phys., 1993, 98, 5648–5652 CrossRef CAS.
  62. P. J. Stephens, F. J. Devlin, C. F. Chabalowski and M. J. Frisch, J. Phys. Chem., 1994, 98, 11623–11627 CrossRef CAS.
  63. C. Lee, W. Yang and R. G. Parr, Phys. Rev. B, 1988, 37, 785–789 CrossRef CAS PubMed.
  64. M. Rudolph, J. Electroanal. Chem., 2003, 543, 23–39 CrossRef CAS.


Electronic supplementary information (ESI) available: Supplemental thermochemical calculations, electrochemical data and simulations and Cartesian coordinates of DFT-optimized intermediates. See DOI: 10.1039/c9sc01766f

This journal is © The Royal Society of Chemistry 2019