Computational NMR of the iron pyrazolylborate complexes [Tp₂Fe]⁺ and Tp₂Fe including solvation and spin-crossover effects

Abstract

Transition metal complexes have important roles in many biological processes as well as applications in fields such as pharmacy, chemistry and materials science. Paramagnetic nuclear magnetic resonance (pNMR) is a valuable tool in understanding such molecules, and theoretical computations are often advantageous or even necessary in the assignment of experimental pNMR signals. We have employed density functional theory (DFT) and the domain-based local pair natural orbital coupled-cluster method with single and double excitations (DLPNO-CCSD), as well as a number of model improvements, to determine the critical hyperfine part of the chemical shifts of the iron pyrazolylborate complexes [Tp₂Fe]⁺ and Tp₂Fe using a modern version of the Kurland–McGarvey theory, which is based on parameterising the hyperfine, electronic Zeeman and zero-field splitting interactions via the parameters of the electron paramagnetic resonance Hamiltonian. In the doublet [Tp₂Fe]⁺ system, the calculations suggest a re-assignment of the ¹³C signal shifts. Consideration of solvent via the conductor-like polarisable continuum model (C-PCM) versus explicit solvent molecules reveals C-PCM alone to be insufficient in capturing the most important solvation effects. Tp₂Fe exhibits a spin-crossover effect between a high-spin quintet (S = 2) and a low-spin singlet (S = 0) state, and its recorded temperature dependence can only be reproduced theoretically by accounting for the thermal Boltzmann distribution of the open-shell excited state and the closed-shell ground-state occupations. In these two cases, DLPNO-CCSD is found, in calculating the hyperfine couplings, to be a viable alternative to DFT, the demonstrated shortcomings of which have been a significant issue in the development of computational pNMR.

---

1 Introduction

Nuclear magnetic resonance (NMR) of open-shell, paramagnetic compounds (pNMR) has become increasingly important in the fields of biochemistry and materials science.^1–5 Theoretical determination of the NMR chemical shifts via electronic structure methods can aid in the assignment and interpretation of experimental spectra,⁶ and sometimes reveal a need for re-assignment of experimental spectra.^7,8 The classic Kurland–McGarvey theory of paramagnetic shielding⁹ expresses the hyperfine part of the shielding tensor in terms of the electron paramagnetic resonance (EPR) parameters of the ground multiplet of 2S + 1 electronic states, where S is the effective electron spin. In modern quantum-chemical implementations^10–12 of this theory, one requires the ability to reliably compute the g-tensor, the zero-field splitting (ZFS, D, in S > 1/2 systems with more than one unpaired electron) and the hyperfine couplings (HFC, A) of the system.

Transition metal ions such as Fe(II) and Fe(III) play important roles in a number of biological processes and have many applications in pharmacy.¹³ Iron systems have also received attention as cost-effective and environmentally friendly catalysts in chemistry and materials science.¹⁴ The present work focuses on the calculation of NMR chemical shifts in isolated paramagnetic iron complexes. Some modern computational pNMR work on iron systems has been previously published.^15–18 Mondal and Kaupp^15,16 computed pNMR chemical shifts of solid-state iron systems as part of studies on a series of lithium-ion battery cathode materials. Martin and Autschbach¹⁷ performed Kohn–Sham (KS) density functional theory (DFT) pNMR calculations on a set of paramagnetic organometallic complexes to investigate the roles of the KS delocalisation error, Gaussian-type versus Slater-type basis sets, relativistic effects, and ZFS. Our earlier work¹⁸ focused on computational pNMR predictions of a series of transition metal pyrazolylborate complexes, including iron systems. In this work, we study the iron pyrazolylborate complexes [Tp₂Fe(III)]⁺ and Tp₂Fe(II) encountered in our previous work¹⁸ more rigorously, exploring various computational methods for calculating their EPR parameters. An image of this type of molecule is shown in Fig. 1. With the former, spin doublet (S = 1/2) complex, we are particularly interested in whether model improvements such as more accurate g-tensor and ZFS calculations, or consideration of solvation effects can explain the observed discrepancy between theory and experiment for carbon shifts, which contrasted with the otherwise generally good theoretical results in the rest of the investigated series of metal complexes.¹⁸ The latter complex is a spin-crossover compound, for which different computational fitting methods are explored in this work. In addition to our earlier work,¹⁸ these systems were studied experimentally by Feher,¹⁹ whose measurements are compared to presently, as well as previously.

Fig. 1 Optimised geometry of Tp₂Fe(II) (d⁶, S = 0 ⇌ 2, charge 0) in the high-spin state with numbering used for NMR signal assignment. [Tp₂Fe(III)]⁺ (d⁵, S = 1/2, charge +1) has a similar structure and numbering. Fe atom in grey, N in yellow, B in green, C in black and H in white.

Some metal complexes exhibit behaviour in which their spin state changes due to an external stimulus, such as temperature, pressure or light, usually involving conversion between low-spin and high-spin configurations.²⁰ This so-called spin-crossover was first observed by Cambi and Szegö²¹ in 1931, and it has various potential applications in, e.g., molecular electronics, data storage and display devices.^22,23 The phenomenon is commonly seen in first-row transition metal complexes with d⁴ to d⁷ electron configuration and octahedral ligand geometry. One system displaying spin-crossover is Tp₂Fe(II), shown in Fig. 1, which switches between a diamagnetic singlet (S = 0) state and a paramagnetic quintet (S = 2) state. Recently, the spin-state energetics of manganese spin-crossover complexes have been studied computationally by ab initio methods by Drosou, Mitsopoulou and Pantazis,²⁴ who found complete active-space self-consistent field (CASSCF)/n-electron valence-state perturbation theory (NEVPT2) calculations unable to accurately describe the spin-state energies of their systems. This suggests that the CASSCF/NEVPT2 method we routinely use for the EPR property calculations of relevance to pNMR shifts may not be well-suited for the determination of energetics. The authors recommend, instead, DLPNO-CCSD(T) for energy calculations, which has also been utilised earlier by Flöser et al.²⁵ Iron spin-crossover compounds, including Tp₂Fe, have been studied several times in the past.^26–30 Recent papers also report investigations of spin-crossover systems and particularly their applications.^31–36 Adjacent to the topic, related spin equilibrium systems have been studied computationally using DFT and Boltzmann statistics by Ke et al.^37,38 with the high-spin and low-spin configurations arising from two distinct paramagnetic centers. In this work, we examine the significance of the spin-crossover effect in predicting the pNMR chemical shifts in an isolated Tp₂Fe(II) molecule, and how different computational fitting methods for handling its equilibrium state succeed in reproducing the experimental results.^18,19

Section 2 gives a brief overview of the theory behind paramagnetic NMR computations including the spin-crossover effect. Section 3 describes the details of the specific computational methods adopted in this work. Section 4.1 discusses the iron complex [Tp₂Fe(III)]⁺ and the various computational methods for prediction its NMR chemical shifts, whereas Section 4.2 focuses on the spin-crossover compound Tp₂Fe(II) and ways of handling its equilibrium state computationally.

2 Theoretical background

The present theoretical prediction of NMR chemical shifts of paramagnetic systems is based on the Kurland–McGarvey theory⁹ in its modern form,^10–12 in which magnetic interactions are parameterised in the ground multiplet by EPR parameters: HFCs, ZFS and the g-tensor. The resulting expression for the Cartesian ετ-component of the pNMR shielding tensor for nucleus K iswhere σ^orb_K,ετ, g_εa and A_K,bτ are the Cartesian components of the orbital shielding tensor, the g-tensor and the hyperfine coupling tensor, respectively. The quantities μ_B, γ_K, ħ, k and T are the Bohr magneton, gyromagnetic ratio of nucleus K, reduced Planck constant, Boltzmann constant and absolute temperature, respectively. The dyadic 〈SS〉 is expressed in terms of the effective spin operators S,where Q_nm areand |n〉 and E_n are the eigenfunctions and eigenvalues of the zero-field splitting (D-tensor) Hamiltonian H_ZFS = S·D·S at the limit of a vanishing external magnetic field. The chemical shift δ_K is expressed with the isotropic average of the shielding tensor,relative to a reference shielding σ_ref, as

δ_K = σ_ref − σ_K,iso. (5)

