Linear Scaling Relationships to Predict pKa's and Reduction Potentials for Bioinspired Hydrogenase Catalysis.

Biomimetic catalysts inspired by the active site of the [FeFe] hydrogenase enzyme can convert protons into molecular hydrogen. Minimizing the overpotential of the electrocatalytic process remains a major challenge for practical application of the catalyst. The catalytic cycle of the hydrogen production follows an ECEC mechanism (E represents an electron transfer step, and C refers to a chemical step), in which the electron and proton transfer steps can be either sequential or coupled (PCET). In this study, we have calculated the pKa's and the reduction potentials for a series of commonly used ligands (80 different complexes) using density functional theory. We establish that the required acid strength for protonation at the Fe-Fe site correlates with the standard reduction potential of the di-iron complexes with a linear energy relationship. These linear relationships allow for fast screening of ligands and tuning of the properties of the catalyst. Our study also suggests that bridgehead ligand properties, such as bulkiness and aromaticity, can be exploited to alter or even break the linear scaling relationships.

Introduction

Electrocatalytic reduction of protons to form molecular hydrogen is an area of major interest in electrochemistry.[1] Recent studies have shown that synthetic compounds mimicking the active site of the [FeFe] hydrogenase enzyme can provide a cheap alternate route for proton reduction compared to the currently used expensive platinum-based heterogeneous catalysts.[2] The challenge in designing an efficient catalyst for practical application is to reduce the overpotential, which is related to the excess energy needed to drive the electrocatalytic cycle. Already, several artificial di-iron complexes have been synthesized and studied for their ability to catalyze the proton reduction and hydrogen formation but thus far with limited success.[3] Especially, choosing the catalyst ligands from the large pool of available compounds and thereby tuning the properties of the catalyst remains a painstaking matter of trial and error. A detailed understanding of the factors that determine the properties of the catalyst is crucial for designing better biomimetic catalysts for hydrogen fuel production.

Figure shows the active site of the [FeFe] hydrogenase enzyme and a prototypical catalyst complex synthesized to mimic this structure.[4,5] We have recently characterized the structure and the electrochemical properties of such di-iron–hydrogenase models using density functional theory (DFT) calculations and DFT-based molecular dynamics simulations.[6,7] The properties of such synthetic catalysts can be tuned by changing the bridgehead ligand, one or more terminal ligands, the chalcogen atoms, and even by changing the metal ions. The catalytic center of the natural enzyme contains an azadithiolate ligand in the bridgehead position and CO and CN– groups as terminal ligands.[8] Organophosphorus ligands (e.g., PMe3, PPh3, PPhMe2, etc.) are commonly used in synthetic catalysts instead of the cyanide ligands at the terminal positions.[9] Catalysts have also been synthesized with substituted bridgehead ligands. Indeed thus far, most experimental studies focus on changing the bridgehead and terminal ligands and by keeping the Fe2S2 unit intact. However, synthesizing new catalysts with modified ligands and characterizing their properties is a cumbersome process. A fast approach for screening ligands would help tremendously to speed up the optimization of the catalyst design. Predictive models based on computational methods can play an important role in this aspect. The chemical properties of new di-iron complexes can be computed with quantum chemical methods before these compounds have been synthesized and experimentally characterized. Computational analysis can also initiate new ideas for selecting ligands.

Figure 1

(Left) Chemical motif of the catalytic site in the [FeFe] hydrogenase enzyme, showing the bridging azadithiolate ligand, the cysteine-linked [Fe4S4] cubane cluster, and the vacant “X” binding site. (Right) Prototypical biomimetic hydrogenase catalyst, with a benzene group as the bridgehead ligand and carbonyl ligands at the terminal positions.

