First-principles paramagnetic NMR of a challenging Fe(V) bis(imido) complex: a case for novel density functionals beyond the zero-sum game

Abstract

We investigate computationally the hyperfine couplings (HFCs) and the consequent paramagnetic nuclear magnetic resonance (pNMR) chemical shifts of a recently synthesised doublet Fe(V) bis(imido) complex. Using conventional global hybrid density-functional theory (DFT) methods with varying exact-exchange admixture, a significant spin contamination problem is observed, leading to a massive spin-density spill-over to the strongly bound imido ligands and to the BH₂ group of the carbene framework. As a result, the computed paramagnetic NMR shifts, which are based on a combination of ab initio calculated g-tensor, DFT-calculated orbital shielding and DFT-based HFCs, disagree strongly with the available experimental ¹H NMR chemical shifts and predict unrealistic ¹³C shifts in the spill-over region. While semi-local functionals like PBE do not suffer from the spin-contamination problem, their inherent delocalisation errors also lead to a too spread-out spin-density distribution. We demonstrate that, by applying novel local hybrid and range-separated local hybrid functionals with correction terms for strong-correlation and/or delocalisation errors to the HFC computations, the spin contamination problem is significantly reduced, while keeping delocalisation errors small. This results in more realistic pNMR shifts obtained for this system, also when compared to data obtained using correlated ab initio HFCs. Problems with the use of a global hybrid functional are already observed at the level of the underlying optimised structure, and employing a semi-local functional in structure optimisation is found to improve the situation. With that, the combined ab initio/DFT method with the latest (range-separated) local hybrid functionals is seen to reasonably reproduce the ¹H experiment and enable plausible predictions for the hitherto unreported ¹³C shifts in this challenging iron complex.

---

1 Introduction

Nuclear magnetic resonance (NMR) spectroscopy is a widely used chemical analysis tool, for example in identifying novel chemical species. Besides the usual practice of applying NMR to closed-shell molecules, the method can also be employed for paramagnetic species (pNMR) having unpaired electrons.¹ Contemporary pNMR research is almost without exception supported by computational electronic-structure studies, which help in assigning, predicting and analysing the observed pNMR signals.² A challenge in calculating pNMR spectral parameters is posed by the fact that the underlying physics is more complicated than in conventional NMR, necessitating in most cases the calculation of several distinct magnetic property tensors, which jointly determine, e.g., the pNMR chemical shift values.

A routine modus operandi was formulated for practical pNMR shift calculations a few years ago^3–5 based on a modern quantum-chemical implementation^3,6–9 of the classic Kurland–McGarvey chemical shift theory.¹⁰ This method is based on parameterising the ground spin multiplet in terms of the electron paramagnetic resonance (EPR) tensors,¹¹ the g-tensor for the Zeeman interaction of the unpaired electron(s) with the magnetic field, the zero-field splitting (ZFS) in systems with more than one unpaired electrons and the hyperfine coupling (HFC) tensor for the interaction of the spin of the unpaired electron(s) with the magnetic nucleus, jointly constituting the hyperfine shielding tensor. To this one should add the orbital shielding tensor involving the corresponding orbital hyperfine and Zeeman interactions, analogous to the theory of the chemical shift in closed-shell systems.¹² While this method is typically used in the limit of small spin–orbit coupling, a more comprehensive theory has been put forward, employing generalised EPR tensors and in which the orbital contribution appears on equal footing with the hyperfine part.¹³ Furthermore, a formulation directly in terms of the ground and thermally available excited electronic states and, hence, avoiding altogether the use of EPR parameters, has been published.¹⁴ These advances are not yet in widespread use in pNMR calculations.

In the standard EPR parameter-based pNMR shift method the property tensors involved can be calculated at different levels of theory according to the predictive power and the computational cost of the different electronic-structure methods in each case, in a mixed approach combining ab initio and density-functional theory (DFT)-based calculations.⁴ In the case of the g- and ZFS tensors, often multireference ab initio approaches have been employed,^{3–5,15–21} either at the state-average complete-active space self-consistent field (CASSCF) level²² or also employing N-electron valence-state perturbation theory (NEVPT2^23–25) on top of the CASSCF optimisation, to capture both static and dynamic correlation effects reasonably at the same time. In these calculations, one-component state-average CASSCF has been performed to construct a basis in which the spin–orbit coupling (SOC) Hamiltonian has been diagonalised to gain access to the g- and ZFS tensors.^26,27

Computations of the HFC interaction can be deemed more problematic, as the genuine multireference ab initio approaches are computationally involved for typical system sizes of pNMR interest. The findings for small main-group radicals were discussed in ref. 28. In large (metal-containing) systems it is difficult to extend the active orbital space all the way down in energy to orbitals that continue to be relevant for gaining HFC of predictive quality. This requires incorporating hand-picked, deep orbitals to sufficiently capture spin polarization in a system-specific manner.²⁹ Recent electron-correlation methodology in the form of domain-based local pair natural orbital coupled-cluster method with single and double excitations (DLPNO-CCSD^30,31) is capable of treating fairly large molecular systems, as well as correlating a large number of orbitals, and has been applied to obtain pNMR shifts.^19,20 It remains a single-reference method, however.

This situation has rendered hybrid DFT the practical method of choice for most HFC computations in the pNMR context. While being applicable to large, metal-ion containing systems and certainly capable of producing useful data for pNMR signal assignment and analysis, DFT-based HFC tensors and the corresponding pNMR shifts show a substantial dependence on the choice of the exchange–correlation functional.^18,20 In standard global hybrid functionals (GHs), the amount of exact-exchange admixture (EXX) is fixed everywhere in the molecule, and the obtained HFC results vary significantly and somewhat unsystematically upon the choice of EXX. Large EXX admixtures are known to be required to describe core–shell spin polarisation, e.g., for metal HFCs.^32,33 However, GHs tend to suffer more from spin contamination than semi-local functionals in certain electronic-structure situations, increasingly so with larger EXX admixtures, due to the overestimation of valence-shell spin polarisation.^32,33 This difficult balance constitutes the “spin-polarisation/spin-contamination dilemma”,^32–34 which may be viewed as a special case of the well-known “zero-sum game” of DFT functionals between the proper description of static correlation effects and a minimisation of delocalisation errors.^35–38

In local hybrid functionals (LHs), EXX admixture is made position-dependent, such that a lower EXX admixture can be used in the valence space, and a larger one in the core region.³⁹ Modern range-separated local hybrids (RSLHs) feature additionally a variation of EXX in different regions of the interelectronic distance.^40,41 The first evaluations of LHs for transition-metal HFCs indicated improvement over GHs but could not yet fully escape the above-mentioned dilemma.^34,42 The introduction of strong-correlation (sc) corrections into the local-mixing functions (LMFs), which determine the position dependence of EXX admixture in LHs,^38,43–45 enabled a more notable departure from the zero-sum game. In this case of scLHs we have to currently distinguish between “damped” sc-factors that are designed to retain the good performance of the underlying LH for weakly correlated systems, and “undamped” sc-factors that provide more notable corrections but also may deteriorate the situation for certain weakly correlated systems.

Most recently, the idea of sc-corrections has been transferred to RSLHs, where the long-range EXX admixture provides the correct asymptotic behavior of the exchange–correlation potential.^38,46 sc-Correction terms in scRSLHs also affect the range separation, but only in such real-space regions where appreciable static correlations have been detected. Notably, the LMFs of some of these scRSLHs have been extended further by including delocalization-error corrections (DEC term) in certain abnormal open-shell spatial regions.^38,46 Indeed, such scRSLHs with DEC terms currently provide the most practical way of abandoning the zero-sum game of DFT functionals^38,46 (we may also view the deep-neural-network functional DM21 as an scRSLH⁴⁷). In addition to improving upon spin-restricted bond-dissociation curves, scLHs and scRSLHs have very recently been investigated for the metal HFCs of a selected set of manganese complexes with variable challenges regarding spin contamination.⁴⁸ Both sc- and DEC-corrections have been found to reduce spin contamination in critical complexes such as MnO₃ or [Mn(CN)₅NO]²⁻, while still providing a good description of core–shell spin polarisation. Notably, only “undamped” sc-corrections in scLHs provided sufficiently large reductions of spin contamination (and associated improvements in dipolar HFCs). In fact, DEC corrections in scRSLHs were found to be even more effective in this context. Interestingly, the improvements by sc-corrections arise from a local reduction of EXX admixture on certain ligand atoms, while the DEC terms in the LMF actually enhance locally the EXX admixture in the metal valence shell. Both effects are operative for the most successful scRSLHs.⁴⁸ Recently, scLHs have also been implemented for the calculation of NMR shifts in diamagnetic systems⁴⁹ (and for magnetisabilities⁵⁰) and were evaluated both for main-group and transition-metal shifts. Here it was also found that “undamped” scLHs outperformed their “damped” counterparts for the systems with the largest static correlations, in spite of some deterioration in other cases.⁴⁹