With sufficiently large S, the assumption on the form of the ZFS Hamiltonian (S·D·S) is approximate and, strictly speaking, a more general formalism³⁹ should be used. It is also possible to directly use optimised ground and excited electronic states and the couplings between them in a multiconfigurational framework^40,41 instead of EPR parameters to calculate the shieldings. This method is not yet in widespread use, however, due to the computational limitations related to the size and the difficulty of selecting the active space.

The spin-crossover system Tp₂Fe undergoes rapid spin interconversion near room temperature between a low-spin, S = 0 electronic ground state and a high-spin, S = 2 excited state. The populations of these states follow the Boltzmann distribution, with the probability of state i being

where is the partition function. The (2S + 1)-fold degeneracy of the electronic states is included implicitly through the entropy term in the exponent. In order to calculate the NMR chemical shift of the equilibrium state, the enthalpy (energy) gap ΔH = H₂ − H₀ and the entropy gap ΔS = S₂ − S₀ must be determined, and the NMR shieldings of both the high-spin and low-spin states must be computed.

3 Computational methods

3.1 Geometry

Throughout the present discussion, unless otherwise specified, geometry optimisations were carried out with the TURBOMOLE programme⁴² using unrestricted Kohn–Sham (UKS) DFT and the PBE0 functional⁴³ with 25% exact Hartree–Fock exchange energy (the PBE-25 notation is later used synonymously). Dispersion effects were accounted for with the DFT-D3 BJ dispersion correction.^44,45 Relativistic effects were approximated by placing the Stuttgart-type scalar relativistic effective core potential ECP10MDF⁴⁶ on the metal center, along with a corresponding 6s5p3d2f1g/8s7p6d2f1g (contracted/uncontracted notation) valence basis set.⁴⁷ The def2-TZVP basis set⁴⁸ was employed for ligand atoms. Coordinates for the optimised structures are provided in the ESI† (Tables S1–S6 and Fig. S1–S4).

3.2 g- and ZFS tensors

The “standard method” adopted previously^49,50 for computing the EPR parameters necessary to determine the NMR shieldings has been successfully applied to a number of isolated molecules in previous works.^7,8,18 It involves calculating the g-tensor and ZFS at the CASSCF⁵¹/NEVPT2^52–54 level of theory. The basic NEVPT2 calculation for the g-tensor and ZFS uses a state-average (SA) cas(n,5) CASSCF wavefunction that correlates n d-electrons of the metal center in the five metal 3d-orbitals. The Douglas–Kroll–Hess (DKH2)^55,56 approximation has been used in this context since ref. 50 to account for scalar relativistic effects in the optimisation of the SA-CASSCF states, and quasidegenerate perturbation theory (QDPT)^57,58 is utilised to incorporate spin–orbit effects using the picture-changed spin–orbit operator. Only the spin–orbit contribution to the ZFS is included, since the spin–spin contribution can be assumed to be small, e.g., as we demonstrated in our earlier work.¹⁸ This calculation is carried out with the ORCA software.^59,60 A locally dense basis set is used, in which the DKH-def2-TZVP basis set⁶¹ is placed on the metal center and atoms directly bonded to it, and a the smaller DKH-def2-SVP basis set is placed on the more distant ligand atoms.

3.3 HFC tensors

The DFT-based, fully relativistic four-component matrix Dirac–Kohn–Sham (mDKS) method^62,63 was mainly used to evaluate the hyperfine couplings. HFC calculations require polarisation of the core electrons in order to obtain reasonable results; for this reason unrestricted DFT (UKS) is necessary as, e.g., restricted open-shell DFT does not allow core spin polarisation. The mDKS calculation of hyperfine couplings is performed with the ReSpect program package⁶⁴ using mostly the PBE0 functional with 25% exact exchange, which again has been the “standard choice” for such calculations since the work in ref. 50. Both scalar relativistic and spin–orbit effects are considered variationally within this method. The decontracted Dyall-VTZ basis set⁶⁵ is placed on the metal atom, and a decontracted DKH-def2-TZVP basis set on the ligand atoms. This is the large wavefunction component basis set, and the small-component basis is generated by kinetic balance.

3.4 Orbital shielding

The orbital shielding term, which is often less critical than the hyperfine term, is calculated at the UKS/PBE0 level of theory employing the Gaussian software⁶⁶ and the DKH-def2-TZVP basis set. The NMR shieldings of the diamagnetic, closed-shell reference compounds, TMS for ¹H and ¹³C, BF₃ (O(C₂H₅)₂) for ¹¹B, and CH₃NO₂ for ¹⁴N, are computed at the restricted Kohn–Sham (RKS)/DKH-def2-TZVP level with the same functional as used in the corresponding HFC calculation. The reference shieldings are provided in Table S7 in the ESI.† In the present work, the computation of the diamagnetic (S = 0) low-spin state of Tp₂Fe follows the same formula.

3.5 Model improvements

In this work, four distinct modifications to the standard method are tested:

(1) The CASSCF active space in the NEVPT2 calculation of the g-tensor and ZFS is expanded from cas(n,5) to cas(n,12) by including two occupied metal–ligand bonding orbitals and five empty metal 4d-orbitals.

(2) Solvation effects⁶⁷ are included in the HFC computations through the conductor-like polarisable continuum model C-PCM,⁶⁸ approximating electrostatic solute–solvent interactions. Importantly, non-electrostatic effects such as hydrogen bonding are not represented by the C-PCM model. Additional calculations with explicit solvent molecules are performed for [Tp₂Fe]⁺ for comparison. An image of the [Tp₂Fe]⁺ structure with twelve acetone molecules is shown in Fig. S3 in the ESI.† Furthermore, the presence of the [PF₆]⁻ counterion is considered for [Tp₂Fe]⁺.

(3) The amount of exact Hartree–Fock exchange energy in the mDKS DFT hybrid functional in HFC calculations is varied (0%, 25%, 40% and 100%).

(4) The HFCs are calculated using the domain-based local pair natural orbital coupled-cluster method with single and double excitations (DLPNO-CCSD) employing lambda equations for unrelaxed coupled-cluster density,⁶⁹ now implemented in ORCA. This is an approximation to CCSD that scales nearly linearly with system size but achieves good accuracy, recovering almost all of the CCSD correlation energy. However, similar to CCSD, it is a single-reference method, which may prove limiting with systems that have strong multireference character. This method has recently been applied with good results to HFC computations in pNMR by Jaworski and Hedin.⁷⁰

3.6 Spin-crossover

The equilibrium-state NMR chemical shifts of the spin-crossover system Tp₂Fe are expressed as Boltzmann averages between the high-spin S = 2 state and the low-spin S = 0 state by fitting the enthalpy gap ΔH and the entropy gap ΔS to experimental NMR data. The sum of squared errors is minimised to find optimal values for ΔH and ΔS. The available experimental data, reported by Feher,^18,19 cover ¹H, ¹³C, ¹¹B and ¹⁴N chemical shifts, each at one temperature (excluding the N2 center) and the NMR temperature series of all four ¹H nuclei from 290 K to 455 K. The fitting is performed in three different ways.

(1) Global fit to the ¹H temperature series of all four protons, resulting in two fitting parameters ΔH(H) and ΔS(H).

(2) Individual fits to each of the four ¹H temperature series with each equivalent group assigned its own fitting parameter, resulting in eight fitting parameters ΔH(H3), ΔH(H4), ΔH(H5), ΔH(BH), ΔS(H3), ΔS(H4), ΔS(H5) and ΔS(BH). For nuclei other than protons, the average values ΔH_avg and ΔS_avg are used.

(3) Global fit to all experimental data including the ¹H temperature series, as well as ¹³C, ¹¹B and ¹⁴N shifts, resulting in two fitting parameters ΔH(all) and ΔS(all).

4 Results

4.1 Doublet iron complex [Tp₂Fe]⁺

The measurement and computational temperatures with this molecule were 305 K (¹H), 300 K (¹³C) and 298 K (¹¹B and ¹⁴N). The full set of results is tabulated in Table S8 in the ESI,† and a computational ¹H temperature series is displayed in Fig. S5 (ESI†).

4.1.1 CASSCF active space. For the doublet [Tp₂Fe]⁺ system, an expanded, 12-orbital CASSCF active space in the NEVPT2 computation of the g-tensor has only a minor impact on the final NMR shielding, with changes from the 5-orbital calculation of up to roughly 2 ppm for hydrogen and 6 ppm for carbon (Table 1). Therefore, at least in the case of this specific system, the increased size of the CASSCF active space is not reflected in significantly improved NMR chemical shifts. Both levels of theory are in very good agreement with ¹H and ¹¹B experiments, except in the case of H5, where the agreement is weaker. The cas(5,5)-level calculation is in fact slightly closer to experimental data, suggesting it is benefiting from error cancellation with, e.g., errors in HFC calculations, such as neglect of solvation, counterion and molecular dynamics effects.