In heterogeneous catalysis, molecular hydrogen formation via proton reduction, aka the hydrogen evolution reaction (HER), has been thoroughly studied. The HER on a metal surface takes place in two steps. In the first step, the proton is adsorbed on the surface (H+ + e– + * → Hads, where “*” denotes a vacant adsorption site), which is called the Volmer reaction.[10] The second reaction step can proceed through two distinct pathways, depending on the metal surface: (1) via a homolytic pathway named the Tafel reaction,[10] in which the second proton is also adsorbed on the surface and the two adsorbed species react to form hydrogen (2 Hads → H2) or (2) via a heterolytic pathway referred to as the Heyrovski reaction,[10] in which the adsorbed species reacts with a proton from the solvent or gas phase to form hydrogen (H+ + e– + Hads → H2). The Volmer reaction is usually fast, while the second (Tafel or Heyrovski) step is usually rate determining.[10] Studies have established linear scaling relationships between the energies of the key intermediates,[11] and computational studies have helped significantly in screening catalytic surfaces for the HER. Such studies for homogeneous molecular catalysts are somehow still rare.[12]

In homogeneous catalysis, Hammett-type linear free energy relationships between the reaction rates or equilibrium constants and the basicities of the reactants or intermediates are used to study the mechanisms and understand the trends in reactivity of the complexes.[13] The Hammett-type linear free energy relationship assumes that there exists a linear correlation between the activation energy and the reaction energy. In other words, the kinetics of the reaction depends linearly on the thermodynamic driving force of the reaction.[14] For example, the study from Jablonskyte et al. established a linear scaling relation between the first oxidation potential (E1) and the rate of first protonation.[15] The basicity of the metal hydride formed in catalytic hydrogenation reactions markedly influences the mechanism of these reactions. Using the Hammett equations, the basicity of the intermediate hydride species can be predicted.

In this study, we have developed a Hammett-type linear scaling relation between the reduction/oxidation potentials and the pKa values for an extended series of artificial hydrogenase catalysts based on the scaffold shown in Figure . The established linear scaling relationship can be used to predict the pKa’s and reduction potentials for the complete reaction cycle. These results can help in the screening of potential candidates from a large pool of contestants to focus on the most promising candidates for a further detailed study. We show that the bridgehead ligands have a more prominent effect on the electrochemical properties of the catalyst than the terminal ligands.

Methods

Computation of the Acidity Constant and the Reduction Potential in Solution

The free energy difference between the protonated and deprotonated forms of a complex is proportional to the pKa of the complex:where R is the ideal gas constant, T is the absolute temperature, and ΔGsolPT is the proton transfer free energy in solution. The latter is computed using a thermodynamic cycle, as schematically drawn for a generic acid, AH, in the left panel of Figure .

Figure 2

Thermodynamic cycles used to incorporate the solvent effects in a reaction free energy calculation of a deprotonation reaction (a) and a reduction reaction (b).

The ΔGsolPT can thus be obtained as a sum of the gas phase reaction free energy, ΔGgasPT, and the difference of the solvation free energies of the product and reactant species involved in the deprotonation reaction:

For the standard molar free energy of a proton in acetonitrile solution, we take the value −252 kcal/mol of Surawatanawong et al. from 2010.[16] To compare our results with previous work in which ΔGsol(H+) = −260.2 kcal/mol, from 2007 by Kelly et al.,[17] is used, our pKa values can simply be lowered by 6.0 units.

The reduction potential of a half reaction, A + e– → A–, is given by eq , where ΔGsolET is the electron transfer free energy of complex A upon reduction in the solvent, and the Faraday constant, F, constitutes the charge of 1 mole of electrons, which is about 96 485 C/mol.Analogous to the pKa calculation, we use a thermodynamic cycle to compute the reduction free energy in solution, ΔGsolET, as shown in the right panel of Figure , so that it can again be calculated as the sum of the reduction free energy in the gas phase, ΔGgasET, and the difference of the solvation free energies of the product and reactant species involved in the reduction reaction:The solvation free energy of the electron is difficult to determine by experiments and computational methods. By using the relative reduction potential with respect to the ferrocene (Fc0/Fc+) electrode in acetonitrile, calculation of the ΔGsol(e–) value can be eliminated.