While the first evaluations of scLHs and scRSLHs for HFCs focused on the transition-metal centre, here we will evaluate them for the ligand ¹H and ¹³C HFCs and the related pNMR shifts of a doublet (S = 1/2) formal Fe(V) bis(imido) complex (Fig. 1), which was synthesised as an intermediate to obtaining a novel closed-shell formal Fe(VI) compound.⁵¹ In this context, the paramagnetic ¹H NMR chemical shifts were recorded from the Fe(V) species featuring distorted tetrahedral coordination of Fe(V) with 3d³ configuration, two strongly coupled imido ligands, and a chelate bis-N-heterocyclic carbene ligand framework. From the point-of-view of computational pNMR, this is a challenging molecule due to its large, 100-atom size. In the present work we tackle the pNMR shift modelling of this system by applying our mixed ab initio/DFT workflow⁴ and focus particularly on the methods for calculating the ¹H and ¹³C HFC tensors in this molecule. Semi-local and standard GH functionals, DLPNO-CCSD, as well as novel scLHs and scRSLHs are applied. It turns out that standard GHs (as well as range-separated hybrids and LHs or RSLHs without proper corrections) lead to massive spin-contamination problems due to the strongly coupled ligands, rendering the calculated pNMR shifts entirely unreliable. Semi-local functionals do not produce large spin contamination but suffer from inherent delocalisation errors. These “zero-sum-game” problems are remedied to varying extents by employing either “undamped” scLHs or scRSLHs with DEC terms for the HFC calculations. It is also found that due to possible spin contamination the selection of the DFT functional is critical also in providing the underlying structure used in the property calculations, and using a semi-local functional rather than a GH in the optimisation helps reduce further the unphysical spill-over of spin density and the concomitant deterioration of the pNMR shift predictions.

Fig. 1 Schematic illustration of the Fe bisimido complex including labeling of the ¹H NMR centres in red and the additional ¹³C centres in black letters. Iron atom in yellow, boron in purple, nitrogen in blue, carbon in orange and hydrogen in white.

2 Methods

2.1 Theory

We employ the standard EPR parameter-based theory, in which the pNMR shielding tensor of nucleus K is calculated for the doublet (S = 1/2) system as^3,6where σ^orb_K is the orbital shielding tensor analogous to that found in closed-shell systems.⁵² g and A_K in the second, hyperfine shielding term are the g-tensor of the complex and the HFC tensor of nucleus K therein. The rest of the symbols in eqn (1) have their usual meanings. The shielding constant is defined from the tensor components as

σ_K,iso = (σ_K,xx + σ_K,yy + σ_K,zz)/3, (2)

and the observable chemical shift in isotropic liquid-phase experiments as

δ_K = σ_ref − σ_K,iso, (3)

where σ_ref is the shielding constant of nucleus K in a diamagnetic reference molecule, in this context tetramethylsilane (TMS).

2.2 Structure, orbital shieldings and g-tensor

Structure optimisation was carried out for the open-shell Fe(V) and closed-shell Fe(VI) complexes with DFT using the global hybrid functional PBE0⁵³ and, in the case of the open-shell complex, also the generalised-gradient approximation (GGA) functional PBE,⁵⁴ applying in both cases DFT D3(BJ) dispersion corrections,^55,56 using the Turbomole programme.^57–59 For the iron centre we applied the Stuttgart-type scalar relativistic effective core potential ECP10MDF,⁶⁰ along with the corresponding 6s5p3d2f1g/8s7p6d2f1g valence basis set.⁶¹ def2-TZVP basis sets⁶² were employed for the ligand atoms. The unrestricted Kohn–Sham (UKS) formalism was used for the Fe(V) complex. The resulting molecular structures are given in the SI.

The ¹H and ¹³C orbital shielding tensors σ^orb were obtained using restricted and unrestricted DFT for the closed- and open-shell complexes, respectively, employing the PBE0 functional and the all-electron DKH-def2-TZVP⁶³ basis in the Gaussian 16 programme.⁶⁴ Similar methodology was employed for the TMS reference, obtaining σ_¹H,ref = 31.8 and σ_¹³C,ref = 188.4 ppm. For comparison, the orbital shielding constants of the Fe(V) system were also calculated with the PBE functional.

The g-tensor was calculated with the ORCA programme⁶⁵ using strongly contracted DLPNO-NEVPT2⁶⁶ based on state-average CASSCF for up to 7 correlated electrons in 12 orbitals, the latter including the five 3d metal orbitals “doubled” with the next five metal d-orbitals, and two metal–ligand hybrid orbitals below the 3d shell. Calculation of g was performed in the basis of one-component states built from orbitals obtained at scalar relativistic second-order Douglas–Kroll–Hess (DKH2) level of relativistic theory.^67,68 The locally dense DKH-def2-TZVP/SVP basis set was employed for the metal and directly bonded atoms/other atoms.

2.3 Hyperfine coupling

In order to permit spin polarisation of core–shell electrons, an important factor in core properties such as HFCs, the UKS formalism was used in the DFT computations. Our basic HFC calculations were carried out at the scalar-relativistic DKH2 level with the ORCA software, using the PBE GGA and the derived GHs PBE0 and PBE0-40, increasing the exact-exchange admixture from 0% to 25% to 40%, respectively. The range-separated CAM-B3LYP functional⁶⁹ was also employed. Ab initio DLPNO-CCSD calculations were carried out for comparison, employing the HFC2 truncation level and frozen core with 10 and 2 electrons for Fe and second-row atoms, respectively.

Subsequent calculations with LHs and scLHs from the Berlin group, LH20t,⁷⁰ scLH22t,⁴⁴ scLH22ta,⁴⁴ scLH23t-mBR,⁴⁵ scLH23t-mBR-P,⁴⁵ and scLH21ct-SVWN-m,⁴³ were done at the 1-component X2C level implemented in Turbomole. Additionally, calculations with the RSLH ωLH22t,⁴¹ a modification ωLH23td with only DEC corrections,⁴⁶ scRSLHs without DEC terms (ωLH23tB, ωLH23tE, ωLH23tP), and the corresponding scRSLHs with DEC terms (ωLH23tdE, ωLH23tdP, and ωLH23tdB) were done with the same implementation.⁴⁶ We will in the following denote the sc-corrected LMF in these functionals as “tq-LMF”, the DEC-corrected LMF as “td-LMF”, and the sc- and DEC-corrected LMF as “tdq-LMF”. For comparison, the previously derived sc-corrected functional KP16 of Kong and Proynov,⁷¹ which had been an inspiration for the sc-corrections of scLHs,^43,44 has also been evaluated.

The HFC tensor calculations were performed employing the DKH-def2-TZVP basis set in their completely uncontracted form, apart from the case of DLPNO-CCSD, for which the contracted basis set was used due to computational limitations. Some uncorrected LHs, particularly LH20t, were overall fairly successfully used for investigating the pNMR shifts in a range of metal complexes in ref. 72, in both the mixed ab initio/DFT context and in an all-DFT implementation.

3 Results

3.1 NMR shifts of the closed-shell complex

The calculated ¹H chemical shifts of the closed-shell formal Fe(VI) bis(imido) d² complex are listed in Table 1. The computational results at the DFT/PBE0 level are in line with measurements and suggest the assignment of the signals b/e and c/f (Fig. 1) which were not resolved experimentally [see the “Exp. (alt)” row of the Table for the alternative experimental assignment]. In the assignment implied by the computations, H-b is to be associated to the higher measured chemical shift, 6.81 ppm, whereas the 6.67-ppm signal belongs to H-e. Similarly, the 1.63-ppm and 1.38-ppm signals appear to belong to centres c and f, respectively. In contrast, fixing the assignment of the signals a/d is not aided by the computations as the experimental shifts are very close to each other.

Table 1 Computational ¹H NMR chemical shifts (in ppm) for the closed-shell Fe(VI) bis(imido) complex, obtained at the DKH2 PBE0 level with different environmental models, using a PBE0-optimised structure

	H-a	H-b	H-c	H-i	H-d	H-e	H-f	H-g	H-h
a Alternative assignment of the experimental signals.b Molecule + counterion at optimal position in vacuo.c Molecule + counterion at optimal position, placed in the implicit C-PCM solvation model (ε = 8.93 appropriate to dichloromethane solution was applied).
Exp.^51	2.23	6.81	1.63	4.01	2.24	6.67	1.38	7.25	7.58
Exp. (alt)^51 ^a	2.24	6.67	1.38		2.23	6.81	1.63
Comp. in vacuo	2.56	7.27	1.75	4.85	2.39	7.00	1.37	7.46	7.99
Comp. with BF4^− counterion^b	2.50	7.30	1.99	4.79	2.79	6.91	1.23	7.45	7.83
Comp. with BF4^− and C-PCM^c	2.56	7.40	1.91	4.69	2.65	6.96	1.28	7.68	7.95

---

Inclusion of the explicit BF₄⁻ counterion, as well as (additionally) implicit solvation via the conductor-like polarizable continuum model⁷³ were found to only have a minor impact on the results. In particular, the signal assignment is not altered by employing these models. Predictions for the shifts of the hitherto unmeasured ¹³C signals are included in the SI.

3.2 Spin contamination in the open-shell compound

Ab initio ligand-field theory analysis⁷⁴ was carried out using the ORCA code on the distorted tetrahedral structure of the doublet Fe(V) complex with d³ configuration. The results based on the CASSCF(3,5) wavefunction, placed in the SI, indicate that the structural distortion is sufficient to cause sufficient splitting between the two lowest, e-symmetry orbitals, to remove any orbital degeneracy issues. With the underlying structure optimisation performed using the PBE functional, the splitting between the two lowest orbitals as well as that to the upper three orbitals become even more pronounced than based on the PBE0 structure.

