|Home | About | Journals | Submit | Contact Us | Français|
Quantum chemical calculations of metal clusters in proteins for redox studies require both computational feasibility as well as accuracies of at least ~50 mV for redox energies but only ~0.05 Å for bond lengths. Thus, optimization of spin-unrestricted density functional theory (DFT) methods, especially the hybrid generalized gradient approximation functionals, for energies while maintaining good geometries is essential. Here, different DFT functionals with effective core potential (ECP) and full core basis sets for [Fe(SCH3)4]2-/1- and [Fe(SCH3)3]1-/0, which are analogs of the iron-sulfur protein rubredoxin, are investigated in comparison to experiment as well as other more computationally intensive electron correlation methods. In particular, redox energies are calibrated against gas-phase photoelectron spectroscopy data so no approximations for the environment are needed. B3LYP gives the best balance of accuracy in energy and geometry compared B97gga1 and BHandH and is better for energies than Møller-Plesset perturbation theory series (MP2, MP3, MP4SDQ) and comparable to coupled cluster [CCSD, CCSD(T)] methods. Of the full core basis sets tested, the 6-31G** basis sets give good geometries, and addition of diffuse functions to only the sulfur significantly improves the energies. Moreover, a basis set with an ECP on only the iron gives a less accurate but still reasonable geometries and energies.
Proteins containing iron-sulfur clusters, which consist of one to four irons tetrahedrally coordinated by cysteine residues and inorganic sulfur atoms,1 are an important class of electron transfer proteins in a wide variety of processes such as photosynthesis, respiration, nitrogen fixation, and hydrogen metabolism. They also play critical roles in enzymatic processes, regulation of gene expression, generation of radicals, and delivery of sulfur and iron for the synthesis of other proteins.2-5 Since redox chemistry is important for many of these proteins, identification of the molecular determinants of the reduction potentials of iron-sulfur redox sites is important for a fundamental understanding of their function as well as for engineering reduction potentials by rational design of mutations.
The physical origins of differences in reduction potentials of proteins are often difficult to identify due to the complexity of the protein environment. The reduction potential E° can be obtained from the free energy ΔG° of the reduction reaction
A common approach is to divide ΔG0 into the intrinsic reduction free energy (ΔGint) of the redox site independent of the protein, the extrinsic reduction free energy (ΔGenv) produced by surrounding protein and solvent at the redox site, and the perturbative interaction between the two (ΔGint/env)6
where n is the number of electrons and F is Faraday’s constant. If ΔGint/env is small, then ΔGint is approximately the reduction free energy of the clusters in the gas phase. Given the complexity of the protein environment, computational methods can play an important role and it is essential that quantum chemical methods reproduce ΔGint for analogs in the gas phase before attempting to apply them in the protein environment. The key issue is to obtain ΔGint with accuracies of at least ~50 mV. In addition, while geometry is obviously an important indicator of the quality of a calculation, accuracies of only ~0.05Å in bond lengths are more than sufficient for modeling redox effects on the proteins. Moreover, since for instance the change in Fe—S bond length upon reduction indicates the increase in ferrous character of the bond, systematic errors in geometry are preferable to random errors.
Significant progress in quantum chemical calculations of metalloprotein redox sites has been made. Density functional theory (DFT)7 is now recognized as both more efficient and more accurate for treating transition-metal systems than some conventional ab initio molecular orbital (MO) methods such as Hartree-Fock (HF) and Møller-Plesset perturbation (MPn) methods, especially for systems containing first-row transition-metals.8 However, since Mössbauer, electron paramagnetic resonance (EPR), and PES experiments have revealed unusual electronic structures of iron-sulfur clusters with spin-coupled interactions between iron sites,2-4 spin-unrestricted methods are required. The spin-spin interactions may be treated by the broken-symmetry (BS) DFT approach,9 in which the ground state is based on a single determinant wavefunction. Noodleman and co-workers have investigated iron-sulfur clusters constrained to experimental geometries using local density approximation (LDA) DFT with nonlocal (NL) corrections for exchange and BS-DFT and spin projection methods, which give good electronic structures and Heisenberg spin coupling contributions to magnetic properties.10,11 However, LDA and non-hybrid generalized gradient approximation (GGA) methods usually give good geometries but poor energies7,12 and systematic studies of the accuracy of the hybrid GGA functionals in BS-DFT for both structural and energetic properties are lacking despite many DFT studies of iron-sulfur proteins and analogs.13-15
Recently, understanding of the intrinsic energetics has been greatly enhanced by gas-phase photoelectron spectroscopy (PES), which is a powerful experimental technique for determining the vertical and adiabatic detachment energies (VDE and ADE, respectively)16,17 of oxidation. Since the entropic contribution is small, ΔGint is approximately equal to the intrinsic reduction energy ΔEint,
where λoxd is the intramolecular relaxation energy of the reaction. The calibration of gas phase calculations against gas phase PES data rather than electrochemical redox potentials in solution means no approximations are needed for the environment. This type of calibration of DFT calculations for iron-sulfur clusters is leading to several important types of studies. First, DFT calculations can be used to study a wide range of clusters and their redox couples18-21 since PES experiments are limited to studies of oxidation of analogs that are sufficiently stable in the gas phase, which are not always the relevant redox couples in proteins. Moreover, our calculations indicate that absolute reduction potentials obtained from gas phase calculations of ΔGint and continuum electrostatic calculations of ΔGext are in good agreement with experiment.22 Second, combined DFT/PES studies of analogs with specially designed ligands can be used to study different physical effects individually such as hydrogen bonds on reduction potentials20 and ligand differences on fission. Finally, the calibration is important for QM/MM simulations of iron-sulfur proteins.
The simplest iron-sulfur redox site is [Fe(SCys)4]1-/2-, which has a single iron core referred to as [1Fe-0S], or [1Fe] for short, and is found in the small proteins rubredoxin (Rd) and desulforedoxin as well as in larger proteins such as desulfoferrodoxin, rubrerythrin, and nigerythrin.23 Calibration of DFT methods for analogs of this site are important because their BS ground state is a pure high-spin state without any spin contamination. Here, a systematic study is reported of the performance of hybrid GGA methods with different basis sets in comparison to experimental measurements for structural and energetic properties of [Fe(SCH3)4]n- and [Fe(SCH3)3]n-. The energies were calibrated against recent gas phase PES measurements of the electronic structure, VDE, and ADE by Wang and co-workers.16,17 In addition, the energies were compared with calculations at the Møller-Plesset perturbation theory series and coupled cluster levels of theory.
The calculations were performed for the [Fe(SCH3)4]2-/1-0 and [Fe(SCH3)3]1-/0 analogs of the rubredoxin redox site. The replacement of the cysteine ligands by thiolate SCH3 groups has been shown to be a good model for geometries.9 Here, the “D2d” ML4 structure, which is an idealization of the geometry of the Rd redox site,24 was utilized for the initial structures of [Fe(SCH3)4]2-/1-/0. In addition, the “D3h” ML3 structure, which is an idealization of the experimental geometry of [Fe(SC6H2-2,4,6-tBu3)3]1-,25 was utilized for the initial structures of [Fe(SCH3)3]1-/0.
Three density functionals with better performance in geometry and energetics for small molecules12 were examined: Becke’s original three parameter fit26 using Lee-Yang-Parr correlation functional (B3LYP),27 Becke’s “half and half” functionals (BHandH),26 and the modified Becke’s 10-parameter functional (B97gga1).12,28 In addition, seven different basis sets were examined. In the first, referred to here as ECPDZ, the iron is represented by the double-ζ (DZ) effective core potential (ECP) basis set by Couty and Hall29 modified from the Hay and Wadt basis set.30 The carbons and sulfurs have the ECPs and double-ζ basis sets with polarization functions (ECPDZP) of Stevens, Basch, and Krauss31 and hydrogens have the Dunning-Huzinaga (31) double-ζ basis set.32,33 In the second, referred to here as ECPDZ+DZVP, the iron has the same basis set as in ECPDZ, and the basis sets for the hydrogen, carbon, and sulfur atoms are the all electron Dgauss DZVP polarized DFT orbital basis sets.34 The third and fourth, DZVP2, and 6-31G**, respectively, are the all electron Dgauss DZVP2 polarized DFT orbital basis sets34 and the standard 6-31G basis sets,35,36 respectively. In the fifth and sixth, referred to here respectively as DZVP2(++)S and 6-31(++)SG**, sp-type diffuse functions from the 6-31++G** basis sets35,37 are added to the sulfur atoms of DZVP2 and 6-31G**, respectively. In the seventh, referred to here as 6-31G(3df)S**, three d-type and one f-type polarization functions from the 6-31G(3df,3pd) basis set35,36 are added to the sulfuratoms of 6-31G**.
The calculated VDE and ADE were refined at the different levels of theory with different basis sets. Since the energies were recalculated for certain geometries, the notation Method(E)/Basis(E)//Method(G)/Basis(G) will be used to denote an energy calculation at the Method(E) level with the Basis(E) sets for the geometry optimized at the Method(G) level with the Basis(G) sets. In addition to the hybrid GGA functionals, the energies were also recalculated at the frozen core MP2, MP3, MP4SDQ, coupled cluster singles and doubles (CCSD), and CCSD with perturbative corrections for triples [CCSD(T)] levels.38 The zero point energy (ZPE) corrections were taken into account only in the final results because they are small and do not significantly affect the comparisons of calculations reported here.
A detailed analysis of the electronic structures of the analogues studied here has been presented elsewhere,18 the results are summarized briefly here to aid the discussion. In a roughly Td weak ligand field41 of the tetrahedral [Fe(SCH3)4]n- site, the t2 and e orbitals have a small splitting energy so that the iron high-spin state is favorable. Taking exchange interaction between α and β spins into account, the spin polarization splits the occupied α-spin orbitals below the β-spin orbitals. Consequently, the molecular orbital interaction between the high-spin iron and the thiolates leads to a less common inverted level pattern, in which the S(3p) ligand orbitals lie below the Fe(3d) minority β-spin orbitals and above the Fe(3d) majority α-spin orbitals.10 The oxidation from [Fe(SCH3)4]2- to [Fe(SCH3)4]1- involves a high-lying Fe(3d) minority β-spin electron orbital, while the oxidation from [Fe(SCH3)4]1- to [Fe(SCH3)4]0 involves a ligand orbital with S(3p) α-spin character. Overall, n = 2, 1, and 0 oxidation states have high-spin ground states with S = 2, 5/2, and 2, respectively, and the major unpaired electron densities localize on the metal orbitals regardless of oxidation state. Thus, both [Fe(SCH3)4]1- and [Fe(SCH3)4]2- retain tetrahedral character, while [Fe(SCH3)4]0 becomes a distorted tetrahedral structure (C2v) because of strong Jahn-Teller distortion arising from the oxidation of the degenerate S(3p) orbitals.
Similarly, in the roughly D3h weak ligand field41 of the triangular [Fe(SCH3)3]n- site, the Fe 3d orbitals transform into e“, a’1, and e’, and spin polarization splits the occupied α-spin orbitals below the β-spin orbitals. The n = 1 and 0 oxidation states have high-spin ground-states with S = 2 and 5/2, respectively,18 and the major unpaired electron densities localize on the metal orbitals in both oxidation states. The oxidation from [Fe(SCH3)3]1- to [Fe(SCH3)3]0 involves a low-lying Fe dxz minority β-spin electron, which has Fe—S π* antibonding character, so that [Fe(SCH3)3]0 is planar.
The dependence of the geometry on the functional and the basis set was examined and compared to experimental crystal structures. Although PES data is available for [Fe(SCH3)4]1-/0 and [Fe(SCH3)3]1-/0,16,17 experimental structures are not available for the oxidized state of the former or either state of the latter so the experimental structures of [Fe(S2-o-xyl)2]2-/1-,42 [Fe(SCH3)4]1-,43 [Fe(1,2-benzenedithiolate)2(PMe3)]0,25b and [Fe(SC6H2-2,4,6-t-Bu3)3]1- 24 were used. However, [Fe(1,2-benzenedithiolate)2(PMe3)]0 is only moderately representative of [Fe(SCH3)4]0 because the FeS4 is nearly planar in the former with a S = 1 spin state whereas it is tetragonal with a high spin (S = 2) ground state in the latter. The calculated geometry of [Fe(SCH3)4]1- and [Fe(SCH3)3]1- are shown in Figure 1 for reference. Overall, the BHandH, B3LYP, and B97gga1 optimized geometries using either ECP or full basis sets of [Fe(SCH3)4]2-/1-/0 (Table S1-5) and [Fe(SCH3)3]1-/0 (Table S1-5, no BHandH) were closer to experiment than previous calculations at the UHF and MP2 levels,18 although all still slightly overestimated the bond lengths.
Focusing first on the functionals, B3LYP and B97gga1 were general superior to BHandH for geometries regardless of the basis set. An important geometric feature that the calculations must reproduce is the Fe—S bond length (Figure 2, Table S1). B97gga1 generally gave the best Fe-S bond lengths while B3LYP tended to be about 0.02Å longer. However, the reduced sites, [Fe(SCH3)4]2- and [Fe(SCH3)3]1-, appeared to be the most sensitive to the functional so that B3LYP was better at reproducing the decrease in bond length upon oxidation. The other internal coordinates generally showed less variation with the type of hybrid GGA functional (Table S2-5). The largest variation was for the thiolate ligand torsion angle, ϕC—S—Fe—S, of [Fe(SCH3)3]1- . Only B3LYP/ECPDZ+DZVP, B97gga1/ECPDZ+DZVP, and B97gga1/6-31G(3df)s** gave values of ~4°, similar to the experimental value of 2°24 whereas the other methods gave values between 30 and 41°, apparently due to stabilization by decreasing the Fe—S π antibonding interaction.18 On the other hand, the small torsion in the experimental complex, [Fe(SC6H2-2,4,6-tBu3)3]1-, appears to caused by steric interference of the ligands so that it is questionable that [Fe(SCH3)3]1- would have a small torsion also. In addition, the CCSD(T)/6-31G** energy of the B97gga1/6-31G** geometry was higher than the CCSD(T)/6-31G** energy of the B3LYP/6-31G** structure (Table S7), indicating that CCSD(T) found the B3LYP geometries more favorable. Moreover, B97gga1 is shown below to give poor energetics.
Now focusing on the basis sets, the larger basis sets generally gave superior geometries, as expected. The Fe—S bond lengths (Figure 2, Table S1) tended to improve as the basis set size increased especially on going to the full core basis sets, but again the reduced [Fe(SCH3)4]2- and [Fe(SCH3)3]1- sites were slightly more sensitive. For B3LYP, the full core 6-31G** and DZVP2 basis sets gave comparable results for the Fe—S bond but 6-31G** gave slightly better decrease in the Fe—S bond length upon oxidation mainly because the DZVP2 Fe-S bond for [Fe(SCH3)4]2- is too long by ~0.07 Å. Adding sp-type diffuse functions to the sulfur atoms changed the Fe—S bond lengths only slightly, while adding more polarization functions caused them to shorten, closer to experiment. Apparently, the radial flexibility obtained from the diffuse functions was less important than angular flexibility from the additional polarization functions for obtaining the correct geometries. Moreover, of the other internal angles, the S—C bond (Table S2) shortens by ~0.015 Å between ECPDZ and ECPDZ+DZVP where the full core replaces the ECP on just the ligand atoms but changes insignificantly on going to either the full core DZVP2 or 6-31G** basis sets; thus, ECPDZ+DZVP is a reasonable compromise for the optimization of larger molecular systems. Also, the B3LYP/DZVP2 S—C bond for [Fe(SCH3)4]2- is about 0.015 Å long, which is a greater error than for B3LYP/6-31G** S—C bonds of any of the analogs.
Overall, the B3LYP method using ECPDZ+DZVP, 6-31G**, and 6-31G(3df)S** as small, medium, and slightly larger basis sets, respectively, give good geometries with respect to the size of calculation. The geometries for these calculations are summarized in Table 1 with results for the others in Table S1-S5. The B3LYP method using 6-31(++)SG**, DZVP2, and DZVP2(++)S also gave good geometries. However, the addition of the diffuse functions for both 6-31G** and DZVP2 basis sets increased the size of the calculation without significantly changing the results. Also, although the 6-31G** and DZVP2 basis sets behaved similarly and are of similar size, DZVP2 basis sets gave slightly poorer geometries for [Fe(SCH3)4]2- and the 6-31G** basis sets are widely used and have been extended to include sp-type diffuse functions (6-31++G**) and additional polarization functions [6-31G(3df,3pd)].
To assess functional and basis set effects on energetics, the ADE and VDE of [Fe(SCH3)4]1- and [Fe(SCH3)3]1- were calculated using energies without ZPE corrections obtained at the B3LYP and B97gga1 levels with the ECPDZ+DZVP, 6-31G**, 6-31(++)SG**, and 6-31G(3df)S** basis sets (Table S6 and S8), using geometries calculated at various levels. The energies were compared to PES results for ADE and VDE of [Fe(SCH3)4]1-/0 and [Fe(SCH3)3]1-/0,16,17 and also to energies calculated at the MPn and coupled cluster levels.
First, the general performance of the B3LYP and B97gga1 functionals for calculating ADE and VDE against experiment and the MP2, MP3, MP4SDQ, CCSD, and CCSD(T) levels of theory were evaluated for the B3LYP/6-31G** and B97gga1/6-31G** geometries (Figure 3, Table S6). The energies were evaluated using the 6-31G** basis so that a consistent basis set could be used even for the very large coupled cluster calculations. The energies were also calculated for the B3LYP/ECPDZ+DZVP geometries (Table S6), but the results showed similar trends as the B3LYP/6-31G** geometries. Although the B97gga1 optimized geometries were in good agreement with experimental values, the ADE and VDE calculated from B97gga1 energies were underestimated by ~0.7 to ~1 eV compared to experiment, while they were generally much better at the B3LYP level where the energies were underestimated by up to ~0.3 eV compared to experiment. Since B97gga1 energy calculations using either B3LYP/6-31G** or B97gga1/6-31G** geometries gave very similar results, this is not due to the poor changes in Fe—S bond length upon reduction in B97gga1 compared to the B3LYP but rather a poorer treatment of the energy. Examining the other methods, the MP series failed as an appropriate method to include correlation energy especially for [Fe(SCH3)3]1- while the coupled cluster methods generally gave better predictions; however, the CCSD(T) ADE and VDE of [Fe(SCH3)3]1- were overestimated by about 0.4 eV apparently due to degeneracy or near-degeneracy problems. Overall, B3LYP/6-31G** gave much better oxidation energies than B97gga1, which equaled the accuracy of the energy calculations at the very costly CCSD level.
Next, the performance of different basis sets for calculating ADE and VDE against experiment were evaluated using the B3LYP functional with both the B3LYP/ECPDZ+DZVP and B3LYP/6-31G** geometries (Figure 4, Table S8). Results for B3LYP/DZVP2 and B3LYP/DZVP2(++)S calculations of geometry and energy also shown for comparison. The performance was also evaluated for the B97gga1 functional with the B97gga1/ECPDZ+DZVP and B97gga1/6-31G** geometries (Table S8); however, the overall energies were underestimated by up to ~1 eV as indicated above and the trends found for B3LYP with increasing basis set size were the same. The calculated ADE and VDE from B3LYP tends to decrease as the basis set size increases from the ECPDZ+DZVP to full core basis sets since larger basis sets can more accurately approximate the orbitals by imposing fewer restrictions on the spatial distribution of the electrons and by describing the defuse nature of anions.44,45 However, relative to the 6-31G** basis sets, adding diffuse functions to the sulfurs increases and thus improves the calculated ADE and VDE with respect to experiment, while adding more polarization functions to the sulfurs decreases and thus worsens these values. In particular, while addition of diffuse functions or additional polarization functions to the sulfurs lowers the energy for both the oxidized and reduced analogs by ~0.1 to 0.2 eV and ~1.2 to 1.5 eV, respectively, the diffuse functions lower the energy of the anionic reduced analog by ~0.1 eV more than that the neutral oxidized analog, while the additional polarization functions lower the energy of the neutral analog by ~0.1 eV more than the anionic analog (Table S9). Other PES and DFT studies on the tetrahedral ferric complexes FeIIIX4- (X = Cl, Br) and the three-coordinate complexes MIIX3- (M = Mn, Fe, Co, Ni; X = Cl, Br) also indicate that increasing the basis set size for all atoms, for example using triple-ζ basis sets, does not significantly change the calculated geometric parameters and redox energies.46
Overall, the B3LYP method using ECPDZ+DZVP, 6-31G**, 6-31(++)SG**, as small, medium, and slightly larger basis sets, respectively, give good ADE and VDE using the B3LYP ECPDZ+DZVP and 6-31G** geometries. The ADE and VDE for these calculations are summarized in Table 2 with results for the other calculations in Table S6 and S8. The B3LYP method using the DZVP2, and DZVP2(++)S basis sets also gave good ADE and VDE, but as mentioned previously, the 6-31G** are widely used and have been extended to include sp-type diffuse functions and additional polarization functions. However, energies using 6-31G(3df)S** are in even poorer agreement than energies using ECPDZ+DZVP at a much greater computational cost so adding diffuse functions is preferred to adding more polarization functions. Also, the energy calculations at the B3LYP level are not very sensitive to the geometries tested here so that reasonable reduction potential energies can be obtained based on optimized geometries with smaller basis sets. In particular, the B3LYP/6-31(++)SG** energies for the 6-31G** and 6-31(++)SG** geometries (i.e., B3LYP/6-31(++)SG**//B3LYP/6-31G** versus B3LYP/6-31(++)SG** in Table S8) are within less than 0.01 eV so that optimization of the geometry can be done using the 6-31G** basis sets.
To find efficient computational methods for structural and energetic properties of iron-sulfur protein redox sites, spin-unrestricted B3LYP, BHandH, B97gga1, MP theory series, CCSD, and CCSD(T) calculations using a variety of ECP and full core basis sets have been performed on analogs of the rubredoxin redox site and compared to the experimental geometries, ADE, and VDE. Overall, the B3LYP method gives the best results in comparison with experiment and are comparable for energies with those at the very costly CCSD or CCSD(T) level of theory. For geometries, B3LYP using ECPDZ+DZVP basis sets overestimated the Fe-S bonds by ~0.06 to ~0.08 Å, using 6-31G** by ~0.03 to ~0.06 Å, and using 6-31G(3df)S** by ~0.02 to ~0.04 Å; but all gave changes in Fe—S bond length upon oxidation within ~0.01 Å of experiment and much more accurate S—C bonds. The energies are relatively insensitive to geometries, giving similar results for the B3LYP/ECPDZ+DZVP and B3LYP/6-31G** geometries. For energy calculations, the B3LYP functional using ECPDZ+DZVP basis sets overestimate by no more than ~0.2 eV of experiment, using 6-31G** underestimate by no more than ~0.3 eV, and using 6-31(++)SG** underestimate by no more than ~0.15 eV. Adding zero point energy corrections improves the 6-31(++)SG** energies so that the ADE are underestimated by no more than 0.05 eV although the VDE for [Fe(SCH3)3]1- is underestimated by 0.18 eV.
In summary, a central issue in the reduction potential calculations of iron-sulfur systems is to obtain relatively accurate energies for both the oxidized and the reduced states. The recommended calculation for iron-sulfur systems is to optimize the geometry using B3LYP with at least a double zeta basis set with polarization functions such as 6-31G** and DZVP2 and then to calculate energies using B3LYP with additional diffuse functions such as from 6-31++G** added to the sulfurs [6-31(++)SG** and DZVP2(++)S]. However, for larger iron-sulfur systems in QM/MM calculations, a reasonable calculation would be to optimize the geometry using B3LYP with an ECP for iron and then to calculate energies calculated using B3LYP/6-31(++)SG** or DZVP2(++)S.
This work was supported by a grant from the National Institutes of Health (GM45303) and by a grant (GC3565 and GC20901) from the Molecular Science Computing Facility (MSCF) in the William R. Wiley Environmental Molecular Sciences Laboratory (EMSL), a national user facility sponsored by the U.S. DOE’s Office of Biological and Environmental Research and located at Pacific Northwest National Laboratory, operated for DOE by Battelle. We thank Dr. David Dixon, Prof. Lai-Sheng Wang, and Dr. Xubin Wang for helpful discussions. TI thanks Dr. Bernard R. Brooks at NIH for his hospitality during part of these studies.