Prediction of 57Fe Mössbauer Nuclear Quadrupole Splittings with Hybrid and Double-Hybrid Density Functionals

TL;DR

This paper compares hybrid and double-hybrid density functionals for predicting nuclear quadrupole splittings in iron-containing molecules, identifying the most accurate and efficient methods.

Contribution

The study introduces a systematic evaluation of newly developed and widely used density functionals for Mössbauer spectroscopy predictions.

Findings

PBE-0DH double-hybrid functional achieves the lowest mean absolute error (0.20 mm/s) for 57Fe NQS predictions.

rSCAN38 hybrid functional offers a balance of accuracy (0.25 mm/s MAE) and computational efficiency.

Double-hybrid functionals outperform hybrid ones for strongly correlated systems like ferrocene.

Abstract

As a crucial parameter in Mössbauer spectroscopy, nuclear quadrupole splitting (NQS) exhibits a strong dependence on quantum chemistry methods, which makes accurate theoretical predictions challenging. Meanwhile, the continuous emergence of new density functionals presents opportunities to advance current NQS research. In this study, we evaluate the performance of eleven hybrid density functionals and twelve double-hybrid density functionals, selected from widely used functionals and newly developed functionals, in predicting the NQS values of the 57Fe nuclide for 32 iron-containing molecules within about 70 atoms. The calculations have incorporated scalar relativistic effects using the exact two-component (X2C) Hamiltonian. In general, the double-hybrid functional PBE-0DH demonstrates superior performance compared to the experimental values, achieving a mean absolute error (MAE) of 0.20 mm/s. Meanwhile, rSCAN38 is the best hybrid functional for our database with an MAE = 0.25 mm/s, and it offers a significant advantage in computational efficiency over PBE-0DH. The +/− sign of NQS has also been considered in our error statistics when it has a clear physical meaning; if neglected, the errors of many functionals decrease, but PBE-0DH and rSCAN38 remain unaffected. Notably, when calculating ferrocene [Fe(C5H5)2], which involves strong static correlations, all hybrid functionals that incorporate more than 10% exact exchange fail, while several double-hybrid functionals continue to deliver reliable results. In addition, we encountered two particularly challenging species characterized by strong static correlations: [Fe(H2O)5NO]2+ and FeO2−-porphyrin. Unfortunately, none of the density functionals tested in our study yielded satisfactory results for the two cases since the density functional theory (DFT) is a single-determinant approach, and it is imperative to explore large-scale multi-configurational methods for these species. This research offers valuable guidance for selecting density functionals in Mössbauer NQS calculations and serves as a reference point for the future development of new density functionals.

Figures

Full text

1. Introduction

Any atomic nucleus with a nuclear spin quantum number (I) greater than 1/2 displays an ellipsoidal distribution of nuclear charge, leading to a quadrupole electric field surrounding the nucleus, referred to as the nuclear quadrupole moment (NQM). It has been found that more than 40% of stable atomic nuclei possess NQMs. NQM can interact with the electric field gradient (EFG) produced by surrounding charges, known as the nuclear quadrupole interaction (NQI). NQI can be investigated using a range of spectroscopic techniques, each yielding a distinct observable. For instance, in Mössbauer spectroscopy, the NQI is quantified through nuclear quadrupole splitting (NQS), while, in gas-phase microwave spectroscopy, NQI is related to the nuclear quadrupole coupling constant (NQCC).

In Mössbauer spectroscopy, three distinct types of hyperfine interactions can be identified [1]:

• The electric monopole interactions between protons in atomic nuclei and electrons (primarily s-electrons and, to a lesser extent, p- or $[eqn]$ -electrons caused by relativistic effects) extend into the nuclear region, which can be measured through the isomer shifts ( $[eqn]$ ).
• The NQM of a nucleus with $[eqn]$ 1/2 interacts with the EFG generated by asymmetric environmental charges (mainly composed of p- and d-electrons from the target atom as well as nuclei and electrons from neighboring atoms), leading to the aforementioned NQS.
• The interactions between the atomic nuclear magnetic dipole moment and the surrounding magnetic field contribute to the Mössbauer magnetic hyperfine Zeeman splitting.

Among these, two essential Mössbauer parameters, namely $[eqn]$ and NQS, can be theoretically examined using modern quantum chemical methods by analyzing the calculated contact density (or effective contact density) [2] and EFG [3], respectively. Of the more than 80 nuclides exhibiting the Mössbauer effects, the majority of theoretical investigations have concentrated on the ^57^Fe nuclide, which serves as a valuable indicator for deducing the bonding characteristics, valence states, and local spins of iron atoms within various compounds. In the literature, nearly all the theoretical studies on the Mössbauer spectroscopy of iron compounds have been performed using density functional theory (DFT), but only a limited number of recent studies among them have taken scalar relativistic effects into account (cf. the summaries in References [2,3]). Given the strong electronic correlations frequently arising in iron compounds, traditional density functionals often fall short in effectively addressing these systems. Consequently, it is essential to identify more appropriate density functionals for accurately calculating the Mössbauer parameters of iron compounds.