The 〈S²〉 expectation values (Fig. 2 and Table 2) obtained with standard GHs, here represented by the PBE-based series of functionals, indicate spin contamination that increases rapidly with EXX admixture. The underlying PBE GGA is relatively free from spin contamination, with a value close to the formal 〈S²〉 expectation value 0.75 appropriate for a doublet state. Similar to the GHs, the CAM-B3LYP range-separated hybrid also features significant spin contamination (1.68 at the PBE0 structure), as do the uncorrected LH20t LH and the ωLH22t RSLH.

Fig. 2 Calculated expectation values 〈S²〉 for the different DFT functionals used in the HFC calculations. The horizontal line indicates the nominal value of 0.75 for a doublet state.

Table 2 Expectation value 〈S²〉 in the DFT HFC calculations with different exchange–correlation functionals

Functional	〈S^2〉 value
	PBE structure	PBE0 structure
a Second-order, scalar relativistic Douglas–Kroll–Hess level of theory.b Scalar relativistic X2C theory.
PBE^a	0.77	0.78
PBE0^a	1.19	1.50
PBE0-40^a	1.97	2.21
LH20t^a	0.93	1.14
scLH22t^b	0.90	1.09
scLH22ta^b	0.77	0.79
scLH23t-mBR^b	0.89	1.05
scLH23t-mBR-P^b	0.90	1.08
scLH21ct-SVWN-m^b	0.77	0.78
KP16^b	0.77	0.77
ωLH22t^b	1.16	1.50
ωLH23tB^b	1.10	1.46
ωLH23tE^b	1.04	1.42
ωLH23tP^b	1.02	1.40
ωLH23td^b	0.79	1.12
ωLH23tdB^b	0.78	0.83
ωLH23tdE^b	0.78	0.83
ωLH23tdP^b	0.78	0.81

---

Moving to the scLHs derived from LH20t, we see that the “damped” versions scLH22t, scLH23t-mBR, and scLH23t-mBR-P reduce 〈S²〉 somewhat compared to LH20t, but significant spin contamination remains. Only the “undamped” scLH22ta and the also undamped, simpler scLH21ct-SVWN-m correct 〈S²〉 almost completely, as does the KP16 functional, which features similar sc-factors.

The scRSLHs are all derived from ωLH22t. Those that exhibit only sc-corrections and no DEC term (ωLH23tB, ωLH23tE, ωLH23tP) reduce the spin contamination only very little compared to ωLH22t. Interestingly, the ωLH23td RSLH with DEC term but without sc-corrections performs already much better at the PBE structure, while some spin contamination remains at the PBE0 structure. Finally, the scRSLHs with DEC term (ωLH23tdB, ωLH23tdE, ωLH23tdP) approach closely the low 〈S²〉 values of the “undamped” scLHs while not reaching them completely. These observations are consistent with recent calculations on metal HFCs and spin-symmetry breaking in critical Mn complexes.⁴⁸ In that work closer analyses showed that the sc-corrections reduce valence-shell spin polarisation by lowering EXX admixture on the relevant, strongly coupled ligand atoms, and that the “undamped” sc-corrections are more effective than the “damped” ones. The current scRSLHs all have damped sc terms, explaining why they seem relatively ineffective. On the other hand, the DEC term enhances EXX admixture locally in the metal valence shell and thereby reduces unphysical spill-over of spin density onto the strongly coupled ligand atoms.⁴⁸ We can confirm here that the latter effect alone is more effective than the former when damping is employed. Overall the combination of having both sc- and DEC-corrections in the LMF has a rather large effect for the scRSLHs, in reducing spin contamination.

Notably, the 〈S²〉 values are generally reduced when using the PBE-optimised rather than the PBE0-optimised structure (Table 2), showing that the shorter metal–ligand distances obtained with the semi-local functional (Table 3) reduce the overall spin-contamination problem. The larger metal–ligand distances with PBE0 are also a consequence of spin contamination with the GH (even larger effects have been found for Fe–N distance and vibrational frequency of a high-oxidation-state iron nitrido complex⁷⁵), which causes an exaggerated population of metal–ligand antibonding orbitals, thereby lengthening the bonds to the metal centre. The latter effect seems to be larger for the Fe–C bonds to the N-heterocyclic carbene ligands than for the imido Fe–N bonds, likely reflecting the steeper potential curve for the latter, formal double bonds.

Table 3 Interatomic Fe(V)–X distances (in Å) in the coordination shell of the iron centre using PBE0 or PBE for structure optimisation

Atom number (element)	Distance from Fe
	PBE structure	PBE0 structure
2 (C)	1.90	1.95
18 (C)	1.90	1.95
19 (N)	1.68	1.70
20 (N)	1.68	1.70

---

The spin densities obtained with example functionals at PBE0 and PBE structures are illustrated in Fig. 3. We show PBE0 as a conventional GH, PBE as an example of a semi-local (GGA) functional and ωLH23tdP as an example of the most advanced RSLHs, with both strong-correlation and delocalisation-error corrections. For both structures a significant spin transfer to the strongly bound imido ligands (and to some extent to the carbene ligands) is seen with the GH, more pronouncedly at the PBE0 structure. The PBE structure reduces this spill-over, and either using PBE or ωLH23tdP to compute the spin density reduces the transfer and the concomitant spin polarisation within the ligand even further. Similar observations have been made recently for a series of Mn complexes.⁴⁸ We note, however, that the PBE spin density is more delocalized than that of ωLH23tdP, both near the BH₂ handle and on some of the mesityl carbon atoms of the bis-carbene ligand. This arises from the spin-polarisation contribution, while the SOMO is in fact less delocalised with PBE than with PBE0.

Fig. 3 Calculated spin density isosurfaces (positive red/negative blue isovalues ±0.001 a.u.) as obtained in the (a) DKH2-PBE0, (b) DKH2-PBE and (c) 1c-X2C-ωLH23tdP calculations at the PBE0 structure, as well as (d) DKH-PBE0, (e) DKH2-PBE and (f) 1c-X2C-ωLH23tdP calculations at the PBE structure. Panels from left to right indicate total spin density, the density corresponding to the singly occupied MO and the difference of the two, representing the spin polarisation contribution.

3.3 g-Tensor calculations

The g-tensor results obtained at different multireference ab initio levels are collected in Table 4. Different aspects of the performed calculations affect the results. The active space initially, in the CASSCF(3,5) calculations, comprised the five metal 3d orbitals. They were supplemented by two occupied metal–ligand bonding orbitals in the CASSCF(7,7) wave functions and, then, by five empty metal 4d-orbitals, in CASSCF(7,12). The description of static correlation is expected to improve in this series of CASSCF wave functions, whereas the incorporation of NEVPT2 introduces an enhanced account of the dynamic correlation on top of each of these CASSCF levels. Finally, the calculations were carried out at both the PBE0 and PBE-optimised structures.

Table 4 g-Tensor eigenvalues and isotropic g-factors from the CASSCF and DLPNO-NEVPT2 calculations

Structure	Method	Eigenvalues
		g11	g22	g33	giso
a Full NEVPT2 calculation instead of DLPNO-NEVPT2.
PBE0	CASSCF(3,5)	1.93	1.97	1.99	1.97
	CASSCF(3,5)+NEVPT2	2.00	2.00	2.00	2.00
	CASSCF(3,5)+NEVPT2^a	2.00	2.00	2.00	2.00
	CASSCF(7,7)	1.97	1.99	1.99	1.98
	CASSCF(7,7)+NEVPT2	1.97	1.99	2.00	1.98
	CASSCF(7,12)	1.96	2.00	2.00	1.99
	CASSCF(7,12)+NEVPT2	1.97	2.00	2.00	1.99
PBE	CASSCF(3,5)	1.93	1.98	1.99	1.97
	CASSCF(3,5)+NEVPT2	2.00	2.00	2.00	2.00
	CASSCF(7,7)	1.97	1.99	2.00	1.99
	CASSCF(7,7)+NEVPT2	1.95	1.99	2.00	1.98
	CASSCF(7,12)	1.96	1.99	2.00	1.98
	CASSCF(7,12)+NEVPT2	1.95	1.98	1.99	1.98
	Experiment^51	1.910	1.995	2.010	1.97

---

At the CASSCF level the g-tensor is close to being cylindrically symmetric, with one unique (g_‖) and two fairly similar, further eigenvalues (g_⊥). The g-anisotropy, g_‖−g_⊥ is consistently negative when the PBE structure is used. The two largest CASSCF wave functions provide fairly similar g-tensors at both PBE0 and PBE structures, implying relative convergence with respect to the choice of the active space choice. Hence, the active space of the CASSCF(7,12) calculation should suffice.

Incorporating the dynamic correlation corrections at the NEVPT2 level changes the g-tensor qualitatively when the smallest, 5-orbital active space is used, whereas only a minor modification due to switching on NEVPT2 is seen with the 7- or 12-orbital active spaces. This, again, points at sufficient convergence of the 12-orbital active space. The overall impression is, hence, that the choice of active space is a more significant factor than incorporation of the NEVPT2 correction, for the g-tensor of this Fe(V) complex. As a side note, we find from Table 4 that very similar data for the g-tensor are obtained from the approximate DLPNO-NEVPT2 calculations as from standard NEVPT2, when tested with the CASSCF(3,5) wave function at the PBE0 structure. This implies adequacy of the DLPNO procedure, which is necessary with the larger active spaces for the present system.