Previous studies of di-iron complexes have used the above methods to calculate the pKa’s and the reduction potentials.[16,18,19] In the current study, DFT calculations are performed using the BP86[20] exchange correlation functional augmented with Grimme’s D3[21] dispersion correction and a triple-ζ basis set[22] using the TIGHTOPT and TIGHTSCF settings as implemented in the ORCA software.[23] All geometries are first optimized in the gas phase, and analytical frequency calculations are carried out to calculate the zero point energy. The structure is further optimized using the COSMO[24] implicit solvation model to calculate the solvent correction with acetonitrile as the solvent. Electrode potentials are calculated relative to the Fc0/Fc+ electrode in acetonitrile.[25,26]

H2 production catalyzed by the hydrogenase enzyme and by most of the synthesized compounds follows an ECEC mechanism, in which “E” represents an electron transfer step leading to reduction of the H-cluster/catalyst, and “C” refers to a chemical step, which is here a proton transfer reaction. The electron and proton transfer steps may take place sequentially, but they can also occur in a coupled and simultaneous manner. In this study, the consensus ECEC mechanism is considered for all complexes. The most common ECEC mechanism starts from the neutral catalyst, as illustrated in the right-hand-side panel of Figure , and consists of the following reaction steps:

Figure 3

Catalytic reaction cycles for H2 production shown in a schematic phase diagram representation. From left to right, complex A becomes more reduced, and from top to bottom, the complex becomes protonated. Figure inspired by Figure 10 in ref (6). (a) Sequence of four reduction and protonation steps for a catalyst starting from a cationic (A+) resting state. (b) Mechanism for a catalyst starting from the neutral resting state (A).

reaction 1:

reaction 2:

reaction 3:

reaction 4:

The first protonation step (i.e., reaction 2) in this reaction mechanism is similar to the Volmer reaction in heterogeneous catalysis, where the proton binds to the catalyst. There can be multiple protonation sites on the catalyst, somewhat similar to a metal surface. The protonation can occur either at the Fe–Fe metal center forming a bridging hydride or at a single Fe site forming a terminal hydride or at one of the sulfur atoms. It is believed that the key intermediate during the hydrogen production by the natural enzyme is a terminal hydride.[27,28] However, previous computational studies of the proton reduction by artificial [FeFe] hydrogenase mimics have explored the different protonation sites and concluded that the bridging hydride is the thermodynamically most favorable intermediate.[16,29−31] Experimental studies of the protonation of the Fe–Fe metal center have also validated the bridging hydride intermediate as most likely.[15] Hence, in this study, we consider the first protonation reaction to result in the bridging hydride intermediate, the chemical structure of which is shown in Figure .

Figure 4

Chemical structure of the bridging hydride intermediate formed after the first protonation step in the catalytic cycle.

Similar as in heterogeneous HER catalysis, the second protonation reaction can occur via two distinct reaction pathways. Unraveling the preferred reaction mechanism for each complex considered in this study would be too demanding. Moreover, this study focuses on calculating the pKa, which is a linear function of the protonation free energy, a state function that is thus independent of the path. Therefore, we consider henceforth a direct reaction mechanism for the second protonation step, where the proton reacts with the adsorbed hydride to form H2 and the catalyst is regenerated. Note the similarity to the Heyrovski type of mechanism on a catalytic surface.

To investigate the influence of the bridgehead ligand on the redox potentials and acidity constants of the di-iron complex, we use the following 10 bridgehead ligands: benzene dithiolate (bdt), dichlorobenzene dithiolate (cl2-bdt), ethane dithiolate (edt), propane dithiolate (pdt), oxa dithiolate (odt), methylpropane dithiolate (Me-pdt and Me2-pdt), isopropylpropane dithiolate (iPr-pdt and iPr2-pdt), and propylpropane dithiolate (Pr2-pdt). The chemical structures of the hexacarbonyl di-iron complexes with these bridgehead ligands are compiled in Figure . To assess the influence of the terminal ligands, we compare these hexacarbonyl complexes with complexes in which one or more CO ligands are substituted by the commonly used organophosphorus ligands PMe3 and PPhMe2. In this manner, we have created 80 different catalyst complexes for analysis; the chemical structures of these complexes are shown in the Supporting Information, Figure S1.