Table 1 Experimental^18,19 and computational ¹H, ¹³C and ¹¹B NMR chemical shifts (ppm) for [Tp₂Fe]⁺ with the g-tensor calculated at the NEVPT2 level based on cas(5,5) and cas(9,12) wavefunctions. HFCs calculated at the mDKS/PBE0 level in vacuo

	H3	H4	H5	BH
Exptl^18,19	−52.9	−13.2	−8.2	39.4
cas(5,5)	−51.3	−10.0	−2.0	41.5
cas(9,12)	−50.9	−9.4	0.3	42.0

---

	C3	C4	C5	B
Exptl^18,19	2.5	175.8	19.6	95.2
cas(5,5)	34.7	−3.6	93.7	91.2
cas(9,12)	34.6	2.5	91.6	92.5

---

In stark contrast to computed proton and boron shifts, there is a severe discrepancy in carbon shifts, as was discussed previously.¹⁸ The potential causes are elaborated on in Section 4.1.5.

Expanding the active space from cas(5,5) to cas(9,12) in CASSCF/NEVPT2 calculations has a noteworthy impact on the g-tensor (Table 2), especially visible in the g₁₁ eigenvalue, and comparable to the differences between the base CASSCF and the NEVPT2 calculation based on it. It is seen that the improved level of theory does not translate into a notable effect on NMR shieldings, which demonstrates their insensitivity to the quality of the g-tensor due to the relatively small range of variation of the g-tensor components compared to the HFC. The improvement of the g-tensor calculations has no meaningful impact on the question of agreement between theory and experiment, or signal shift assignment in this molecule.

Table 2 Computed ZFS (spin–orbit contribution only) parameters as well as g-tensor eigenvalues and the isotropic g-factor at both the CASSCF level and the NEVPT2 level based on CASSCF wavefunctions. The ZFS parameters D and E/D are explained in the ESI. Geometry optimised in vacuo has been used

System/level of theory	S	D (cm^−1)	E/D	g 11	g 22	g 33	g iso
a Doublet systems have no ZFS. b High-spin multiplet.
[Tp2Fe]^+^a	1/2
cas(5,5), CASSCF				0.074	0.897	4.594	1.855
cas(5,5), NEVPT2				0.507	0.990	4.378	1.958
cas(9,12), CASSCF				0.225	0.880	4.289	1.798
cas(9,12), NEVPT2				0.340	0.871	4.290	1.834
Tp2Fe^b	2
cas(6,5), CASSCF		9.68	0.0089	1.994	2.140	2.143	2.092
cas(6,5), NEVPT2		8.69	0.0088	1.996	2.111	2.113	2.074
cas(10,12), CASSCF		8.44	0.0088	1.995	2.120	2.122	2.079
cas(10,12), NEVPT2		8.25	0.0086	1.997	2.106	2.108	2.070

---

4.1.2 Solvation model in HFC calculations. In Tables 3 and 4, NMR shifts have been calculated using the C-PCM continuum solvation model for the HFCs. Due to software limitations, these calculations were performed using the scalar relativistic DKH2 approximation on the ORCA software, instead of the fully relativistic mDKS method; in vacuo mDKS and DKH results are included for comparison. The HFC tensors extracted from the DKH calculations do not include spin–orbit contributions. The same basis set was employed as for the mDKS calculations, with a decontracted Dyall-VTZ basis set placed on the metal center and a decontracted DKH-def2-TZVP basis set used for the ligand atoms.

Table 3 Experimental^18,19 and computational ¹H NMR chemical shifts (ppm) for [Tp₂Fe]⁺ with HFCs calculated at the in vacuo mDKS, DKH and DKH + solvent(C-PCM) levels, as well as with 12 explicit acetone solvent molecules using the PBE0 functional. The g-tensor was calculated in vacuo at the NEVPT2 level based on the cas(9,12) wavefunction. Geometry optimised in vacuo was used in g-tensor calculations. In counterion and explicit solvent HFC calculations, the relevant geometry was used; otherwise in vacuo unless differently indicated. Temperature for ¹H: 305 K

HFC level of theory	H3	H4	H5	BH
a Acetone-d6 (ε = 20.7). b COSMO-optimised geometry used in HFC calculation.
Exptl^18,19	−52.9	−13.2	−8.2	39.4
mDKS in vacuo	−50.9	−9.4	0.3	42.0
DKH in vacuo	−43.9	−9.1	−2.6	36.6
DKH + counterion [PF6]^−	−41.9	−11.5	−3.7	34.4
DKH + solvent (C-PCM)^a	−44.3	−8.9	−2.5	36.6
DKH + solvent (C-PCM)^a^b	−43.8	−9.3	−2.3	36.4
DKH + 12 acetones	−41.8	−7.9	−4.7	34.6

---

Table 4 As Table 3 but for ¹³C and ¹¹B chemical shifts (ppm). Temperature for ¹³C: 300 K; ¹¹B: 298 K

HFC level of theory	C3	C4	C5	B
a Chloroform (ε = 4.81). b Dichloromethane (ε = 8.93). c Acetone-d6 (ε = 20.7). d COSMO-optimised geometry used in HFC calculation.
Exptl^18,19	2.5	175.8	19.6	95.2
mDKS in vacuo	34.6	2.5	91.6	92.5
DKH in vacuo	61.2	2.9	88.5	78.5
DKH + counterion [PF6]^−	58.2	8.5	96.4	74.1
DKH + solvent (C-PCM)	61.7^a	4.3^a	89.9^a	78.9^b
DKH + solvent (C-PCM)^c	60.3	2.5	89.0	78.9
DKH + solvent (C-PCM)^c^d	70.3	3.8	102.3	78.5
DKH + 12 acetones	68.1	−2.0	106.0	75.1

---

Calculations with the presence of twelve explicit acetone molecules in the hyperfine computations have been included to illustrate the difference to the C-PCM continuum solvation model. This structure was optimised using the PBE functional to reduce the computational cost. For carbon and boron, an additional C-PCM calculation using acetone as solvent is shown for valid comparison between the different computational methods.

The scalar relativistic DKH versus fully relativistic mDKS hyperfine calculations are in greater contrast for this molecule than solvent effects are to in vacuo computations, with differences of several ppm for hydrogen, over 26 ppm for C3 and 14 ppm for boron. For H3 and B, mDKS is clearly in better agreement with experiment, while for H5, DKH is closer to measurements.

The C-PCM model results in very small deviations from the non-solvated DKH model, of no more than 1 ppm even for carbon. Non-specific electrostatic effects approximated by the C-PCM model thus appear to be unimportant for this system. In contrast, moving to an explicit solvation model has a notably larger impact of a couple ppm for hydrogen, and up to 17 ppm for carbons C3 and C5. Therefore, in this case, the continuum solvation model is not sufficient to capture the most important solvent effects. If inclusion of the solvent is desired, a more accurate model is necessary. In some cases, such as H5, explicit solvation improves compatibility with experiment, whereas in others compatibility is worsened. In the latter cases, the non-solvated calculations likely again benefit from error cancellation.

Furthermore, a calculation in which HFCs were computed using a COSMO-optimised geometry is included for completeness. The impact on proton and boron chemical shifts is small, with results changing by less than 1 ppm. For carbon, we see somewhat more substantial changes of the order of 10 ppm for C3 and C5. For these nuclei, the discrepancies between the C-PCM model and explicit solvation are alleviated, suggesting that the exact geometry may have some importance, although such an alleviation is not consistent across all the other nuclei.

Finally, hyperfine couplings have been calculated otherwise in vacuo but in the presence of the counterion [PF₆]⁻ used in the measurements. The counterion changes chemical shifts roughly by 2 ppm for hydrogen, and by a few ppm for carbon and boron. The effect is stronger than that of the C-PCM model, although explicit solvation by 12 solvent molecules remains overall more impactful. In calculations requiring great accuracy, including the counterion effect may be meaningful.

Overall, computations underestimate the proton H3, H4 and H5 shifts across the board, and solvent considerations provide no systematic improvement.