There have been many recent studies in the literature that evaluate density functionals (for example, see References [4,5,6,7,8,9,10,11,12,13,14,15]), but the conclusions drawn from these investigations often differ significantly, influenced by the specific properties being examined and the molecular systems under consideration. Notably, some newly developed hybrid functionals and double-hybrid functionals exhibit better accuracy than others [6,7,8,9,12,14]. In contrast, traditional functionals frequently inaccurately describe the electronic structures and properties of strongly correlated systems. In the most extreme cases, such as diatomic CuX (X = H, F, Cl, Br, and I) molecules, certain commonly used functionals can produce qualitative errors regarding the $[eqn]$ signs of EFG for the copper nuclide [16,17,18,19,20]. Over the past decade, many new hybrid and double-hybrid density functionals have been developed, but, to our knowledge, there is currently no comprehensive database dedicated to the functional development that accounts for the electron distribution around the nucleus. Consequently, the accuracy of these functionals in calculating Mössbauer parameters, particularly for iron compounds, remains uncertain.

In a recent paper [3], we employed the spin-free X2C (exact two-component) [21,22,23,24] relativistic DFT method to compute the NQS values of ^57^Fe for seventeen iron compounds, each containing up to 60 atoms. We utilized two randomly selected density functionals: the hybrid functional B3LYP and the double-hybrid functional DSD-PBEP86. Our findings indicated that the latter functional yielded better results than the former. In this paper, we extend our investigation to include 32 iron compounds, incorporating a range of hybrid and double-hybrid functionals, with a particular emphasis on the latest developments in this area.

This paper is organized as follows: The results are presented in Section 2 and discussed in Section 3, respectively. Section 4 presents the fundamental formulas, defines relevant symbols, and outlines the computational details. Finally, some conclusions are drawn in the last section.

2. Results

2.1. Preliminary Evaluation of Density Functionals

Given the vast number of hybrid and double-hybrid functionals available for testing (see Section 4.2), a preliminary screening is essential. Small molecules containing iron often exhibit open-shell degenerate states with significant multi-configurational characteristics or lack reliable experimental data [25], and, therefore, we chose the closed-shell CuF molecule as our initial testing system, which was originally proposed by Bast and Schwerdtfeger for evaluating new functionals [17]. While CuF may appear straightforward for calculations at first glance, it presents considerable challenges when approached with conventional functionals [17,18,19,20], as well as with advanced ab initio methods [26]. Compounds containing 3d metals are widely recognized for their strong correlations, which make them particularly challenging to calculate [27]. Several explanations have been proposed to account for this phenomenon from the perspective of static–dynamic correlations [28], whereas the phenomenological explanation in condensed matter physics seems not helpful, especially given that DFT with a Hubbard U correction usually performs worse than hybrid DFT methods [29]. A prevalent explanation suggests that these strong correlations stem from the competition between the 3 $[eqn]$ and 3 $[eqn]$ occupations, along with the large discrepancies in their correlation errors caused by the different spatial distributions of the two types of 3d orbitals. Another explanation claims that significant $[eqn]$ local excitations contribute to the emergence of strong correlations. Anyhow, the strong correlations present in 3d metal systems can be effectively handled using multi-reference methods that incorporate the so-called “double-d shell effect” [30,31] (that is, a 3 $[eqn]$ -shell has to be included in the active space, which is essentially the 4d-shell while allowing for mixing with 4p when symmetry permits). Additionally, in cases where the Hartree–Fock (HF) reference wavefunctions are qualitatively correct, single-reference high-order coupled-cluster methods are also applicable. Given that both copper and iron compounds share similar sources of error, it is reasonable to anticipate comparable performance across various functionals applied to these molecules.

CuF has been studied using microwave rotational spectroscopy. From the experimental NQCC value of $[eqn]$ = 21.9562 MHz for ^63^CuF [32] and NMQ value of Q = −0.220 barn for ^63^Cu [33], the experimental principal value $[eqn]$ of EFG is determined to be −0.425 a.u. using Equation (1). Table S1 in the Supplementary Materials compiles the theoretical $[eqn]$ results for the Cu nucleus in CuF ( $[eqn]$ = 1.7449 Å [32]) computed by hybrid functionals, double-hybrid functionals, and advanced ab initio methods. As seen in the table, $[eqn]$ = −0.449 a.u. by CCSD(T) closely aligns with the experimental value, and this result improves to −0.433 a.u. when accounting for the correlations from the innermost core electrons. Further enhancement of the results necessitates the corrections from the complete basis set limit and high-order excitations beyond perturbative triples [26]. The EFG results obtained from three multi-configurational methods utilizing nine active orbitals (corresponding to Cu 3d4s and F 2p) are significantly inaccurate. However, after incorporating the Cu 3 $[eqn]$ -shell into the active space, leading to a total of fourteen active orbitals, the new EFG results are qualitatively accurate. Notably, the MRAQCC result of −0.367 a.u. is only 0.058 a.u. larger than the experimental value. To elucidate the discrepancies in the results, we conduct a natural bond orbital (NBO) analysis [34] on the canonical or natural orbitals produced by the aforementioned methods, with the natural electron configurations [35] of Cu detailed in Table S2. For the three types of reference wavefunctions, the copper atom in the CuF molecule demonstrates an occupation of approximately 3 $[eqn]$ 4 $[eqn]$ as predicted by HF and CASSCF(14o), while CASSCF(9o) supports the 3 $[eqn]$ 4 $[eqn]$ occupation. It is evident that CASSCF(9o) markedly overestimates the energy associated with the 3 $[eqn]$ occupation due to a substantial correlation energy loss. Therefore, the strong correlations observed in CuF arise from the competition among different 3d orbital occupations, rather than from the 3 $[eqn]$ 4d correlations (CASSCF(14o) has contained all the $[eqn]$ excitations but its $[eqn]$ result is still not satisfactory), which is likely applicable to other compounds involving 3d metals.