Comparison with the X-band EPR experiment carried out in toluene at 20 K⁵¹ is only partially justified, as the present g-tensor calculations were carried out for the complex in vacuo, i.e., omitting any matrix or intermolecular interaction effects. The experimental tensor is, similarly as in the computational results very nearly cylindrically symmetric. The measured g_‖–g_⊥ equals −0.09, i.e., a negative value, as in our supposedly best calculation, CASSCF(7,12)+NEVPT2 at the PBE structure. This wave function seems to, however, underestimate the anisotropy as compared to the experiment. The experimental isotropic g-value (1.97) is well-reproduced by the CASSCF(7,12)+NEVPT2 model, particularly when the PBE structure is employed. In contrast, the experimental range of the g-tensor eigenvalues (1.91…2.01) is underestimated in our best calculation (1.95…1.99). This may be caused in part by the omission of intramolecular molecular dynamics effects and (most presumably) the fact that we treated an isolated molecule in these calculations.

3.4 ¹H pNMR shifts with different methods for the hyperfine coupling tensors

3.4.1 HFCs from conventional functionals. Table 5 shows the ¹H chemical shifts computed for the paramagnetic Fe(V) bis(imido) complex based on HFCs obtained with conventional DFT functionals, i.e., the PBE GGA and the related GHs PBE0 and PBE0-40, with 25% and 40% EXX admixture, respectively. For the PBE0-based structure, also the CAM-B3LYP range-separated hybrid was evaluated.

Table 5 Computational ¹H NMR chemical shifts (in ppm) for the open-shell Fe(V) bis(imido) complex, obtained using hyperfine coupling tensors from PBE and PBE0 functionals, the latter with both 25% (PBE0) and 40% (PBE0-40) exact-exchange admixtures. Also CAM-B3LYP was employed. Calculations at the DKH2 level of relativity, unless otherwise noted. Orbital shieldings were obtained using the PBE0 functional and the g-tensor is adopted from CASSCF(7,12)+NEVPT2 calculations

Structure	HFC method	H-a	H-b	H-c	H-i	H-d	H-e	H-f	H-g	H-h
a Ref. 51.b Ref. 51. Alternative assignment.c Fully relativistic matrix-Dirac-Kohn–Sham calculation.
	Exp.^a	−1.74		4.50	73.16	2.04	6.69	1.44	−7.35	−0.84
	Exp. (alt)^a^b								−0.84	−7.35
PBE0	PBE	−39.03	31.30	−14.77	176.47	2.03	7.16	0.43	−14.17	−2.94
	PBE0	−281.87	−102.52	−199.02	149.86	2.68	7.74	−2.81	−16.04	−15.01
	PBE0 (mDKS^c)	−281.26	−102.09	−198.12	150.29	2.67	7.72	−2.79	−16.29	−14.95
	PBE0-40	−433.76	−245.93	−331.03	147.14	2.30	7.43	−3.80	−9.68	−14.95
	CAM-B3LYP	−337.54	−132.58	−252.70	130.30	2.59	7.42	−2.84	−4.47	−10.30
PBE	PBE	−16.63	40.47	−1.25	189.52	2.38	7.54	0.68	−12.97	−1.28
	PBE0	−195.10	−58.81	−135.49	151.70	2.73	7.84	−2.60	−18.01	−12.22
	PBE0-40	−380.93	−219.30	−293.16	163.50	2.01	7.18	−3.55	−13.07	−15.98

---

For the centres H-d, H-e and H-f in the trimethylphenyl ligands, the computations with the standard functionals are in fairly good agreement with experiment. The correct order of the signals is reproduced and semiquantitative agreement with experiment is obtained for the magnitude of the H-d and H-e shifts for all the functionals. H-f shows a somewhat larger dependence on the chosen functional and the best results are obtained using PBE for the HFC tensors, regardless of whether the PBE or PBE0 structure is used.

Moving still closer to the metal centre, the computed H-g and H-h signals of the imidazole fragment reside somewhat further from the experiment, where ambiguity has remained on the mutual assignment. Data based on the PBE functional suggest assigning the −7.35 ppm and −0.84 ppm shifts to H-g and H-h, respectively, whereas the other functionals in Table 5 either suggest the opposite or remain undecided. The H-i signal of the BH₂ group has experimentally a large positive shift, which is overestimated by 100% or more by all the more traditional functionals. Interestingly, in this case the deviations from experiment are even larger with PBE than with the GHs or CAM-B3LYP. This may suggest that delocalisation errors with PBE increase the shifts for these two protons.

In the imido ligand the signals of the H-a and H-c sites have been recorded, whereas H-b has not been observed experimentally. The conventional functionals of Table 5 completely fail for H-a and H-c, by producing chemical shifts one to two orders of magnitude larger than the experimental data. This holds even for the PBE GGA, while EXX admixture in the GHs or CAM-B3LYP lead to dramatically large deviations, increasingly so with larger EXX admixture. The large negative shifts for these three nuclei with GHs or CAM-B3LYP reflect very large negative (i.e., β-spin) spin density at these imido-ligand nuclei, consistent with Fig. 3: large, both negative and positive spin density lobes are seen at this part of the molecule in the PBE0-level calculations depicted in the right-most panels of Fig. 3(a) and (d). As the data in these panels stem from the total spin density from which the density of the singly occupied molecular orbital (SOMO, representing spin delocalisation) has been subtracted, they indicate the spin polarisation contribution. Comparison with the much less pronounced spin polarisation in the corresponding, right-most panel of Fig. 3(b), obtained using the PBE functional, points at exaggerated valence-shell spin polarisation by the GH functionals, in this ‘hard’ part of the molecule, as the culprit. On the other hand, considering the relative similarity of the experimental values for H-a and H-c with those for the diamagnetic complex (Table 1) suggests that relatively little spin density should be present.

The apparent failure of the traditional hybrid functionals is clearly correlated with the spin contamination data presented above (Table 2): 〈S²〉 is 1.50 (1.19) and 2.21 (1.97) for PBE0 and PBE0-40, respectively, at the PBE0 (PBE) structure, deviating significantly from the nominal value 0.75. PBE does not exhibit appreciable spin contamination yet shows smaller but still significant deviations of the shifts for these three nuclei from experiment.

As a side remark, Table 5 includes also a comparison of results obtained using fully relativistic HFCs calculated at the mDKS-level using the ReSpect code⁷⁶ (the present HFC basis set has been used in a fully uncontracted manner in the large component) and the present main level of theory, scalar relativistic DKH2, performed at the PBE0-optimised structure. There is very little difference in the obtained ¹H chemical shifts, a result that eliminates speculation of the role of the omitted spin–orbit interaction in the scalar relativistic calculations of the HFCs.

3.4.2 HFCs from LHs, RSLHs, scLHs, scRSLHs, KP16, and DLPNO-CCSD.
3.4.2.1 The LH20t local hybrid functional and its strong-correlation corrected variants. Table 6 compares the ¹H shifts obtained using the more advanced functionals, as well as DLPNO-CCSD, with the experimental data. The mean average errors (MAEs) at these levels are visualised in Fig. 4, which also includes a comparison with PBE, PBE0, and PBE0-40. Fig. S1 and S2 in the SI illustrate the signed deviations of the calculated ¹H shifts from experiment.

Table 6 Computational ¹H NMR chemical shifts (in ppm) for the open-shell Fe(V) bis(imido) complex, using hyperfine coupling tensors obtained at the 1c-X2C level from LHs, RSLHs, scLHs, scRSLHs, and KP16, as well as at 1c-DKH2 level with DLPNO-CCSD. Abbreviations for RSLHs and scRSLHs start with an “ω ”, those of scLHs with an “sc”

