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

Abstract

The free energy landscape of catalytic intermediates in the two-electron reduction of proton donors and/or CO₂ to H₂, CO and HCO₂⁻ 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 HCO₂⁻, 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 d⁶ 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, pK_a'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 CO₂ reduction at the predicted redox potentials with no trapping by CO. However, contrary to the in silico predictions, we found electrochemical evidence of CO₂ 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 CO₂ intermediates, further refinement and validation of the DFT methods are needed.

---

Introduction

The multi-electron reduction of abundant feedstocks such as protons and CO₂ to fuels and value-added products^1–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 CO₂ 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 CO₂ to methanol or C₂ 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 H₂ and of CO₂ to HCO₂⁻ 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 H₂ and CO production^8–11 than for HCO₂⁻ production,^12–17 and even the few examples of the latter involve precious metals and require high overpotentials.

HA + CO₂ + 2e⁻ → HCO₂⁻ + A⁻ (1)

2HA + 2e⁻ → H₂ + 2A⁻ (2)

HA + 2CO₂ + 2e⁻ → CO + HCO₃⁻ + 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.

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, pK_a'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 CO₂ and protons to produce HCO₂⁻ and H₂ 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 CO₂ insertion into a ruthenium hydride complex.²⁴ The hydride complex we were studying despite being an active transfer hydrogenation catalyst for ketones³² was not found to be a catalyst for CO₂ reduction due to product inhibition. In another case, a cyclopentadienyl Ru complex³³ was predicted by using DFT to have near ergoneutral CO₂ insertion into the corresponding Ru-H intermediate. However, CO₂ 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,³⁴ 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 CO₂ 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 CO₂ insertion to produce HCO₂⁻, under electrocatalytic conditions, CO is the major product.³⁸ 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 d⁶ 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 [Fe^II(bpy)₂(py)(CH₃CN)]²⁺ complex would not be experimentally stable due to ligand scrambling at ambient temperatures, a pentadentate iron complex, [Fe(bpy2PYMe)S]²⁺ (bpy2PYMe = 1-(2-pyridyl)-1-bis-(6-2,2′-bipyridyl)ethane, S = CH₃CN), reported previously by Long et al.⁴⁰ was synthesized. [Fe(bpy2PYMe)S]²⁺ was indeed found to electrocatalytically reduce protons as well as CO₂ 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 CO₂ 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 HCO₂⁻. 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–CO₂ 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 CO₂ to HCO₂⁻ and of protons to H₂ in the standard states of reactants and products, the relative driving force in acetonitrile as a function of the pK_a of the stoichiometric proton donor, HA, is shown in eqn (4). They are based on the thermochemical analyses shown in Table S1.†

Given the difference in proton stoichiometry between eqn (1) and (2), plotting versus the pK_a of the proton donor HA in acetonitrile yields a straight line that crosses zero at a pK_a of ca. 24 (Fig. S1†). With weaker acids in acetonitrile (pK_a > 24) there will therefore be a thermodynamic bias for CO₂ reduction to HCO₂⁻versus proton reduction.⁴¹ For example, phenol (pK_a = 29.14)⁴² should disfavor H₂ while acetic acid (pK_a = 23)⁴² should favor H₂.

Under weakly acidic conditions, the two-electron, one-proton reduction of CO₂ 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 CO₂ to CO and water (eqn S1†). We estimate this latter process (eqn S1†) to be roughly uphill by ca. 4.8 kcal mol⁻¹ relative to H₂ production under water levels of ∼1 mM in acetonitrile (eqn S2†), in agreement with reported values.⁴³

Catalytic intermediates in the two-electron reductions