In Table S1, only half of the hybrid functionals and half of the double-hybrid functionals exhibit relatively small errors within 0.22 a.u. Additionally, the double-hybrid functional DSD-PBEP86 also performed well in our previous study involving seventeen iron compounds [3], and therefore eleven hybrid functionals and fourteen double-hybrid functionals are initially selected for further analysis. Table S2 also details the natural electron configurations of copper determined by two well-performing functionals (BH&HLYP and r2SCAN-CIDH) and two less effective functionals (B3LYP and B2PLYP), revealing that these configurations are quite similar irrespective of the EFG results. Therefore, it can be concluded that the EFG error in DFT calculations primarily originates from the inner electrons, which, however, are not addressed in the NBO analysis. Since the correlations in core orbitals are relatively weak, the core orbitals at the HF level may be regarded as nearly exact and can be used to cure the poor functionals. For instance, when examining the B3LYP functional, substituting each of its core orbitals with the corresponding HF orbital reveals a significant change in the $[eqn]$ value for the three Cu 3p-like orbitals, from the original 0.350 a.u. to −0.375 a.u. Moreover, increasing the proportion of the HF term in the functional recipe yields a similar effect [17], as evidenced by the hybrid functional results presented in Table S1, including those from the PBEn and TPSSn families. Therefore, the $[eqn]$ error associated with the copper nucleus is influenced not only by the treatment of valence electrons but also by the inner 3p electrons. This conclusion aligns with findings from the study of CuCl as well [19], and has been attributed to the defects in the exchange part of the functional as noted by Schwerdtfeger et al. [16,17].

The second testing system examined is [Fe=O(tmc)(NCCH_3_)]^2^+ (i.e., the molecule 24 in the test set), which exhibits an experimental NQS value of $[eqn]$ = −1.24 mm/s for the ^57^Fe nuclide [36]. Our previous research indicates that its $[eqn]$ is difficult to predict using B3LYP [3], and the complete results by different functionals have been provided in the next subsection. Following initial screening, all eleven hybrid functionals demonstrate satisfactory performance, and the predicted $[eqn]$ values of ^57^Fe are close to the experimental one. For the double-hybrid functionals, however, two of them yield inaccurate predictions: 2.19 mm/s by $[eqn]$ SCAN-QIDH and 2.39 mm/s by P $[eqn]$ SCAN69. Consequently, these two functionals have to be excluded from further calculations.

In conclusion, our formal evaluation includes a total of eleven hybrid functionals and twelve double-hybrid functionals (cf. Tables S3 and S4 for the name list).

2.2. Comparison of Density Functionals

The thirty-two medium-sized iron compounds are illustrated in Figure 1, and their optimized Cartesian coordinates are taken from the literature [36,37,38,39]. Considering the limitations of the computational efficiency of double-hybrid functionals, only the molecules containing approximately 70 atoms or fewer are selected, with the largest being Fe(PyO)I(ArS)(CO)2_PPh_3 (19) which comprises 71 atoms. Notably, the ferrocene molecule [Fe(C_5_H_5_)2; 26] has two distinct conformations: the eclipsed global minimum in $[eqn]$ symmetry and the staggered saddle point in $[eqn]$ symmetry with a tiny energy difference between them [40,41]. Our tests indicate that there is a negligible difference in $[eqn]$ between these two conformations, and therefore the coordinates of the staggered conformation from Reference [38] are adopted.

Figure 2 shows the error distribution of theoretical $[eqn]$ values for the selected 32 iron compounds, and the raw data along with experimental $[eqn]$ values are provided in Table S3 for eleven hybrid functionals and Table S4 for twelve double-hybrid functionals, respectively, in the Supplementary Materials. The signs of the experimental $[eqn]$ values are taken from the literature when available, including the experimental measurements for [Fe(H_2_O)5_NO]^2^+ (molecule 4) [42], [FeS_4_C_8_O_4]^2^− (7) [43], [Fe(SPh)4]^2^− (10) [43], and the molecules collected in References [36,39]. In cases where experimental data were not identified, the signs are determined based on the theoretical results from this study. As discussed in Section 4.1, predicting the signs of $[eqn]$ (^57^Fe) in theoretical calculations can be difficult when $[eqn]$ 0.4 mm/s or $[eqn]$ ; these specific iron compounds are indicated with an asterisk in the tables.

In contrast to the calculations of nuclear contact density, which remains largely unaffected by the choice of functionals [13], the values of $[eqn]$ are significantly influenced by different functionals. As seen in Tables S3 and S4, all the selected hybrid and double-hybrid functionals generally provide reasonable predictions of $[eqn]$ values for most molecules. However, for the remaining molecules, the $[eqn]$ results exhibit considerable variation. Notably, the case of [Fe(por)(O_2_)]^−^ (molecule 15) stands out, as all functionals tend to overestimate its $[eqn]$ by 1.5 mm/s or more. It is widely known that the metal–porphyrin systems are characterized by strong static correlations. In the literature, accurate calculations of the iron porphyrin systems can be achieved using full-CI quantum Monte Carlo (FCIQMC) [44,45], generalized active space self-consistent field (GASSCF) [46], density matrix renormalization group (DMRG) [47], and selected-CI (sCI) [48,49] algorithms that can handle very large active spaces. Given that this study is limited to the DFT method, [Fe(por)(O_2_)]^−^ has to be excluded from the error analysis.