Struct.^a	HFC method	H-a	H-b	H-c	H-i	MAE^b	H-d	H-e	H-f	H-g	H-h	MAE^c
a Underlying level of structure optimisation.b Mean absolute error (in ppm) with respect to the experimental shifts in the calculations for protons H-a, H-c, and H-i.c Mean absolute error (in ppm) with respect to the experimental shifts in the calculations for protons H-d, H-e, H-f, H-g, and H-h.d Ref. 51.e Ref. 51. Alternative assignment.
PBE0	Exp.^d	−1.74		4.50	73.16		2.04	6.69	1.44	−7.35	−0.84
	Exp. (alt)^d^e									−0.84	−7.35
	LH20t	−170.51	−23.61	−111.07	112.85	109.17	2.33	7.28	−0.50	−10.94	−10.68	3.59
	scLH22t	−161.84	−17.62	−103.59	118.34	105.62	2.31	7.30	−0.62	−11.45	−10.34	3.64
	scLH22ta	−47.91	17.35	−23.12	132.36	45.49	2.15	7.22	1.95	−8.93	−2.60	1.24
	scLH23t-mBR	−144.83	−11.49	−90.62	108.13	92.22	2.28	7.23	−0.13	−11.86	−9.23	3.38
	scLH23t-mBR-P	−154.52	−15.42	−98.33	110.36	98.76	2.30	7.24	−0.37	−11.34	−9.88	3.46
	scLH21ct-SVWN-m	−39.52	20.65	−16.29	168.57	52.49	2.20	7.33	0.49	−8.81	−1.65	1.14
	KP16	−31.12	16.05	−9.43	165.47	46.36	−0.08	7.35	−0.47	3.89	−0.83	3.52
	ωLH22t	−263.30	−74.83	−192.32	109.53	166.08	2.33	7.11	−1.64	−2.85	−10.27	3.88
	ωLH23tB	−259.49	−67.37	−188.34	115.88	165.59	2.33	7.16	−1.95	−3.64	−10.83	3.91
	ωLH23tE	−242.02	−60.61	−174.55	105.71	151.78	2.34	7.11	−1.32	−3.82	−9.99	3.57
	ωLH23tP	−233.85	−53.71	−167.45	105.14	146.50	2.33	7.09	−1.43	−3.51	−9.88	3.62
	ωLH23td	−118.93	42.74	−68.89	63.28	67.98	2.03	6.88	−2.46	12.83	−0.86	5.19
	ωLH23tdB	−44.25	11.22	−21.87	51.27	31.41	2.08	6.88	−0.92	9.83	5.13	4.81
	ωLH23tdE	−48.39	11.63	−23.52	56.33	31.66	2.05	6.86	−0.78	9.05	4.40	4.47
	ωLH23tdP	−33.85	15.30	−12.89	52.45	24.56	2.09	6.91	−0.37	4.98	2.65	3.24
	DLPNO-CCSD	−96.68	−36.18	−74.83	49.63	67.10	2.07	6.73	0.62	−1.44	0.32	1.47
PBE	LH20t	−102.05	0.52	−62.91	110.82	68.46	2.45	7.64	−2.19	−9.97	−7.23	2.80
	scLH22t	−93.83	5.57	−56.07	116.85	65.45	2.43	7.64	−1.99	−9.86	−6.73	2.63
	scLH22ta	−23.76	27.11	−6.98	140.62	33.65	2.51	7.57	1.01	−7.46	−0.93	0.39
	scLH23t-mBR	−85.03	8.27	−48.95	108.38	57.32	2.43	7.57	−1.42	−10.21	−5.94	2.42
	scLH23t-mBR-P	−90.43	6.30	−53.48	109.42	60.98	2.43	7.59	−1.74	−9.98	−6.44	2.54
	scLH21ct-SVWN-m	−18.87	29.30	−3.33	181.19	44.33	2.57	7.72	0.75	−7.66	−0.19	0.64
	KP16	−19.28	21.24	−0.90	182.35	44.04	3.27	6.72	0.68	4.62	0.09	2.98
	ωLH22t	−173.85	−39.80	−125.91	108.66	112.67	2.46	7.57	−2.55	−5.18	−7.94	2.91
	ωLH23tB	−161.33	−29.72	−115.24	114.04	106.74	2.46	7.59	−2.65	−5.85	−7.81	2.77
	ωLH23tE	−139.28	−20.96	−97.36	103.07	89.77	2.47	7.52	−1.68	−6.59	−6.45	2.15
	ωLH23tP	−131.61	−15.31	−90.91	102.03	84.72	2.46	7.49	−1.65	−6.16	−6.23	2.18
	ωLH23td	−28.62	21.74	−11.47	68.20	15.93	2.48	7.40	−0.09	11.64	5.56	5.61
	ωLH23tdB	−16.99	18.15	−4.26	54.92	14.08	2.56	7.41	0.57	11.41	7.55	5.85
	ωLH23tdE	−20.63	19.83	−5.71	60.50	13.92	2.54	7.39	0.44	10.93	6.94	5.65
	ωLH23tdP	−16.94	21.78	−2.46	55.82	13.16	2.58	7.43	0.56	7.30	4.47	4.43
	DLPNO-CCSD	−61.47	−16.76	−47.30	55.03	43.22	2.38	7.14	1.31	−2.29	1.43	1.65

---

Fig. 4 Comparison of mean absolute errors in the calculated ¹H NMR chemical shifts with different functionals relative to experiment. Open bars pertain to PBE0 structures, filled bars to PBE structures. DLPNO-CCSD results are also given. Panel (a): ‘hard’ part of the complex. Panel (b): ‘easy’ part of the complex.

The results are divided into two groups. The ‘hard’ nuclei H-a, H-c and H-i reside in the parts of the molecule that receive too large (negative or positive) spin density in the GH calculations. For the BH₂ group protons in the H-i case, this leads to overestimated positive shifts by several tens of ppm. For the methyl protons H-a and H-c in the mesitylimido ligands, overestimated negative shifts that increase very rapidly with EXX admixture are found, as discussed above and as seen in the leftmost bars of the MAE image 3(a). We omit the H-b centre from the plot as there is no experimental datum for it. On the other hand, the ‘easy’ centres from H-d to H-h imply far smaller errors by the standard DFT calculations, albeit the GHs PBE0 and PBE0-40 fare worse than PBE even in these cases [Table 3(b)].

The first general finding is that using the PBE structure improves the shifts significantly compared to the PBE0 structure, consistent with the lower spin contamination (see above). The bars representing the MAEs of the ‘hard’ cases in Fig. 4(a) become unequivocally shorter, reflecting also the smaller leakage of spin polarisation to the strongly bound imido ligands (and to the BH₂-group) at the PBE structure (Fig. 3). The MAEs are almost in all cases also correspondingly reduced for the more well-behaving cases depicted in Fig. 4(b). This underlines the importance of the quality of the input structure for the pNMR predictions.

Comparing the performance of the different functionals in the HFC calculations, as reflected in the predicted pNMR shifts, the scLHs clearly improve performance compared to LH20t or the GHs. For H-a and H-c the improvement is clearly most pronounced for the “undamped” scLH22ta and scLH21ct-SVWN-m, while a damping of the sc-factors in the LMF (scLH22t, scLH23t-mBR, scLH23t-mBR-P) obviously prevents the scLHs from achieving larger corrections. On the other hand the “undamped” scLHs overshoot for H-i, while the damped ones remain closer to the experimental shift. Overall the MAE for this subset is smaller for the undamped functionals. The related rung-4 functional KP16 also uses “undamped” sc-factors and improves even more for H-a and H-b while also overshooting for H-i. We note in passing that this functional gave large deviations for manganese HFCs in ref. 48, which have been traced back to artefacts in the core shells of the metal centre. Overall the improvements compared to the GHs clearly correlate with reduced spin contamination (Table 2) for the ‘hard’ nuclei. The remaining deviations are of similar magnitude as for the PBE or DLPNO-CCSD data, albeit distributed somewhat differently (Table 3). The “undamped” scLHs also fare best for the ‘easy’ signals, in this case clearly better than either DLPNO-CCSD or PBE.

The remaining significant deviations for the H-a, H-b, and H-i shifts, regardless of the employed GH, LH or scLH employed for the HFCs, remind of the situation found in ref. 72. In there, both larger discrepancies and larger dependence on EXX were found in the proton shifts of pyridine ligands of a Re complex than in the other parts (phosphane ligands). No spin expectation values were reported in that work, so we cannot be sure that the origin is the same as in the present case. In ref. 72, the structure was also optimised at the GGA level.

3.4.2.2 Range-separated local hybrid functionals. We move on to the ωLH22t RSLH and the related scRSLHs (Table 6 and Fig. 4, Fig. S1, S2) and start again with the ‘hard’ nuclei (H-a, H-c, H-i). We had seen above that adding only the (damped) sc-corrections in ωLH23tB, ωLH23tE or ωLH23tP did not reduce the spin contamination very much. This is reflected in the ¹H HFCs and pNMR shifts, which exhibit MAEs on the order of around 85–110 ppm (using the PBE structure). In each case the H-a and H-c signals are predicted much too upfield, by around 100 ppm or more, in comparison with experiment. These errors actually exceed those of the “damped” scLHs (see above), whereas the overestimation for H-i is comparable.

Interestingly, the DEC correction term in the “td”-LMF is more effective than the damped sc-corrections, not only in reducing the spin contamination but also in improving the pNMR shifts (Table 6). This is already the case for the ωLH23td RSLH, which does not feature sc-corrections. While the H-a and H-c signals remain much too negative at the PBE0 structure, consistent with the still elevated 〈S²〉 value (Table 2), the results are much improved at the PBE structure. The H-i shift from the –BH₂ group is now already slightly too low rather than too high with GHs, LHs or scLHs (see above). Adding various variants of sc-corrections on top of the td-LMF (ωLH23tdB, ωLH23tdE, ωLH23tdP) provides further improvements for H-a and H-c, and ωLH23tdP has now the overall lowest MAE for the ‘hard’ cases, clearly outperforming DLPNO-CCSD or PBE. However, the MAE for the ‘easy’ cases increases somewhat compared to DLPNO-CCSD and to scLHs and RSLHs/scRSLHs without the DEC term, due to somewhat larger deviations for the H-g and H-h signals. Here the “undamped” scLHs perform best. They suggest the initial assignment with the more negative value for H-g to hold (Table 6). The same ordering holds when using PBE, where the deviations are somewhat larger (Table 5). We note that those functionals (scLH22ta, scLH21ct-SVWN-m, ωLH23tdB, ωLH23tdE, ωLH23tdP, ωLH23td, KP16) that provide the lowest MAE for the “hard” protons H-a, H-c, and H-i, agree in predicting a value near +20 ppm for the unobserved H-b.