The reactants, products and catalytic intermediates involved in the two-electron reduction of CO₂ and protons as expected from chemical precedent are shown in Scheme 2, mediated by a generic transition metal complex [M–S]ⁿ⁺, where M is a transition metal in a six-coordinate ligand environment in which S is a labile solvent ligand (CH₃CN 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 CO₂ and HCO₂⁻ ergoneutral in the standard state, H₂ at a relative free energy of 6.7 kcal mol⁻¹ (eqn (4)), and CO with an upper bound of ca. 11.5 kcal mol⁻¹ (see the ESI†).⁴⁴

Scheme 2 Two-electron reduction pathways for a generic transition metal electrocatalyst [M–S]ⁿ⁺, HA = proton donor. The free energies of catalytic intermediates are chosen arbitrarily.

Two-electron reduction of the metal complex [M–S]ⁿ⁺ with the dissociation of S is followed by protonation at the metal center to form the metal hydride intermediate. CO₂ insertion into the metal hydride bond leads to HCO₂⁻ while direct protonation of the metal hydride yields H₂. Additionally, the singly or doubly reduced metal complex, [M]⁽ⁿ⁻¹⁾⁺ and [M]⁽ⁿ⁻²⁾⁺ can directly bind CO₂ at the metal center³³ subsequently leading to the potential release of CO via an [M–CO]ⁿ⁺ 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 and the metal carbonyl intermediates relative to the thermodynamically favorably product with weak proton donors, HCO₂⁻. Given the same proton stoichiometry in eqn (1) and (3), the comparison of the two thermodynamic quantities and relative to HCO₂⁻ 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 d⁶ 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.

Scheme 3 Families of d⁶ transition metal complexes screened (S = CH₃CN, bpy = 2,2′-bipyridine, dppe = bis(diphenylphosphino) ethane, and en = 1,2-ethylene diamine).

The two thermodynamic quantities and relative to HCO₂⁻ 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 HCO₂⁻ production, the optimal value of is zero and is greater than zero, and vice versa for CO production. values greater than 10 kcal mol⁻¹, for example, would inhibit the formation of [M–CO]ⁿ⁺ at room temperature.⁴⁵

Fig. 1 2D plot of the calculated values of and (in kcal mol⁻¹ units) for the d⁶ 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]ⁿ⁺ decreases and the free energy of [M–H]⁽ⁿ⁻¹⁾⁺ goes up. This is reasonable as an increase in electron density at the metal center is expected to increase the [M–H]⁽ⁿ⁻¹⁾⁺ free energy but lower the [M–CO]ⁿ⁺ 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 HCO₂⁻ and therefore unreactive towards CO₂. Conversely, in the Mn(I) systems, the hydrides lie higher than HCO₂⁻ in relative energy, and therefore will react with CO₂, 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 and 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⁻¹) suggesting that [M–CO]ⁿ⁺ 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⁻¹ (Fig. 1C), resulting in a ca. 10 kcal mol⁻¹ spread, which is greater than the ca. 5 kcal mol⁻¹ error associated with such calculations. The bis-bipyridine iron complex with L = pyridine (, ) 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⁻¹ 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 , and the latter for destabilizing [M–CO]ⁿ⁺ 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 and for the [Fe(bpy)₂(py)S]²⁺ 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, for the [Fe(bpy)(tpy)]²⁺ system increases by ca. 3 kcal mol⁻¹ with almost no effect on . Secondly, substitution at the para-position of the pyridine ring with a nitro group has a similar effect on and . Consistent with the known trapping of the [Ru(bpy)(tpy)]²⁺ system by CO,³⁸ the Ru analogs have much lower values of 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.²⁴

Fig. 2 Effect of conjugation and substitution of the pyridine ring on and in the [M(bpy)₂(py)]²⁺ system (M = Fe, Ru).

Electrocatalytic activity of [Fe(bpy2PYMe)(CH₃CN)]²⁺

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).⁴⁰ 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)(CH₃CN)]²⁺.

Scheme 4 Structure of [Fe(bpy2PYMe)S](CF₃SO₃)₂, where bpy2PYMe = 1-(2-pyridyl)-1-(6-2,2′- bipyridyl)ethane and S = CH₃CN.

The calculated free energies of the catalytic intermediates for the different two-electron reduction processes (eqn (1)–(3)) mediated by the complex [Fe(bpy2PYMe)(CH₃CN)]²⁺ are shown in Fig. 3, with phenol as the stoichiometric proton donor. In addition to possessing low corrugations with respect to and , viz. 1.7 kcal mol⁻¹ and 7.4 kcal mol⁻¹ respectively, the most endergonic on-path intermediates for this system are within 5 kcal mol⁻¹. 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

Fig. 3 DFT-calculated energy level diagram (standard state) for CO₂ reduction to HCO₂⁻ (in black) and CO (in grey), and proton reduction to H₂ (in grey), mediated by [Fe(bpy2PYMe)(CH₃CN)]²⁺, with phenol as the proton donor.