Based on the maximum error (MaxE) and mean absolute error (MAE), the leading three hybrid functionals are PBE50 ≺ TPSS38 ≺ rSCAN38, with MaxE and MAE being within 1.60 mm/s and 0.31 mm/s, respectively. The top three double-hybrid functionals are $[eqn]$ SCAN-CIDH ≺ $[eqn]$ SCAN-0DH ≺ PBE-0DH, which exhibit slightly smaller errors. Table 1 provides a comprehensive overview of the results for these six hybrid and double-hybrid functionals. Notably, PBE-0DH achieves the best performance with an MAE of 0.20 mm/s, while PBE50 ranks lowest with an MAE of 0.31 mm/s. The remaining four functionals fall within a middle range, exhibiting MAE values between 0.23 and 0.27 mm/s. When analyzing the absolute values of $[eqn]$ , as commonly reported in the literature, the MAE values of the majority of the double-hybrid functionals listed in Table S4, including the top three performers, remain unchanged. Conversely, among the hybrid functionals, only the MAE of rSCAN38 is unaffected, while the MAE values of all other hybrid functionals show a decrease. Especially, the MAE of SCAN38 is halved, allowing it to surpass PBE50 and secure the position of the third-ranked hybrid functional. This suggests that most double-hybrid functionals excel in qualitatively predicting the signs of $[eqn]$ compared to the hybrid functionals, with the exception of rSCAN38. Interestingly, the performance of many double-hybrid functionals is on par with or even inferior to that of the hybrid functional rSCAN38, highlighting significant opportunities for further optimization in the development of double-hybrid functionals.

3. Discussion

Among the six functionals presented in Table 1, the most significant discrepancies with MaxE exceeding 1.0 mm/s are observed in three molecules (excluding molecule 15), namely [Fe(H_2_O)_5_NO]^2^+ (molecule 4), Fe(PyO)I(ArS)(CO)2_PPh_3 (19), and ferrocene (26). The MaxE value for molecule 19 arises from an incorrect sign of $[eqn]$ , as previously discussed, which occurs only in the hybrid functionals PBE50 and TPSS38, whereas the central radical FeNO^2^+ in 4, according to the DMRG calculations by Boguslawski et al. [50], is a difficult system with strong static correlations, which should be also true for molecule 4. Regarding ferrocene, while it is typically classified as a single-reference system, its ground state exhibits significant correlation effects [51]. In our results, all hybrid functionals and the truncated hybrid components of double-hybrid functionals consistently overestimate its $[eqn]$ by more than 1.0 mm/s. Only the TPSSh hybrid functional that incorporates a small fraction of exact exchange ( $[eqn]$ = 0.1) and the pure functional TPSS are able to accurately calculate its $[eqn]$ , being 2.69 and 2.36 mm/s, respectively. It is evident that the inclusion of the exact exchange term contributes to the deterioration of these results. A closer examination through NBO analysis indicates that the iron atom in ferrocene exhibits the configuration 3 $[eqn]$ 4 $[eqn]$ according to HF calculations, while TPSS suggests 3 $[eqn]$ 4 $[eqn]$ (that is, Fe^0^ with two 4s electrons back donated to 3d). So, this is the competition between the 3 $[eqn]$ and 3 $[eqn]$ occupations (n = 8 now), as seen in the case of CuF. However, perhaps due to the fact that iron possesses fewer 3d electrons compared to copper, HF incorrectly predict the 3 $[eqn]$ occupation for ferrocene. Ideally, the performance of the exchange functional should be consistent with the exact exchange limit by HF. However, our calculations reveal that the TPSS exchange-only functional by removing the correlation functional continues to support the 3 $[eqn]$ 4 $[eqn]$ configuration, suggesting that this configuration is caused by the exchange functional instead of the correlation one. This implies that the favorable $[eqn]$ results seen in TPSS and TPSSh arise from a coincidentally correct 3d occupation, which is a consequence of inherent limitations in the approximate exchange functional. Similar trends have also been observed in other hybrid functional families, such as PBEn, and the truncated hybrid components of double-hybrid functionals. Upon incorporating correlations from virtual orbitals through second-order perturbation (PT2) corrections, the $[eqn]$ values for all double-hybrid functionals significantly decrease, with most errors falling below 0.5 mm/s (refer to Table S4). So, it is clear that ferrocene is distinctly marked by strong static correlations, but it also demonstrates specific dynamic correlation behaviors since it can be partially corrected through the PT2 corrections.