4.1.3 Exact exchange. Paramagnetic NMR shieldings within the Kurland–McGarvey theory are significantly more sensitive to the quality of the HFC computations than to that of the ZFS and g-tensor computations;⁶³ their insensitivity to the g-tensor was already seen in Section 4.1.1. To demonstrate the issue, the hyperfine couplings were calculated at the fully relativistic mDKS level with different exact exchange admixtures in the DFT functional (0%, 25%, 40% and 100%), with results displayed in Fig. 2. Full numerical results can be found in Table S8 in the ESI.† In Fig. 2, the black bars represent the range of NMR shifts obtained using different exact exchange admixture parameters. Results calculated with the more conventional choices of 25% and 40% have been explicitly marked. In addition, shifts using hyperfine calculations at the DKH + C-PCM solvation level are shown for comparison.

Fig. 2 Experimental^18,19 and computational NMR chemical shifts for [Tp₂Fe]⁺. Bars indicate the range of in vacuo mDKS results for hyperfine coupling tensors obtained with hybrid, PBE0-based exchange–correlation functionals with different exact exchange admixtures (0%, 25%, 40% and 100%).

With the functional PBE-100 with 100% exact exchange it should be noted that spin contamination in hyperfine DFT calculations may be significant and may distort results, although DFT computations at the DKH level do not indicate major spin contamination for either [Tp₂Fe]⁺ or Tp₂Fe based on deviation from the theoretical expectation value 〈S²〉. Here the full 0–100% range has been included for illustrative purposes, but of greater practical interest is the 25–40% range.

The mDKS results with varying proportions of exact exchange cover a range of up to 40 ppm for hydrogen and up to 300 ppm for carbon, although PBE-25 and PBE-40 are closer together. This kind of large range of results constitutes a major problem with employing DFT in the computation of HFCs, since the amount of exact exchange is based on empirical considerations – there is no theoretical justification for, e.g., choosing PBE-25 over the 40% exact exchange admixture. The latter was pragmatically recommended by Kaupp et al.⁶³

Pure DFT functionals are known to overestimate the covalent character of ligand–metal bonding especially for open-shell metal complexes with large ligand contact shifts, causing spin density to extend too far from the paramagnetic metal center to the ligand atoms.^71,72 This spin delocalisation error can be reduced by hybrid functionals containing exact Hartree–Fock exchange energy. Conversely, however, an overly localised electronic structure may underestimate the covalent nature of metal–ligand interactions and likewise lead to an unrealistic picture of the electronic structure. In addition, hybrid functionals containing a high admixture of exact exchange tend to overestimate valence-shell spin polarisation, which may lead to significant spin contamination.⁷³ On the other hand, hybrid functionals are also needed to satisfactorily capture core–shell spin polarisation, which is important in the computation of core properties such as hyperfine couplings. We see the expected trend in spin delocalisation in Fig. 3, where we have plotted the spin density isosurfaces of [Tp₂Fe]⁺ at the UKS/DKH level of theory with 0%, 25%, 40% and 100% exact exchange. Increasing the exact exchange admixture causes positive spin density isosurfaces to slightly shrink, indicating decreasing spin delocalisation. Simultaneously, the negative spin density isosurfaces clearly increase, possibly reflecting growing spin polarisation and eventual spin contamination problems.

Fig. 3 Spin density isosurfaces for [Tp₂Fe]⁺ at the UKS/DKH level with 0%, 25%, 40% and 100% exact exchange admixture in the hybrid, PBE0-based exchange–correlation functional. Positive isosurfaces are drawn in red and negative in blue. Isovalues of ±0.0001 used throughout.

In our computed NMR chemical shift results, 25% or 40% exact exchange typically reproduce experimental shifts better than 0% or 100%. This is in line with the expectation that a mid-range portion of exact exchange should provide a more balanced description of spin delocalisation and spin polarisation effects for a global hybrid functional than either of the two extremes.

4.1.4 DLPNO-CCSD. DLPNO-CCSD presents itself as a potential alternative ab initio method for computing HFCs that does not suffer from the ambiguity of selecting the exact exchange admixture, present in the DFT calculations. For the proton shifts, the DLPNO-CCSD-based HFCs lead to an overall satisfactory agreement with experiment, sufficient, e.g., for unambiguous signal assignment (Fig. 2). For signals other than H4, the DLPNO-CCSD data appear within the DFT range of results, closest to the shifts obtained with the pure PBE functional. For H4, DLPNO-CCSD produces a shift between the DFT data points corresponding to 40% and 100% exact exchange admixtures. For the ¹³C and ¹¹B data, DLPNO-CCSD again results in shifts mostly within the DFT range, and at different positions within that range depending on the particular nucleus. The comparison with experimental ¹³C shifts is discussed in the next subsection.

DLPNO-CCSD produces overall slightly worse shifts in comparison to PBE-25 or PBE-40 for this system. Only the BH signal is reproduced better by DLPNO-CCSD. In this comparison one must bear in mind that other factors, e.g., the presently omitted rovibrational effects⁷⁰ may disfavour the ab initio data in this case. The fact that the DLPNO-CCSD data are quasirelativistic, as opposed to the fully relativistic nature of the DFT data, also contributes to the differences.

4.1.5 Discrepancy in ¹³C shifts. The shift results calculated for [Tp₂Fe]⁺ are in reasonably good agreement with experiment for hydrogen and boron. However, the calculations fail to reproduce experimental carbon shifts, as discussed in our previous work.¹⁸ Neither improving the CASSCF active space, nor accounting for either the solvent effects nor the counterion can explain these discrepancies. The most obvious possible explanation for this would be an error in signal assignment – re-assigning the experimental signals cyclically as C3 → C4, C4 → C5, C5 → C3 would greatly improve compatibility of the DFT range of calculated results with experiment. Alternatively, C4 and C5 could simply be interchanged, particularly if one focuses on the shifts calculated with DLPNO-CCSD-based HFC data. However, this is not supported by the experimental linewidth data, as also previously noted¹⁸ (though this is not a definite refutation of the re-assignment argument). Even after re-assignment, however, agreement with experiment would be unsatisfactory in the quantitative sense for C3 and C5, although compatibility for C3 would likely be improved if an mDKS + explicit solvation calculation were possible.

Another plausible explanation is in contributions from a low-lying excited state 601.1 cm⁻¹ above the ground state, possessing the same doublet multiplicity as in the ground state multiplet, predicted by cas(5,5)/NEVPT2 calculations in our earlier work.¹⁸ Our newest cas(9,12)/NEVPT2 computations predict the excited doublet 538.8 cm⁻¹ above the ground state (see below). This explanation was speculated upon in ref. 18 by very approximate calculations based on a fitting procedure to obtain hypothetical NMR shieldings for the excited state and Boltzmann-weighted thermal occupations of the ground and excited multiplets. The possibility was further explored in the present work by estimating the NMR shieldings of the excited doublet and fitting a thermal equilibrium of the ground and excited doublets to experimental data. An improved agreement with experiment was not achieved, however, possibly due to the inadequacy of the employed approximations. See the section “Excited state of [Tp₂Fe]⁺” in the ESI† for more details. A more likely explanation is that a re-assignment of the experimental ¹³C signals is appropriate. In particular, both DFT and DLPNO-CCSD calculations seem to support re-assignment.

The physical contributions to the nuclear shieldings of [Tp₂Fe]⁺ in the ground doublet state, following an explanation of the various terms in Table S11 in the ESI,† have been included in Table S12 (and correspondingly for the excited quintet of Tp₂Fe in Table S13, ESI†). Across the series of calculations, the shieldings in [Tp₂Fe]⁺ are dominated by the contact and pseudocontact terms (1 and 9 in the numbering of the table), with smaller contributions by the “con,3” term number 6. The prominence of the pseudocontact and “con,3” terms arises from the significant g-tensor anisotropy of the system.

4.2 Spin-crossover compound Tp₂Fe

Measurements and calculations were carried out at the temperatures 290 K (¹H), 305 K (¹³C) and 298 K (¹¹B and ¹⁴N). In addition, a temperature series from 290 K to 455 K has been recorded for the proton shifts. The full set of results is tabulated in Tables S14 and S15 in the ESI.†

It should be noted that our results correspond to the equilibrium between the S = 0 ground and S = 2 excited multiplets. The possibility of an excited triplet (S = 1) multiplet receiving significant thermal occupation in addition to the ground singlet and excited quintet was considered, but ruled out on the basis that the triplet was found computationally to lie far higher above the ground state than the quintet (see Section 4.2.7).