The cyclic voltammogram of [Fe(bpy2PYMe)(CH₃CN)]²⁺ in acetonitrile is shown in Fig. 4. The complex undergoes two one-electron reductions at −1.55 V and −1.66 V vs. Cp₂Fe^+/0 (Cp₂Fe = 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 CD₃CN solution of [Fe(bpy2PYMe)(CH₃CN)]²⁺ in a J-Young NMR tube for ca. 30 minutes. No changes in the ¹H NMR spectrum of [Fe(bpy2PYMe)(CH₃CN)]²⁺ were observed, suggesting that the corresponding [Fe(bpy2PYMe)(CO)]²⁺ complex is in negligible concentration, which is in agreement with our calculations (Fig. 3).

Fig. 4 Cyclic voltammograms of 2 mM [Fe(bpy2PYMe)(CH₃CN)]²⁺ at 100 mV s. (A) [Fe(bpy2PYMe)(CH₃CN)]²⁺ only (black solid), with the addition of 0.3 M phenol (red dashed) followed by satd. CO₂ (ca. 0.28 M, blue dotted). (B) [Fe(bpy2PYMe)(CH₃CN)]²⁺ 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 (i_cat) to the peak current in the absence of any substrate (i_pc) 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)(CH₃CN)]²⁺

	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 (pK_a = 29.14 in acetonitrile⁴²) to the CV solution containing 2 mM [Fe(bpy2PYMe)(CH₃CN)]²⁺, there is a slight enhancement in the current at the second reduction potential suggesting that H₂ evolution is taking place through the formation of a metal hydride and subsequent protonation. Saturating this solution with CO₂ (ca. 0.28 M)⁴⁷ results in a higher (ca. 2.5-fold) current enhancement near the second reduction potential (ca. −1.66 V vs. Cp₂Fe^+/0), indicating that CO₂ is a substrate for an electrocatalytic process mediated by [Fe(bpy2PYMe)(CH₃CN)]²⁺. Assuming a 2 mM concentration of phenoxide and a 300 mM concentration of phenol⁴⁸ we calculate that this electrocatalytic process is occurring at an overpotential of ca. 200 mV, from eqn (4) with pK_a = 29.14 for phenol. Scan rate dependent ratios of the current in the presence of CO₂ 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⁻¹ for the two-electron electrocatalytic process involving the reduction of CO₂. 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 CO₂ 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 CO₂ binding accompanying the first electrochemical reduction of [Fe(bpy2PYMe)(CH₃CN)]²⁺, giving rise to the positive shift in the onset in the presence of CO₂. On this basis, we assign the oxidative feature around ca. −1.35 V vs. Cp₂Fe^+/0 to the [Fe(bpy2PYMe)(CO₂)]^2+/+ couple. Fig. S4† shows that the experimental data can be approximately simulated by CO₂ binding and a small amount of reductive deoxygenation at extremely negative potentials, presumably due to carbonate formation.³³ This unexpected result of favorable CO₂ 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)(CH₃CN)]²⁺ was explored with acetic acid (pK_a = 23.51 in acetonitrile)⁴² as the proton donor instead of phenol. In the cyclic voltammograms (Fig. 4B), the electrocatalytic currents near the second reduction potential of [Fe(bpy2PYMe)(CH₃CN)]²⁺ are over three times greater than with phenol. However, the overpotential, as compared to the phenol case, is about 200 mV higher, as in the standard state (Table S1†). With 0.5 M acetic acid, a pseudo-first order rate constant greater than 150 s⁻¹ (the actual rate is limited by the cell resistance) was estimated from scan-rate dependent CVs (Fig. S2†). Upon the addition of CO₂ (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)(CH₃CN)]²⁺ was first performed with 0.3 M phenol and 0.1 M KPF₆ as the supporting electrolyte in CO₂-saturated CD₃CN (3 mL solution) at a potential of −1.65 V vs. Cp₂Fe^+/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⁻⁶ moles of CO and 2.6 × 10⁻⁵ moles of H₂. 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 ¹H NMR spectrum of the electrolyte solution, the peak corresponding to HCO₂⁻ at ca. 8.2 ppm⁴⁹ 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, H₂ and HCO₂⁻ 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 H₂ 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)(CH₃CN)]²⁺, −1.66 V vs. Cp₂Fe^+/0 (and the electrolysis potential), being ca. 400 mV more negative than the potential for H₂ production under the electrolysis conditions (vide supra), together explain the observed H₂ in the electrolysis.

We hypothesize that the lack of the thermodynamically favored product, HCO₂⁻, under the electrolysis conditions, is due to the rate of CO₂ insertion into the metal hydride being significantly slower than the rate of protonation of the metal hydride. CO₂ insertion into hydrides is typically a slow process (ca. 1.82 × 10⁻² M⁻¹ s⁻¹ for [Ru(tpy)(bpy)H]⁺ at 298 K in acetonitrile⁵⁰). The incorporation of proton-directing groups into the ligand backbone to accelerate the rate of CO₂ insertion into the metal hydride, as we have explored in previous work,²⁴ would favor the HCO₂⁻ pathway.