In double-hybrid functional calculations, three types of density are identified: the self-consistent density produced by the truncated hybrid functional (HFun), the unrelaxed PT2 density with fixed molecular orbitals (UnRlx), and the relaxed PT2 density that incorporates orbital response (Rlx) by solving the coupled-perturbed equations. Table S4 presents the $[eqn]$ results and associated errors for all three density types, which can assist in identifying the sources of these errors [13]. It is important to highlight that, for the top three double-hybrid functionals PBE-0DH, $[eqn]$ SCAN-0DH, and $[eqn]$ SCAN-CIDH, their HFun and UnRlx densities exhibit the most favorable $[eqn]$ results, which not only surpass those of all other double-hybrid functionals but also demonstrate errors that are comparable to those associated with the PBE50 functional. In most instances, the HFun and UnRlx densities yield qualitatively accurate $[eqn]$ results, while the Rlx densities provide additional minor corrections. However, for molecules 1, 6, 13, 16, and 29, it is evident that the Rlx densities from certain functionals can overcorrect $[eqn]$ , resulting in significantly poorer outcomes. This observation suggests that the proportion of PT2 term is the primary source of error for these molecules. For molecules exhibiting strong static correlations, such as 4 and 15, the HFun and UnRlx densities are already qualitatively inaccurate, rendering any corrections from the Rlx densities ineffective.

The discrepancies in the signs of $[eqn]$ presented in Tables S3 and S4 can be attributed to several key factors:

• The parameter $[eqn]$ approaches the critical threshold of 0.75, particularly illustrated by molecule 19. Notably, the $[eqn]$ values predicted by certain standard and truncated hybrid functionals even exceed the critical threshold.
• An improper proportion of the PT2 term contributes to overcorrection, as evidenced by the results of molecule 29 obtained from several double-hybrid functionals.
• The self-consistent field (SCF) iterations utilizing some (truncated) hybrid functionals converge to distinct occupation patterns within the Fe 3d-shell. For instance, in molecule 7, the contributions of $[eqn]$ electrons in 3d to the EFG tensor are minimal in the elements $[eqn]$ and $[eqn]$ (depending on the coordinate orientation employed in our calculations); however, these contributions are erroneously calculated as −1.2 a.u. by SCAN38, SCAN50, and $[eqn]$ SCAN38, leading to an incorrect sign reversal upon diagonalization. A comparable case is also observed in molecule 25, where the contributions of $[eqn]$ electrons in 3d to $[eqn]$ and $[eqn]$ are −1.0 a.u. but are significantly underestimated to be 0.1 a.u. by the truncated hybrid functional components in $[eqn]$ SCAN0-2 and DSD-PBEP86.

There may be other reasons, such as SCF converging to the configurations with varying 3d occupations, a phenomenon observed in CuF and ferrocene, but in our testing we have not encountered any instances where this leads to a reversal in the sign of $[eqn]$ (^57^Fe).

Based on the aforementioned findings, it can be concluded that there is no universally applicable functional for EFG and Mössbauer NQS calculations of the selected iron compounds. Nevertheless, by excluding the most tricky molecule 15 characterized by strong static correlations, some valuable insights can still be obtained. Among the various hybrid and double-hybrid functionals evaluated, PBE-0DH stands out as the most effective, and the calculated $[eqn]$ parameters for the tested compounds show good agreement with the experimental ones with an MAE of only 0.20 mm/s. Among the less optimal functionals, rSCAN38, $[eqn]$ SCAN-CIDH, and $[eqn]$ SCAN-0DH demonstrate comparable performance, with MAE values ranging from 0.23 to 0.25 mm/s. However, the hybrid functional rSCAN38 shows greater promise since it is superior to the two double-hybrid functionals $[eqn]$ SCAN-CIDH and $[eqn]$ SCAN-0DH in computational efficiency and has the ability to predict the sign of $[eqn]$ . Unfortunately, these newly developed functionals, along with the RIJCOSX approximation, have not been implemented in most quantum chemistry programs. Therefore, if a slight reduction in the tolerance for $[eqn]$ accuracy is acceptable and the signs of $[eqn]$ are not critical, the hybrid functionals BH&HLYP and M06-2X are both good options for predicting satisfactory $[eqn]$ results in most cases.

In the literature over the past decade, numerous studies have explored $[eqn]$ (^57^Fe) using various functionals [13,52,53,54]. However, readers may notice that the conclusions drawn in these studies differ from those presented in our research. For instance, some studies suggest that pure functionals outperform hybrid functionals (our results show it is correct only in some cases), or double-hybrid functionals do not perform as well as hybrid ones. The reasons for the varied conclusions can be attributed to several factors, including the following:

• 1.The exclusion of scalar relativistic effects in DFT calculations.
• 2.The electronic correlations present in the dataset, which may be relatively straightforward to manage or distinctly unique.
• 3.The systematic exclusion of pure functionals alongside the inclusion of newly developed hybrid and, particularly, double-hybrid functionals in the evaluation.
• 4.The optimized basis sets for contact density calculations may be inadequate for EFG.
• 5.The neglect of the sign of $[eqn]$ .
• 6.The SCF iterations converge towards a specific excited state.
• 7.The molecular structures optimized by different methods may affect the errors, among other factors.

Some of these factors may also offset one another, resulting in improved outcomes. Consequently, it is advisable to consider these aspects in future researches, and especially the “Two Wrongs Make a Right” phenomenon in terms of functionals: the incorrect behavior of electronic density predicted by the exchange functional in dealing with transition metals [16,17] might surprisingly result in the correct electronic configuration of Fe whereas the HF method yields the contrary, as evidenced in our results of ferrocene (26). Moreover, the effectiveness of the implicit solvent model in simulating the molecular crystal environment has not been comprehensively assessed in the existing literature and in this study, while it often demonstrates superior performance compared to gas phase calculations [36,52], but lacks essential information regarding crystal packing and anisotropic interactions. Alternative approaches for modeling solid environments are worth trying, including the explicit solvent model combined with QM/MM [55], embedded many-body expansion [56], and the more precise (relativistic) DFT method with periodic boundary conditions [57,58], which can significantly enhance the description of the lattice environment and thereby improve the accuracy of theoretical simulations.