4.2.1 CASSCF active space. The impact of expanding the CASSCF active space in calculating the g-tensor and the ZFS for the Tp₂Fe system is shown in Table 5. The differences between these two levels of theory, 5- and 12-orbital CASSCF, are even smaller than in the case of [Tp₂Fe]⁺, with changes of the order of 0.1 ppm for hydrogen and 1 ppm for carbon. Both are reasonably well in line with experimental measurements, although the carbon results are not stellar. We only show the results using the global ¹H Boltzmann fit, but the same conclusions can be drawn using the other fitting methods (numerical results available in Tables S14 and S15 in the ESI†).

Table 5 Experimental^18,19 and computational ¹H, ¹³C, ¹¹B and ¹⁴N NMR chemical shifts (ppm) for Tp₂Fe with g-tensor and ZFS calculated at the NEVPT2 level based on cas(6,5) and cas(10,12) wavefunctions. HFCs calculated at the mDKS/PBE0 level in vacuo. Global Boltzmann fit performed using the ¹H temperature series

	H3	H4	H5	BH
Exptl^18,19	13.0	12.9	8.1	−2.7
cas(6,5)	15.3	15.5	9.1	−2.0
cas(10,12)	15.1	15.5	9.2	−1.8

---

	C3	C4	C5	B	N1
Exptl^18,19	266	168	229	−26.4	−126
cas(6,5)	341.3	231.9	332.5	−37.5	−143.9
cas(10,12)	342.4	232.8	333.9	−37.0	−143.4

---

As the NMR results suggest, the g-tensor in the high-spin (S = 2) state is only a little changed by the cas(6,5) to cas(10,12) expansion (Table 2), unlike in the case of [Tp₂Fe]⁺. Likewise, the ZFS D-parameters produced by the two active spaces are quite similar at the NEVPT2 level, although at the CASSCF level a somewhat larger change is seen. The ZFS of commonly occurring magnitude has previously been reported to have negligible impact on NMR shifts at ambient temperatures – its importance is significant at very low temperatures and in cases where pseudocontact mechanisms are dominant in the total shift.^12,17,49 The ZFS and the g-tensor determine the magnetic susceptibility,^74,75 which can be used to formulate pseudocontact shifts.¹

4.2.2 Solvation model in HFC calculations. Application of the C-PCM solvation model has a very minor effect on the NMR shifts in comparison to moving from fully relativistic to scalar relativistic calculations (Tables 6 and 7), which itself changes the results only on the order of 1 ppm for hydrogen, and a few ppm for carbon. In this case, there is no systematic difference between the fully relativistic mDKS and scalar relativistic DKH calculations.

Table 6 Experimental^18,19 and computational ¹H NMR chemical shifts (ppm) for Tp₂Fe with HFCs calculated at the in vacuo mDKS, DKH and DKH + solvent(C-PCM) levels using the PBE0 functional. The ZFS and g-tensor were calculated in vacuo at the NEVPT2 level based on the cas(10,12) wavefunction. Geometry optimised in vacuo was used in all calculations. Global Boltzmann fit was performed using the ¹H temperature series. Temperature for ¹H: 290 K

HFC level of theory	H3	H4	H5	BH
a Chloroform-d (ε = 4.81).
Exptl^18,19	13.0	12.9	8.1	−2.7
mDKS in vacuo	15.1	15.5	9.2	−1.8
DKH in vacuo	14.1	15.3	9.5	−0.8
DKH + solvent (C-PCM)^a	14.0	15.2	9.5	−0.8

---

Table 7 As Table 6 but for ¹³C, ¹¹B and ¹⁴N chemical shifts (ppm). Temperature for ¹³C: 305 K; ¹¹B and ¹⁴N: 298 K

HFC level of theory	C3	C4	C5	B	N1
a Toluene (ε = 2.38).
Exptl^18,19	266	168	229	−26.4	−126
mDKS in vacuo	342.4	232.8	333.9	−37.0	−143.4
DKH in vacuo	345.5	235.8	340.8	−34.3	−137.9
DKH + solvent (C-PCM)^a	346.6	234.1	336.4	−34.2	−137.4

---

4.2.3 Signature of spin-crossover in the ¹H shifts. In Fig. 4(a), the ¹H NMR temperature series of Tp₂Fe has been calculated in vacuo at the cas(10,12) + NEVPT2/mDKS level for the pure paramagnetic quintet (S = 2) state, the pure diamagnetic singlet (S = 0) state, and the high-spin/low-spin equilibrium state resulting from Boltzmann averaging (eqn (6)) with fitted ΔH and ΔS parameters. The singlet-state chemical shifts do not depend explicitly on temperature (temperature dependence would be introduced through molecular dynamics, which we currently neglect), whereas the quintet shifts exhibit an explicit Curie-type 1/T temperature dependence. The experimental temperature dependence is non-monotonic, resulting from the excited paramagnetic S = 2 state having an increasingly large influence as the paramagnetic state occupation increases with temperature. At lower temperatures, the occupation of the paramagnetic state diminishes, and the experimental shifts approach the diamagnetic S = 0 shifts. We see this behaviour in our results, where the pure quintet state produces approximately correct results at high temperatures, but diverges from the experimental data points at lower temperatures, whereas the equilibrium state reproduces a qualitatively correct temperature dependence throughout the temperature series. This demonstrates that it is crucial to take the spin-crossover effect into account, especially at temperatures at which both spin states are expected to have significant occupation.

Fig. 4 Experimental^18,19 and computational ¹H temperature series for Tp₂Fe at the cas(10,12) + NEVPT2/mDKS/PBE0 level with (a) a global fit to ¹H temperature series (b) individual proton fits to ¹H temperature series and (c) a global fit to all experimental data. Dashed lines: S = 2 state; dotted lines: S = 0 state; solid lines: high-spin/low-spin equilibrium state. Symbols indicate experimental data.

4.2.4 Exact exchange. Variation of the proportion of exact exchange in the DFT functional with a global fit to the proton temperature series (Fig. 5(a)) produces quite a wide range of NMR shifts, with proton shifts covering a range of up to almost 10 ppm and carbon shifts a range of approximately 100 ppm. We see that neither PBE-25 nor PBE-40 produces systematically better results than the other. In C3 and C5, both are quite far from the experimental shift.

Fig. 5 Experimental^18,19 and computational NMR chemical shifts for Tp₂Fe with (a) a global fit to ¹H temperature series, (b) proton fits to individual ¹H temperature series, (c) a global fit to all experimental data and (d) experimental ΔH and ΔS (ref. 18). Bars indicate the range of in vacuo mDKS results for hyperfine coupling tensors obtained with hybrid, PBE0-based exchange–correlation functionals with different exact exchange admixtures. Temperature for ¹H: 290 K; ¹³C: 305 K; ¹¹B and ¹⁴N: 298 K.

In the spin density isosurface plots as a function of exact exchange for the high-spin (S = 2) state of Tp₂Fe in Fig. 6, we see more clearly the same trends regarding spin delocalisation as in the case of [Tp₂Fe]⁺. Here, however, the positive spin density extends much further from the metal center than the negative, with a more balanced distribution appearing only with 100% exact exchange. This is explained by the lack, in this d⁶ system, of vacant spin-up acceptor 3d orbitals on the metal center, causing an unbalanced ligand–metal back-bonding, similarly to what was discussed by Pritchard and Autschbach⁷¹ in the cases of Cr(acac)₃ and Fe(acac)₃. In contrast, the d⁵ system [Tp₂Fe]⁺ has two empty metal 3d orbitals, which allows a more balanced ligand–metal back-bonding as shown by the spin density plots of that system (Fig. 3).

Fig. 6 Spin density isosurfaces for the high-spin (S = 2) state of Tp₂Fe at the UKS/DKH level with 0%, 25%, 40% and 100% exact exchange admixture in the hybrid, PBE0-based exchange–correlation functional. Positive isosurfaces are drawn in red and negative in blue. Isovalues of ±0.0001 used throughout.

The physical contributions to the nuclear shieldings for the high-spin (S = 2) state of Tp₂Fe are listed in Table S13 in the ESI.† We have reported these at the temperature of 455 K, where the high-spin state is the primary contributor to the equilibrium state. In all calculations, we see the most significant contributions arising from the contact, dipolar and pseudocontact terms 1, 2 and 9 (see ESI† for numbering). The latter two are a consequence of the appreciable ZFS (D-parameter of around 9 cm⁻¹) and g-tensor anisotropy, respectively. In ¹³C shieldings, the contact term dominates despite notable additional contributions by the “con,3” term 6, which originates from an isotropic g-shift and contact hyperfine coupling.