Figure 5

Chemical structures of biomimetic hexacarbonyl di-iron complexes showing the 10 bridgehead ligands analyzed in this study. See Figure S1 in the Supporting Information for the complete listing of the 80 structures investigated in this work.

Results and Discussion

For catalytic cycle I, which starts from the neutral resting state (see the right-hand-side panel in Figure ), the relevant reduction potentials and pKa’s that determine the thermodynamic driving force for the four reaction steps are E2, pKa2, E4, and pKa4, respectively. In Figure , we plot for each of the 80 di-iron complexes pKa2 versus E2, relevant for the first reduction and protonation steps, and pKa4 versus E4, connected to the second electron and proton transfer reactions. The computed values are shown for the complexes in acetonitrile solvent, color-coded by the bridgehead ligand as indicated in the figure. Two noticeable linear correlations are seen. Generally, complexes that are easily reduced are subsequently more difficult to protonate and vice versa. The red line in the figure is the linear fitting curve through the computed data points, which has the following function:and a squared correlation coefficient, R2, of 0.94. The blue dashed line is the fitting curve with function:and an R2 of 0.95. On closer inspection, two outlier data points may be discerned for the iPr2-pdt bridgehead ligand at (E = −3.5, pKa = 50) and at (E = −2.4, pKa = 42). This deviation from the linear relationship is due to a significant distortion in the molecular structure, which we address in more detail hereafter.

Figure 6

Acidity constant versus reduction potential computed for 80 di-iron complexes for catalytic cycle I. Filled circles, pKa2 versus E2; open squares, pKa4 versus E4. Linear fit functions to the data are shown by a red solid line for pKa2 versus E2 and a blue dashed line for pKa4 versus E4. Data points are color-coded by bridgehead ligand: bdt (black), cl2-bdt (orange), edt (green), pdt (blue), odt (yellow), Me-pdt (red), Me2-pdt (brown), iPr-pdt (pink), iPr2-pdt (purple), and Pr2-pdt (light blue).

The plots show that for each bridgehead ligand, the data points are scattered over a broad range of the pKa and E scales when the terminal ligands are substituted. In particular, replacing more CO groups by PMe3 or PPhMe2 groups at the terminal positions results in more electron donation at the metal centers, which leads to an increase in the required pKa of the proton donating acid in the solution. In other words, less strong acids are needed to protonate the complex during the catalytic process when more phosphine terminal ligands are used. Consequently, due to the linear scaling relationships, more negative reduction potentials are needed for the electron transfer. Previous work has shown that the basicity of the complex indeed increases with the number of phosphine ligands.[32]

Next we proceed with catalytic cycle II, i.e., the sequence of four reaction steps starting from the cationic resting state shown in the left panel of Figure . In Figure , we show the plots of E1 versus pKa1 associated with the first reduction and protonation steps, and E3 versus pKa3 related to the second pair of charge transfer reactions. The data points clearly show more scatter than we saw for the data of catalytic cycle I. In particular, when drawing a linear fit through the data points of E1 versus pKa1 (shown in Figure S2 in the Supporting Information), the data from the three bridgehead ligands Me2-pdt, iPr2-pdt, and Pr2-pdt cluster below the fitting line, while the data from the other bridgehead ligands, bdt, cl2-bdt, edt, pdt, odt, Me-pdt, and iPr-pdt, are grouped above this line. For the data points of E3 versus pKa3, we see the same grouping but with the opposite trend. Not surprisingly, dividing the data in these two subsets results in much better linear fits, which are shown by the solid and dashed lines in Figure . The main similarity among the latter three bridgeheads groups seems to be that they are aliphatic and bulky near the dithiolate bridge compared to the other ligands. Therefore, we will refer to this selection as the “bulky bridgehead ligands” and the group containing all other complexes as the “less bulky bridgehead ligands” (although one might argue that some ligands in this selection, such as cl2-bdt, are rather bulky).