4. Materials and Methods

4.1. Electric Field Gradients and Nuclear Quadrupole Interactions

Any nucleus with a nuclear spin $[eqn]$ exhibits a non-zero NQM tensor $[eqn]$ , in addition to a scalar NQM value Q as collected in the literature [33,59]. The EFG tensor $[eqn]$ provides insight into the asymmetric charge distribution surrounding the nucleus. Both $[eqn]$ and $[eqn]$ are symmetric and traceless 3 × 3 matrices. The interactions between $[eqn]$ and $[eqn]$ can lead to various NQI quantities that can be measured experimentally. Therefore, a primary objective of theoretical research is to calculate EFG as accurately as possible.

By convention, $[eqn]$ is transformed into $[eqn]$ in a special principal axis system, resulting in $[eqn]$ being represented as a diagonal matrix with the eigenvalues $[eqn]$ . Because of its traceless nature (i.e., $[eqn]$ ), the EFG tensor can be described using just two coordinate-independent parameters: the principal value $[eqn]$ and the asymmetry parameter $[eqn]$ [1]. When $[eqn]$ (i.e., $[eqn]$ = $[eqn]$ = $[eqn]$ ), the nucleus is positioned at the center of an axisymmetric molecule along the $[eqn]$ axis ( $[eqn]$ 2); in other words, there is only one independent parameter $[eqn]$ in this special case. Examples include the linear molecule HCN (for all the three nuclei), PCl_5_ in $[eqn]$ symmetry (P nucleus), and Ge(CN)4 in $[eqn]$ symmetry (Ge nucleus). It is important to highlight that, for molecules exhibiting spatial degeneracy (e.g., FeO has a doubly degenerate ground state $[eqn]$ ), it is essential to compute the sum of the $[eqn]$ tensors for all the degenerate components, and failing to do so may lead to a breakdown of symmetry.

After $[eqn]$ being available, NQCC (referred to as $[eqn]$ ) can be calculated by the definition

$$[eqn]$$ $$[eqn]$$

where Q and $[eqn]$ are expressed in barns (1 barn = $[eqn]$ $[eqn]$ ) and atomic units (a.u.), respectively, $[eqn]$ is a conversion factor, $[eqn]$ is the Rydberg constant (2 $[eqn]$ = 1 Hartree = 219,474.63137 $[eqn]$ ), $[eqn]$ = 0.5291772083 Å is the Bohr radius, and $[eqn]$ cm/s is the speed of light in the CGS unit system (cf. the notation of Equation (47) in Reference [60]).

Another important NQI quantity is the Mössbauer NQS (represented by the symbol $[eqn]$ ), which also depends on $[eqn]$ and I. The electric quadrupole interaction energy associated with the magnetic quantum number $[eqn]$ ( $[eqn]$ = $[eqn]$ , $[eqn]$ , …, I) is [1]

$$[eqn]$$

For the first excited state of ^57^Fe nuclide, I is 3/2; thus, the energy difference between $[eqn]$ and $[eqn]$ is

$$[eqn]$$ $$[eqn]$$

where Q is 0.160 barn for ^57^Fe in the I = 3/2 excited state [33], $[eqn]$ = 34.84924 × $[eqn]$ MHz (=14.413 KeV [1]) is the $[eqn]$ -radiation energy of ^57^Fe thus $[eqn]$ = 11.6248 MHz·s/mm, and the meanings of the other symbols have been clarified previously.

Due to the difficulty in measuring the sign of $[eqn]$ experimentally and the frequent inconsistent signs of $[eqn]$ between theoretical calculations and experimental measurements, only the absolute values of $[eqn]$ are usually compared in the literature. However, in the context of NQCC investigation, the signs of $[eqn]$ are crucial (see the examples of CuX [16,17,18,19,20]), implying that it is not appropriate to take the absolute values of $[eqn]$ in a general way [52]. Furthermore, the intricate electronic structures of iron compounds mean that overlooking the signs of $[eqn]$ could obscure fundamental shortcomings of theoretical approaches. There are only two exceptions where determining the sign of $[eqn]$ becomes challenging, even at advanced theoretical levels, which consequently diminishes its physical meaning.

If $[eqn]$ is very small or even zero (for example, the sulfur nucleus in SF_6_ with $[eqn]$ symmetry), some minor theoretical errors may result in a swap between $[eqn]$ and $[eqn]$ , so the absolute value of the new $[eqn]$ (i.e., the old $[eqn]$ ) remains nearly unchanged but with an opposite sign. The suggested effective range for $[eqn]$ is $[eqn]$ > 0.25 a.u., which approximately corresponds to $[eqn]$ > 0.4 mm/s for ^57^Fe, as indicated by Equation (4). $[eqn]$ approaches one, i.e., $[eqn]$ 0 and $[eqn]$ , which leads to an uncertainty regarding the sign of $[eqn]$ since $[eqn]$ and $[eqn]$ may be interchanged by theoretical errors (cf. page 95 of Reference [1]). The schematic structure for the case of $[eqn]$ 1 is the Ge nucleus in the model molecule “GeHe_2_F_4_” with $[eqn]$ symmetry, as illustrated in Figure 3. In this model system, the central Ge nucleus experiences a symmetrical charge distribution along the z-axis while the x- and y-directions show asymmetric charge distributions of equal magnitude. Consequently, $[eqn]$ 0. It has been found in real systems that even a minor adjustment in the dihedral angle can cause a reversal of the sign of $[eqn]$ when the $[eqn]$ value surpasses a specific critical point [39]. In this work, $[eqn]$ is suggested to make the sign of $[eqn]$ valid (that is, $[eqn]$ or equivalently $[eqn]$ ).