Along with the effect of the underlying structure, in our case the choice between optimisation using the PBE0 or PBE functionals, a related question may be posed regarding the method choice for the orbital shielding tensor. Table S7 in the SI compares the calculated orbital ¹H chemical shifts obtained using these two functionals. The shift differences for the hard centres (H-a, H-b, H-c and H-i) are all tiny as compared to the overall deviations from the experiment—in a fairly complicated molecule the interplay between the structure and functional used in the shielding calculation easily gives changes of this magnitude. For the easy nuclei (H-d, H-e, H-f, H-g, H-h), for which the changes could be more important, changing the functional used for shielding calculations gives the same, negative, direction of change (PBE gives smaller shift than PBE0) at both the PBE and PBE0 structures, for three of the five nuclei. At the apparently better structure (PBE) the shift difference PBE–PBE0 is negative for all the five centres, with the values ranging from −0.29 to −0.96 ppm. Functional-dependent changes of up to, say, 1 ppm in σ_orb could play a role in the relative performance of the strong correlation-corrected LHs and the RSLHs, which do not have the delocalisation correction, as the differences between the MAEs [Fig. 4(b)] of functionals belonging to these two groups have about this magnitude. However, the differences between the delocalisation-corrected and -uncorrected functionals is larger, in the few-ppm range, which implies that the dependence of σ_orb on the functional may not explain away the apparent, slightly inferior performance of the delocalisation-corrected RSLHs for the ‘easy’ centres.

A point to note is that we have not presently attempted to include solvation or intramolecular dynamics effects in the calculations of the paramagnetic Fe(V) complex. This will eventually become necessary with the development of the electronic-structure methodology in treating NMR shifts in challenging open-shell materials, for a detailed comparison with experiment.

3.5 Predicted ¹³C pNMR shifts

As so far no ¹³C pNMR shifts for this paramagnetic Fe(V) complex have been measured, we can use the present calculations for predicting them. We use a selection of methods based on the above-documented performance for the ¹H shifts and low spin contamination. We have selected the standard PBE functional, the “undamped” scLH22ta, the ωLH23tdP scRSLH, as well as DLPNO-CCSD methods for the computation of the HFCs. The same orbital shielding and g-tensor are used as above for the ¹H results. We present results obtained at the PBE structure in Table 7.

Table 7 Predicted ¹³C chemical shifts (in ppm) for the open-shell Fe(V) bis(imido) complex, with HFCs obtained at the 1c-X2C level with selected DFT functionals, or at the DKH2 level with DLPNO-CCSD. Results provided for the PBE structure, with other computational choices as for the ¹H results. The labeling of the ¹³C centres follows Fig. 1

Method	Type	C-a	C-b	C-c	C-j	C-k	C-l	C-p	C-d	C-e	C-f	C-g	C-h	C-m	C-n	C-o
PBE	GGA	47	361	15	66	1179	2144	−2183	21	141	11	114	−10	142	160	151
scLH22ta	scLH	51	338	50	37	959	2329	−2444	25	143	20	120	32	147	160	159
ωLH23tdP	scRSLH	38	338	57	515	1102	2064	−624	25	139	14	−56	−1	151	153	143
DLPNO-CCSD	Ab initio	244	924	202	−925	−76	3221	−3745	26	141	−1	155	−36	143	148	144

---

Focusing first on the ‘hard’ part of the molecule, the predicted ¹³C shifts in the imido ligands appear mostly systematically at common values across the three different-type DFT functionals. The signals for the centres C-a and C-c (the methyl group carbons) can, based on these calculations, be expected in the range 10…60 ppm, slightly downfield from the TMS reference, as the functionals agree sufficiently in their predictions. The CH-group carbon, C-b, signal is similarly anticipated at around 340…360 ppm, i.e., significantly downfield. The C-k and C-l centres gain significant positive spin density due to delocalisation of the SOMO, and the calculations place these centres at roughly 1000 and 2000 ppm, respectively. It is likely that observing these signals will present experimental challenges due to the presumably rapid relaxation. This last remark also applies to the diazole centre C-p, which is directly bonded to the metal atom and features a large negative predicted shift. Whereas PBE and scLH22ta place the signal at around −2000 ppm, ωLH23tdP and DLPNO-CCSD produce HFCs that result in much smaller and larger upfield shifts, respectively, of this signal. The DFT functionals are, hence, not in mutual agreement for either the C-p, nor the C-j signals, the latter of which is trans to the imido linkage. In these two cases the delocalisation error-corrected ωLH23tdP scRSLH predicts much smaller and larger shifts, respectively, as compared to the example GGA and scLH functionals in the table. Clearly, an experimental observation of these shift values could help further distinguish between functionals.

In contrast to the above-discussed situation of the proton signals, for the ¹³C shifts in the ‘hard’ part of the molecule, the shifts obtained using ab initio HFCs computed by the DLPNO-CCSD method do not even qualitatively agree with those in which the best present DFT functionals have been used. Instead, for each of these carbons the predicted shifts are very different from the DFT-based data, and in four of the total of seven cases the DLPNO-CCSD-based signals are found in significantly downfield positions as compared to the range of DFT-based results. The qualitatively different behaviour of the DLPNO-CCSD HFC in the proton and carbon cases suggests that the use of the frozen carbon 1s core in these calculations may play a role in causing the inferior performance of the method for ¹³C, as the potentially significant, negative spin-polarisation contribution from these orbitals is neglected. We do not presently have the computational capacity for testing this hypothesis for a model of the size of the present complex. On the other hand, two of the DLPNO-CCSD-based carbon shifts in the imido ligand (C-j and C-k), as well as C-p, deviate from the DFT ‘consensus’ to the opposite direction, upfield from the TMS reference. Clarification of the situation calls for further investigation.

Moving to the ¹³C signals within the less critical part of the molecule, C-d…C-h and C-m…C-o, we note that the results fall into two categories. The outermost 2,4,6-trimethylphenyl group comprises the carbon centres C-d…C-f and C-m…C-o, and the methods in Table 7 are quite unanimous concerning their shifts, suggesting a rather high predictive value for these nuclei. Whereas the centres in the phenyl ring, C-e and C-m…C-o are found in the range from 140 to 150 ppm, downfield from the reference, the methyl group signals C-f (the 2,6-positions) and C-d (the 4-position) are found at 10…20 ppm and 20…25 ppm. For C-f, however, the DLPNO-CCSD-based result already deviates somewhat from the trend of the DFT-based values. The mutual agreement of the computational predictions is less convincing for the ¹³C signals of the C-g and C-h nuclei in the CH-groups. Three of the four methods in Table 7 agree on the C-g shift residing in the 110…150 ppm range, while ωLH23tdP surprisingly predicts a negative shift. It was observed in the case of the ‘easy’ ¹H signals that this scRSLH in fact produced larger average errors in the ‘easy’ part than either the scLHs or the PBE GGA. This may suggest that, while the delocalisation-error correction in functionals such as ωLH23tdP alleviates the gross problems of spin-density spill-over and exaggerated valence-shell spin polarisation, as in the imido groups of the present complex, it may do so at the expense of precision in the less critical parts of the molecule. More comprehensive studies will be needed to shed light on this issue. The last signal, C-f, is predicted at 0…20 ppm, with the DLPNO-CCSD-based datum somewhat deviating from the bunch of DFT results. All in all, despite the occasional scatter, the best present calculations should provide a meaningful starting point for the experimental search of the hitherto unavailable ¹³C signals in the present paramagnetic complex.

4 Conclusions

We have applied first-principles paramagnetic NMR chemical shift theory to a challenging, novel Fe(V) bis(imido) compound, a spin doublet system for which experimental ¹H chemical shifts have recently been reported. In the calculations we applied methodology combining different electronic-structure models for the various quantities appearing in the expression of the shielding tensor, the g- and orbital shielding tensors and, in particular, the critically important hyperfine coupling tensors. As a fully ab initio treatment of the latter for a system like this is still a challenge, we have in particular evaluated state-of-the-art density functionals for the hyperfine couplings.

Standard global hybrid functionals lead to severe spin-density delocalisation and spin-polarisation problems in the present complex, causing excessive spillage of spin density to the strongly bound imido ligands and to the BH₂ group of the bis-carbene framework. Hence, unrealistically large positive or (due to the concomitant spin polarisation) negative proton chemical shifts are computed for these parts of the molecule, with no resemblance of the experimental signal pattern whatsoever. The computed 〈S²〉 expectation values indicate that these calculations are significantly affected by spin contamination. The semi-local PBE functional does not produce spin contamination but nevertheless apparently still features somewhat too much spin delocalisation. Notably, using PBE rather than spin-contaminated PBE0 to optimise the structures leads to shorter key metal–ligand distances, which overall diminishes the effects of spin contamination on the hyperfine couplings.

As these observations point to the crucial zero-sum game between self-interaction (delocalisation) and static-correlation errors in DFT, we subsequently evaluated recently developed families of local hybrid (LH) and range-separated local hybrid (RSLH) functionals based on the exact-exchange energy density. Several of these include strong-correlation and/or delocalisation-error corrections and have been shown to escape to some extent the zero-sum game. We also compared ab initio DLPNO-CCSD results for the hyperfine couplings.