The production of CO as a minor product indicates direct binding of CO₂ during the electrochemical reduction of [Fe(bpy2PYMe)(CH₃CN)]²⁺. The direct CO₂ 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 CO₂ after the first reduction is uphill by ca. 7 kcal mol⁻¹ (Fig. 3). The same calculation using the hybrid functional B3LYP instead of BP86 suggests that the CO₂ binding step after the first reduction is downhill by about 4 kcal mol⁻¹. 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 CO₂ 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 CO₂ binding to the metal center could lead to CO production, as evidenced in the case of [Fe(bpy2PYMe)(CH₃CN)]²⁺. This motivates further refinement and experimental validation of DFT methods to adequately model the thermodynamics of CO₂ binding to reduced metal complexes, in order to map the free energies of the whole range of intermediates in a catalytic pathway.

Conclusions

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 CO₂ or protons at low driving forces. We applied this catalyst screening approach on a library of common ligands around three earth-abundant d⁶ 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 CO₂ reduction as well as H₂ production at the predicted potentials. The model, however, underestimated the free energy of CO₂ 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.⁴⁰¹H NMR of bpy2PYMe (CDCl₃, 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(CF₃SO₃)₂ 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 CH₃CN–diethyl ether. The resulting dark-red/black crystals were filtered and dried in a vacuum oven at 50 °C for three days. ¹H NMR (CD₃CN, 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 N₂ glovebox fitted with a CO₂ 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/AgNO₃ (10 mM)/TBAPF₆ (100 mM) reference electrode and a Pt wire counter electrode in a sealed cell. The pseudo-first-order rate constant k_cat for the catalytic waves in Fig. 4 was determined using eqn (5)¹² (ν is the scan rate).

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 N₂-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 N₂ carrier gas. Percent H₂ 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 09⁵¹ (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 LANL2DZ⁵⁷ 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.⁵⁸ 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

Oxⁿ⁺ + Cp₂Fe ⇌ Red⁽ⁿ⁻¹⁾⁺ + Cp₂Fe⁺ (6)

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.²⁰ The B3LYP^52,61–63 functional was employed in one instance for comparison of the CO₂ binding energies to the singly reduced Fe complex. Simulations of the electrochemical data were performed with Digielch⁶⁴ using a planar semi-infinite 1D diffusion model (see the ESI†).

Conflicts of interest

The authors declare no competing financial interests.

Acknowledgements

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.

References

1. M. Mikkelsen, M. Jørgensen and F. C. Krebs, Energy Environ. Sci., 2010, 3, 43–81.
2. D. L. Dubois, Inorg. Chem., 2014, 53, 3935–3960.
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.
4. C. W. Li and M. W. Kanan, J. Am. Chem. Soc., 2012, 134, 7231–7234.
5. Y. Chen, C. W. Li and M. W. Kanan, J. Am. Chem. Soc., 2012, 134, 19969–19972.
6. Y. Hori, in Modern Aspects of Electrochemistry, Springer New York, New York, NY, 2008, pp. 89–189.
7. K. P. Kuhl, E. R. Cave, D. N. Abram and T. F. Jaramillo, Energy Environ. Sci., 2012, 5, 7050–7059.
8. K. A. Grice and C. P. Kubiak, Adv. Inorg. Chem., 2014, 66, 163–188.
9. C. Costentin, M. Robert and J.-M. Savéant, Chem. Soc. Rev., 2013, 42, 2423–2436.
10. M. L. Helm, M. P. Stewart, R. M. Bullock, M. R. DuBois and D. L. DuBois, Science, 2011, 333, 863–866.
11. J. R. McKone, S. C. Marinescu, B. S. Brunschwig, J. R. Winkler and H. B. Gray, Chem. Sci., 2014, 5, 865–878.
12. P. Kang, C. Cheng, Z. Chen, C. K. Schauer, T. J. Meyer and M. Brookhart, J. Am. Chem. Soc., 2012, 134, 5500–5503.
13. J. R. Pugh, M. R. M. Bruce, B. P. Sullivan and T. J. Meyer, Inorg. Chem., 1991, 30, 86–91.
14. C. Caix, S. Chardon-Noblat and A. Deronzier, J. Electroanal. Chem., 1997, 434, 163–170.
15. A. Taheri and L. A. Berben, Chem. Commun., 2016, 52, 1768–1777.
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.
17. B. M. Ceballos and J. Y. Yang, Proc. Natl. Acad. Sci. U. S. A., 2018, 115, 12686–12691.
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.
19. C. J. Cramer and D. G. Truhlar, Phys. Chem. Chem. Phys., 2009, 11, 10757–10816.
20. F. Neese, Coord. Chem. Rev., 2009, 253, 526–563.
21. L. E. Roy, E. Jakubikova, M. Graham Guthrie and E. R. Batista, J. Phys. Chem. A, 2009, 113, 6745–6750.
22. B. H. Solis and S. Hammes-Schiffer, Inorg. Chem., 2011, 50, 11252–11262.
23. Y. Zhu, Y. Fan and K. Burgess, J. Am. Chem. Soc., 2010, 132, 6249–6253.
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.
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.
26. A. J. Göttle and M. T. M. Koper, J. Am. Chem. Soc., 2018, 140, 4826–4834.
27. B. H. Solis and S. Hammes-Schiffer, Inorg. Chem., 2014, 53, 6427–6443.
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.
29. S. I. Johnson, R. J. Nielsen and W. A. Goddard, ACS Catal., 2016, 6, 6362–6371.
30. C. Tsay, B. N. Livesay, S. Ruelas and J. Y. Yang, J. Am. Chem. Soc., 2015, 137, 14114–14121.
31. M. S. Jeletic, M. T. Mock, A. M. Appel and J. C. Linehan, J. Am. Chem. Soc., 2013, 135, 11533–11536.
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.
33. S. Ramakrishnan and C. E. D. Chidsey, Inorg. Chem., 2017, 56, 8326–8333.
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.
35. C. Creutz and M. H. Chou, J. Am. Chem. Soc., 2009, 131, 2794–2795.
36. C. Creutz, M. H. Chou, H. Hou and J. T. Muckerman, Inorg. Chem., 2010, 49, 9809–9822.
37. Y. Matsubara, E. Fujita, M. D. Doherty, J. T. Muckerman and C. Creutz, J. Am. Chem. Soc., 2012, 134, 15743–15757.
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.
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.
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.
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.
42. B. D. McCarthy, D. J. Martin, E. S. Rountree, A. C. Ullman and J. L. Dempsey, Inorg. Chem., 2014, 53, 8350–8361.
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.
44. These values neglect the homoconjugation of phenol with phenoxide which is discussed in greater detail later.
45. The addition of 43 kcal mol⁻¹ (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 H₂ 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.
48. Assuming a 1 mM concentration of the conjugate base, the proton activity is ca. 300 K_a, 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.
50. H. Konno, A. Kobayashi, K. Sakamoto, F. Fagalde, N. E. Katz, H. Saitoh and O. Ishitani, Inorg. Chim. Acta, 2000, 299, 155–163.
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.
52. A. D. Becke, Phys. Rev. A, 1988, 38, 3098–3100.
53. J. P. Perdew, Phys. Rev. B, 1986, 33, 8822–8824.
54. R. Ditchfield, W. J. Hehre and J. A. Pople, J. Chem. Phys., 1971, 54, 724–728.
55. P. C. Hariharan and J. A. Pople, Theor. Chim. Acta, 1973, 28, 213–222.
56. V. A. Rassolov, M. A. Ratner, J. A. Pople, P. C. Redfern and L. A. Curtiss, J. Comput. Chem., 2001, 22, 976–984.
57. P. J. Hay and W. R. Wadt, J. Chem. Phys., 1985, 82, 270–283.
58. A. V. Marenich, C. J. Cramer and D. G. Truhlar, J. Phys. Chem. B, 2009, 113, 6378–6396.
59. D. P. Estes, A. K. Vannucci, A. R. Hall, D. L. Lichtenberger and J. R. Norton, Organometallics, 2011, 30, 3444–3447.
60. K. M. Waldie, A. L. Ostericher, M. H. Reineke, A. F. Sasayama and C. P. Kubiak, ACS Catal., 2018, 8, 1313–1324.
61. A. D. Becke, J. Chem. Phys., 1993, 98, 5648–5652.
62. P. J. Stephens, F. J. Devlin, C. F. Chabalowski and M. J. Frisch, J. Phys. Chem., 1994, 98, 11623–11627.
63. C. Lee, W. Yang and R. G. Parr, Phys. Rev. B, 1988, 37, 785–789.
64. M. Rudolph, J. Electroanal. Chem., 2003, 543, 23–39.

---

Footnote

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

---