|Home | About | Journals | Submit | Contact Us | Français|
The EPR spectral parameters, i.e. g-tensors and molybdenum hyperfine couplings, for several d1 systems of the general formula [MoVEX4]n−, [MoVOX5]2−, and [MoVOX4(H2O)]− (E = O, N; X = F, Cl, Br; n = 1 or 2) were calculated using Density Functional Theory. The influence of basis sets, their contraction scheme, the type of exchange-correlation functional, the amount of Hartree-Fock exchange, molecular geometry, and relativistic effects on the calculated EPR spectra parameters have been discussed. The g-tensors and molybdenum hyperfine coupling parameters were calculated using a relativistic Hamiltonian coupled with several LDA, GGA, and ‘hybrid’ exchange-correlation functionals and uncontracted full-electron DGauss DZVP basis sets. The calculated EPR parameters are found to be sensitive to the Mo=E distance and E=Mo–Cl angle, and thus the choice of starting molecular geometry should be considered as an important factor in predicting the g-tensors and hyperfine coupling constants in oxo-molybdenum compounds. In the present case, the GGA exchange-correlation functionals provide a better agreement between the theory and the experiment.
Mononuclear molybdoenzymes are involved in the global cycling of nitrogen, sulfur, and arsenic.1 During catalysis, the Mo-center cycles through MoVI, MoV, and MoIV oxidation states and participates in a variety of oxygen-atom transfer reactions. The nature of the active center in MoIV or MoVI oxidation states of various enzymes has been characterized by X-ray crystallography.2 However, information about the transient MoV state relies on spectroscopic techniques such as electron paramagnetic resonance (EPR), UV-visible spectroscopy, and magnetic circular dichroism (MCD). Specifically, EPR spectroscopy has proven to be a very valuable tool in investigating metal-ligand interactions and details concerning the first coordination sphere in the MoV state.3 In many cases, several forms of the same enzyme, such as sulfite oxidase and xanthine oxidase, have been defined based on the EPR parameters. These forms have provided mechanistic details about the functioning of the enzymes. Experimental EPR spectral parameters, such as g-tensors, hyperfine coupling constants (HFC), and their orientations in Euler space, can be extracted from the raw data using modern EPR simulation approaches. These parameters, in most cases, can be directly correlated with the metal–ligand interactions in molybdoenzymes. 1 In some cases, however, the relationship between the observed g-tensors and HFC and the enzymatic MoV center is not obvious and must be confirmed by different experimental methods or a reliable theoretical approach.4,5
An early theoretical approach for predicting the EPR spectral parameters in molybdenum systems was introduced by Sunil et al.6 and this approach was later modified by Westmoreland and coworkers.7 The method is based on ligand-field theory and requires an analysis of the ground-state wavefunction coupled with the determination of vertical excitation energies in MoV complexes using the Slater method along with empirically evaluated spin-orbit coupling constants. In spite of good agreements between the predicted and experimentally observed EPR parameters for small, high-symmetry complexes, this method is difficult to use in the modeling of EPR parameters in low-symmetry MoV active site model complexes and in inherently low-symmetry molybdoenzymes. Later, semi-empirical INDO/S and INDO/S-CI-Stone methods coupled with a second-order perturbation theory methodology were used with some success in predicting the g-tensors in several simple MoV complexes.8 In 1999, Patchkovskii and Ziegler published a very detailed study on DFT predicted g-tensors in [MEX4]n− d1 transition-metal compounds, including several MoV complexes.9 This report also discussed the influence of molecular geometry as well as local density (LDA) and gradient-corrected approximation (GGA) exchange-correlation (XC) functionals on the calculated g-values and provided a basis for applying modern DFT methods for the accurate prediction of EPR g-tensors in 4d transition-metal complexes. A similar approach has been reported by Neese et al. who used GGA (BP86) and hybrid (B3LYP) XC functionals for evaluating the EPR spectral parameters of several MoV complexes coordinated by sulfur donors,10 with both XC functionals providing good agreements between theory and experiment. Similarly, in 2001, Neese computed g-values in transition-metal compounds using coupled perturbed Hartree-Fock (HF) and Kohn-Sham theories.11 This approach was further modified by Kaupp’s group using sophisticated, but more computationally demanding, treatments of relativistic effects and a more accurate calculation of spin-orbit constants.12–14 The latter group also suggested that for accurate predictions of EPR spectral parameters in transition-metal complexes, the XC functionals should have 30–40% Hartree-Fock exchange.12–14 However, this suggestion was disputed by Neese.10 While the approach suggested by Kaupp provides a more sophisticated means of calculating g-values and HFCs and includes higher-order spin-orbit contributions, it is also more computationally demanding. In contrast, at least for computationally simple d1 systems (e.g. MoV) where no electron–electron repulsion need to be considered, the method outlined by Neese11 can provide a fast and reliable outcome.
In this manuscript, a coupled perturbed method similar to that proposed by Neese is used for testing the effects of a variety of parameters on calculated g-values and HFCs in MoV d1 complexes. The aim of this article is to systematically compute and compare EPR spectral parameters, specifically g-tensors and HFCs, in discrete MoV complexes. The following aspects will be discussed: (i) can the current “standard” XC functionals be used for accurate prediction of both the g-values and the HFC parameters in MoV complexes; (ii) can a relatively small full-electron basis set (DZVP) be used for the accurate prediction of EPR spectral parameters in MoV complexes; and (iii) what is the influence of molecular geometry on the calculated values of g-tensors and HFCs. A clear understanding of these factors should allow a relatively fast, yet accurate, prediction of EPR parameters in representative MoV d1 complexes.15
The following small, well-known MoV d1 systems were tested: [MoVOF4]−, [MoVOCl4]−, [MoVOBr4]−, [MoVOF5]2−, [MoVOCl5]2−, [MoVOBr5]2−, [MoVNCl4]2−, [MoVOCl4(H2O)]−, and [MoVOBr4(H2O)]− (Fig. 1). These systems were used because of their well-described EPR properties and molecular structures.
Four sets of molecular geometries were used in the present investigation. The first one represents crystallographically determined geometries obtained either from the CCDC database16 or original publications.17 In the second set, molecular geometries were optimized using the 1997 hybrid functional of Perdew, Burke and Ernzerhof.18 This functional, referred to as PBE1PBE, uses ~25% of Hartree-Fock exchange. In the third set, molecular geometries were optimized using Becke’s three-parameter hybrid exchange functional19 (~20% of Hartree-Fock exchange) and a Lee-Yang-Parr non-local correlation functional20 (B3LYP). Finally, in the fourth set, molecular geometries were optimized using Becke’s exchange functional21 and Perdew’s non-local correlation functional22 (BP86, 0% of Hartree-Fock exchange). In all cases, full-electron DGauss DZVP21 and 6-311G(d)22 basis sets were used for molybdenum and all other atoms, respectively. Such a combination of basis sets has proved to be reliable in the prediction of molecular geometries in numerous molybdenum complexes.23 For all optimized structures, frequency calculations were carried out to ensure the optimized geometries represented minima on their respective potential energy surfaces. Optimized geometries of [MoVEX4]n− and [MoVOX5]2− (where X = F, Cl or Br) have the expected C4v symmetries. The global minima for [MoVOX4(H2O)]− compounds can only be achieved within C1 symmetries, while one or two imaginary frequencies were always observed in the cases of all possible higher (Cs and C2v) molecular symmetries.
In the present study, the relatively small (18s12p9d) DGauss full-electron DZVP basis set has been tailored for the accurate prediction of EPR g-tensors and A values by substituting the original p and d basis functions in the DZVP basis set with either Stuttgart/Dresden effective core potentials or basis functions from Ahlrich’s basis set. These modified basis sets are presented in the ESI and Table S1.†
The following XC functionals were used for the calculation of EPR parameters: the local density approximation (LDA) Xα,24 SVWN,25 SVWN5,25 and HFS;26 the generalized gradient approximation (GGA) BP86,21,22 BPW91,19,27 PW91PW91,28 HFB;31 hybrid functionals with ~20% Hartree-Fock exchange, B3P86,19,22 B3LYP,19,20 B3PW91,19,27 B98;29 hybrid functionals with ~25% Hartree-Fock exchange PBE1PBE,30, MPW1PW91;28 and hybrid functionals with 50% Hartree-Fock exchange, BHandH,31 and BHandHLYP.31 In addition, user-defined BLYP-based hybrid functionals of the general formula:
were used, where a is the amount of Hartree-Fock exchange constructed in the Gaussian 03 program, with values of 0, 10, 20, 30, 40, 50, 60, and 70% of exchange, and used to investigate the influence of Hartree-Fock exchange.
Initial calculations on selected MoV d1 systems suggest that relativistic effects should be considered. Therefore relativistic effects, using a relativistic elimination of small components (RESC) approach, were included where possible.32 In addition, in selected cases, first- and second-order Douglas-Kroll-Hess (DKH)33 approaches were also tested.
Because the Fermi contact term is sensitive to the quality of the numerical integration, a relatively large integral grid with 128 radial shells and 770 angular points per shell were utilized in all calculations.
All EPR calculations were performed using the Gaussian 98 or Gaussian 03 program family running on either a Windows or UNIX operating system34 using methodology implemented into the Gaussian 03 software (see ESI for details†). Mulliken35 charges for all atoms of interest were calculated using the standard procedures that are implemented into the Gaussian software package as reported elsewhere.36 When necessary, the percentage contributions of atomic orbitals to molecular orbitals were calculated using the VMOdes program.37 In all cases, a tight energy (10−8 au) SCF convergence criterion was used.
To date, EPR parameters of molybdenum(V) d1 systems have been calculated primarily using fully uncontracted large triple-ζ quality Slater-9 or Gaussian-type5,10,12 basis sets to facilitate basis set flexibility. The calculated Fermi contact term, AF, depends on the s-electron density present at the Mo center. Thus, it is sensitive to an accurate description of the s-part of the basis set.38 In contrast, the anisotropic component of the hyperfine coupling constant requires an accurate description of the molybdenum d orbitals. Thus, it is imperative to have a basis set that accurately reflects both s and d parts.
The medium-size (18s12p9d) DGauss full-electron DZVP basis set has been tailored for [MoOCl4]−, [MoNCl4]2−, [MoOF5]−, [MoOCl5]2− and [MoOBr5]2− complexes using a hybrid B3P86 XC functional. The calculations utilized crystallographically determined geometries (ESI Tables S2–S5†). The original DGauss contraction scheme for this basis set allows for the accurate prediction of g-tensors and anisotropic HFCs, while the calculated Fermi contact term was found to be much lower than the experimental values. As expected, the Fermi contact term was calculated more accurately with an uncontracted s-part. A small improvement (~8%) on the calculated AF term was achieved by the uncontraction of the d-part of the basis set as well as the addition of polarization f-functions (BS 2). On the other hand, stepwise uncontraction of the p-part of the basis set and addition of several s-functions with larger exponents showed39 little effect on the calculated values of AF. Thus, it seems that uncontraction of the s- and d-parts of the medium size DZVP basis set results in an acceptable agreement between theoretical and experimental g-values and anisotropic A-tensors in MoV d1 systems. Taking this into consideration, the following basis sets were tested in detail: (i) the original DZVP basis set with completely uncontracted s- and d-parts (referred as BS1); (ii) the more flexible in the d-part BS2, described in the basis set section.
In order to investigate the influence of different relativistic scalar effects on the calculated EPR spectral parameters, first- as well as second-order Douglass-Kroll-Hess (DKH) along with RESC calculations were conducted on complexes with C4v symmetry (i.e. [MoOF4]−, [MoOCl4]−, [MoOBr4]−, [MoOF5]2−, [MoOCl5]2−, [MoOBr5]2−, and [MoNCl4]2−) using either X-ray derived or optimized geometries, and a B3P86 exchange correlation functional (ESI Table S6†). Overall, the calculated g-tensors and anisotropic contact terms, calculated using either of the first- or second-order DKH as well as RESC levels of theory, are very close to those obtained from non-relativistic calculations. Because of this, we used the RESC level of theory for the remaining calculations presented in this study. As expected, the only significant difference is in Fermi contact terms calculated with or without relativistic corrections. Specifically, in the case of all non-relativistic calculations, the calculated values of the AF term were significantly underestimated, similar to previous non-relativistic calculations on MoV systems.7 On the other hand, when BS1 and BS2 were used on Mo, with relativistic corrections, the results were in good agreement with experimental values of AF (ESI Table S6†).
The deviation from the free electron g-value in the case of Mo-enzymes and model complexes follows a general trend of ge > g||(gz) > g(gx,gy), and this trend has been attributed predominantly to metal–ligand covalency and large values of ligand spin-orbit coupling, which can be either directly or indirectly geometry dependent. Thus, an accurate description of geometry is of significant importance. Although it is expected that the crystal packing forces would affect the geometries of the complexes of interest as compared to those in solution, when both EPR single crystal and solution data are available (i.e. [MoOCl4]− and [MoOCl4(H2O)]−), only very minor differences were observed. The influence of molecular geometry on the DFT predicted EPR parameters is exemplified in Fig. 2, which shows the dependence of the calculated g-tensors on Mo=O bond distance and the O–Mo–Cl bond angle. Similarly, Fig. 3 shows the dependence with regards to the calculated AMo HFCs. Clearly, the magnitude of DFT predicted EPR parameters vary greatly with subtle changes in bond lengths and angles. Previous studies on small model MoV complexes included geometries optimized with BP86 and SVWN XC functionals. These resulted in Mo=O bond lengths in [MoOCl4]− of 1.692 and 1.706 Å, respectively, and 1.690 and 1.703 Å, respectively, in [MoOBr4]−.9 In the present work, three different XC functionals (BP86, B3P86 and PBE1PBE) were used for geometry optimization. These functionals contain 0, 20, and 25% Hartree-Fock exchange, respectively. The bond distances and angles, which are presented in Table 1, are clearly suggestive of an XC functional geometry dependence. The bond distances in the structures optimized with pure GGA XC functionals are longer than the majority of available X-ray crystal structures due to the overall increase in the covalency of metal–ligand bonds. Fig. S1 in the ESI† presents a histogram of Mo=O bond distances for reported crystal structures available in the CCDC. These high-quality structures have Mo=O groups in which Mo is bonded to at least two halogens. The majority of these structures have a Mo=O bond distance of approximately 1.64 to 1.69 Å, which can be directly compared to those obtained from the geometry optimization. This histogram clearly indicates that Mo=O bond distances are slightly overestimated by GGA XC functionals as compared to those obtained using hybrid XC functionals. Taking into consideration the clear geometry dependence of the predicted EPR spectral parameters, it is important to compare calculated EPR g-tensors and HFCs for different optimized, as well as X-ray determined, geometries.
The accuracy of calculated spectroscopic parameters (e.g., Mössbauer, UV-vis, EPR parameters)39,40 depends on the type of XC functional, and so far no “universal” XC functional has been found for predicting different spectroscopic properties of transition metal complexes. Thus, a clear understanding of the type of DFT method, whether it is based on pure (i.e. LDA and GGA), hybrid, or “half-and-half” hybrid XC functionals, that consistently model experimentally observed EPR parameters of transition metal complexes (e.g., oxomolybdenum compounds) is very important. Results of correlation analyses for the calculated g-tensors, calculated Fermi contact term, and the largest anisotropic contribution (AMo3) to the A33 tensor for 16 different exchange correlation functionals, two different basis sets and four different geometries are summarized in Tables S7–S14 in the ESI† and graphically presented in Fig. 4 and Fig. S2–S9 in the ESI.† The results clearly display a delicate interplay between the amount of Hartree-Fock exchange used in the exchange correlation functional and the accuracy of the computed EPR parameters. In order to evaluate the quality of the results, MAD criterion (Δg||<0.03, Δg<0.03, ΔAF < 10 cm−1, ΔAMo3 < 3 cm−1, ESI Tables S7–S14†) and ‘border conditions’ (absolute deviations presented as dotted lines in Fig. 4 and ESI Fig. S2–S9†) were used.
In general, g|| in test systems is more accurately predicted by GGA-based XC functionals (i.e., BP86, BPW91, HFB, and PW91PW91) when all geometries and both BS1 and BS2 were used. An exception was the computed value of g|| for the [MoOBr4(H2O)]− complex, which was significantly over-estimated. It should be noted, however, that the reported value of g|| in this complex (1.98) is far below that expected for [MoOBr4X]n− complexes (X=ligand trans to the Mo=O bond), for which g|| was suggested to be ~2.1 based on the small influence of the trans-ligand X. The calculated g values for all four geometries and both basis sets are slightly under-estimated, with the largest error observed for the [MoOF5]2− complex. Again, GGA-based XC functionals (i.e., BP86, BPW91, HFB, and PW91PW91) provide a better agreement between theory and experiment. The DFT predicted values of AF are the most dependent on the amount of Hartree-Fock exchange present in the XC functional as well as the type of XC functional. Indeed, LDA-based SVWN and SVWN5 XC functionals slightly underestimate the AF term, while other LDA (Xα and HFS) as well as GGA- and hybrid (~20% of Hartree-Fock exchange) XC functionals (i.e. BP86, BPW91, HFB, PW91PW91, B3P86, B3LYP, and B3PW91) provide a good agreement between theory and experiment. A further increase in Hartree-Fock exchange leads to the overestimation of AF values in the test systems. Finally, the anisotropic components of HFCs are less dependent on the XC functional, as shown in Fig. 3 and Fig. S8 and S9 in the ESI,† with the more flexible BS2 providing a slightly better agreement between theory and experiment. Overall, BPW91, PW91PW91, and HFB XC functionals provide the best agreement between theory and experiment for all geometries when BS1 is used. In addition, hybrid B3P86, B3PW91, and B3LYP XC functionals can be used with X-ray geometries. In the case when BS2 was used, GGA-based BP86, BPW91, PW91PW91, and HFB XC functionals again can be used for the calculation of the spectral parameters of MoV d1 systems.
An admixture of 30–40% HF exchange to the exchange correlation functional has been suggested to be a prerequisite in obtaining a satisfactory agreement between predicted and experimental EPR spectral parameters. To test whether this is applicable in the present case, we adjusted the Hartree-Fock exchange in a BLYP-xx XC functional in a stepwise manner. The results of these calculations for the [MoOCl4]− complex are presented in Fig. 5. Indeed, the calculated g-factors and A-tensors change with increasing amounts of Hartree-Fock exchange, with 10–20% of Hartree-Fock exchange providing the best agreement between theory and experiment. Of course, such a Hartree-Fock dependency varies with computational method (i.e., the inclusion of higher-order spin-orbit contributions12–14 can lead to the requirement of incorporating 30–40% of Hartree-Fock exchange in calculations).
The g-tensors and molybdenum hyperfine coupling constants for a set of d1 systems of the general formulae [MoVEX4]n−, [MoVOX5]2−, and [MoVOX4(H2O)]− (E = O, N; X = F, Cl, Br; n = 1 or 2) were calculated using Density Functional Theory. The influence of the basis set, basis set contraction scheme, type of XC functional, amount of Hartree-Fock exchange, molecular geometry, and relativistic effects on the calculated EPR spectral parameters have been discussed in detail. The EPR g-tensors and molybdenum hyperfine coupling parameters calculated using a relativistic Hamiltonian coupled with several GGA and hybrid XC functionals and specifically tailored medium-size DZVP basis sets were found to be in excellent agreement with the experimental data. The calculated EPR parameters were found to be very sensitive to the Mo=E distance and E=Mo–Cl angle. Taken together, the accurate prediction of the EPR parameters of MoV compounds reflects a complex interplay between molecular geometry, XC functional, basis set, and relativistic effects. An important finding of this investigation is that the geometry of the system should be defined first which ultimately controls the parameters. A small change in the bond distance and angle changes the orbital interaction, leading to a change in the EPR parameters. Overall, BPW91, PW91PW91, and HFB XC functionals provide the best agreement between theory and experiment. The gradient corrected methods provide a better agreement than the local density approximation. In addition, the inclusion of HF exchange correlation has a negative impact on the results. While we have used only oxo-Mo(V) centers, which are supposedly simpler due to a lack of inter-electron repulsion, as representatives of d1 systems, we anticipate the same may hold true for other d1 systems.
We thank the National Science Foundation (Grant CHE-0809203) as well as the National Institute of Health (GM 650155 to PB) for financial support and the Minnesota Supercomputing Institute for the generous support of computer time, as well as Undergraduate Research opportunity grants to RH and JO.
†Electronic supplementary information (ESI) available: Calculation of EPR spectral parameters; basis sets tested for the prediction of EPR parameters.