4.2. Hybrid and Double-Hybrid Density Functionals

Hybrid density functionals incorporate an exact exchange term in the HF form. These functionals primarily consist of three-parameter (3P) hybrid functionals [61,62] and one-parameter (1P) hybrid functionals [63], characterized respectively by the following energy formulas:

$$[eqn]$$

In this study, the following hybrid functionals are selected: BH&HLYP [61], B3LYP [62], CAM-B3LYP [64], PW6B95 [65], M06-2X [66], $[eqn]$ B97XD [67], PBEn family (PBE0 [63], PBE38 [68], and PBE50 [69]), TPSSn family (TPSSh [70], TPSS0, TPSS38, and TPSS50) [8], SCANn family (SCAN0 [71], SCAN38, and SCAN50) [8], rSCANn family (rSCAN0, rSCAN38, and rSCAN50) [9], and $[eqn]$ SCANn family ( $[eqn]$ SCAN0 [72], $[eqn]$ SCAN38 [9], and $[eqn]$ SCAN50 [72]). This study excludes pure functionals due to their poor performance in computing the electronic structures of 3d metal-containing systems, as found in the studies of CuX [16,17,18,19,20].

Double-hybrid density functionals enhance this approach by incorporating a PT2 term, typically without the need for orbital optimization. There are mainly two categories of double-hybrid functionals: the B2PLYP type introduced by Grimme [73] and the XYG3 type developed by Zhang, Xu, and Goddard [74]. Their general energy formulas are [75].

$$[eqn]$$

with $[eqn]$ , where subscript “SCF” means orbital optimization through SCF iterations, “H” denotes a full hybrid functional for B2PLYP or a standard hybrid functional (like B3LYP and PBE0) for XYG3, while “h” indicates the hybrid component of the double-hybrid functional in the XYG3 formulation, which is a reparameterized version of the standard hybrid functional.

Since the natural orbital and expectation value algorithms associated with the XYG3 type have not been extensively implemented, this study focuses solely on the B2PLYP type of double-hybrid functionals, including B2PLYP [73], mPW2PLYP [76], Martin’s reparameterizations of B2PLYP (B2GP-PLYP, B2K-PLYP, and B2T-PLYP) [77,78], $[eqn]$ B97X-2 [79], DSD-BLYP [80], DSD-PBEP86 [81], PBE-0DH [82] and its range-separated version RSX-0DH [83], PWPB95 [84], DSD-PBEB95 [85], PBE-QIDH [86] and its range-separated version RSX-QIDH [87], $[eqn]$ B2PLYP and $[eqn]$ B2GP-PLYP [88], $[eqn]$ PBEPP86 and $[eqn]$ B88PP86 [7], 2019 version of DSD/DOD family [6] (DOD-SCAN-D3(BJ), noDispSD-SCAN69, revDSD-PBEP86-D3(BJ), revDSD-BLYP-D3(BJ), and revDOD-PBEP86-D3(BJ)), and 2023 version of $[eqn]$ SCAN family [89] ( $[eqn]$ SCAN-0DH, $[eqn]$ SCAN-CIDH, $[eqn]$ SCAN-QIDH, $[eqn]$ SCAN0-2, P $[eqn]$ SCAN50, and P $[eqn]$ SCAN69). Certain spin-component scaled (SCS) or spin-opposite scaled (SOS) double-hybrid functionals are specifically defined for excitation energy calculations without affecting the ground state results of the parent functionals (for example, SCS/SOS-B2PLYP21 vs. B2PLYP), and, therefore, these functionals are not pertinent to the current study.

4.3. Computational Methods

Firstly, both ORCA 6.0 [90,91] and BDF 2024A [92] program packages are utilized to perform unrestricted BH&HLYP [61] calculations for each molecule independently, ensuring consistency in total energy and $[eqn]$ values. Scalar relativistic effects are addressed using the spin-free X2C relativistic Hamiltonian [21,22,23,24] along with the Gaussian-type finite-size nuclear charge distribution [93]. Compared with the approximate relativistic Hamiltonians, X2C offers significant advantages in terms of “simplicity, accuracy, and efficiency” [22] and therefore has been utilized in both this study and our previous Mössbauer parameter investigations [2,3]. The integration grid level is set to defgrid3 in ORCA and ultrafine in BDF, respectively. Since Mössbauer spectroscopy is measured in the solid phase, methanol serves as a solvent to implicitly simulate the crystalline environment using the polarizable continuum model (PCM) [94]. This approach, in conjunction with the conductor-like screening solvation model (COSMO), is deemed appropriate for organometallic systems [36,37,38,52,53,95,96,97]. To enhance the efficiency of integral calculations, the RIJCOSX [98] and aMPEC+aCOSX [99] approximations have been enabled in ORCA and BDF, respectively. For the majority of molecules, both programs yield very close total energies and $[eqn]$ values, but, for two molecules, [Fe(H_2_O)_5_NO]^2^+ (4) and [Fe(PyS)I(CO)2_PPh_3] (20), BDF produces slightly lower energies. As a result, the utilities molden2fch and fch2mkl in MOKIT 1.2.5 [100] are utilized to transfer molecular orbitals from BDF to ORCA.