4.2.5 DLPNO-CCSD. In Fig. 7, the proton temperature series are compared with a global fit to the ¹H temperature series with HFCs computed at the mDKS/PBE0 level (a) and at the DLPNO-CCSD level (b). The DLPNO-CCSD results disagree on the order of H3 and H4 with both mDKS and experiment. This is evidently caused by a worse reproduction of H4 by DLPNO-CCSD at high temperatures. On the other hand, the DLPNO-CCSD curve of H3 is further from experiment than mDKS at low temperatures. Conversely, for BH, DLPNO-CCSD reproduces the measured temperature series better than mDKS at low temperatures, and for H5, throughout the temperature range. Looking at the situation at the main measurement temperature (Fig. 5(a)), DLPNO-CCSD reproduces the experimental signal shifts for most nuclei better than or as well as DFT. A notable exception is C3, for which the measured shift falls outside the range of all computational results.

Fig. 7 Experimental^18,19 and computational ¹H temperature series for Tp₂Fe with a global fit to ¹H temperature series (a) at the cas(10,12) + NEVPT2/mDKS/PBE0 level and (b) at the cas(10,12) + NEVPT2/DLPNO-CCSD level. Dashed lines: S = 2 state; dotted lines: S = 0 state; solid lines: high-spin/low-spin equilibrium state. Symbols indicate experimental data.

4.2.6 Boltzmann averaging. The above-discussed global proton temperature series fitted with data based on DFT-calculated HFCs produce the enthalphy and entropy parameters ΔH in the range of approximately 20–30 kJ mol⁻¹ and ΔS in the range 40–90 J K⁻¹ mol⁻¹ (Table S14, ESI†), respectively, depending on the level of theory. These are to be compared to the experimental values of 16.1 kJ mol⁻¹ and 47.7 J K⁻¹ mol⁻¹, respectively, in ref. 26 obtained from temperature-dependent measurements of the magnetic susceptibility in solution. An improved match is obtained with the newer experiments of ref. 18, based on the temperature series of H4 and yielding the values 24 kJ mol⁻¹ and 70 J K⁻¹ mol⁻¹, respectively. Using HFCs from DLPNO-CCSD calculations produces the higher values of 48.1 kJ mol⁻¹ and 154.8 J K⁻¹ mol⁻¹. It should be noted that these experimental ΔH and ΔS are not as accurate as the primary NMR measurements. The Gibbs free energy gaps ΔG = ΔH − TΔS have also been calculated in Table S14 (ESI†).

In, instead, fitting individual equivalent groups of protons to ¹H temperature series (Fig. 5(b)), we see a narrower range of DFT results with varying exact exchange for hydrogen, but not for other nuclei, as expected due to the proton-focused fitting method. In the cases of C3 and C5, the computations performed for HFCs at the PBE-25 and PBE-40 levels yield NMR shifts nearly at the opposite ends of the range of DFT results. The experimental data point at C5 again falls rather far outside the range of computational results. In this case, the obtained DLPNO-CCSD fitting parameters were dependent on the initial values of the fitting algorithm, and no meaningful shifts could be obtained. These results have therefore been omitted in Tables S14 and S15 in the ESI† as well as Fig. 5.

In the individual proton fits to ¹H temperature series, the fitted average ΔH and ΔS fall in roughly the same range of 20–30 kJ mol⁻¹ and 40–90 J K⁻¹ mol⁻¹, respectively, as in the global ¹H temperature series fit (Table S14, ESI†).

By fitting computational parameters globally to all experimental data, i.e. besides the ¹H temperature series, also the ¹³C, ¹¹B and ¹⁴N signal shift measurements (Fig. 5(c)), we get the best agreement with experiment for carbon shifts, but at the expense of worse performance for proton shifts. This is a natural outcome, since we now move away from a proton-focused fitting technique and take experimental carbon results into account as well. DLPNO-CCSD is again overall competitive with DFT. DFT produces clearly better results only for H4.

With this all-data fit, ΔH and ΔS become larger than with the other two fitting methods, with typical ranges of 30–40 kJ mol⁻¹ and 100–120 J K⁻¹ mol⁻¹, respectively (Table S14, ESI†), although again it should be remembered that here the experimental data is not necessarily of the greatest accuracy. The DLPNO-CCSD fit again stands as an outlier with ΔH of 63.2 kJ mol⁻¹ and ΔS of 196.5 J K⁻¹ mol⁻¹.

Finally, the NMR shifts of Tp₂Fe have been thermally averaged without a fitting procedure, by using the experimental values of ΔH and ΔS from ref. 18, in Fig. 5(d). Here we, again, see a quite large range of the mDKS results as a function of exact exchange admixture. Agreement with experiment is overall comparable to the ¹H temperature series fit (a). DLPNO-CCSD is better than or comparable with DFT for most nuclei, with the exceptions of H4 and BH. The NMR shifts have additionally been calculated using the experimental ΔH and ΔS of ref. 26 (Fig. S6 in the ESI†). With these alternative parameters, agreement with experiment is notably worse in most cases, particularly for the ¹³C shifts.

The proton temperature series calculated with each of the three fitting methods are shown side by side in Fig. 4. The best agreement with measurements is gained with individual proton fits (b), though this also results in weaker performance for other nuclei. A global fit to all experimental data (c) yields the poorest ¹H temperature series in comparison to experiment but reproduces carbon shifts better, as mentioned before. Fitting globally to the proton temperature series (a) produces an overall good agreement with experiment as an intermediate form between (b) and (c). All in all, each of the three fitting methods is quite well suited for qualitatively understanding NMR shifts, and the choice of the method can be guided by the availability of experimental data and the primary information targeted.

4.2.7 Energetics. We have performed additional calculations to evaluate the energy gap between the high-spin and low-spin states of the spin-crossover complex Tp₂Fe at the state-specific CASSCF/NEVPT2 and DLPNO-CCSD(T) levels of theory (Table S16 in the ESI†), following the work done by Drosou, Mitsopoulou and Pantazis²⁴ for manganese spin-crossover complexes. The CASSCF wavefunctions were calculated at both cas(6,5) and cas(10,12) levels. The DLPNO-CCSD(T) calculations were performed with the “semi-canonical” perturbative triples correction approximation (T0) instead of the more accurate, iterative perturbative triple excitations (T1) recommended by Drosou et al.²⁴ due to computational limitations. A locally dense def2-SVP/def2-TZVP basis set was used, as well as a larger locally dense def2-SVP/def2-QZVPP basis set, for comparison. The geometry used in each calculation was optimised (as described above) for the appropriate spin state.

It is clear that CASSCF yields very poor energies in all cases, with energy gaps ΔH = H₂ − H₀ of the wrong sign, implying S = 2 ground state, on the order of around −200 kJ mol⁻¹ compared to the experimental values of 16.1 kJ mol^−1 26 or 24 kJ mol⁻¹.¹⁸ Introducing dynamic electron correlation via NEVPT2 significantly improves the results, with the calculation based on a cas(10,12) wavefunction with the SVP/QZVPP basis set yielding an energy gap of 13.1 kJ mol⁻¹, approaching the range of experimental results. The best results are obtained with DLPNO-CCSD(T), giving a result of 14.6 kJ mol⁻¹ with the SVP/QZVPP basis set. The choice of basis set is clearly significant, as enlarging it around the metal center from triple-zeta to quadruple-zeta is required at the cas(10,12) + NEVPT2 and DLPNO-CCSD(T) levels to obtain an energy gap of the correct sign. Overall, DLPNO-CCSD(T) performs better than NEVPT2, although the agreement with experiment (ref. 18 and 26) cannot be considered quantitative. Ref. 24 found CASSCF to be a poor method for calculating energy gaps in their spin-crossover molecules even with a NEVPT2 correction, preferring DLPNO-CCSD(T). Our findings are largely in line with this, although we get a fairly satisfactory NEVPT2 result with the larger cas(10,12) CASSCF active space.