While results with the “uncorrected” LH20t local hybrid functional also suffer from spin-contamination, albeit less so than, e.g., with PBE0, use of its strong-correlation corrected variants provide significant improvement both in the strongly bound imido ligands (as well as the BH₂ group) and in the less dramatically affected diazole and phenyl groups of the complex. This holds in particular for scLH22ta and the earlier scLH21ct-SVWN-m functional, which both exhibit “undamped” strong-correlation factors. Indeed, these two functionals provide the best performance overall, better than the ab initio DLPNO-CCSD method, at much lower computational cost. While the achieved precision is not better than to within tens of ppm from the experiment for the ‘hard’ part of the molecule, the calculations using these scLH-based hyperfine couplings are much more useful for NMR assignment and analysis than those employing standard functionals that can be hundreds of ppm off. The reduced spin contamination problem is also reflected in the corresponding 〈S²〉 values that reside much closer to the formal value of 0.75 for a doublet system, in contrast to global hybrid functionals. The related KP16 functional exhibits “undamped” strong-correlation terms and also performs well for the ‘hard’ positions.

Results with the “uncorrected” ωLH22t range-separated local hybrid suffer from appreciable spin contamination. Several variants with strong-correlation and/or delocalization error corrections have been evaluated. Functionals with only (damped) strong-correlation terms (ωLH23tX; X = B, E, P) reduce spin contamination only moderately and therefore still give large errors for the ‘hard’ ¹H nuclei. Inclusion of only the delocalisation-error corrections in ωLH23td is more effective in reducing spin contamination and the resulting artefacts. Range-separated local hybrid functionals with both types of contributions (ωLH23tdX; X = B, E, P) provide the best performance for the ‘hard’ ¹H nuclei, with deviations of more than an order of magnitude lower than for uncorrected (global, range-separated, local, and local range-separated) hybrids, but they give somewhat larger deviations for the H-g and H-h positions in the ‘easy’ part than the best scLHs or PBE. The DLPNO-CCSD method performs better than global hybrids for the more difficult positions but cannot compete with the best scLHs and scRSLHs or even with PBE.

The methods that fared best in different aspects of the proton shift calculation, the standard PBE, one selected representative of both scLH and scRSLH families, as well as the DLPNO-CCSD method, were selected to provide the hyperfine coupling tensors for predicting the ¹³C shifts, which have not been experimentally reported for this system. It is notable that, for most of the nuclei, the three quite different DFT methods resulted in relatively narrow ranges of predicted ¹³C shifts, implying useful predictive value of the calculations. The scRSLH ωLH23tdP, which provided the best ¹H data in the region affected by the spin-density spill-over, did result in outlier predictions for a couple of the ¹³C centres. While the ultimate criterion for the DFT performance will be provided by the eventual ¹³C experiment, it appears clear that the present DLPNO-CCSD-based hyperfine couplings do not furnish a suitable basis for predicting the ¹³C shifts in the present system. Particularly for the ‘hard’ part of the molecule, the data bear no resemblance to the relative DFT consensus, possibly because for computational reasons we had to employ the frozen core approximation encompassing, e.g., the carbon 1s orbitals, in these calculations.

Overall it appears that functionals based on the exact-exchange energy density, that provide at least a partial escape from the so-called zero-sum game between delocalisation and static-correlation errors, provide an important tool to predict the hyperfine couplings and the related pNMR shifts in challenging transition-metal complexes. The fact that some of the scRSLHs performing best for the ‘hard’ part of the complex but have somewhat larger deviations for some of the nuclei in the ‘easy’ part suggests room for further improvement in the development of modern exchange–correlation functionals.

Author contributions

A. Pyykkönen: general investigation, visualisation of methods and results, data curation except for the HFC calculations using the novel local hybrid functionals. A. Wodyński: HFC calculations using novel hybrid functionals. M. Kaupp: conceptualisation and writing of the original draft. J. Vaara: conceptualisation, writing of the original draft, supervision. All authors contributed together to reviewing/editing of the manuscript.

Conflicts of interest

There are no conflicts to declare.

Data availability

The data supporting this article have been included as part of the SI. Optimised structures of the paramagnetic complex at the PBE and PBE0 levels, optimised structure of the diamagnetic complex at the PBE0 level, predicted ¹³C shifts of the diamagnetic complex, results of ab initio ligand-field theory analysis for the paramagnetic complex and figures illustrating the deviation of the predicted ¹H shifts from the experiment and comparison of the orbital shifts of the paramagnetic complex obtained by using the PBE and PBE0 functionals. See DOI: https://doi.org/10.1039/d5cp02544c

Acknowledgements

We acknowledge funding from the Academy of Finland (grants 331008 and 361326) and U. Oulu (Kvantum Institute). Computations were carried out at CSC – the Finnish IT Centre for Science. Work in Berlin has been supported by Deutsche Forschungsgemeinschaft under project KA1187/14-2.

Notes and references