Next, the ORCA program package is used to conduct a range of hybrid and double-hybrid functional calculations. To maximize convergence to the same electronic states, initial orbitals are read from the previously performed BH&HLYP calculations. Sample ORCA input files for undefined functionals can be found in the original literature associated with these functionals. The default setting of frozen core in ORCA is adopted in the double-hybrid functional calculations, that is, in addition to valence electrons, the semi-core electrons in 3s3p of Fe and (n-1)d of post-3d elements are also correlated.

All the calculations are carried out using the x2c-TZVPPall-f (for pre-3d atoms) and x2c-TZVPPall [101] relativistic basis sets. As noted by Santra et al. [13], the standard x2c-TZVPPall basis set of the iron atom performs poorly in the EFG calculations and therefore has to be slightly modified in our study: the functions with angular quantum number $[eqn]$ are decontracted to enhance flexibility in describing the electron density distribution near the nucleus; in contrast, the spherical s-functions, which do not explicitly contribute to EFG, are retained in their contracted form. Furthermore, EFG calculations often require the use of some very tight functions with $[eqn]$ to achieve saturation of the basis set, which is essential for precise EFG results [26]. Our test calculations at the HF level suggest that the x2c-TZVPPall basis set for Fe is deficient in tight p- and d-functions. To achieve acceptable convergence with the $[eqn]$ error remaining below 0.01 a.u. (or a $[eqn]$ (^57^Fe) error of less than about 0.016 mm/s; see Table S5), it is recommended to add one additional tight p-function and three tight d-functions, and their Gaussian exponents generated by the even-tempered series can be found in Reference [3].

In addition to the 32 iron compounds, the strongly correlated gas-phase molecule CuF is also calculated using a variety of hybrid and double-hybrid density functionals, as well as several ab initio methods implemented in the Molpro 2015.1 program package [102]. The single-reference ab initio methods employed include CCSD (coupled-cluster with single and double excitations) and CCSD(T) (CCSD with perturbative triple excitations) [103], while the multi-configurational methods are CASSCF (complete active space self-consistent field), internally contracted MRCI (multi-reference configuration interaction with single and double excitations) [104], and internally contracted MRAQCC (multi-reference averaged quadratic coupled-cluster with single and double excitations) [105]. The calculation process of EFG at these theoretical levels is detailed in Reference [3]. The scalar relativistic effects are calculated again at the spin-free X2C level, and the core–valence correlations from the Cu 3s3p electrons are counted in the double-hybrid functional and advanced ab initio calculations. The basis sets employed are x2c-TZVPPall-f for the fluorine atom and x2c-TZVPPall for the copper atom [101], respectively, but the latter requires some modifications, similar to the adjustments made for the iron atom, by supplementing four additional Gaussian exponents with $[eqn]$ = 5813.6650 and $[eqn]$ = 1428.9269, 491.72241, and 169.21155.

5. Conclusions

While high-precision EFG calculations for small molecules have been successfully addressed, as demonstrated in References [20,26], predicting the Mössbauer NQS for medium-sized molecules continues to pose significant challenges. In this context, the DFT method is indispensable. Consequently, identifying the “best” density functionals tailored for EFG calculations has emerged as a critical objective in this field.

In this study, the $[eqn]$ values of the ^57^Fe nuclide for 32 iron-containing molecules have been calculated using selected hybrid and double-hybrid functionals, and are compared with experimental values. In error statistics, the signs of $[eqn]$ with the considering of two exceptional cases, namely $[eqn]$ (^57^Fe) < 0.4 mm/s and $[eqn]$ , have also been examined. Our results lead us to recommend the following functionals for NQS calculations of ^57^Fe nuclide.

The double-hybrid functional PBE-0DH demonstrates strong agreement with experimental results, outperforming other functionals with an MAE of 0.20 mm/s.If computational cost is a primary concern, the hybrid functional rSCAN38 is recommended, as it exhibits a slightly larger MAE of 0.25 mm/s, while still delivering satisfactory results for most molecules.In cases where the quantum chemistry program does not support the aforementioned functionals, the older hybrid functionals BH&HLYP and M06-2X can be utilized, albeit with a greater MAE of 0.33 mm/s.

The most challenging systems to study are those characterized by strong static correlations. While individual hybrid, double-hybrid, or even pure functionals can provide reasonable results, they exceed the capabilities of DFT and the results often hinge on chance, showing the complexities involved in accurately modeling such systems. In contrast, advanced multi-configurational methods, such as DMRG and sCI, offer significant advantages as demonstrated in our previous research on Mössbauer parameters [2,3].

The pursuit of identifying optimal functionals for Mössbauer NQS calculations of iron compounds remains an active area of research, which, however, has not garnered significant attention in the realm of functional development. Future advancements in the development of new functionals may incorporate $[eqn]$ data into their evaluation, whereas the contact density data are still not enough. This approach may represent a promising avenue for enhancing the accuracy of DFT calculations.