In addition, the entropy gap from the harmonic vibrational frequencies was computed with TURBOMOLE at room temperature (298 K) and normal pressure (0.1 MPa) to be ΔS_vib = S_2,vib − S_0,vib = 73.9 J K⁻¹ mol⁻¹, while the electronic contribution arising from the five-fold degenerate quintet is 13.4 J K⁻¹ mol⁻¹.⁷⁶ The sum of these contributions is 87.3 J K⁻¹ mol⁻¹, which is somewhat higher than the experimental result in ref. 18, but consistent with the higher end of the range of our fitted results.

At the DLPNO-CCSD(T0) level of theory, with a locally dense def2-SVP/def2/QZVPP basis set, the lowest triplet multiplet was computed to be 93.0 kJ mol⁻¹ above the ground state, compared to 14.6 kJ mol⁻¹ for the quintet. NEVPT2 based on a cas(10,12) wavefunction yields a similar result (Table S16, ESI†). The difference is considerable enough that the triplet is not expected to have significant thermal occupation, which justifies its omission from the present consideration of the thermal equilibrium.

5 Conclusions

We have conducted a study of the challenging iron complexes [Tp₂Fe]⁺ and Tp₂Fe, with attention paid to the computational factors influencing ¹³C signal assignment in the former and the spin-crossover behavior of the latter compound. The impact of a number of model improvements and computational modifications to the “standard method” of pNMR calculations have been considered: more accurate g-tensor and ZFS calculations by expanding the CASSCF active space, solvation effects via the implicit C-PCM method and, for [Tp₂Fe]⁺, explicit solvent molecules, use of the ab initio DLPNO-CCSD method to compute HFCs, as well as DFT calculation of HFCs with different proportions of exact exchange.

Of interest has especially been the discrepancy between theory and measurements of carbon shifts in [Tp₂Fe]⁺. For the spin-crossover compound Tp₂Fe, three fitting methods for the Boltzmann averaging between the diamagnetic ground state and the thermally occupied quintet excited state have been employed, the fitting parameters being the enthalpy and entropy gaps between the spin states: a global fit to ¹H temperature series, individual proton fits to ¹H temperature series, and global fit to all ¹H, ¹³C, ¹¹B and ¹⁴N experimental data.

Expanding the CASSCF active space from five to twelve active orbitals in the g-tensor and ZFS calculations has little effect on the NMR shieldings of the [Tp₂Fe]⁺ and Tp₂Fe systems due to the relatively small range of variation of the g-tensor components and the negligible impact of the ZFS at ambient temperatures. This does not appear to be a promising way to improve the accuracy of paramagnetic NMR calculations except in some special cases.

Comparison between C-PCM and explicit solvent calculations for [Tp₂Fe]⁺ reveals that the continuum solvation model is not sufficient to account for the important solvent effects in this system. Non-electrostatic effects that cannot be reproduced by continuum solvation evidently play a significant role. In accurate calculations accounting for solvent contributions, a more complex model is therefore needed.

A major issue with using DFT in the computation of hyperfine couplings is the resulting wide range of NMR shifts with different proportions of exact exchange. Even though it is possible to tune the exact exchange to improve computational results for certain nuclei or systems, no theoretical argument can be made for why, e.g., 25% should be preferred over 40%. The ab initio DLPNO-CCSD method was explored as a possible alternative, and the results seem promising, overall producing shifts comparable with DFT, albeit with significantly greater demands on computational resources.

The disagreement between theoretical and experimental ¹³C NMR shifts of the [Tp₂Fe]⁺ system encountered previously¹⁸ has been reviewed in the present work. The assignment of measured signal shifts appears to be the most likely cause despite the lack of support by linewidth data. However, the theoretical results are not fully satisfactory even with re-assignment and consideration of explicit solvation or counterion effects. Still, the fact that both DFT and DLPNO-CCSD calculations support re-assignment strengthens this point of view. A low-lying excited state is another candidate, which was tentatively supported by earlier rough calculations. We could not verify this in the present slightly more advanced, but still very approximate calculations of the contributions of the excited multiplet. There is currently no practical way to reliably compute the excited-state contributions.

With the spin-crossover system Tp₂Fe, of note is the observation that the paramagnetic quintet state alone cannot reproduce a qualitatively correct temperature behavior of the NMR chemical shifts – the singlet state occupation must be taken into account to obtain reasonable results at the lower end of the temperature range, including room temperature. All three fitting techniques used produce results qualitatively in line with signal shift measurements, although which data is used in the fitting procedure results in better or worse compatibility with experiment for different groups of nuclei.

Author contributions

A. Pyykkönen was responsible for data curation, general investigation, visualisation of methods and results, and writing of the original draft; J. Vaara supervised the entire work. Both authors contributed together to conceptualisation of the project and reviewing/editing of the manuscript.

Conflicts of interest

There are no conflicts to declare.

Acknowledgements

We thank Prof. Dr Frank H. Köhler (Munich) for useful discussions and comments. We are also thankful for the suggestions of one of the anonymous reviewers of the first submitted version of this paper. This work was funded by the Academy of Finland (project no. 331008), the University of Oulu (Kvantum Institute) and the Finnish Cultural Foundation. Computational resources were provided by the CSC – IT Center for Science (Espoo, Finland) and the Finnish Grid and Cloud Infrastructure project (persistent identifier urn:nbn:fi:research-infras-2016072533).

Notes and references

