Golokesh
Santra
,
Frank
Neese
* and
Dimitrios A.
Pantazis
*
Max-Planck-Institut für Kohlenforschung, Kaiser-Wilhelm-Platz 1, 45470 Mülheim an der Ruhr, Germany. E-mail: dimitrios.pantazis@kofo.mpg.de
First published on 23rd August 2024
Mössbauer spectroscopy is a powerful technique for probing the local electronic structure of iron compounds, because it reports in an element-selective manner on both the oxidation state and coordination environment of the Fe ion. Computational prediction of the two main Mössbauer parameters, isomer shift (δ) and quadrupole splitting (ΔEQ), has long been targeted by quantum chemical studies, and useful protocols based on density functional theory have been proposed. Here we present an extensive curated reference set of Fe compounds that is considerably larger and more diverse than literature precedents. We make a distinction between low-temperature and high-temperature experimental subgroups. This set is employed for optimizing a refined computational protocol utilizing the scalar version of the exact 2-component (X2C) Hamiltonian with the finite nucleus approximation. Attention is devoted to having an accurate and flexible all-electron basis set for Fe. We assess the performance of several DFT methods that cover all representative families and rungs of functionals and find that hybrid functionals with ca. 25–30% exact exchange offer the best accuracy for isomer shifts. The work establishes a refined general protocol of wide applicability that achieves good performance for the prediction of isomer shifts in a wider variety of systems than before, but the limitations of DFT for quadrupole splittings are also highlighted. Finally, comparison of calculated values with high-temperature experimental results shows that the use of an empirical correction factor is required to account for the second-order Doppler shift and to achieve the same quality of correlation as with the low-temperature data.
The two most important parameters encoded in a Mössbauer spectrum are the isomer shift (δ) and the quadrupole splitting (ΔEQ). The isomer shift of a metal center is directly related to the electron density at the nucleus, while ΔEQ is proportional to the electric field gradient (EFG), which originates from a nonsymmetrical distribution of electrons in the valence shell as well as charges on the neighboring ligands.
The isomer shift of 57Fe measures the shift in the energy of the γ-ray absorption relative to a standard, usually iron foil. The isomer shift is sensitive to the electron density at the nucleus, and indirectly probes changes in iron–ligand bond lengths, covalency and nature of its bonds, and shielding due to the 3d orbital occupation pattern. As a result, it can successfully probe oxidation and spin states, and the coordination environment of Fe. This shift between an absorber (A) and a source (S) comes from the difference in the electrostatic interactions between electronic and nuclear charge distributions, which originates from the difference in their electron densities as well as the change in the nuclear radius upon gamma transition. Considering the nucleus to be a uniformly charged sphere, the mathematical expression for the Mössbauer isomer shift is
δ = α(![]() ![]() | (1) |
![]() | (2) |
δ = a[ρ(0) − C] + b | (3) |
Over the years, linear regression analysis has been extensively applied for the calculation of isomer shifts using semi-empirical, Hartree–Fock (HF), density functional theory (DFT), and wave function based ab initio methods like the domain based local pair-natural orbital coupled-cluster theory (DLPNO-CCSD).11,13–34 All these studies have demonstrated good correlation between theory and experiment. Despite its simplicity, this approach is fairly reliable and efficient and is known to predict isomer shifts with an accuracy of up to ∼0.1 mm s−1.27 Accurate computation of the contact density is challenging and depends on factors like the choice of the quantum chemical method, the basis set, and the proper treatment of relativistic effects. An ideal basis set which can adequately describe the contact density of Fe has to be sufficiently large in the region where a cusp in the electron density will occur. The known nonrelativistic HF limit for ρ(0) is ∼11903.987 a.u.−3, whereas with a good basis set one can only obtain up to 11
820 a.u.−3.16,19 Although in terms of absolute contact density this error is negligible, compared to the variation of the electron density over the chemical range (∼10 a.u.−3) it is nontrivial.
The relativistic effects on the electron density at the nucleus are large for iron, increasing ρ(0) by a factor of 1.3 compared to the nonrelativistic electron density of 57Fe.17,19 On the other hand, Saue and co-workers have shown that for the calculation of Mössbauer isomer shifts spin–orbit coupling can be safely ignored while using the eXact 2-Component (X2C) Hamiltonian35–38 and a finite nucleus model,39 therefore the consideration of scalar relativistic effects is sufficient.40 As already demonstrated,11 the contributions from the core 1s and 2s orbitals remain nearly constant to the absolute contact density at the iron nucleus. The major contribution to the variation of the contact density due to its electronic configuration and ligand environment arise from the valence and subvalence regions. Hence, these inner valence and outer valence orbitals of iron matter the most for obtaining a correct calibration.11,19
The other important parameter in a Mössbauer spectrum is the quadrupole splitting (ΔEQ), which arises from the interaction of the nuclear quadrupole moment of the excited state with the EFG at the nucleus. Quadrupole splittings can also be used as a sensitive probe for the coordination environment of iron centers. The quadrupole splitting is obtained from the EFGs using the expression
![]() | (4) |
Another important factor is the absence of an accurate value of the 57Fe nuclear quadrupole moment (NQM). As the NQM is impossible to determine experimentally, the only possible way to obtain those values is via linear regression analysis using the experimental ΔEQ values and the theoretical EFGs. Therefore, a wide range of values between 0.1 to 0.3b can be found in the literature. In the present study, we have used Q = 0.160b for the calculation of quadrupole splittings.17
Over the years, several benchmark studies on 57Fe Mössbauer parameters have been performed by employing scalar relativistic Hamiltonians like ZORA (zero-order regular approximation),41 DKH2 (second-order Douglas–Kroll–Hess),42–48 and X2C (exact two-component).18,19,21,22,40,49 Although changing the Hamiltonian had little effect on the overall correlation of experimental δ and calculated contact densities of iron,18,22 a detailed evaluation of a variety of density functional methods and basis sets with the scalar relativistic X2C Hamiltonian against a large and diverse dataset is still absent. Another potential constraint of previous studies is related to the set of iron complexes considered for benchmarking quantum-chemistry methods. These datasets are either composed of limited number of compounds, a mixture of molecular and solid-state systems, a specific type of ligand, or only a limited number of spin and oxidation states of iron. The main objectives of the present study are:
(a) To construct a complete database of Fe complexes with well-established Mössbauer parameters, which is curated so that it includes most of the known spin and oxidation state of iron and is representative of the wide range of chemical types encountered in iron coordination chemistry;
(b) To conduct a thorough benchmark study of basis sets and DFT functionals against our new database in order to develop a refined computational protocol based on the popular X2C Hamiltonian.
In the present study the performance of 10 different basis sets for iron is evaluated: CP(PPP),11 x2c-TZVPPall,56 x2c-TZVPPall-s,62 aug-cc-pVTZ-J (or aVTZ-J),63 DKH-def2-TZVPP (exponents from def2-TZVPP64 were recontracted for scalar-relativistic DKH243–47,65 Hamiltonian), ANO-RCC-VTZP,66 s-decontracted x2c-TZVPPall, s-decontracted aug-cc-pVTZ-J, s-decontracted aug-cc-pVTZ-J(-dfg), and aug-cc-pVTZ-Jmod.67 The last basis set was proposed by Gómez-Piñeiro et al. for the calculations of Cu(II) core properties, where the aug-cc-pVTZ-J was modified by decontracting the s functions and removing the three innermost primitives.67 On the other hand, the outermost d-, f- and g-primitives are removed from s-decontracted aug-cc-pVTZ-J to obtain the s-decontracted aug-cc-pVTZ-J(-dfg) basis set. For the ligand atoms, x2c-TZVPPall56 was used throughout. 15 different density functionals from all five rungs of Jacob's ladder68 were calibrated: SVWN5,69,70 BP86,71,72 PBE,73 BLYP,72,74 TPSS,53 PBE0,75 B1LYP,76 B3LYP,74,77,78 TPSSh,52,53 TPSS0,79 M06,80 LC-BLYP,81 CAM-B3LYP,82 ωB97X,83 and B2PLYP.84 For the PT2 part of the double hybrid functionals, we correlated all core electrons and employed both relaxed and unrelaxed densities. The ORCA sample input files for the calculation of Mössbauer properties and the modified basis sets are provided in the ESI.† In a previous study, the accuracy of selected DFT functionals was evaluated on a small set of 20 iron-containing compounds.18 The structures were optimized there using the TPSS functional and def2-TZVP85 basis set. Following the same protocol for purposes of comparison, the complexes of MPMIC80 were also reoptimized and the nonrelativistic Mössbauer parameters were calculated by employing the B3LYP functional, CP(PPP)11 basis set for iron, and def2-TZVP64 for the ligand atoms.
These results are worse than literature expectations,18,22 owing to the considerably expanded reference set of compounds in the present study. The mean absolute error and standard deviation of calculated isomer shift (δcal) with respect to δexp are 0.07 and 0.09 mm s−1, respectively. Closer inspection reveals that the calculated isomer shifts of all nine Fe(ii) complexes with quintet spin multiplicity (2S + 1 = 5) deviate significantly from the experimental data. This stresses the necessity of reconsidering the linear correlation parameters obtained from more restricted reference sets.
Basis set for Fe | a | b | C | R 2 | MADc (mm s−1) | Max. dev.d (mm s−1) | St. dev.e (mm s−1) |
---|---|---|---|---|---|---|---|
a With first-order picture change effect and finite nucleus model. For the results obtained with the point nucleus model and picture change effect, see Table S4 in the ESI. b a and b are the fitting coefficients obtained from the linear fit of eqn (3). C is a constant, which is very close to the calculated ρ(0) value. The units of a, b, and C are a.u.3 mm s−1, mm s−1, and a.u.−3, respectively. R2 is the coefficient of determination from the linear fit. c Mean absolute deviations of the calculated ISs with respect to the experimental ISs, where the former ones are obtained by using eqn (3). d Maximum deviation of the calculated ISs from the experimental ISs. e Standard deviation of the calculated ISs. | |||||||
CP(PPP) | −0.30 | 0.9704 | 14![]() |
0.956 | 0.054 | 0.191 | 0.070 |
x2c-TZVPPall | −0.22 | 1.0294 | 13![]() |
0.820 | 0.108 | 0.346 | 0.142 |
x2c-TZVPPall-s | −0.22 | 0.8917 | 13![]() |
0.806 | 0.113 | 0.356 | 0.147 |
aug-cc-pVTZ-J | −0.41 | 1.3456 | 10![]() |
0.956 | 0.054 | 0.181 | 0.070 |
ANO-RCC-VTZP | −0.18 | 0.8903 | 15![]() |
0.414 | 0.184 | 0.752 | 0.255 |
DKH-def2-TZVPP | −0.30 | 1.0315 | 14![]() |
0.956 | 0.055 | 0.190 | 0.070 |
s-decontracted x2c-TZVPPall | −0.31 | 1.1020 | 13![]() |
0.955 | 0.054 | 0.197 | 0.071 |
s-decontracted aug-cc-pVTZ-J | −0.29 | 1.1073 | 14![]() |
0.958 | 0.053 | 0.180 | 0.069 |
s-decontracted aug-cc-pVTZ-J (-dfg) | −0.29 | 1.1051 | 14![]() |
0.958 | 0.053 | 0.181 | 0.069 |
aug-cc-pVTZ-Jmod67 | −0.29 | 1.1015 | 14![]() |
0.958 | 0.053 | 0.180 | 0.069 |
Exp.12 | −0.31 ± 0.04 |
Our evaluation included first some standard basis sets. With the default x2c-TZVPPall basis set the calculated ρ(0) values, which are reflected in the very large value of C in eqn (3), are significantly smaller than the fully relativistic electron density 15070 a.u.−3.86 The coefficient of determination (R2) obtained from linear regression analysis is very low and the value of “a” is far from the experimentally determined value (see Table 1). As a result, the MAD of the calculated isomer shifts from the experimental ones is also high (0.108 mm s−1). This is to be expected because of the lack of enough tight s basis functions in the x2c-TZVPPall basis set of Fe, which has been proven to be critical to obtain correct ρ(0) values.11 For the same reason, shifting from the default x2c-TZVPPall to Weigend's segmented contracted relativistic basis set for NMR shielding constants (x2c-TZVPPall-s),62 which has no additional s-space flexibility, does not result in any improvement.
The CP(PPP) basis set was originally proposed precisely for the prediction of Mössbauer parameters, albeit in a non-relativistic context.11 Nevertheless, we test it and is a clear improvement also using the X2C Hamiltonian, underlining the leading importance of core flexibility for these properties. However, even though the calculated contact densities are better than those obtained with x2c-TZVPPall, those are still significantly smaller than 15070 a.u.−3. Next, we test the aug-cc-pVTZ-J basis set, which was optimized by Sauer and coworkers for the calculation of electron paramagnetic resonance (EPR) hyperfine coupling constants.63 Interestingly, we get a good R2 value from the linear fit and the MAD is close to what was obtained with CP(PPP). However, the calculated ρ(0) values are very small and consequently the “a” value obtained from linear fit significantly deviates from the experimental calibration constant α = −0.31 ± 0.04. The probable reason behind this unusually small contact densities could be the presence of contracted s functions in aug-cc-pVTZ-J, which restricts the core flexibility of Fe. Due to the lack of sufficiently tight s-primitives, the ρ(0) values calculated using the relativistically contracted ANO-family basis set for iron, ANO-RCC-VTZP,66 has very poor linear correlation (R2 = 0.414) with the experimental isomer shifts. As a result, the fitted parameter “a” deviates significantly from the experimentally determined value (see Table 1).
Another standard basis set we have tested for Fe is DKH-def2-TZVPP, where the exponents of def2-TZVPP64 were recontracted for the DKH2 Hamiltonian with a looser contraction.43–47,65 Although the calculated ρ(0) values are still not close to the fully relativistic density, the R2 and “a” obtained from the linear regression analysis are considerably improved. Additionally, the mean absolute error and standard deviation of δcal relative to δexp are significantly better than those obtained with the x2c-TZVPPall basis set. The importance of flexibility in the s-functions is highlighted when we look at the correlation between the calculated ρ(0) and experimental isomer shifts with the s-decontracted x2c-TZVPPall, which is clearly better than standard x2c-TZVPPall. By using s-decontracted aug-cc-pVTZ-J the calculated contact densities of the iron centre are very close to the fully relativistic value 15070 a.u.−3. Moreover, we also achieve a very good correlation between ρ(0) and δexp (R2 = 0.958).
Finally, we test two more modifications of aug-cc-pVTZ-J: (a) fully decontracted s functions and removal of the three innermost s-primitives (i.e., aug-cc-pVTZ-Jmod),67 (b) fully decontracted s-functions and removal of the outermost d-, f-, and g-primitives (denoted s-decontracted aug-cc-pVTZ-J(-dfg)). The first one was recommended for the prediction of Cu(II) hyperfine coupling constants in a scalar relativistic approach,67 whereas the second modification might be useful for large Fe-complexes where linear dependencies in the basis set may arise if diffuse functions are included. These modified basis sets yield effectively indistinguishable results. However, the calculated ρ(0) values with s-decontracted aug-cc-pVTZ-J(-dfg) are closer to the fully relativistic electron density (i.e., 15070 a.u.−3) than those obtained with aug-cc-pVTZ-Jmod (see Table 1). With s-decontracted aug-cc-pVTZ-J and its two modifications, the mean absolute errors and the R2 and “a” values obtained from the linear regression analysis are the same. We note that using a set of 12 iron clusters, Kurian and Filatov also found that the computed isomer shift values are only marginally affected by the addition (or subtraction) of the tightest primitive functions to a sizable, uncontracted basis set.88 Unlike that study, however, here we find that the performance of contracted basis sets is noticeably worse compared to the large decontracted ones for the MPMIC80 set.
For each basis set, the mean absolute error of the isomer shifts calculated using the point nucleus model is slightly higher than those obtained with the Gaussian finite nucleus model (see Table 1 and Table S4 in the ESI†). Except for the aug-cc-pVTZ-J and its modifications, the linear fitting parameters of other basis sets are quite similar regardless of whether a point nucleus or finite nucleus model is used.
Our conclusion regarding basis set selection for the calculation of the Mössbauer isomer shifts is that within the X2C approach with the Gaussian finite nucleus model a basis set with very tight s functions is necessary. Among the standard and non-standard basis sets evaluated here, s-decontracted versions of aug-cc-pVTZ-J and aug-cc-pVTZ-J(-dfg) are equally good and are among the best choices. For the next step, we benchmark different DFT methods in combination with s-decontracted aug-cc-pVTZ-J, CP(PPP), and DKH-def2-TZVPP.
Basis set for Fe | Methods | a | b | C | R 2 | MADb (mm s−1) | Max. dev.c (mm s−1) | St. dev.d (mm s−1) |
---|---|---|---|---|---|---|---|---|
a a and b are the fitting coefficients obtained from the linear fit of eqn (3). C is a constant, which is very close to the calculated ρ(0) value. The units of a, b, and C are a.u.3 mm s−1, mm s−1, and a.u.−3, respectively. R2 is the coefficient of determination from the linear fit. b Mean absolute deviations of the calculated ISs with respect to the experimental ISs, where the former ones are obtained by using eqn (3). c Maximum deviation of the calculated ISs from the experimental ISs. d Standard deviation of the calculated ISs. | ||||||||
s-decontracted aug-cc-pVTZ-J | SVWN5 | −0.32 | 0.8946 | 14![]() |
0.913 | 0.071 | 0.296 | 0.098 |
BP86 | −0.32 | 0.9957 | 14![]() |
0.922 | 0.067 | 0.270 | 0.093 | |
PBE | −0.32 | 0.9970 | 14![]() |
0.921 | 0.067 | 0.267 | 0.094 | |
BLYP | −0.32 | 1.0974 | 14![]() |
0.917 | 0.070 | 0.275 | 0.096 | |
TPSS | −0.32 | 1.0045 | 14![]() |
0.935 | 0.060 | 0.243 | 0.085 | |
PBE0 | −0.28 | 1.1225 |
14![]() |
0.965 | 0.050 | 0.148 | 0.063 | |
B1LYP | −0.28 | 1.1358 | 14![]() |
0.961 | 0.054 | 0.150 | 0.066 | |
B3LYP | −0.29 | 1.1073 | 14![]() |
0.958 | 0.053 | 0.180 | 0.069 | |
TPSSh | −0.30 | 1.1976 | 14![]() |
0.954 | 0.052 | 0.227 | 0.072 | |
TPSS0 | −0.28 | 1.1425 |
14![]() |
0.967 | 0.049 | 0.122 | 0.060 | |
M06 | −0.31 | 0.9446 | 15![]() |
0.936 | 0.065 | 0.277 | 0.084 | |
LC-BLYP | −0.29 | 1.0845 | 14![]() |
0.950 | 0.062 | 0.168 | 0.075 | |
CAM-B3LYP | −0.27 | 1.2024 | 14![]() |
0.960 | 0.055 | 0.137 | 0.067 | |
ωB97X | −0.27 | 1.1689 | 15![]() |
0.958 | 0.055 | 0.147 | 0.068 | |
B2PLYP | −0.24 | 1.1367 | 14![]() |
0.939 | 0.066 | 0.183 | 0.083 | |
CP(PPP) | SVWN5 | −0.33 | 1.1214 | 14![]() |
0.910 | 0.072 | 0.295 | 0.100 |
BP86 | −0.33 | 0.8681 | 14![]() |
0.919 | 0.068 | 0.270 | 0.095 | |
PBE | −0.33 | 1.1551 | 14![]() |
0.917 | 0.068 | 0.266 | 0.096 | |
BLYP | −0.33 | 1.1620 | 14![]() |
0.914 | 0.070 | 0.275 | 0.098 | |
TPSS | −0.33 | 1.1673 | 14![]() |
0.932 | 0.061 | 0.253 | 0.087 | |
PBE0 | −0.29 | 1.0344 | 14![]() |
0.965 | 0.049 | 0.157 | 0.063 | |
B1LYP | −0.29 | 1.1515 | 14![]() |
0.963 | 0.051 | 0.159 | 0.064 | |
B3LYP | −0.30 | 0.9704 | 14![]() |
0.956 | 0.054 | 0.191 | 0.070 | |
TPSSh | −0.32 | 0.9775 | 14![]() |
0.952 | 0.053 | 0.237 | 0.073 | |
TPSS0 | −0.29 | 1.1373 | 14![]() |
0.966 | 0.049 | 0.127 | 0.062 | |
M06 | −0.32 | 1.0447 | 14![]() |
0.934 | 0.065 | 0.277 | 0.085 | |
LC-BLYP | −0.30 | 1.2454 | 14![]() |
0.951 | 0.061 | 0.177 | 0.074 | |
CAM-B3LYP | −0.29 | 1.0531 | 14![]() |
0.961 | 0.054 | 0.134 | 0.066 | |
ωB97X | −0.28 | 1.1687 | 14![]() |
0.959 | 0.055 | 0.147 | 0.068 | |
B2PLYP | −0.25 | 1.1856 | 14![]() |
0.940 | 0.066 | 0.176 | 0.082 | |
DKH-def2-TZVPP | SVWN5 | −0.33 | 0.9272 | 14![]() |
0.909 | 0.074 | 0.305 | 0.101 |
BP86 | −0.33 | 0.9640 | 14![]() |
0.919 | 0.069 | 0.277 | 0.095 | |
PBE | −0.33 | 1.1343 | 14![]() |
0.919 | 0.069 | 0.273 | 0.095 | |
BLYP | −0.33 | 1.0651 | 14![]() |
0.914 | 0.072 | 0.280 | 0.098 | |
TPSS | −0.33 | 1.1151 | 14![]() |
0.934 | 0.061 | 0.247 | 0.086 | |
PBE0 | −0.29 | 1.0853 | 14![]() |
0.964 | 0.051 | 0.156 | 0.063 | |
B1LYP | −0.29 | 1.1600 | 14![]() |
0.962 | 0.054 | 0.161 | 0.065 | |
B3LYP | −0.30 | 1.0315 | 14![]() |
0.956 | 0.055 | 0.190 | 0.070 | |
TPSSh | −0.32 | 0.9873 | 14![]() |
0.953 | 0.053 | 0.220 | 0.072 | |
TPSS0 | −0.29 | 1.0877 | 14![]() |
0.967 | 0.049 | 0.127 | 0.060 | |
M06 | −0.32 | 1.2105 | 14![]() |
0.934 | 0.067 | 0.282 | 0.086 | |
LC-BLYP | −0.30 | 1.0128 | 14![]() |
0.945 | 0.065 | 0.163 | 0.078 | |
CAM-B3LYP | −0.29 | 1.2614 | 14![]() |
0.957 | 0.058 | 0.138 | 0.069 | |
ωB97X | −0.28 | 1.2380 | 14![]() |
0.957 | 0.056 | 0.146 | 0.069 | |
B2PLYP | −0.25 | 1.1934 | 14![]() |
0.938 | 0.66 | 0.178 | 0.083 |
Unlike the different basis sets, using different density functionals does not have much influence on the calculated contact densities. Based on the R2 and mean absolute errors, climbing the rungs of “Jacob's ladder” improves accuracy gradually from the 1st to the 4th rung. As in earlier findings, hybrid functionals demonstrate a significant improvement over the performance of pure GGA and meta-GGA approaches.26,27,33,88,89 The inferior correlations obtained with pure density functionals can be attributed to their incorrect behavior near the nucleus.90 The only exception is the extensively parametrized hybrid functional M06, which offers accuracy similar to the pure meta-GGA functionals. Among the hybrid functionals, range separation offers no benefit over the global hybrid variants, and it does more harm than good when there is no short-range HF-exchange involved (e.g., in LC-BLYP). Contrary to what Römelt et al.18 found with the scalar-relativistic ZORA Hamiltonian41,91,92 on a much smaller dataset, the performance of the double hybrid functional for the extensive MPMIC80 set is inferior to standard hybrids like PBE0 and TPSS0.
For the TPSS exchange- and correlation-based hybrid functionals, increasing the percentage of exact exchange from 10% to 25% improves its accuracy. Therefore the question may arise, what is the optimum percentage of exact exchange for predicting the isomer shifts in MPMIC80? To investigate this, we chose the BLYP and TPSS functionals and their hybrid counterparts with varying percentage of exact (Hartree–Fock, HF) exchange (i.e., %HFx) while employing the s-decontracted aug-cc-pVTZ-J basis set for iron. Polynomial fitting of MADs with respect to the %HFx suggest minima near 30% for the BLYP-based and near 25% for the TPSS-based functionals (see Fig. 3 and Table S6 in the ESI†). The MAD gap between the BLYP- and TPSS-based functionals decreases with the increase of exact exchange. We obtained the lowest mean absolute error with TPSS0 (25% HFx) among all functionals, and it starts to increase rapidly with higher percentage of exact exchange. For the present dataset the fitting parameter “a” from linear regression has strong dependence on the exact exchange, but it is not sensitive to specific exchange and correlation combinations (Fig. S1, ESI†). Finally, we note that despite having 27% HF exchange, the M06 functional gives inferior electron densities and hence its performance for isomer shifts is inferior to other hybrid functionals (Table 2).
As the GGA exchange and correlation parts of B2PLYP (i.e., B88 and LYP) work fine for our dataset, the origin of the poor performance of the double hybrid functional could be either the PT2 correlation part or the amount of exact exchange used. We note that for the calculation of Cu(II) hyperfine coupling constants, the use of relaxed densities for the PT2 part was recommended.67 Following the same protocol and using the s-decontracted aug-cc-pVTZ-J basis set for Fe improves the “a” value of the linear fit from −0.24 to −0.27. However, R2 and mean absolute error are worse than the unrelaxed density calculations, allowing us to rule out the first possible source of error (see Table S7, ESI†). Hence, the high fraction of exact exchange (53%) used in the B2PLYP seems to be the principal reason behind its poorer performance compared to standard hybrid DFT functionals.
Earlier studies have shown that spin-component scaled double hybrids are more accurate than simple double hybrids like B2PLYP for calculating the energetics and spectral properties of transition metal complexes.93–96 The contact densities computed using the DSD-PBEP86 functional and the s-decontracted aug-cc-pVTZ-J basis set for iron, along with the relaxed PT2 density, were fitted to the experimental isomer shifts. This resulted in the fitting parameter a = −0.26 with R2 = 0.895 (see Fig. S2, ESI†). The relatively high percentage (68%) of HF exchange used in DSD-PBEP86 could be the reason behind its marginally poorer performance than B2PLYP. However, for most of the quintet FeII and sextet FeIII complexes of our dataset, the isomer shifts calculated using DSD-PBEP86 are closer to the experimental values than δcal using TPSS0 or B2PLYP (see Table S11, ESI†). We note in passing that a double hybrid calculation with relaxed density is approximately three times more expensive and requires more memory allocation than a similar calculation with unrelaxed density. The “unrelaxed” PT2 density corresponds to simple PT2 expectation value density, whereas the “relaxed” one incorporates orbital relaxation.
In conclusion, hybrid functionals incorporating ca. 25–30% exact exchange offer the best linear correlation and accuracy for predicting isomer shifts. Among all methods tested, our best picks are TPSS0 and PBE0, combined with the s-decontracted aug-cc-pVTZ-J basis set for Fe.
With respect to the literature, we note that using 15 iron-containing compounds and the CP(PPP) basis set for Fe, an R2 = 0.972 was obtained from the linear regression analysis of the nonrelativistic electron densities calculated with B3LYP against experimental isomer shifts.11 The reported standard deviation for the theoretical prediction of isomer shifts was 0.09 mm s−1. Later, with a slightly larger and more diverse set, B3LYP offered a marginally better correlation (R2 = 0.980) and a standard deviation of 0.09 mm s−1, while using a larger dataset, Pápai and Vankó found R2 = 0.975 and MAD = 0.06 mm s−1 for B3LYP.30 However, Kurian and Filatov found a better correlation with experimental results using BH&HLYP instead of B3LYP, which might be a result of employing a small dataset.88 Using a dataset comprising 69 iron compounds, Comas-Vilà and Salvador showed that replacing the density at the iron nucleus with the density integrated in a sphere of radius 0.06 au surrounding the iron center can provide excellent correlation (R2 = 0.976) when using the conventional def2-TZVP basis set.89 In a recent study, using 20 molecular Fe complexes, Gallenkamp et al. found the best performance with TPSSh and PBE0 (R2= 0.978 and 0.976) with mean absolute errors of 0.05 and 0.06 mm s−1, respectively.33 With a set of 21 iron-complexes, using the scalar-relativistic ZORA Hamiltonian and ZORA-def2-TZVP basis set for Fe, hybrid, and double hybrid functionals offered similar accuracy (R2 = 0.970 and Standard deviation = 0.08 mm s−1).18 Employing the DKH2 Hamiltonian and custom-defined basis sets for Fe did not offer any further improvement in the performance of hybrid functionals.22 Although our best pick, TPSS0 with the s-decontracted aug-cc-pVTZ-J basis set for Fe, has a slightly smaller R2 value, the standard deviation is better than what was reported previously. The smaller R2 can be attributed to using a dataset that is more than three times larger than what was used in ref. 18 and 22, Perdew and coworkers argue that highly parameterized density functional methods are often significantly inferior to the functionals developed by constraint satisfaction while calculating the electron densities.97 Among the 128 functionals tested in their work, PBE0 is one of the best for calculating electron density distributions compared to the all-electron coupled cluster singles and doubles (CCSD-full) densities. PBE0 and TPSS0 should therefore yield reliable electron densities and, hence, good performance for Mössbauer isomer shift computations.
CS = δ + SODS | (5) |
As we already have a curated subset of molecules (21 Fe-complexes) with high-T experimental isomer shirts, we can examine whether the solution proposed by Noodleman is transferable (see Table S2 in the ESI†). For this purpose, we chose the TPSS0 functional with the s-decontracted aug-cc-pVTZ-J basis set for iron. Considering the 80 low-temperature isomer shifts, linear regression analysis of ρ(0) vs. δexp gives R2 = 0.967. Adding the 21 high-T isomer shifts to the MPMIC80 and refitting the linear equation we obtain a significantly lower R2 value of 0.922. However, if Noodleman's correction is applied to the 21 high-T isomer shifts, linear regression analysis for the mixed set of 101 compounds almost recovers the correlation obtained with the original MPMIC80 set (see Fig. 4a–c). The same trend is observed with PBE0/s-decontracted aug-cc-pVTZ-J(Fe) and B3LYP/CP(PPP)(Fe) (see Fig. S3 and S4, ESI†).
Now, refitting a high-T isomer shift correction factor against 21 complexes we get the value 0.16 mm s−1. However, using a larger IS correction factor only marginally improves correlation compared to Noodleman's correction (see Fig. 4c and d).
The mean absolute and root-mean-square errors of calculated |ΔEQ| with respect to the experimental |ΔEQ| for ten different basis sets evaluated using B3LYP functional, finite nucleus model, and first-order picture change effect are listed in Table 3. Among the six standard basis sets, DKH-def2-TZVPP and x2c-TZVPPall are the best and the worst performer, respectively. Unlike what we observed for isomer shifts, x2c-TZVPPall-s basis set, developed for NMR shielding constants, offers a noticeably better accuracy than x2c-TZVPPall. Decontraction of the s primitives does more harm than good for x2c-TZVPPall.
Basis set for Fe | MAD (mm s−1) | RMSD (mm s−1) |
---|---|---|
a For the results with point nucleus model, see Table S8 in the ESI. b With CP(PPP) basis set for Fe and def2-TZVP for the ligand atoms. | ||
CP(PPP) | 0.312 | 0.451 |
x2c-TZVPPall | 0.357 | 0.504 |
x2c-TZVPPall-s | 0.299 | 0.436 |
aug-cc-pVTZ-J | 0.352 | 0.495 |
DKH-def2-TZVPP | 0.241 | 0.381 |
ANO-RCC-VTZP | 0.283 | 0.421 |
s-decontracted x2c-TZVPPall | 0.361 | 0.514 |
s-decontracted aug-cc-pVTZ-J | 0.347 | 0.487 |
s-decontracted aug-cc-pVTZ-J (-dfg) | 0.349 | 0.489 |
aug-cc-pVTZ-Jmod67 | 0.347 | 0.487 |
Non-relativisticb | 0.331 | 0.443 |
Next, for the evaluation of DFT methods to calculate |ΔEQ| values three basis sets were chosen: CP(PPP), DKH-def2-TZVPP, and s-decontracted aug-cc-pVTZ-J. The MAD and RMSD statistics for each the method and basis set combinations are listed in Table 4. (For the error statistics obtained using the s-decontracted aug-cc-pVTZ-J(-dfg) basis set for Fe, see Table S9 in the ESI†).
Basis set for Fe | Methods | MAD (mm s−1) | RMSD (mm s−1) |
---|---|---|---|
s-decontracted aug-cc-pVTZ-J | SVWN5 | 0.319 | 0.495 |
BP86 | 0.274 | 0.450 | |
PBE | 0.280 | 0.457 | |
BLYP | 0.272 | 0.450 | |
TPSS | 0.269 | 0.437 | |
PBE0 | 0.374 | 0.523 | |
B1LYP | 0.404 | 0.557 | |
B3LYP | 0.347 | 0.487 | |
TPSSh | 0.261 | 0.410 | |
TPSS0 | 0.383 | 0.525 | |
M06 | 0.544 | 0.949 | |
LC-BLYP | 0.432 | 0.562 | |
CAM-B3LYP | 0.474 | 0.632 | |
ωB97X | 0.493 | 0.649 | |
B2PLYP | 0.648 | 0.857 | |
CP(PPP) | SVWN5 | 0.348 | 0.516 |
BP86 | 0.308 | 0.473 | |
PBE | 0.319 | 0.484 | |
BLYP | 0.303 | 0.475 | |
TPSS | 0.305 | 0.463 | |
PBE0 | 0.326 | 0.457 | |
B1LYP | 0.346 | 0.480 | |
B3LYP | 0.312 | 0.451 | |
TPSSh | 0.251 | 0.405 | |
TPSS0 | 0.353 | 0.489 | |
M06 | 0.500 | 0.891 | |
LC-BLYP | 0.371 | 0.487 | |
CAM-B3LYP | 0.409 | 0.553 | |
ωB97X | 0.419 | 0.546 | |
B2PLYP | 0.587 | 0.779 | |
DKH-def2-TZVPP | SVWN5 | 0.451 | 0.602 |
BP86 | 0.398 | 0.545 | |
PBE | 0.408 | 0.557 | |
BLYP | 0.396 | 0.542 | |
TPSS | 0.396 | 0.534 | |
PBE0 | 0.258 | 0.399 | |
B1LYP | 0.258 | 0.392 | |
B3LYP | 0.241 | 0.381 | |
TPSSh | 0.291 | 0.431 | |
TPSS0 | 0.270 | 0.412 | |
M06 | 0.423 | 0.813 | |
LC-BLYP | 0.270 | 0.381 | |
CAM-B3LYP | 0.301 | 0.440 | |
ωB97X | 0.305 | 0.423 | |
B2PLYP | 0.466 | 0.635 |
For all these basis sets, GGA and meta-GGA functionals offer better accuracy than the LDA functional. Importantly, and unlike what we observed for the isomer shifts, hybrid functionals produce larger errors than GGA and meta-GGA functionals. Range separation leads to deterioration rather than improvement. The only exception is when the DKH-def2-TZVPP functional is used, in which case hybrid functionals are better performers than the 2nd and 3rd rung functionals. Interestingly, using either the CP(PPP) or the s-decontracted aug-cc-pVTZ-J basis set for Fe combined with TPSSh functional provides noticeably lower mean absolute error than TPSS0.
Once again, B2PLYP is the worst performer among all functionals tested. While using the s-decontracted aug-cc-pVTZ-J basis set for Fe, employing relaxed density for PT2 correlation improve its accuracy (MAD goes down from 0.648 to 0.485 mm s−1). A large percentage of exact exchange in B2PLYP also contributes a significant share to its large mean absolute error. Analyzing the results obtained from varying percentage of exact exchange in PBE-, BLYP-, and TPSS-based hybrid functionals, it is evident that beyond 10% the mean absolute error in calculated |ΔEQ| compared to experiment increases rapidly (see Table S10 in the ESI†). However, pure GGA and meta-GGA functionals are also good and cheaper alternatives. This observation agrees with the recommendation of Nemykin and Hadt for using the pure GGA functional BPW91 over B3LYP to predict quadrupole splittings in ferrocenes accurately.26 Overall, the quadrupole splittings are harder to reproduce systematically with equally high level of reliability as the isomer shifts and remain a challenge that will have to be addressed more satisfactorily in future studies.
So far, we have used the 57Fe nuclear quadrupole moment Q = 0.16b to calculate the quadrupole splitting values. A reasonable question is whether there is anything to be gained by refitting the quadrupole moment of the 57Fe nucleus against the MPMIC80 set. Using the s-decontracted aug-cc-pVTZ-J basis set for iron, we selected nine representative functionals for this purpose and plotted the experimental QSs following eqn (4). The new Q(57Fe) value is determined from the slope of the linear fit to the plot (see Table 5). With TPSS-based hybrid functionals, increasing the percentage of HF-exchange gradually decreases the fitted Q(57Fe).
Functionals | Calculated Q(57Fe)a | MAD (mm s−1) | RMSD (mm s−1) |
---|---|---|---|
a Slope of the linear fit to the plot of experimental QS vs.![]() |
|||
BP86 | 0.165 | 0.268 | 0.446 |
PBE | 0.166 | 0.271 | 0.452 |
TPSS | 0.167 | 0.254 | 0.429 |
TPSSh | 0.155 | 0.248 | 0.405 |
B3LYP | 0.141 | 0.272 | 0.408 |
PBE0 | 0.139 | 0.285 | 0.424 |
TPSS0 | 0.141 | 0.303 | 0.435 |
TPSS30 | 0.137 | 0.328 | 0.456 |
TPSS50 | 0.125 | 0.429 | 0.569 |
Using the nuclear quadrupole moment thus obtained, we reevaluated the quadrupole splittings and calculated the MAD and RMSD errors against the experimental ΔEQ values (see Table 5). It can be seen that refitting the 57Fe NQMs only marginally improves the performance of pure GGA and mGGA functionals. This improvement is negligible for the hybrid functionals at a small fraction of HF exchange, and increasing the %HFx also increases the performance gap.
• Irrespective of the choice of basis sets and density functional methods, use of finite nucleus model and picture change effect is recommended with the scalar relativistic X2C-Hamiltonian.
• Tight s-primitives are necessary for the calculation of contact densities close to the fully relativistic limit (ρ = 15070 a.u.−3). Consequently, standard basis sets that are either general-purpose or have been optimized for other properties (e.g., for NMR chemical shifts) are inappropriate poor choices for these calculations, whereas the s-decontracted version of aug-cc-pVTZ-J is a solid choice.
• For isomer shifts, 25–30% is the optimum percentage of exact (HF) exchange for hybrid functionals with scalar relativistic X2C Hamiltonian. Best accuracy was found with the TPSS0 and PBE0 functionals with s-decontracted aug-cc-pVTZ-J. Range separation has no benefit over global hybrid functionals.
• Double hybrid functionals are not recommended. The performance of B2PLYP is worse than any hybrid functional, which is most likely the result of the high percentage of exact exchange used.
• Refitting the original 80 isomer shifts and 21 high-T isomer shifts adjusted to 4.2 K by a linear correction against the calculated contact densities almost fully recovers the R2 value obtained with the low-temperature MPMIC80.
• Quadruple splitting values are less systematically predicted with DFT. Unlike isomer shifts, pure GGA, meta-GGA, and hybrid functionals with small fraction of HF exchange (∼10%) are preferred for these calculations, but there is considerable room for improvement.
Footnote |
† Electronic supplementary information (ESI) available. See DOI: https://doi.org/10.1039/d4cp00431k |
This journal is © the Owner Societies 2024 |