Figure 7

Acidity constant versus reduction potential computed for 80 di-iron complexes for catalytic cycle II. Circles and black linear fit lines show pKa1 versus E1; squares and red fit lines show pKa3 versus E3. Filled symbols and solid lines show the data associated with the bulky bridgehead ligands and open symbols and dashed lines show that of the less bulky bridgehead ligands. Data points are color-coded by bridgehead ligand, with the less bulky ligands being bdt (black), cl2-bdt (orange), edt (green), pdt (blue), odt (yellow), Me-pdt (red), iPr-pdt (pink); and the bulky bridgehead ligands are Me2-pdt (brown), iPr2-pdt (purple), and Pr2-pdt (light blue).

The linear fit functions for E1 versus pKa1 for complexes with the bulky and the less bulky bridgehead ligands are respectively:For E3 versus pK, the fit functions are given by

A comparison of some of the computed redox potentials and acidity constants with available experimental numbers is shown in Table S1 in the Supporting Information. Previous experimental work on di-iron dithiolate complexes by Jablonskyte and co-workers established a relationship between the first oxidation potential and the rate of protonation. They also reported that the stereoelectronic effect of bulky ligands increases the electronic energy level of the highest occupied molecular orbital (HOMO), which enhances the rate of protonation.[15] The slope of the correlation between the protonation rate constant and the oxidation potential for the first step was estimated to be 11.7.[15] In our study the slope between pKa and the oxidation potential (over all data points) is 14.02.

Previous work by Bordwell et al.[33,34] established a linear scaling relation between the pKa and the redox potential for a series of carbon bases and fluorenide anions. A so-called Brønsted type plot of Eox vs pKa is linear with a Brønsted slope near unity when the x and y axes are in the same units. Also, their work established a linear relationship between the rate of the electron transfer and the pKa of the complexes. Similar to the Eox vs pKa plot, the Brønsted coefficient (β) is close to one for the plot of the rate of the electron transfer versus the pKa. Moreover, another study observed that the enthalpy of the protonation of the metal complex correlates linearly with the pKa for single electron transfer substitution reactions.[32] Similar to these results, we also establish such a linear relationship for the complexes considered in this study. The magnitude of the dimensionless Brønsted coefficient can be calculated by multiplying the slope of the pKa versus E curve by 2.303RT/F. The Brønsted coefficients for the pKa versus E plots range between 0.83 and 1.2. The near unity magnitude of the Brønsted coefficient indicates that the activation energy for the protonation of the metal–metal bond correlates with the basicity of the bond.

These relationships allow for the prediction of the pKa’s of the intermediate species before actually synthesizing the complex experimentally and to choose the ligands for synthesizing a complex with desirable pKa values. A catalytic mechanism in which the resting state of the catalyst is in its cationic oxidized form (see Figure a) can be more influenced by the choice of bridgehead ligands than the catalytic cycle with a neutral resting state. We note that some catalysts have a one electron reduced resting state;[5] these cases are not considered in this study.

For the final part of our investigation, we go beyond correlating reduction potentials with the pKa’s of the subsequent protonation steps and try to establish further interesting relationships, for example, among all reduction potentials shown in Figure . Figure shows the correlation between the reduction potentials, E1 to E4, for all complexes. For the less bulky bridgehead ligands, clear linear scaling relationships are seen between E1 vs E2, E1 vs E3, and E2 vs E3. However, the plots of E4 with the other potentials show more scatter, and there is not a generic linear correlation. Significantly better fits for these correlations involving E4 can be obtained if we make a second division into a group containing the aromatic bdt (black dots) and cl2-bdt (orange) bridgehead ligands and a group containing all other ligands (data not shown). For the bulky bridgehead ligands, the correlation is rather poor (see Figure S4 of the Supporting Information). Nevertheless, this means that to estimate the redox properties along reaction mechanism 2 for the less bulky ligands, only one (e.g., the first) reduction potential needs to be known. All other potentials can then be predicted to reasonable accuracy using our linear relationships. This is very helpful, as it allows for a quick screening before having to perform a more accurate study of the catalyst in order to determine the exact mechanism and the intermediate structures.

