Article pubs.acs.org/JPCC
On the Accuracy of Calculated Reduction Potentials of Selected Group 8 (Fe, Ru, and Os) Octahedral Complexes Lubomír Rulíšek* Institute of Organic Chemistry and Biochemistry, Gilead Sciences Research Center & IOCB, Academy of Sciences of the Czech Republic, Flemingovo nám. 2, 166 10 Praha 6, Czech Republic S Supporting Information *
ABSTRACT: The theoretical calculations of reduction potentials for the [M(H2O)6]2+/3+, [M(NH3)6]2+/3+, [M(en)3]2+/3+, [M(bipy)3]2+/3+, [M(CN)6]4−/3−, and [MCl6]4−/3− systems (M = Fe, Os, Ru) were carried out. The DFT(PBE)/ def2-TZVP//DFT(PBE)/def2-SVP quantum chemical method was employed to obtain presumably accurate ionization energies, whereas the conductor-like screening model for real solvents (COSMO-RS) was selected as the most suitable method for calculations of solvation energies of the oxidized and reduced forms of the studied species. It has been shown that COSMO-RS may overcome problems related to directionality of hydrogen bonds in the second solvation sphere that previously led to errors of ∼1 V for the [Ru(H2O)6]2+ complex employing PCM-like models. Thus, most of the values for (2+) → (3+) oxidations are now within 0.1−0.2 V from the experimental data, once the anticipated spin−orbit coupling effects in Os complexes (downshifting the calculated reduction potentials by ∼0.3 V) are taken into account. The robustness of the DFT(PBE)/COSMO-RS computational protocol is further verified by showing that reduction potentials obtained for selected 2+/3+ redox pairs with and without the inclusion of explicit second-sphere water molecules are almost identical. At the same time, it must be admitted that the calculated values of reduction potentials for systems involving quadruple charged species, exemplified here by [M(CN)6]4−/3− and [MCl6]4−/3− redox pairs, might still not be within the grasp of contemporary solvation models, possibly due to the large values of solvation energies of their reduced (4−) forms that are in the range of 700−750 kcal mol−1 (30−33 eV) and possibly larger errors inherent in their calculations. Finally, a comparison is made with M06-L//SMD computational protocol, which is also shown to correct some of the deficiencies of previous PCM models. E0abs(NHE) values are 4.26 V as suggested by Cramer and Truhlar,1 and 4.281 V as advocated by Isse and Gennaro.8 The Gibbs free energies in eq 1 can be calculated by using the following formula:
1. INTRODUCTION Theoretical calculations of reduction potentials have attracted the attention of computational chemists for many years (see, e.g., refs 1−3 and the references therein). These efforts were mostly driven by the exciting possibility of connecting the experimental data with one (or several) plausible structural representations of the studied system.4 A successful structural interpretation of the experimental electrochemical data may then assist in the elucidation of the reaction mechanisms of redox-active bioinorganic systems5,6 or electrochemical reactions.7 Also, the theoretical calculations can be potentially useful in the prediction of these values for unknown compounds (possibly saving some synthetic efforts). For one-electron processes, reduction potentials can be expressed as:
G = Eel + ΔGsolv + EZPE − RT ln(qtransqrotqvib)
Combining eqs 1 and 2, it can be seen that one essentially computes the in vacuo energy difference between the reduced and oxidized structures (the ionization energy or electron affinity for oxidations and reductions, respectively), the difference in their solvation energies, and the enthalpic and entropic contributions to gas-phase Gibbs free energies, mostly calculated at the harmonic oscillator/rigid rotor level. In recent years, it has been shown that fairly encouraging results can be obtained for smaller molecules (see, e.g., refs 9−11), whereas larger systems, such as metalloproteins, are probably still elusive to computational efforts, and it is difficult
0 E 0 [V] = 27.21(Gox [au] − Gred [au]) − Eabs (NHE) [V]
(1)
where Gox/red is the Gibbs free energy of the oxidized (reduced) form and E0abs(NHE) is the absolute potential of the normal hydrogen electrode (NHE). The latest recommended © 2013 American Chemical Society
(2)
Received: July 9, 2013 Revised: July 25, 2013 Published: July 26, 2013 16871
dx.doi.org/10.1021/jp406772u | J. Phys. Chem. C 2013, 117, 16871−16877
The Journal of Physical Chemistry C
Article
shown to be systematic and approximately equal to −0.07 V shifts in the calculated reduction potentials for the Ru2+/3+ and −0.3 V for the Os2+/3+ octahedral systems.17 The aim of this work is the application of two modern solvation methods: “semiexplicit” COSMO-RS15 and densitybased SMD model20 to calculate presumably more accurate solvation energies of doubly and triply charged species and thus greatly extend the applicability of the protocol described above (e.g., to various bioinorganic systems). The COSMO-RS can be considered as an improvement over the standard PCM-like models, because it treats both the solute and the solvent molecules as the interacting entities by matching their so-called σ-profiles. Therefore, it might be expected that the COSMORS solvation method may account for a part of the secondsphere directionality in the hope that the calculated solvation energies of [MX6]2+/3+ systems might be more accurate; the same might hold true for the SMD. It can be noted that similar efforts were reported21,22 for theoretical calculations of stability constants of miscellaneous inorganic systems and a broad range of metal ions (Mn2+, Fe2+, Co2+, Ni2+, Cu2+, Zn2+, Cd2+, and Hg2+), and it was shown that the COSMO-RS protocol yields satisfactory results. The complexes selected represent a common, but reasonably diverse, set of systems, the [M(H2O)6]2+/3+, [M(NH3)6]2+/3+, [M(en)3]2+/3+, and [M(bipy)3]2+/3+ complexes, to which two systems with higher formal charges, [M(CN)6]4−/3− and [MCl6]4−/3−, have been added to examine the possible extensions of the computational protocol used in the study.
(and most likely rather accidental or fortuitous) to obtain quantitative agreement between the experimental and calculated data from the first principles, mostly because of the insufficient treatment of the protein surroundings (see ref 12 and the discussion therein). Many computational and synthetic chemists are particularly interested in the complexes of group 8 elements (Fe, Ru, Os) as they belong to common redox systems and can also be used in synthetic chemistry, for example, as single electron transfer (SET) oxidants.13 Recently, we have calculated reduction potentials corresponding to oxidations of a large series of ferrocene derivatives14 and have shown that the agreement between the calculated and experimental data was better than 0.03 V (mean unsigned error) using the protocol based on DFT(PBE)/def2-TZVP calculations and the conductor-like screening model for real solvents (COSMO-RS) method for solvation energy calculations.15 Moreover, from the limited comparison of experimental and calculated ionization energies for the parent system, ferrocene/ferrocenium couple (IEexp = 6.71 ± 0.08 eV16 vs IEcalc = 6.75 eV), it seems likely that the above remarkable agreement obtained is not caused by the cancellation of errors, but both major components (gas-phase ionization and solvation energy) of the free energy change upon oxidation might be theoretically predicted to a very good accuracy. Such quantitative agreement implies that quantum chemical calculations can be used as a predictive tool and save considerable synthetic efforts. Less satisfactory agreement was obtained in our previous study dealing with the reduction potentials of a series of Ru2+/3+ and Os2+/3+ octahedral complexes.17 Therein, we used the COSMO (one of the standard PCM models) method and indicated that the calculated solvation energies differences, ΔΔGsolv(red/ox), are likely to be inaccurate (up to 1 V for the [Ru(H2O)6]2+/3+ complex) unless the second solvation sphere comprising 12−14 water molecules is taken explicitly into account.17 The same phenomenon was observed earlier by Jaque et al.18 and revisited in a recent work19 from the same group dealing with Pourbaix diagrams for ruthenium-based water-oxidation catalysts. In the latter study, authors employed two new density functionals, M11-L and M11, and the SMD implicit solvation model. This combination of DFT functional and solvation model yielded the reduction potential of 0.2 V for the Ru2+|Ru3+ redox pair, in almost quantitative agreement with the experiment, but, again, the second shell of water molecules had to be used (i.e., Ru(H2O)182+/3+ system). In our work,17 we also showed that the “second-sphere solvation effect” vanishes with the growing size and presumably lower polarity of the ligands (on the example of the series of H2O, NH3, en, and bipy ligands, for which the “second-sphere” effects decreased from above ∼1 V to almost zero). Therefore, it appears that the directionality of the hydrogen bonding between the first and second spheres, mostly exemplified in [M(H2O)6]2+/3+ complexes, prevents PCM-like models from being satisfactorily accurate for the calculations of the reduction potentials of these types of systems unless an elaborate treatment (such as the inclusion of the second-sphere water molecules concomitantly invoking a question of the conformational sampling of these less tightly bound and structurally less clearly defined solvent species) is applied. At this point, it should be noted that there are also problems related to the spin−orbit coupling in both Ru(II/III) and Os(II/III) systems that systematically influence the ionization energies (and reduction potentials) of these systems. However, these were
2. COMPUTATIONAL DETAILS Quantum chemical calculations were performed using the TURBOMOLE 6.423 and Gaussian 09 (revision D1)24 programs. Geometry optimizations were carried out at the DFT level, using the Perdew−Burke−Ernzerhof (PBE) functional.25 DFT (PBE) calculations were expedited by expanding the Coulomb integrals in an auxiliary basis set, the resolutionof-identity (RI-J) approximation.26 For the geometry optimization, the def2-SVP basis set was employed on all of the atoms,27 whereas the def2-TZVP basis set was used for the final single-point calculations to obtain presumably accurate molecular energies.28 The single point calculations were carried out using the PBE, M06-L,29 and M11-L30 functionals. For all open-shell systems, the unrestricted Kohn−Sham (UKS) formalism has been used. Solvation effects were taken into account by using the SMD20 and COSMO-RS (conductor-like screening model for real solvents)15,31 methods. The SMD calculations were carried out using the recommended protocol in Gaussian 09. COSMORS calculations were done using TURBOMOLE 6.4 for the COSMO calculation32 with εr = ∞ (ideal screening) and the COSMOtherm33 program for the subsequent COSMO-RS calculation. The recommended protocol involving the Becke− Perdew (B−P) functional34 for the in vacuo and the εr = ∞ calculations along with the def-TZVP basis set was used. Gibbs free energy and reduction potentials were calculated according to eqs 1 and 2 (vide supra), where Eel is the in vacuo energy of the system and Gsolv is the solvation free energy (calculated using the RI-BP/def-TZVP(COSMO-RS, ε = 1, ε = ∞) or SMD method as described above). EZPVE is the zeropoint energy, and −RT ln(qtransqrotqvib) accounts for thermal corrections to the enthalpy and entropic terms. EZPVE and −RT ln(qtransqrotqvib) are obtained from a frequency calculation using the same method and software as for the geometry 16872
dx.doi.org/10.1021/jp406772u | J. Phys. Chem. C 2013, 117, 16871−16877
The Journal of Physical Chemistry C
Article
Table 1. Calculated Reduction Potentials of Iron(II/III), Ruthenium(II/III), and Osmium(II/III) Complexes Using the DFT(RI-PBE)/def2-TZVP//COSMO-RS Protocola system
total spin (S) {red; ox} 2+/3+
[Ru(H2O)6] [Ru(NH3)6]2+/3+ [Ru(en)3]2+/3+ [Ru(bipy)3]2+/3+ [Ru(CN)6]4−/3− [RuCl6]4−/3− [Os(H2O)6]2+/3+ [Os(NH3)6]2+/3+ [Os(en)3]2+/3+ [Os(bipy)3]2+/3+ [Os(CN)6]4−/3− [OsCl6]4−/3− [Fe(H2O)6]2+/3+ [Fe(NH3)6]2+/3+ [Fe(en)3]2+/3+ [Fe(bipy)3]2+/3+ [Fe(CN)6]4−/3−
{0, {0, {0, {0, {0, {0, {0, {0, {0, {0, {0, {0, {2, {2, {2, {0, {0,
1/2} 1/2} 1/2} 1/2} 1/2} 1/2} 1/2} 1/2} 1/2} 1/2} 1/2} 1/2} 5/2} 5/2} 5/2} 1/2} 1/2}
Δ(EZPVE − RT ln Q) [eV]d
IE [eV]b,c 15.56 14.47 13.59 12.15 −7.50 −9.40 14.94 14.02 13.21 11.97 −7.61 −10.52 16.11 14.96 13.86 11.94 −8.25
ΔGsolv(red/ox) [eV]e
−0.10 0.14 0.01 −0.01 −0.05 0.21 −0.19 0.18 0.02 0.01 −0.05 0.23 0.24 0.15 0.04 0.01 −0.03
g
(15.76) (14.56) (13.66) (11.95) (−7.56) (−9.37) (14.95) (14.00) (13.16) (11.72) (−7.68) (−10.04) (16.63) (15.25) (14.25) (11.69) (−8.50)
−11.11 −10.13 −9.04 −6.27 13.77 14.00 −11.19 −10.08 −9.02 −6.26 13.72 14.05 −11.12 −10.11 −8.98 −6.31 14.09
E0(calc) [V]
E0(exp) [V]f
0.06 0.19 0.28 1.59 1.94 −0.01 −0.73 −0.16 −0.07 1.44 1.78 −0.52 0.94 0.71 0.63 1.35 1.54
0.23 0.10 0.21 1.24 0.86
(−10.02) (−9.83) (−8.87) (−6.32) (12.30) (12.62) (−9.96) (−9.80) (−8.84) (−6.38) (12.26) (12.60)
0.80
0.77
1.03 0.54
a The computational protocol described above was used (eqs 1 and 2). bFor the sake of convenience, all energy terms are expressed in eV (V). cThe adiabatic ionization energies calculated in vacuo at the PBE/def2-TZVP//PBE-def2-SVP level. dThe difference in the gas-phase chemical potential (including ZPVE) calculated at the PBE/def2-SVP level. eThe difference in the solvation energies calculated using the COSMO-RS method; the PCM-like (COSMO) values are in italics in parentheses (taken from ref 17). fThe experimental data are taken from ref 35. gThe B3LYP/def2-TZVP values of ionization energies, taken from ref 17, for the Ru and Os systems and calculated in this work for the Fe systems (in parentheses).
Table 2. Calculated Reduction Potentials of Iron(II/III), Ruthenium(II/III), and Osmium(II/III) Complexes Using the DFT(M06-L)/def2-TZVP//SMD Protocola system
IE [eV]b
Δ(EZPVE − RT ln Q) [eV]c
ΔΔGsolv(red/ox) [eV]d
E0(calc) [V]
E0(exp) [V]e
[Ru(H2O)6]2+/3+ [Ru(NH3)6]2+/3+ [Ru(en)3]2+/3+ [Ru(bipy)3]2+/3+ [Ru(CN)6]4−/3− [RuCl6]4−/3− [Os(H2O)6]2+/3+ [Os(NH3)6]2+/3+ [Os(en)3]2+/3+ [Os(bipy)3]2+/3+ [Os(CN)6]4−/3− [OsCl6]4−/3− [Fe(H2O)6]2+/3+ [Fe(NH3)6]2+/3+ [Fe(en)3]2+/3+ [Fe(bipy)3]2+/3+ [Fe(CN)6]4−/3−
15.60 14.34 13.48 11.96 −7.64 −9.57 14.94 13.88 13.08 11.78 −7.73 −10.08 15.92 14.76 13.82 11.98 −8.32
−0.10 0.14 0.01 −0.01 −0.05 0.21 −0.19 0.18 0.02 0.01 −0.05 0.23 0.24 0.15 0.04 0.01 −0.03
−11.14 −10.22 −9.41 −6.86 12.10 11.62 −11.07 −10.18 −9.43 −6.89 12.06 11.61 −10.99 −10.19 −9.42 −6.86 12.30
0.07 −0.01 −0.20 0.81 0.12 −2.02 −0.60 −0.40 −0.61 0.61 0.00 −2.51 0.89 0.43 0.16 0.84 −0.33
0.23 0.10 0.21 1.24 0.86
0.80
0.77
1.03 0.54
a The computational protocol described above was used (eqs 1 and 2); the spin states were the same as in Table 1. bThe adiabatic ionization energies calculated in vacuo at the M06-L/def2-TZVP//PBE-def2-SVP level. cThe difference in the gas-phase chemical potential (including ZPVE) calculated at the PBE/def2-SVP level. dThe difference in the solvation energies calculated using the SMD method. eThe experimental data are taken from ref 35.
puted at the B3LYP/def2-TZVP level; both are considered as important for subsequent discussions (for iron(II/III) systems, the B3LYP values were calculated in this work). The equilibrium structures of the selected species studied in this work are depicted in Figure 1. Beginning our discussion with the Ru(II/III) systems and the DFT(PBE)/COSMO-RS computational protocol, it can be seen that the agreement between calculated and experimental values is mostly satisfactory. Even for the “notoriously” problematic [Ru(H2O)6]2+/3+ system, the ΔE0(exp/calc) is 0.17 V (cf., ∼1 V calculated using the PCM-like models).17,18
optimization (RI-PBE/def2-SVP level), 298 K, and 1 atm using the ideal-gas approximation.
3. RESULTS AND DISCUSSION The calculated reduction potentials of the [M(H2O)6]2+/3+, [M(NH 3 ) 6 ] 2+/3+ , [M(en) 3 ] 2+/3+ , [M(bipy) 3 ] 2+/3+ , [M(CN)6]4−/3−, and [MCl6]4−/3− octahedral complexes are summarized in Tables 1 and 2. Besides the data provided by the DFT(PBE)/COSMO-RS calculations, also two sets of data taken from ref 17 were included in Table 1: the PCM (COSMO) solvation energies and ionization energies com16873
dx.doi.org/10.1021/jp406772u | J. Phys. Chem. C 2013, 117, 16871−16877
The Journal of Physical Chemistry C
Article
tested in this work. There is only one experimental value with which to compare; see the [Ru(CN) 6 ]4−/3− complex. Apparently, the calculated value (1.94 V) is more than 1 V higher than the experimental one (0.86 V). This might be caused either by the incapability of the solvation model to describe the solvation of quadruple charged species, ΔGsolv(reduced form) ≈ 700−750 kcal mol−1 (i.e., 30−33 eV), or by the inadequacy of the model system (involving species with 4- charge in vacuum, which leads, e.g., to negative ionization energies, counterbalanced by the solvation energies, as can be seen in Table 1). Unfortunately, a larger deviation between the calculated and experimental values is also observed for the [Ru(bipy)3]2+/3+ system. The value of 1.59 V was obtained in water and 1.55 V in acetonitrile (not shown in Table 1), which can be compared to the experimental values of 1.24 V in the CRC tables35 (the solvent used could not be identified) or 1.302 V obtained for the E0{[Ru(bipy)3]2+/3+} in the equimolar mixture of water and acetonitrile.11 In any case, it does appear that the DFT(PBE)/COSMO-RS values are ∼0.3 V higher than the experimental ones. As has been shown, the spin−orbit coupling lowers this value by 0.04 V.17 The rest might originate in an overestimation of ionization energies by the PBE functional. In fact, the SOC-corrected B3LYP/ COSMO-RS value (∼1.35 V in water; 1.31 V in acetonitrile) would be in almost quantitative agreement with the experiment (ΔE0(exp/calc) < 0.1 V). As mentioned above, the PCM(COSMO) solvation energies (if used instead of COSMO-RS) would not significantly deteriorate this agreement. Concerning the studied osmium(II/III) complexes, it was possible to find only one experimental value available, for [Os(bipy)3]2+/3+ complex(es). A part of a fairly large overestimation of its reduction potential computed at the DFT(PBE)+COSMO-RS level (ΔE0(exp/calc) = +0.64 V) is readily understood as the effect of spin−orbit coupling that has been shown to lower its value by 0.25 V.17 Keeping this relativistic effect in mind, we obtain almost the same overestimation of the [Os(bipy)3]2+/3+ reduction potential as for the [Ru(bipy)3]2+/3+ system (0.39 V) using the PBE functional for the ionization energy, whereas the SOCcorrected B3LYP/def2-TZVP +COSMO-RS value would point to the value of 0.94 V, reasonably close to the experimental value. The reduction potentials of the corresponding iron(II/III) systems have been calculated as well (Table 1). Here, the agreement is somewhat less satisfactory. While the SOC effects are expected to play a minor role in the calculated reduction potentials, the proximity of the spin states for many of these complexes makes these calculations nontrivial. It can be mentioned that the calculations presented in Table 1 were carried out for high-spin states (quintets and sextets for the Fe2+ and Fe3+, respectively) for the [Fe(H2O)6]2+/3+, [Fe(NH3)6]2+/3+, and [Fe(en)3]2+/3+ systems and closed-shell lowspin states for the [Fe(CN)6]4−/3− and [Fe(bipy)3]2+/3+ systems (singlets and doublets), which are presumably the ground electronic states of these molecules. The same as discussed above holds true for the [Fe(CN)6]4−/3− system, obviously poorly described by the computational protocol used, which leads to an overestimation of the reduction potential by 1 V for the [Fe(CN)6]4−/3− pair. Again, whether this is attributable solely to the incapability of the current models to describe the solvation of quadruplecharged species or whether it is a failure of the model itself is a nontrivial question (partially addressed in this study, vide infra).
Figure 1. The selected molecular structures (equilibrium geometries) studied in this work: (a) [Ru(NH3)6]2+/3+, (b) [Ru(H2O)6]2+/3+, (c) [Ru(en)3]2+/3+, (d) [Ru(bipy)3]2+/3+, (e) [Ru(CN)6]4−/3−, and (f) [RuCl6]4−/3−. The Ru−L distances, in Å, are depicted for both the reduced and the oxidized (in italics in parentheses) forms. The coordinate files for all of the structures studied in this work can be found in the Supporting Information.
Using B3LYP/def-TZVP ionization energy (Table 1, the values in italics), it would be even smaller (0.26 V calculated vs 0.23 V experimental). The same holds true for the [Ru(NH3)6]2+/3+ and [Ru(en)3]2+/3+ systems, where the agreement is well below 0.1 V (and would be also quite satisfactory using B3LYP ionization energies, which differ by