1. A. J. Pell, G. Pintacuda and C. P. Grey, Progr. NMR Spectrosc., 2019, 111, 1–271.
2. J. Vaara, Ann. Acad. Sci. Fenn., 2023, 2, 96–114.
3. J. Vaara, S. A. Rouf and J. Mareš, J. Chem. Theory Comput., 2015, 11, 4840–4849.
4. S. A. Rouf, PhD thesis, University of Oulu, Oulu, Finland, 2017.
5. S. A. Rouf, J. Mareš and J. Vaara, J. Chem. Theory Comput., 2017, 13, 3721–3745.
6. S. Moon and S. Patchkovskii, in Calculation of NMR and EPR Parameters: Theory and Applications, ed. M. Kaupp, M. Bühl and V. G. Malkin, Wiley-VCH, 2004, pp. 325–338.
7. T. O. Pennanen and J. Vaara, J. Chem. Phys., 2005, 123, 174102.
8. W. Van den Heuvel and A. Soncini, J. Chem. Phys., 2013, 138, 054113.
9. B. Martin and J. Autschbach, J. Chem. Phys., 2015, 142, 054108.
10. R. J. Kurland and B. R. McGarvey, J. Magn. Reson., 1970, 2, 286–301.
11. J. E. Harriman, Theoretical Foundations of Electron Spin Resonance, Academic Press, 1978.
12. N. F. Ramsey, Phys. Rev., 1952, 86, 243–246.
13. W. Van den Heuvel and A. Soncini, Phys. Rev. Lett., 2012, 109, 073001.
14. F. Gendron, K. Sharkas and J. Autschbach, J. Phys. Chem. Lett., 2015, 6, 2183–2188.
15. S. A. Rouf, J. Mareš and J. Vaara, J. Chem. Theory Comput., 2015, 11, 1683–1691.
16. S. A. Rouf, V. B. Jakobsen, J. Mareš, N. D. Jensen, C. J. McKenzie, J. Vaara and U. G. Nielsen, Solid State Nucl. Magn. Reson., 2017, 87, 29–37.
17. A. B. A. Andersen, A. Pyykkönen, H. J. A. Jensen, V. McKee, J. Vaara and U. G. Nielsen, Phys. Chem. Chem. Phys., 2020, 22, 8048–8059.
18. A. Pyykkönen, R. Feher, F. H. Köhler and J. Vaara, Inorg. Chem., 2020, 59, 9294–9307.
19. A. Jaworski and N. Hedin, Phys. Chem. Chem. Phys., 2022, 24, 15230–15244.
20. A. Pyykkönen and J. Vaara, Phys. Chem. Chem. Phys., 2023, 25, 3121–3135.
21. M. Mohan, A. B. A. Andersen, J. Mareš, N. D. Jensen, U. G. Nielsen and J. Vaara, Phys. Chem. Chem. Phys., 2023, 25, 24081–24096.
22. B. O. Roos, Chem. Phys. Lett., 1972, 15, 153–159.
23. C. Angeli, R. Cimiraglia, S. Evangelisti, T. Leininger and J.-P. Malrieu, J. Chem. Phys., 2001, 114, 10252–10264.
24. C. Angeli, R. Cimiraglia and J.-P. Malrieu, Chem. Phys. Lett., 2001, 350, 297–305.
25. C. Angeli, R. Cimiraglia and J.-P. Malrieu, J. Chem. Phys., 2002, 117, 9138–9153.
26. D. Ganyushin and F. Neese, J. Chem. Phys., 2013, 138, 104113.
27. D. Ganyushin and F. Neese, J. Chem. Phys., 2006, 125, 024103.
28. B. Engels, in Calculation of NMR and EPR Parameters. Theory and Applications, ed. M. Kaupp, M. Bühl and V. G. Malkin, Wiley-VCH Verlag GmbH & Co. KGaA, 2004, ch. 30, pp. 483–492.
29. K. Sharkas, B. Pritchard and J. Autschbach, J. Chem. Theory Comput., 2015, 11, 538–549.
30. C. Riplinger and F. Neese, J. Chem. Phys., 2013, 138, 034106.
31. C. Riplinger, P. Pinski, U. Becker, E. F. Valeev and F. Neese, J. Chem. Phys., 2016, 144, 024109.
32. M. Munzarová and M. Kaupp, J. Phys. Chem. A, 1999, 103, 9966–9983.
33. M. L. Munzarová and M. Kaupp, J. Am. Chem. Soc., 2000, 122, 11900.
34. C. J. Schattenberg, T. M. Maier and M. Kaupp, J. Chem. Theory Comput., 2018, 14, 5653–5672.
35. A. J. Cohen, P. Mori-Sánchez and W. Yang, J. Chem. Phys., 2008, 129, 121104.
36. A. J. Cohen, P. Mori-Sánchez and W. Yang, Chem. Rev., 2012, 112, 289–320.
37. B. G. Janesko, E. Proynov, J. Kong, G. Scalmani and M. J. Frisch, J. Phys. Chem. Lett., 2017, 8, 4314–4318.
38. M. Kaupp, A. Wodyński, A. V. Arbuznikov, S. Fürst and C. J. Schattenberg, Acc. Chem. Res., 2024, 57, 1815–1826.
39. T. M. Maier, A. V. Arbuznikov and M. Kaupp, Wiley Interdiscip. Rev.: Comput. Mol. Sci., 2019, 9, e1378.
40. R. Haunschild and G. E. Scuseria, J. Chem. Phys., 2010, 132, 224106.
41. S. Fürst, M. Haasler, R. Grotjahn and M. Kaupp, J. Chem. Theory Comput., 2023, 19, 488–502.
42. A. Wodyński and M. Kaupp, J. Chem. Theory Comput., 2020, 16, 314–325.
43. A. Wodyński, A. V. Arbuznikov and M. Kaupp, J. Chem. Phys., 2021, 155, 144101.
44. A. Wodyński and M. Kaupp, J. Chem. Theory Comput., 2022, 18, 6111–6123.
45. A. Wodyński, A. V. Arbuznikov and M. Kaupp, J. Chem. Phys., 2023, 158, 244117.
46. S. Fürst, M. Kaupp and A. Wodyński, J. Chem. Theory Comput., 2023, 19, 8639–8653.
47. J. Kirkpatrick, B. McMorrow, D. H. P. Turban, A. L. Gaunt, J. S. Spencer, A. G. Matthews, A. Obika, L. Thiry, M. Fortunato and D. Pfau, et al., Science, 2021, 374, 1385–1389.
48. A. Wodyński, B. Lauw, M. Reimann and M. Kaupp, J. Chem. Theory Comput., 2024, 20, 2033–2048.
49. C. J. Schattenberg and M. Kaupp, J. Phys. Chem. A, 2024, 128, 2253–2271.
50. C. J. Schattenberg, A. Wodyński, H. Åström, D. Sundholm, M. Kaupp and S. Lehtola, J. Phys. Chem. A, 2023, 127, 10896–10907.
51. J. L. Martinez, S. A. Lutz, H. Yang, J. Xie, J. Telser, B. M. Hoffman, V. Carta, M. Pink, Y. Losovyj and J. M. Smith, Science, 2020, 370, 356–359.
52. N. Ramsey, Phys. Rev., 1950, 78, 303–307.
53. C. Adamo and V. Barone, J. Chem. Phys., 1999, 110, 6158–6170.
54. J. P. Perdew, K. Burke and M. Ernzerhof, Phys. Rev. Lett., 1996, 77, 3865–3868; J. P. Perdew, K. Burke and M. Ernzerhof, Phys. Rev. Lett., 1997, 78, 1396.
55. S. Grimme, J. Antony, S. Ehrlich and H. Krieg, J. Chem. Phys., 2010, 132, 154104.
56. S. Grimme, S. Ehrlich and L. Goerigk, J. Comput. Chem., 2011, 32, 1456–1465.
57. TURBOMOLE V7.2 2017, a development of University of Karlsruhe and Forschungszentrum Karlsruhe GmbH, 1989-2007, TURBOMOLE GmbH, since 2007, Available from https://www.turbomole.com.
58. S. G. Balasubramani, G. P. Chen, S. Coriani, M. Diedenhofen, M. S. Frank, Y. J. Franzke, F. Furche, R. Grotjahn, M. E. Harding, C. Hättig, A. Hellweg, B. Helmich-Paris, C. Holzer, U. Huniar, M. Kaupp, A. Marefat Khah, S. Karbalaei Khani, T. Müller, F. Mack, B. D. Nguyen, S. M. Parker, E. Perlt, D. Rappoport, K. Reiter, S. Roy, M. Ruckert, G. Schmitz, M. Sierka, E. Tapavicza, D. P. Tew, C. van Wüllen, V. K. Voora, F. Weigend, A. Wodynski and J. M. Yu, J. Chem. Phys., 2020, 152, 184107.
59. Y. J. Franzke, C. Holzer, J. H. Andersen, T. Begušić, F. Bruder, S. Coriani, F. Della Sala, E. Fabiano, D. A. Fedotov, S. Fürst, S. Gillhuber, R. Grotjahn, M. Kaupp, M. Kehry, M. Krstić, F. Mack, S. Majumdar, B. D. Nguyen, S. M. Parker, F. Pauly, A. Pausch, E. Perlt, G. S. Phun, A. Rajabi, D. Rappoport, B. Samal, T. Schrader, M. Sharma, E. Tapavicza, R. S. Treß, V. Voora, A. Wodyński, J. M. Yu, B. Zerulla, F. Furche, C. Hättig, M. Sierka, D. P. Tew and F. Weigend, J. Chem. Theory Comput., 2023, 19, 6859–6890.
60. M. Dolg, U. Wedig, H. Stoll and H. Preuss, J. Chem. Phys., 1987, 86, 866–872.
61. J. M. L. Martin and A. Sundermann, J. Chem. Phys., 2001, 114, 3408–3420.
62. F. Weigend and R. Ahlrichs, Phys. Chem. Chem. Phys., 2005, 7, 3297–3305.
63. Ref. 62, recontracted for DKH2 by D. A. Pantazis.
64. M. J. Frisch, G. W. Trucks, H. B. Schlegel, G. E. Scuseria, M. A. Robb, J. R. Cheeseman, G. Scalmani, V. Barone, G. A. Petersson, H. Nakatsuji, X. Li, M. Caricato, A. V. Marenich, J. Bloino, B. G. Janesko, R. Gomperts, B. Mennucci, H. P. Hratchian, J. V. Ortiz, A. F. Izmaylov, J. L. Sonnenberg, D. Williams-Young, F. Ding, F. Lipparini, F. Egidi, J. Goings, B. Peng, A. Petrone, T. Henderson, D. Ranasinghe, V. G. Zakrzewski, J. Gao, N. Rega, G. Zheng, W. Liang, M. Hada, M. Ehara, K. Toyota, R. Fukuda, J. Hasegawa, M. Ishida, T. Nakajima, Y. Honda, O. Kitao, H. Nakai, T. Vreven, K. Throssell, J. A. Montgomery, Jr., J. E. Peralta, F. Ogliaro, M. J. Bearpark, J. J. Heyd, E. N. Brothers, K. N. Kudin, V. N. Staroverov, T. A. Keith, R. Kobayashi, J. Normand, K. Raghavachari, A. P. Rendell, J. C. Burant, S. S. Iyengar, J. Tomasi, M. Cossi, J. M. Millam, M. Klene, C. Adamo, R. Cammi, J. W. Ochterski, R. L. Martin, K. Morokuma, O. Farkas, J. B. Foresman and D. J. Fox, Gaussian 16, Revision C.02, Gaussian, Inc., Wallingford CT, 2019.
65. F. Neese, Wiley Interdiscip. Rev.: Comput. Mol. Sci., 2018, 8, e1327.
66. Y. Guo, K. Sivalingam, E. F. Valeev and F. Neese, J. Chem. Phys., 2016, 144, 094111.
67. M. Douglas and N. M. Kroll, Ann. Phys., 1974, 82, 89–155.
68. B. A. Hess, Phys. Rev. A:At., Mol., Opt. Phys., 1986, 33, 3742–3748.
69. T. Yanai, D. P. Tew and N. C. Handy, Chem. Phys. Lett., 2004, 393, 51.
70. M. Haasler, T. M. Maier, R. Grotjahn, S. Gückel, A. V. Arbuznikov and M. Kaupp, J. Chem. Theory Comput., 2020, 16, 5645–5657.
71. J. Kong and E. Proynov, J. Chem. Theory Comput., 2016, 12, 133–143.
72. Y. J. Franzke, F. Bruder, S. Gillhuber, C. Holzer and F. Weigend, J. Phys. Chem. A, 2024, 128, 670–686.
73. V. Barone and M. Cossi, J. Phys. Chem. A, 1998, 102, 1995–2001.
74. M. Atanasov, D. Ganyushin, K. Sivalingam and F. Neese, in Molecular Electronic Structures of Transition Metal Complexes II, no. 143 in Structure and Bonding, ed. D. M. P. Mingos, P. Day and J. P. Dahl, Springer, Berlin, Heidelberg, 2011, pp. 149–220.
75. K. Theilacker, PhD thesis, Technische Universität Berlin, 2016.
76. M. Repisky, S. Komorovsky, M. Kadek, L. Konecny, U. Ekström, M. Kaupp, K. Ruud, O. Malkina and V. Malkin, J. Chem. Phys., 2020, 152, 184101.

---