1. I. Bertini, C. Luchinat and G. Parigi, Solution NMR of Paramagnetic Molecules, Elsevier, Amsterdam, 2001.
2. C. P. Grey and N. Dupré, Chem. Rev., 2004, 104, 4493–4512.
3. M. Kaupp and F. H. Köhler, Coord. Chem. Rev., 2009, 253, 2376–2386.
4. P. H. Keizers and M. Ubbink, Prog. Nucl. Magn. Reson. Spectrosc., 2011, 58, 88–96.
5. A. J. Pell, G. Pintacuda and C. P. Grey, Prog. Nucl. Magn. Reson. Spectrosc., 2019, 111, 1–271.
6. M. Kaupp, M. Bühl and V. G. Malkin, Calculation of NMR and EPR Parameters: Theory and Applications, Wiley-VCH, Weinheim, 2004.
7. A. B. A. Andersen, A. Pyykkönen, H. J. Aa. Jensen, V. McKee, J. Vaara and U. G. Nielsen, Phys. Chem. Chem. Phys., 2020, 22, 8048–8059.
8. 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.
9. R. J. Kurland and B. R. McGarvey, J. Magn. Reson., 1970, 2, 286–301.
10. W. Van den Heuvel and A. Soncini, J. Chem. Phys., 2013, 138, 054113.
11. B. Martin and J. Autschbach, J. Chem. Phys., 2015, 142, 054108.
12. J. Vaara, S. A. Rouf and J. Mareš, J. Chem. Theory Comput., 2015, 11, 4840–4849.
13. Metal Ions in Biological Systems, ed. H. Sigel, Dekker, Basel, 1982, vol. 14.
14. D. Lopez-Tejedor, R. Benavente and J. M. Palomo, Catal. Sci. Technol., 2018, 8, 1754–1776.
15. A. Mondal and M. Kaupp, J. Phys. Chem. C, 2019, 123, 8387–8405.
16. A. Mondal and M. Kaupp, J. Phys. Chem. Lett., 2018, 9, 1480–1484.
17. B. Martin and J. Autschbach, Phys. Chem. Chem. Phys., 2016, 18, 21051–21068.
18. A. Pyykkönen, R. Feher, F. H. Köhler and J. Vaara, Inorg. Chem., 2020, 59, 9294–9307.
19. R. Feher, Dissertation, Technische Universität München, 1996.
20. P. Gütlich and H. Goodwin, in Spin Crossover in Transition Metal Compounds I, ed. P. Gütlich and H. A. Goodwin, Springer-Verlag Berlin, Heidelberger Platz 3, D 14197 Berlin, Germany, 2004, vol. 233 of Topics in Current Chemistry-Series, pp. 1–47.
21. L. Cambi and L. Szegö, Ber. Dtsch. Chem. Ges. (A and B Series), 1931, 64, 2591–2598.
22. J.-F. Létard, P. Guionneau and L. Goux-Capes, Towards Spin Crossover Applications, Spin Crossover in Transition Metal Compounds III. Topics in Current Chemistry, Springer, Berlin, Heidelberg, 2004, vol. 235.
23. K. S. Kumar and M. Ruben, Angew. Chem., Int. Ed., 2021, 60, 7502–7521.
24. M. Drosou, C. A. Mitsopoulou and D. A. Pantazis, Polyhedron, 2021, 208, 115399.
25. B. M. Flöser, Y. Guo, C. Riplinger, F. Tuczek and F. Neese, J. Chem. Theory Comput., 2020, 16, 2224–2235.
26. J. P. Jesson, S. Trofimenko and D. R. Eaton, J. Am. Chem. Soc., 1967, 89, 3158–3164.
27. J. P. Jesson, J. F. Weiher and S. Trofimenko, J. Chem. Phys., 1968, 48, 2058–2066.
28. B. Hutchinson, L. Daniels, E. Henderson, P. Neill, G. J. Long and L. W. Becker, J. Chem. Soc., Chem. Commun., 1979, 1003–1004.
29. F. Grandjean, G. J. Long, B. B. Hutchinson, L. Ohlhausen, P. Neill and J. D. Holcomb, Inorg. Chem., 1989, 28, 4406–4414.
30. L. Salmon, G. Molnár, S. Cobo, P. Oulié, M. Etienne, T. Mahfoud, P. Demont, A. Eguchi, H. Watanabe, K. Tanaka and A. Bousseksou, New J. Chem., 2009, 33, 1283–1289.
31. A. A. Pavlov, G. L. Denisov, M. A. Kiskin, Y. V. Nelyubina and V. V. Novikov, Inorg. Chem., 2017, 56, 14759–14762.
32. F. Milocco, F. de Vries, I. M. A. Bartels, R. W. A. Havenith, J. Cirera, S. Demeshko, F. Meyer and E. Otten, J. Am. Chem. Soc., 2020, 142, 20170–20181.
33. J. Huang, R. Xie, Y. Hu, S. Lei and Q. Li, Chem. Phys. Lett., 2020, 758, 137925.
34. J. Villalva, A. Develioglu, N. Montenegro, R. Sánchez-de Armas, A. Gamonal, M. Gaarcia, M. Ruiz-González, J. S. Costa, C. Calzado, E. Pérez and E. Burzuri, Nat. Commun., 2021, 12, 1578.
35. S. K. Karuppannan, A. Martín-Rodríguez, E. Ruiz, P. Harding, D. J. Harding, X. Yu, A. Tadich, B. Cowie, D. Qi and C. A. Nijhuis, Chem. Sci., 2021, 12, 2381–2388.
36. H. Zhang, C. Niu, J. Zhang, L. Zou, Z. Zeng and X. Wang, Phys. Chem. Chem. Phys., 2021, 23, 9679–9685.
37. Z. Ke, L. E. Jamieson, D. M. Dawson, S. E. Ashbrook and M. Bühl, Solid State Nucl. Magn. Reson., 2019, 101, 31–37.
38. Z. Ke, D. M. Dawson, S. E. Ashbrook and M. Bühl, Origin of temperature dependence of ¹³C pNMR shifts for copper paddlewheel metal-organic frameworks, 2021, poster, Center for Magnetic Resonance Day, University of St. Andrews, Scotland, June 3rd (2021).
39. W. Van den Heuvel and A. Soncini, Phys. Rev. Lett., 2012, 109, 073001.
40. F. Gendron, K. Sharkas and J. Autschbach, J. Phys. Chem. Lett., 2015, 6, 2183–2188.
41. F. Gendron and J. Autschbach, J. Chem. Theory Comput., 2016, 12, 5309–5321.
42. TURBOMOLE V7.2 2017, a development of University of Karlsruhe and Forschungszentrum Karlsruhe GmbH, 1989–2007, TURBOMOLE GmbH, since 2007, available from http://www.turbomole.com.
43. C. Adamo and V. Barone, J. Chem. Phys., 1999, 110, 6158–6170.
44. S. Grimme, J. Antony, S. Ehrlich and H. Krieg, J. Chem. Phys., 2010, 132, 154104.
45. S. Grimme, S. Ehrlich and L. Goerigk, J. Comput. Chem., 2011, 32, 1456–1465.
46. M. Dolg, U. Wedig, H. Stoll and H. Preuss, J. Chem. Phys., 1987, 86, 866–872.
47. J. M. L. Martin and A. Sundermann, J. Chem. Phys., 2001, 114, 3408–3420.
48. F. Weigend and R. Ahlrichs, Phys. Chem. Chem. Phys., 2005, 7, 3297–3305.
49. S. A. Rouf, J. Mareš and J. Vaara, J. Chem. Theory Comput., 2015, 11, 1683–1691.
50. S. A. Rouf, J. Mareš and J. Vaara, J. Chem. Theory Comput., 2017, 13, 3731–3745.
51. B. Roos, Chem. Phys. Lett., 1972, 15, 153–159.
52. C. Angeli, R. Cimiraglia, S. Evangelisti, T. Leininger and J.-P. Malrieu, J. Chem. Phys., 2001, 114, 10252–10264.
53. C. Angeli, R. Cimiraglia and J.-P. Malrieu, Chem. Phys. Lett., 2001, 350, 297–305.
54. C. Angeli, R. Cimiraglia and J.-P. Malrieu, J. Chem. Phys., 2002, 117, 9138–9153.
55. M. Douglas and N. M. Kroll, Ann. Phys., 1974, 82, 89–155.
56. B. A. Hess, Phys. Rev. A: At., Mol., Opt. Phys., 1986, 33, 3742–3748.
57. D. Ganyushin and F. Neese, J. Chem. Phys., 2006, 125, 024103.
58. D. Ganyushin and F. Neese, J. Chem. Phys., 2013, 138, 104113.
59. F. Neese, Wiley Interdiscip. Rev.: Comput. Mol. Sci., 2012, 2, 73–78.
60. F. Neese, Wiley Interdiscip. Rev.: Comput. Mol. Sci., 2018, 8, e1327.
61. D. A. Pantazis, X.-Y. Chen, C. R. Landis and F. Neese, J. Chem. Theory Comput., 2008, 4, 908–919.
62. E. Malkin, M. Repiský, S. Komorovský, P. Mach, O. L. Malkina and V. G. Malkin, J. Chem. Phys., 2011, 134, 044111.
63. S. Gohr, P. Hrobárik, M. Repiský, S. Komorovský, K. Ruud and M. Kaupp, J. Phys. Chem. A, 2015, 119, 12892–12905.
64. ReSpect, version 5.1.0 Relativistic Spectroscopy DFT program of authors M. Repisky, S. Komorovsky, V. G. Malkin, O. L. Malkina, M. Kaupp, K. Ruud, with contributions from R. Bast, U. Ekstrom, M. Kadek, S. Knecht, L. Konecny, E. Malkin, I. Malkin-Ondik, 2019, (see http://www.respectprogram.org).
65. K. G. Dyall and A. S. P. Gomes, unpublished, Available from the Dirac web site, http://dirac.chem.sdu.dk/basisarchives/dyall/.
66. 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.01, Gaussian, Inc., Wallingford CT, 2019.
67. J. Novotný, M. Sojka, S. Komorovsky, N. Marek and R. Marek, J. Am. Chem. Soc., 2016, 138, 8432–8445.
68. V. Barone and M. Cossi, J. Phys. Chem. A, 1998, 102, 1995–2001.
69. C. Riplinger and F. Neese, J. Chem. Phys., 2013, 138, 034106.
70. A. Jaworski and N. Hedin, Phys. Chem. Chem. Phys., 2022, 24, 15230–15244.
71. B. Pritchard and J. Autschbach, Inorg. Chem., 2012, 51, 8340–8351.
72. A. J. Cohen, P. Mori-Sánchez and W. Yang, Science, 2008, 321, 792–794.
73. C. J. Schattenberg, T. M. Maier and M. Kaupp, J. Chem. Theory Comput., 2018, 14, 5653–5672.
74. L. Benda, J. Mareš, E. Ravera, G. Parigi, C. Luchinat, M. Kaupp and J. Vaara, Angew. Chem., Int. Ed., 2016, 55, 14713–14717.
75. L. Lang, E. Ravera, G. Parigi, C. Luchinat and F. Neese, J. Phys. Chem. Lett., 2020, 11, 8735–8744.
76. G. Privault, J.-Y. Mevellec, M. Lorenc, B. Humbert, E. Janod, N. Daro, G. Chastanet, A. Subedi and E. Collet, Cryst. Growth Des., 2022, 22, 5100–5109.

---

Footnote

† Electronic supplementary information (ESI) available. See DOI: https://doi.org/10.1039/d2cp03721a

---