Figure 8

Correlation between all pairs of reduction potentials, E1 to E4, for all complexes. Solid lines and filled circles show the data associated with the bulky bridgehead ligands and dashed lines and open circles for the less bulky bridgehead ligands.

Interestingly, in addition to the shifted linear trends seen for the bulky bridgehead complexes and the aromatic bridgehead complexes, there are a few other complexes within the less bulky bridgehead ligands that do not show the expected linear scaling behavior. An example is the case of the Fe2(edt)(CO)3(PMe3)3 complex, for which the outliers from the linear trends become particularly evident in the plots of the pKa’s versus the reduction potentials, shown in Figure S5 of the Supporting Information. In particular, all plots involving E4 or pKa4 (except the plot of E4 versus pKa4) show a poor correlation. A closer inspection of the Fe2(edt)(CO)3(PMe3)3 structures reveals that the AH– species has a somewhat distorted geometry involving a broken Fe–S bond, which is different from all other complexes with an edt bridgehead ligand in the AH– state. For comparison, the AH– and AH optimized structures are shown in Figure . To test the effect of this structural change, we have also calculated the single point energy of the AH– species with the same geometry as AH. Now the data points involving E4 or pKa4, shown in Figure S6 of the Supporting Information, are much closer to the linear correlation curves, confirming that indeed the structural changes in the complex during the catalytic cycle can cause large deviations from the linear scaling relationships.

Figure 9

Top view of the spatial geometries of the Fe2(edt)(CO)3(PMe3)3 complex in the AH– (left) and AH (right) states. Note the broken Fe–S bond in the AH– complex, resulting in a somewhat rotated dithiolate ligand and an under-coordinated iron ion (purple balls). Side view of these complexes are shown in the Supporting Information, Figure S9.

Conclusions

In this study, we have established linear scaling relationships between the reduction potentials and the pKa values of bioinspired di-iron hydrogenase catalysts. This type of widely explored Hammett-type linear free energy relationships is effectively utilized in this study with the aim to allow for fast prediction of the catalytic properties of these iron complexes designed for efficient and clean hydrogen fuel production. We find not only such linear scaling rules between the pKa values and the reduction potentials of successive reduction and protonation steps (which in practice may also occur simultaneously as proton-coupled electron transfer reactions) but also between the different reduction potentials and between the different pKa values. However, for the latter scaling behavior, there is no generic linear scaling rule. Within the set of 10 different bridgehead ligands investigated here, we observed different scaling trends for bridgehead ligands that are aliphatic and bulky at the dithiolate bridge in comparison to ligands that are less bulky. Second, we observed a shifted linear scaling of complexes with an aromatic bridgehead ligand for correlations involving the fourth reduction potential, E4. Third, certain bridgehead ligands, such as edt, undergo a relatively large structural change in certain oxidation states, which causes large deviations from the linear scaling relationships. This means that for the classes of ligands for which the linear relationships hold, fast screening studies can be performed, because only one reduction potential has to be measured or calculated to estimate with reasonable accuracy all other reduction potentials and pKa values that are relevant for the molecular hydrogen producing catalytic cycle. On the other hand, however, different classes of bridgehead ligands can provide ways to alter (i.e., shift or tilt) the linear relationships or even break with the linear scaling rules whenever this is desirable. These results serve in the preliminary stage as a screening tool before further detailed study needs to be carried out for the most promising candidates to understand in further detail the reaction mechanisms and the electronic properties of new biomimetic catalysts.