Magnetic exchange and valence delocalization in a mixed valence [Fe²⁺Fe³⁺Te₂]⁺ complex: insights from theory and interpretations of magnetic and spectroscopic data

Abstract

A mixed valence binuclear Fe^2.5+–Fe^2.5+ (Robin–Day Class III) transition metal complex, [Fe^2.5+μTe₂Fe^2.5+]¹⁻, composed of two FeN₂Te₂ pseudo-tetrahedral units with μ-bridging Te²⁻ ligands was reported to exist in an unprecedented S = 3/2 ground state (Nature Chemistry, https://doi.org/10.1038/s41557-021-00853-5). For this and the homologous complexes containing Se²⁻ and S²⁻, the Anderson-Hasegawa double exchange spin-Hamiltonian was broadly used to interpret the corresponding structural, spectroscopic and magnetic data. First principles multireference ab initio calculations are used here to simulate magnetic and spectroscopic EPR data; analysis of the results affords a rationale for the stabilization of the S = 3/2 ground state of the Fe₂ pair. Complete Active Space Self-Consistent Field (CASSCF) calculations and dynamical correlation accounted for by means of N-Electron Valence Perturbation Theory to Second Order (NEVPT2) reproduce well the g-factors determined from simulations of X-band EPR spectra. A crucial technical tool to achieve these results is: (i) use of a localized orbital formulation of the many-particle problem at the scalar-relativistic CASSCF step; (ii) choice of state averaging over states of a given spin (at the CASCI/NEVPT2 step); and (iii) accounting for spin–orbit coupling within the non-relativistic Born–Oppenheimer (BO) many-particle basis using Quasi-Degenerate Perturbation Theory (QDPT). The inclusion of the S = 5/2 spin manifold reproduced the observed increase in the magnetic susceptibility (χT) in the high temperature range (T > 100 K), which is explained by thermal population of the S = 5/2 excited state at energy 160 cm⁻¹ above the S = 3/2 ground state. Theoretical values of χT from experimentally reported data points in the temperature range from 3 to 30 K were further computed and analyzed using a model which takes spin–phonon coupling into account. The model considerations and the computational protocols of this study are generally applicable to any Class I/II mixed valence dimer. The work can potentially stimulate further experimental and theoretical work on bi- and oligonuclear transition metal complexes of importance to bioinorganic chemistry and life sciences.

---

1. Introduction

Mixed-valence transition metal compounds are of great prominence in bioinorganic chemistry. The oxygen evolving complex, a Mn³⁺₃Mn⁴⁺CaO₅ cluster catalyzing water splitting (2H₂O → O₂ + 4H⁺ + 4e⁻) in the cycle of photosystem II¹ of natural photosynthesis and the iron-molybdenum cofactor (FeMoco), a Fe²⁺Fe^2.5+₄Fe³⁺₂Mo³⁺S₉C complex (Fig. 1), facilitating formation of ammonia from atmospheric N₂ in a well-known, but yet not completely understood biological nitrogen fixation process,² are prominent examples, among many others. To achieve reliable information about the electronic and geometric structure of these oligonuclear open-shell spin systems, a combination of spectroscopic techniques – electron paramagnetic resonance (EPR), Mössbauer spectroscopy, element and spin-specific X-ray spectroscopies, and single X-ray crystal structural determinations, all aided by theoretical calculations,^2,3 have been employed with some initial success. For such big molecules, geometry optimizations and spin-spin coupling computed using broken symmetry Density Functional Theory^4,5 (DFT) was (and still is) the only practicable tool.^6,7Ab initio multireference approaches have been applied with some success to mono- and binuclear transition metal dimers,^8,9 but, because of the large system size of real systems, are impractical in computing and analyzing electronic structures of complexes with three or more spin centers and two or more unpaired electrons per site. In order to understand the electronic structure of such large clusters using a bottom up approach, one possible choice is to consider mono- and binuclear subunits of the big cluster by using diamagnetic substitution of all other centers. Taking for example the FeMoco complex (Fig. 1), diamagnetic substitution of Zn²⁺ and Ga³⁺ in place of Fe²⁺ and Fe³⁺, respectively, would allow for the determination of their magnetic tensors and of the spin-Hamiltonian parameters. These centers would then be coupled in a second step into a (composed) model Hamiltonian of the entire cluster. In an attempt to approach this goal experimentally, a series of binuclear complexes [L₂Fe₂Q₂]⁻ with Q₂ bis-μ bridging S²⁻, Se²⁻ and Te²⁻ ligands supported by β-diketiminate terminal ligands (L) have been recently reported and characterized by X-ray structural analysis, X-band EPR, Mössbauer spectroscopy and magnetometry (Fig. 2a).¹⁰ While the S²⁻ (1) and Se²⁻ (2) congeners feature two overlapping quadrupole doublets in a 1 : 1 ratio, reflecting valence localized formally high-spin Fe²⁺ (d⁶) and Fe³⁺(d⁵) centers (Robin–Day Classes I or II),¹¹ the Mössbauer spectrum of the Te²⁻ analogue (3, Fig. 2b) shows a single quadrupole doublet at T = 80 K, suggesting a fully valence delocalized formally Fe^2.5+–Fe^2.5+ complex (Robin–Day Class III).¹¹ This finding is further supported by the X-ray structure with two crystallographically equivalent iron sites. X-Band EPR spectra of frozen solutions of 3 (Fig. 2c) have been simulated employing a S = 3/2 spin Hamiltonian. From the interpretation of the EPR-spectrum using a S′ = 1/2 pseudo-spin, large anisotropic g-factors (g = 5.740, 1.950 and 1.515) have been determined (Fig. 2c).¹⁰ The low field EPR transition (corresponding to g = 5.740) was relatively sharp; the dependence of its Lorentzian linewidth above T = 10 K allowed estimation of the energy separation between the ground and first excited spin state assuming an Orbach relaxation mechanism. This afforded an upper limit for the axial zero-field splitting parameter: D ≤ 15 cm⁻¹.¹⁰ A very large orthorhombic |E/D| ratio (close to the limiting value of 1/3), was extracted from the simulation of the EPR spectrum.¹⁰ Variable-temperature SQUID magnetic susceptibility data of solid powder samples (Fig. 2d) has been interpreted using the Bleaney–Bowers equation for mixed valence complexes considering the magnetic exchange between the high-spin Fe²⁺ (S₁ = 2) and Fe³⁺ (S₂ = 5/2) centers (Heisenberg exchange Hamiltonian H_exch = −2JŜ₁·Ŝ₂) and a double exchange term, H_trans = −BT̂₁₂ accounting for the transfer of the extra electron in the 3d-shell of Fe²⁺ (d⁶) [described by B-the transfer (hopping) energy] to the half field open d-shell of Fe³⁺ (d⁵). Because of the equivalence of the two valence forms (Fe₁²⁺(d⁶)–Fe₂³⁺(d⁵) and Fe₁³⁺(d⁵)–Fe₂²⁺(d⁶)), each spin of the exchange coupled iron centers (S = |S₁ − S₂|, |S₁ − S₂| + 1, … |S₁ + S₂| (S = 1/2, 3/2, 5/2, 7/2 and 9/2) appears in pairs which are split into two sublevels by electron transfer. To describe and interpret the spectra and magnetic properties of these systems, the Anderson–Hasegawa double exchange model yielding eqn (1.1) for the dependence of the energy on the spin S has been extensively used.^12–15

E_±(S) = −JS(S + 1) ± B(S + 1/2) (1.1)

Fig. 1 Model complex for the iron-molybdenum cofactor (FeMoco) from a DFT geometry optimization of the resting S = 3/2 state.³ The Figure has been constructed here after deleting water and other neutral residues used for a correct account of the FeMoco cluster anchored to the protein. The cluster allows one to study the local electronic structures of the Mo and Fe centers and of the binuclear subunits -the three Mo_μS₂Fe, Fe_μS₂Fe and Fe-_μC/SFe dimer units in dependence on the oxidation state and the electron counts on each center − d²(Mo⁴⁺)/d³(Mo³⁺) and d⁵(Fe³⁺) and total cluster charges −1/−2.

Fig. 2 Experimental X-ray structure of the mixed valence complex Fe^2.5+-_μTe₂-Fe^2.5+ (a); its T = 80 K Mössbauer spectrum in zero magnetic field (b); continuous-wave X-band EPR spectrum (9.63 GHz) of a frozen solution at 3.6–8 K, collected under non-saturating microwave power conditions (black line) and simulations (red lines): pseudo-spin S′ = 1/2, g′ = 5.740, 1.950, 1.515 (solid red line) and a S = 3/2, spin-Hamiltonian, |D| ≫ hν, |E/D| = 1/3 (assuming g = g_e = 2.0023) and a value of D < 15 cm⁻¹, estimated from the temperature dependence of the Lorentzian linewidth (c), variable-temperature SQUID magnetic susceptibility of a solid sample (d); data have been reproduced in a modified form from J. T. Henthorn, G. E. Cutsail III, T. Weyhermüller and S. DeBeer, “Stabilization of intermediate spin-states in mixed-valence diiron dichalcogenide complexes”, Nat. Chem., https://doi.org/10.1038/s41557-021-00853-5.

The iron centers in 1, 2 and 3 couple antiferromagnetically (J < 0), leading to a dependence of the lower energy branch E₋(S)/|B| on the ratio |B/J| (Fig. 3) which shows, that S = 3/2 may become the ground state in the 3 < |B/J| < 5 range. From the simulation of χT vs. T, values for J and B of −200(25) and 750(110) cm⁻¹, respectively, have been reported.¹⁰ The model was refined by introducing vibronic coupling.¹⁰ In the linear coupling limit, the effect of vibronic coupling on the spin levels is quantified by the vibronic coupling energy λ = V²/K (eqn (1.2), but notice the presence of an erroneous factor of (1/2) in the square root of eqn (2) and somewhat different notations used in ref. 10), with V- the linear vibronic coupling, K- the harmonic force constant, and q- the displacement along the symmetry breaking vibrations q = (V/K)x₋(q-in Å, x₋-dimensionless, Fig. 4a). This effective normal mode q is not an eigenfunction of the vibrational Hamiltonian but a superposition of several normal modes, approximated in ref. 10 (and also here), by a uniform shift of the two μ-bridging tellurium ligands from the left to the right of the symmetric Fe^2.5+–Fe^2.5+ pair or vice versa,¹⁶ thus assisting localizations into Fe₁²⁺(d⁶)–Fe₂³⁺(d⁵) and Fe₁³⁺(d⁵)–Fe₂²⁺(d⁶) distorted structures, respectively.

E(S, ∓) = −JS(S + 1) + (1/2)λx₋²∓[λ²x₋² + B²(S + 1/2)²]^1/2 (1.2)

Fig. 3 Energies of spin levels E(S, −)/J versus the |B/J| ratio; in a 3 < |B/J| < 5 range the = 3/2 state is stabilized with respect to the other spin states.

Fig. 4 (a) Representation of the anion [L₂Fe_2μTe₂]⁻, the red arrows depict the asymmetric (symmetry breaking) mode (q); (b) adiabatic ground state (S = 3/2) and two lowest (S = 1/2, S = 5/2) excited states along q; The curves were plotted taking the values of J = −200 cm⁻¹ and B = 750 cm⁻¹ (from ref. 10) and the vibronic coupling energy λ = 480 cm⁻¹ (see Table 7).

Using the previously reported values, J = −200 cm⁻¹, B = 750 cm⁻¹ (B/J = 3.7),¹⁰ and λ = 480 cm⁻¹ (see Section 7; notice the 3.5 times larger reported value λ = 1675 cm⁻¹ ref. 10, which in turn leads to a S = 1/2 ground state), we illustrate in Fig. 4b the dependence of the lower energy branches E(S, −) on the nuclear coordinate q for the lowest S = 1/2, 3/2 and 5/2 spin states. This set of parameters correctly reproduce the stabilization of the S = 3/2 ground state in a range −0.95 < q < 0.95 with a single minimum corresponding to a non-distorted (symmetric) structure and an energy separation between the ground S = 3/2 and the lowest S = 1/2 excited state of Δ = [E(1/2, −) − E(3/2, −)]_q=0 = 150 cm⁻¹. This is comparable to the vibrational energy ħω = 140 cm⁻¹ and indicates the break-down of the Born–Oppenheimer approximation (BOA) for 3, which is only fulfilled if Δ ≫ ħω. It is the main goal of the present work to compute and analyze, using first principles multireference ab initio calculations, the electronic structure and magnetic and spectral properties of 3 going beyond the BOA by taking spin–phonon coupling into account.

2. Outline

Mixed valence dimers of transition metals represent a challenge to first principle (ab initio) quantum chemistry studies. Here we apply a computational protocol put forward by Domingo et al.¹⁷ and implemented in the program package ORCA.^18,19 Its application to 1 and 2 allowed for computing and reproducing, with a reasonable accuracy, the Heisenberg exchange coupling parameters.²⁰ In Section 3 the various steps according to the implementation of this protocol in the ORCA code and its application to 3 will be described. In this work we use a bottom-up approach toward the analysis and understanding of the electronic structure and magnetic properties of 3. We first characterize the high-spin Fe²⁺(d⁶) and Fe³⁺(d⁵) coordination centers by artificially eliminating exchange and electron transfer (Section 4). This is achieved using diamagnetic substitution of Fe³⁺ and Fe²⁺ by the closed shell Ga³⁺ and Zn²⁺ ions, respectively. In the same section Ab Initio Ligand Field Theory (AILFT)^21–23 is applied to characterize the iron-ligand bonding with special emphasis on the bridging tellurium ligands and the effect of low-symmetry on the local (single ion) zero-field splitting (ZFS), the D-tensor and the g-factors (g-matrix). In section 5 we employ results from scalar relativistic CASSCF²⁴ calculations to quantify the Heisenberg exchange within the Fe²⁺(d⁶)–Fe³⁺(d⁵) pair, keeping to a description which utilizes localized 3d-orbitals. The effect of dynamical correlation as well as the extension of the CAS(11,10) active space with the valence 5p electrons of tellurium on the Heisenberg exchange (the value of J) will be analyzed. Dynamical correlation has been accounted for using NEVPT2.^25–27 In Section 6 parameters of the spin-Hamiltonian inferred from SOC calculations of the Fe²⁺(d⁶)–Fe³⁺(d⁵) cluster are used to compute anisotropic g-factors and characterize the magnetic response of the thermally populated magnetic sublevels originating from the S = 1/2, 3/2, 5/2, 7/2 and 9/2 spin states. Section 7 is dedicated to the computation of vibronic coupling parameters defining the adiabatic potential of the S = 3/2 ground state under the scalar relativistic and spin–orbit coupling regime. Derivation and application of the formalism used to compute vibronic levels and wave functions and magnetic expectation values of valence delocalized Fe^2.5+–Fe^2.5+ complex (Class III) is the subject of Section 8. In Section 9 a critical evaluation of the double-exchange Hamiltonian in light of the results herein is briefly discussed. Finally, in Section 10 we summarize the conclusions of this work with outlooks for further applications using the same formalism.

3. The computational protocol

All computations were carried out using the recent release of the ORCA code.^18,19 A model for complex 3′ was derived from the X-ray structure of 3 by replacing the bulky substituents (1,3-diisopropylcyclohexa-1,3,5-trienyl groups) at each N-donor with methyl groups. The X-ray geometry of the N₂FeTe₂FeN₂ core was then kept frozen, while Cartesian coordinates of the methyl groups were optimized using density functional theory (DFT). Coordinate axes orientations facilitating ligand field interpretations of the results are shown in Fig. 5.

Fig. 5 Model complex 3′ and Cartesian axes orientation.

(1) CASSCF/NEVPT2 wave functions/energies of the states split out from the d⁶ (Fe²⁺) and d⁵ (Fe³⁺) configurations of iron in the diamagnetically substituted model complexes Fe²⁺(d⁶)Ga[ZnFe³⁺(d⁵)] (keeping to the same geometry as 3′), the five, 5 S = 2(45 S = 1) [1 S = 5/2(24 S = 3/2)], provide the many-particle basis for SOC calculations. These yield fine structure tensors – the ZFS tensor D and g-matrices of the spin-Hamiltonians of the S = 2[S = 5/2] ground states of the Fe²⁺[Fe³⁺] centers. Iron-ligand bonding was characterized using AILFT,^21–23 which allows for the determination of the respective 5 × 5 ligand-field matrices. The latter were further analyzed using the angular overlap model (AOM).^28,29 A best fit to the five diagonal and ten off-diagonal matrix elements affords bonding parameters which quantify the σ and π donor functions of the bridging tellurium and terminal nitrogen ligands.

(2) Energies of the lowest spin-levels of the Fe²⁺–Fe³⁺ exchange coupled pair (the ladder of spin-states) were computed using the non-distorted centrosymmetric structure of 3′. The five S₁ = 2 states on Fe²⁺(left) couple with the S₂ = 5/2 ground state of Fe³⁺ (right) yielding a manifold of five roots of each cluster total spin S (“%casscf, mult 2,4,6,8,10, nroots 5,5,5,5,5”). Already a small artificial shift by 0.01 Å of the two bridging tellurium atoms from Fe²⁺(left) to Fe³⁺(right) was enough to trigger Fe²⁺(d⁶)–Fe³⁺(d⁵) localization. The resulting “geometry, basis, wavefunctions” file, “gbw” provides a proper orbital guess in all calculations, both of the centrosymmetric and the geometrically distorted N₂FeTe₂FeN₂ cores. An important condition for the convergence of these calculations is keeping to the use of local orbitals at each and every step of the CASSCF procedure of the optimization of orbitals and CI coefficients (“%casscf, actorbs unchanged, actconstraints 1”). Different amounts of spin dependent dynamical correlation require computing each spin state of lowest energy separately (state and spin specific CASSCF with an active space consisting of 11 electrons distributed on 10 local 3d-MOs, CAS(11,10) followed by NEVPT2). Because of the two possible localizations of the extra electron, each spin state is, to a first approximation, doubly degenerate. Owing to the non-orthogonality and interactions of the localized wave functions describing the two open-shell iron centers, small splitting for each pair of spins results. This requires including both states in a given run (“%casscf, mult 2S + 1, nroots = 2”). Averaging the energies of the two roots for each spin affords the energies of the spin-ladder and extraction of the Heisenberg exchange and biquadratic exchange parameters, J and J′, respectively, from a best fit to the energies. To study the effect of the tellurium ligands on J and J′, an extended active space was introduced, including six more doubly-occupied valence orbitals (the 5p_x,5p_y,5p_z AOs on each tellurium) and 12 more electrons, CAS(23,16).

(3) Accounting for SOC in the CAS(11,10) scalar relativistic state manifold (configuration state functions (CSFs)), inclusion of 27720, 19 800, 4950, 440 and 10 CSFs for S = 1/2, 3/2, 5/2, 7/2 and 9/2, respectively, is mandatory in order obtain ZFS tensors D and g-matrices for a direct comparison with experimental – EPR and magnetic data. In the ORCA code this is achieved using QDPT^30,31 by diagonalization of the SOC Hamiltonian within the many particle basis of the CASSCF wave functions, large enough in order to capture properly the effect of SOC on the lowest, thermally-accessible magnetic sublevels. It follows that SOC calculations can never be state-specific and of the same quality as the scalar relativistic S states of lowest energy (see above). A compromise between accuracy and feasibility is inevitable. A trial and error analysis has shown that taking the 40 lowest scalar relativistic roots is enough to properly approximate the ZFS tensor D and g-matrices of the spin-Hamiltonians of the lowest levels for each spin. To this end, wave functions from a converged CASSCF calculation resulting from the “%casscf, mult 2,4,6,8,10, nroots 5,5,5,5,5” option were employed as the input for a spin-specific, state average SOC calculation using simple Complete Active Space Configuration Interaction, CASCI (i.e. not iterating orbitals, “%casscf maxiter 1, mult 2S + 1, nroots 40” but accounting for dynamical correlation “nevpt2 true”).

(4) The Adiabatic Potential Energy Surface (APES) resulting from the scalar relativistic calculation of the S = 3/2 ground state of 3′ was derived using a scan of the total CASSCF/NEVPT2 energy along the q mode in the spin-free space. The SOC calculation of the APES of the ground state Kramers quartet was accomplished by a Herzberg-Teller expansion of the full ZFS tensor D (not the traceless one) and the g-matrix up to second order in q.

(5) The formalism employed to calculate vibronic wave functions, energy levels, the ZFS and g-matrices using of 3′ is described in Section 8.

Throughout this work Douglas–Kroll–Hess scalar relativistic basis sets, DKH-Def2-TZVP^32,33 for all elements were employed, except for tellurium where “old-DKH-TZVP” was used. More details and sample input files used in the calculations may be found in the electronic (ESI†).

4. Electronic structure and magnetic anisotropy of the Fe²⁺Te₂N₂ and Fe³⁺Te₂N₂ complex units: ab initio ligand field analysis of the Fe–Te and Fe–N bonding

In tetrahedral complexes and coordinate axes orientations shown in Fig. 5, the 3d-orbitals split into e(d_xy, d_z2) and t₂(d_x2−y2, d_xz and d_yz) MO subsets, which are π- and σ + π antibonding, respectively. Ab initio ligand field diagrams for the Fe²⁺ and Fe³⁺ sites are depicted in Fig. 6. Orbital contour plots of the Fe²⁺ complex (Fig. 6, left) show significant participation of ligand functions in the upper (t₂) subset and almost no such ligand tails in the lower e subset with a tiny low-symmetry splitting. Orbital energy splitting is much larger in the Fe³⁺ complex (Fig. 6, right) where both the t₂ and the e -type orbitals are affected by the symmetry lowering. This symmetry lowering stems for the presence of two different ligands in the coordination sphere of the two centers – two nitrogen and two tellurium ligands with quite different donor properties. While tellurium can be considered as σ and a weaker π-donor, for ketimine nitrogen being an sp²-type ligand both σ and out-of-plane π interactions (in-plane nitrogen orbitals are involved in strong interactions with their neighboring carbon atoms) with iron are expected. This is supported by values of the antibonding energies extracted from a best fit of matrix elements of the 5 × 5 ligand field matrix parametrized using the AOM. When applied to the Fe²⁺N₂Te₂ and Fe³⁺N₂Te₂ complex units the general matrix element of the ligand field matrix is given by eqn (4.1)While the factors “F” in this equation account for the dependence of the ligand field on the angular geometry described by the Euler angles θ, φ and ψ (in the case of N^imin), the pre-factors e^Te_σ, e^Te_π and e^N_σ, e^N_πs quantify the energy of antibonding imposed by the tellurium and nitrogen donors, respectively. Values of these parameters listed in Table 1 reflect the weak ligand field in the case of Fe²⁺ with almost zero contributions from e^Te_π and the much stronger ligand field for Fe³⁺ where for both tellurium and nitrogen, appreciable contributions from both σ and π (out-of-plane in the case of N^imin)-type antibonding to the ligand field splitting are computed. In addition to the effect of the ligands on the low-symmetry splitting of the parent tetrahedral orbitals, the large deviations of the Te–Fe–Te and N–Fe–N angles (104° and 91°, respectively) from the angle of an ideal tetrahedral geometry (109.47°) on both Fe²⁺ and Fe³⁺ sites contribute to the lower symmetry. The rather low effective symmetry on each site is also reflected in the large anisotropy of the S = 2 ground state ZFS tensor for the Fe²⁺ site with D = −4.205 cm⁻¹ and quite large orthorhombic splitting ratio E/D = 0.110; an even larger orthorhombicity (in agreement with the stronger ligand field), E/D = 0.295 for the S = 5/2 ground state at the Fe³⁺ site is computed. In Table 2 we list S = 2 spin-allowed d-d transitions for Fe²⁺ and the lowest spin-forbidden S = 5/2 to S = 3/2 transitions for Fe³⁺, SOC splitting of the S = 2 and S = 5/2 ground states and anisotropic spin-Hamiltonian g-factors for Fe²⁺. As will be shown in the next Section, these findings help rationalize the magnetic anisotropy of the exchange coupled iron pair.

Fig. 6 AILFT 3d-MO energies and orbital contour plots (contour values: 0.030 e Bohr⁻³) from CASSCF/NEVPT2 calculations of the diamagnetically substituted Fe²⁺_μTe₂Ga³⁺ (left) and Zn²⁺_μTe₂Fe³⁺ (right) model complexes with the geometry of 3′. For the Cartesian axes orientation see Fig. 5; orbital notations are based on the leading contribution of the corresponding AILFT MO as follows (in %): 3d_z2(87), 3d_xy(99) 3d_x2−y2(87), 3d_yz(95), 3d_xz(100%) (Fe^II, left) and 3d_z2(69), 3d_xy(100) 3d_x2−y2(69), 3d_yz(98) 3d_xz(99%) (Fe^III, right).

Table 1 AILFT ligand field 3d-MO energies (NEVPT2) and local AOM parameters (in cm⁻¹) for the Fe²⁺μTe₂N^imin₂ and Fe³⁺μTe₂N^imin₂ model complexes derived from corresponding diamagnetically substituted [FeGaTe₂]⁺ and [FeZnTe₂]⁺ model clusters^a

3d-Type MOs(Td)	Fe^2+μTe2N^imin2	Fe^3+μTe2N^imin2
a Symmetry notations for the parent tetrahedral (Td) point group are listed in parenthesis.
3dz2(e)	0	0
3dxy(e)	395	1595
3dx2−y2(t2)	2184	5735
3dyz(t2)	3550	8931
3dxz(t2)	6296	9415
e ^Teσ	2188	7552
e ^Teπ	33	2619
e ^Niminσ	4243	6277
e ^Niminπs	434	850
Rms deviation NEVPT2/comp.	256	426

---

Table 2 Scalar relativistic d–d transitions, ground state SOC zero-field split sublevels, spin-Hamiltonian zero-field-splitting parameters D, E/D and g-factors of Fe²⁺μTe₂N^imin₂ and Fe³⁺μTe₂N^imin₂ from CASSCF/NEVPT2 calculations of Fe²⁺μTe₂Ga³⁺ and Zn²⁺μTe₂Fe³⁺ diamagnetic substitutions of 3′^a

Fe^2+Te2N2	Fe^3+Te2N2
a KD = Kramer's doublet.
d–d Transitions/NEVPT2/cm^−1	d–d Transitions/NEVPT2/cm^−1
S = 2 ^5E 0	S = 5/2 ^6A1 0
S = 2 580	S = 3/2 ^4T1 19 218
S = 2 ^5T2 2490	19 588
S = 2 4503	23 163
S = 2 7521	S = 3/2 ^4T2 24 861
	26 177
	26 279
SOC-QDPT sublevels S = 2 ground state/cm^−1	SOC-QDPT sublevels S = 5/2 ground state/cm^−1
0	KD1 0
0.15	KD2 2.78
11.17	KD3 5.84
14.02
17.30
D = −4.205 cm^−1	D = 0.851 cm^−1
E/D = 0.110	E/D = 0.295
g eff = 2.059, 2.071, 2.162	g eff = 2.002, 2.002, 2.002
	g KD1 = 0.735, 1.083, 9.585
	g KD2 = 4.070, 4.258, 4.498
	g KD3 = 0.491, 0.665, 9.770

---

5. Exchange coupling in the Fe²⁺(d⁶)–Fe³⁺(d⁵) pair: the spin ladder

Ligand field analysis in Section 4 shows five S = 2 states of Fe²⁺ in a range of 7000 cm⁻¹ and a much larger energy gap between the S = 5/2 ground state and the lowest spin forbidden d-d excitations at the Fe³⁺ site. The small low-symmetry splitting of the ⁵E tetrahedral ground state of Fe²⁺ leads one to expect large anisotropy from competing exchange and spin–orbit coupling. The latter is neglected in scalar relativistic calculations, which, in this approximation, allows one to consider S as a good quantum number and focus on the ground state only. Because of the two possible localizations of the surplus electron, each spin state is, to a first approximation, doubly degenerate but actually slightly split due to the electron communication between the two centers. In Table 3 we list energies of lowest pairs for each spin (S = 1/2, 3/2, 5/2, 7/2 and 9/2) obtained from state and spin specific CAS(11,10)/NEVPT2 calculations. The results allow one to derive a mean energy for each pair (Table 3, second column) and a splitting due to the small (because of the localization procedure imposed), yet non-negligible resonance delocalization (Table 3, fourth column). Focusing on the first set of data, least square fitting of the parameters of the spin-Hamiltonian including Heisenberg (J) and biquadratic (J′) exchange (eqn (5.1) and (5.2)) allows for the determination of values for J = −69.04 cm⁻¹ and J′ = 0.65 cm⁻¹. Taking the second set, a spin-dependent delocalization parameter B (eqn (1)) can be extracted from the data using expressions listed in Table 3, column 5. A value of B = 13 cm⁻¹ reproduces the five data points reasonably well. However, given the localization procedure applied here, this value cannot be considered as relevant for the case where the extra electron is delocalized over the two centers; this is the situation expected for complex 3.

Ĥ_exch = −2JŜ₁Ŝ₂ − 2J^'(Ŝ₁Ŝ₂)² (5.1)

E_exch = −JS(S + 1) − J^'[S(S + 1)]² (5.2)

Table 3 Lowest non-relativistic S = 1/2, 3/2, 5/2, 7/2, 9/2 spin- and state- specific energies from CASSCF/NEVPT2 calculations^a

S=	Energy (cm^−1)	E −(E+) pair energy	Resonance splitting	2(S + 1/2)B expression	Spin-dependent delocalization parameter B
a CAS(11,10), “mult 2S + 1”;“nroots 2”. Parameters J = −69.04 cm^−1; J′ = 0.6549 cm^−1 (root mean square deviation 2.77 cm^−1) were deduced from a least square fitting of the Heisenberg Hamiltonian extended with a biquadratic term: Ĥexch= −2JŜ1Ŝ2 − 2J^'(Ŝ1Ŝ2)^2 and the following energy dependence on the total spin: Eexch = −JS(S + 1) − J^'[S(S + 1)]^2.
1/2	0	0	30	2B	15
		30
3/2	201	194	44	4B	11
		238
5/2	505	485	70	6B	12
		555
7/2	870	831	109	8B	14
		940
9/2	1257	1199	147	10B	15
		1346

---

In Table 4 we compare values of J and J′ determined using various combinations of options, including: active spaces CAS(11,10)/CAS(23,16), state average/state and spin specific, neglecting (CASSCF) or accounting (CASSCF + NEVPT2) for dynamical correlation as well as the case of a simple CASCI + NEVPT2(CAS(23,16)) calculation. Based on the reported value of J = −200 cm⁻¹¹⁰ (albeit with some precaution because of large and unspecified experimental error bars) we may conclude that, irrespective of the option employed, the electronic structure method leads to an underestimate of the computed value of J, except in the CAS(23,16) CASCI(NEVPT2 + CASCI) method, which matches closely to the reported value. The largest difference, (one order of magnitude smaller than the experimentally determined value of J) is obtained when using state averaged CASSCF and excluding dynamical correlation (Table 4, second column). Dynamical correlation (NEVPT2) improves only slightly the value of J. It becomes twice larger for the state- and spin-specific option within the CAS(11,10) model space (Table 4, columns 4 and 5) and is exclusively large when going to the extended CAS(23,16) active space (Table 4, columns 7 and 9). One can conclude that dynamical correlation is of extreme importance for magnetic exchange; in the case of 3′ it contributes to larger negative values of the exchange integral J. One can try understanding this result by the fact that the effect of dynamical correlation, which leads to negative corrections to the CASSCF reference energy, decreases with increasing S as reflected already by the CAS(11,10) results (Fig. 7 (left)). Its contribution to the energy separation between the ground S = 1/2 and the S = 3/2, 5/2, 7/2 and 9/2 spin states is illustrated in Fig. 7 (right). Using these data, a best fit of the parameters yields values of −43.4 cm⁻¹ and 0.47 cm⁻¹ for J and J′, respectively. Thus, as much as 62% of the total value of J (J = −69.1 cm⁻¹, Table 4, column 5) comes from dynamical correlation.

Table 4 Energies of the spin ladder S = 1/2, 3/2, 5/2, 7/2 and 9/2 and parameters of the Heisenberg (J) and biquadratic exchange (J′) (in cm⁻¹) from a best fit to state and spin average and state- and spin-specific CAS(11,10) and CAS(23,16) with CASSCF and CASSCF/NEVPT2^a

S	CAS(11,10)				CAS(23,16)
	State average over the five scalar relativistic states for each spin		State and spin specific		State and spin-specific
	CASSCF	CASSCF/NEVPT2	CASSCF	CASSCF/NEVPT2	CASSCF	CASSCF/NEVPT2	CASCI	CASCI/NEVPT2
a Spin-energies reproduced using the best fit J and J′ values are listed in parenthesis.
1/2	0	0	0	0	0	0	0	0
3/2	65(65)	106(107)	75(74)	201(198)	107(107)	348(372)	−14(47)	936(560)
5/2	172(172)	280(280)	191(191)	505(503)	274(274)	915(945)	137(139)	1336(1320)
7/2	320(320)	514(514)	338(338)	870(874)	483(484)	1689(1649)	328(297)	1816(2016)
9/2	505(505)	797(797)	500(500)	1257(1256)	713(713)	2372(2384)	541(552)	2352(2282)
J	−21.81	−36.14	−25.65	−69.05	−36.84	−129.37	−14.05	−206.47
J′	0.030	0.114	0.189	0.655	0.279	1.178	−0.350	4.37
Stand dev.	0.089	0.470	0.388	2.777	0.296	33.24	39.52	258.88

---

Fig. 7 Dynamical energy corrections ΔE_NEVPT2 to the total CASSCF reference scalar relativistic ground state energy from state and spin-specific CASSCF/NEVPT2 calculations in dependence total spin S of the exchange coupled Fe²⁺_μTe₂Fe³⁺ pair (left); dynamical correlation energy contributions to the energy separations between the S = 1/2 ground and the S = 3/2, 5/2,7/2 and 9/2 excited states (right).

We finally mention results from calculations which employ, as is usually done, canonical orbitals as a guess in the CASSCF/NEVPT2 procedure for 3′. The latter imply a full delocalization of the extra electron over the two Fe sites. These studies have shown stabilization of a S = 9/2 ground state and energy spacing to the S = 7/2, 5/2, 3/2 and 1/2 states of 30248, 27657, 25356 and 38537 cm⁻¹ – thus violating the Landé interval rules (eqn (5.2)) and rendering parameters J and J′ ill defined. Dynamic correlation corrections to the CASSCF reference energy computed using such orbitals exceed by as much as 13 167 to 24 140 cm⁻¹ those obtained when imposing a Fe²⁺μTe₂Fe³⁺ electron localization.

6. Magnetic sublevels from spin–orbit coupling calculations: comparison with experiment

Isotropic ground state spin-levels, S = 1/2, 3/2, 5/2,7/2, 9/2 get largely modified under the combined effects of low symmetry (see Section 4) and spin–orbit coupling. States with S > 1/2 undergo splitting in the absence of a magnetic field into two (S = 3/2), three (S = 5/2), four (S = 7/2) and five (S = 9/2) Kramers doublets with quite different anisotropies. This is quantified by their g-factors (g-matrices) and the orientations of their canonical magnetic axes (the axes system in which the ZFS tensor D and/or g·g′ become diagonal). Zero-field splittings (D and E/D parameters) and g-factors for the sublevels of the S states of lowest energy are listed in Table 5. In these calculations, SOC within the manifold of scalar relativistic states was considered. Mixing between S and S ± 1 states is presently neglected (vide infra). In all four S = 3/2–9/2 states, D is computed to be negative and therefore dominated by the relatively large local ZFS at the Fe²⁺ site (D = −4.2 cm⁻¹), a result stemming from the near orbital degeneracy of the ⁵E ground state of the Fe²⁺ site. It is remarkable that the g-factors of the lowest Kramers doublet resulting from the S = 3/2 calculation (Table 5, third row) exactly reproduce the reported experimental ones (Table 5, bottom row); also the computed value of E/D = 0.313 nearly perfectly agrees with the reported value (E/D close to 1/3). Do the same results also reproduce the χT vs. T data? Fig. 8 (left) shows that this is definitely not the case; except for the two lowest temperature data points, computed values for χT are much lower that experimental ones. To this end, we extended the spin–orbit coupling calculation for the S = 3/2 manifold with the S = 5/2 states (“CAS(11,10), mult 6,4, nroots 40,40”). As a side remark, that in contrast to single center transition metal complexes, the S to S ± 1 mixing rule for SOC appears not to hold for exchange coupled spins. Energies of the lowest scalar relativistic S = 3/2 and S = 5/2 states and SOC split sublevels in the range from 0 to 200 cm⁻¹ are listed in Table 6. While this extension does not affect g-factors for the lowest S = 3/2 state, thermally accessible magnetic sublevels originating from the lowest S = 5/2 state in the range 159–194 cm⁻¹ appear. They readily explain the increase in the χT curve at temperatures higher that 100 K. Simulations of the χT vs. T behavior lend support of this proposition (Fig. 8, right). Yet, as indicated in Fig. 8 (right), there is a deviation between experimental and simulated data points in the temperature range 25 < T < 100 K. Below we will try analyze this in terms of a vibronic coupling model.

Table 5 g-Factors and energies of Kramers doublets originating from the S = 1/2, 3/2, 5/2, 7/2 and 9/2 spin-free states of lowest energy from spin-specific SOC-QDPT NEVPT2 calculations on the mixed valence [Fe²⁺μTe₂Fe³⁺]⁻ dimer^a

Lowest spin state S	Energy/cm^−1	g 1	g 2	g 3	D/cm^−1	E/D	E tot eH (cm^−1) −16731.
a CASCI/NEVPT2 calculations based on state average CASSCF (CAS(11,10), “mult 10,8,6,4,2”; “nroots = 5,5,5,5,5”) results and the X-ray geometry of a [Fe^2+μTe2Fe^3+]^− truncated model complex; included are zero-field split energies and values of D and E/D ratios, orientations of the canonical –magnetic axes with respect to the adopted molecular Cartesian frame (Fig. 5).
1/2	0	0.652	0.713	0.876	—	—	−.742988(56380)
3/2	0	1.412(x,y)	1.915(y,x)	5.748(z)	−1.836	0.313	−.742560(56437)
	4.2	1.568(z)	2.131(y,x)	5.411(x,y)
5/2	0	0.019(x,y)	0.020(y,x)	10.794(z)	−2.453	0.044	−.741847(56630)
	9.4	0.948(y,x)	0.963(x,y)	6.467(z)
	14.9	2.131(z)	4.974(y,x)	6.934(x,y)
7/2	0	0.022(y,z)	0.026(x)	14.777(z,x)	−1.827	0.0185	−.740886(56841)
	7.9	1.114(y,z)	1.411(x)	10.208(z,y)
	12.6	4.364(z)	5.760(x)	6.513(y)
	17.6	0.583(z,y)	0.985(x)	13.349(y,z)
9/2	0	0(x,y)	0(y,x)	18.718(z)	−1.348	0.176	−.999893(0)
	14.33	0.000(x,y)	0.000(y,x)	14.560(z)
	25.15	0.057(x,y)	0.058(y,x)	10.394(z)
	32.37	2.151(x,y)	2.170(y,x)	4.141(z)
	36.10	1.989(z)	7.824(y,x)	12.133(x,y)
X-EPR simul.		1.515	1.950	5.740	Max. 15	1/3

---

Fig. 8 Experimental (red points) and simulated χT vs. T dependencies from spin–orbit coupling calculations using the manifold of 40 lowest S = 3/2 scalar relativistic roots (left) and including in addition the manifold of 40 lowest S = 5/2 roots (right). The latter results explain the increase of χT at T > 100 K with Boltzmann populations of excited S = 5/2 levels (see Table 6).

Table 6 Scalar relativistic and SOC energies of the lowest energy levels of 3′ from CAS(11,10) calculations including lowest 40 S = 3/2 and 40 S = 5/2 roots

	Scalar-relativistic	SOC-QDPT	g 1	g 2	g 3
3/2	0	0	1.685(x)	2.004(z)	5.496(y)
		5.59	1.506(y)	2.218(z)	5.500(x)
3/2	48	28.74	1.623(z)	3.476(x)	4.378(y)
		48.87	0.479(x,y)	0.535(y,x)	5.633(z)
5/2	160	158.98	0.043(x,y)	0.048(x,y)	10.691(z)
		181.85	2.442(x)	2.864(y,z)	5.924(z,y)
		194.03	1.652(z)	2.846(y)	8.100(x)

---

7. The adiabatic potential of 3′ in the lowest S = 3/2 state: effect of spin–orbit coupling on the topology of the ground state adiabatic potential surface

Thus far, the geometry of the N₂FeTe₂FeN₂ core in the model complex 3′ was identified with the static X-ray structure of 3 showing two equivalent iron sites; bond angles and bond distances show a centrosymmetric iron dimer (Fig. 9, left). In a mixed valence class III complex, such a geometry may, or may not be at ground state energy minimum – a double well adiabatic ground state potential energy surface may arise; DFT geometry optimizations can potentially show how a geometry of higher symmetry relaxes to a minimum of energy at lower symmetry, while frequency calculations at the centrosymmetric structure can show the shape of the active normal modes (imaginary frequencies) along which the complex will eventually distort. Such calculations are at best scalar relativistic and the role of spin–orbit coupling, expected to counteract vibronic forces (tending to support high-symmetry configurations) cannot be routinely accounted for. Moreover, numerical instabilities may show artificial minima, that depend on the DFT functional, the basis set and the auxiliary basis functions used in the calculations. In order to supply a more realistic model for the adiabatic potential energy surface (APES) which takes SOC into account, a scan of the ground state energy along the symmetry breaking mode q was carried out (Fig. 9 right). Starting with the geometry at q = 0 (corresponding to the one shown in Fig. 9, left) and using a grid of points at positive and negative values of q, the dependence of the total energy on q under the effect of SOC for the S = 3/2 ground Kramers doublet could be traced.³⁴ Data points (Fig. 9, right) could be fitted consistently using quadratic polynomial (eqn (7.1)) thus accounting for superfluous anharmonicity.

Fig. 9 The centrosymmetric geometry of 3′ (left) used as reference (q = 0) to scan the potential surface along the symmetry breaking local mode q (right).

The coefficients of this polynomial can be directly related with the relativistic SOC energy reference E_o = c, the SOC linear vibronic coupling (V = b) and the SOC harmonic force constant K (K = 2c).

aq² + bq + c = (1/2) Kq² + Vq + E_o (7.1)

The resulting energy expression (eqn (7.2)) allows one to relate

E = (1/2) Kq² + Vq + E_o (7.2)

the geometrical distortion q_o, the energy stabilization ΔE_stab and the vibronic energy parameter λ with V and K from one side, and to compare these energies with the results obtained from scalar relativistic calculations from the other side (Table 7). The comparison shows, that accounting for SOC leads to twice smaller distortions (q = 0.055 vs. 0.114 Å) and three times smaller stabilization energies (235 vs. 718 cm⁻¹, see also Fig. S1 in the ESI† facilitating the comparison). It follows from Table 7 listing the set of parameters extracted from the analysis of the two cases – with and without SOC, that one can trace back these changes to lowering V (increase of K) when SOC is considered. Values of the harmonic force constant K allow one to estimate the harmonic frequency corresponding to the effective mode q − ħω = 121 cm⁻¹ (scalar relativistic) vs. 142 cm⁻¹ (SOC). This will be used in section 8. Computed scalar relativistic values of ΔE_stab (718 cm⁻¹) are close to ones extracted from the analysis of experimental χT vs. T data using a Bleaney–Bowers model (ΔE_stab = 839 cm⁻¹)¹⁰ but are much larger than the one given by the SOC calculation (235 cm⁻¹).

Table 7 Vibronic coupling parameters for the lowest S = 3/2 Kramers doublet from QDPT spin–orbit coupling NEVPT2 calculations and from scalar relativistic NEVPT2 calculations of the S = 3/2 ground state of 3′

	Scalar relativistic CASCI/NEVPT2	QDPT-SOC CASCI/NEVPT2	Reported In ref. 1
q o (Å)	0.114	0.055	0.20
V (cm^−1 Å^−1)	12 605	8369	8375
K (cm^−1 Å^−2)	110 500	151 057	41 873
ħω (cm^−1)	121	142	140
ΔEstab	718	235	838
	1438	464	1675

---

8. Spin-phonon coupling in Fe^2.5+μTe₂Fe^2.5+

The APES of a Class III mixed valence binuclear complex is spin-dependent (Section 1, eqn (1.1) and (1.2)). Focusing on the S = 3/2 ground state of the Fe^2.5+μTe₂Fe^2.5+ pair, and taking localized valence bond structures Fe²⁺(_μTe₂)Fe³⁺ (left) and Fe³⁺(_μTe₂)Fe²⁺(right), the APES takes the analytical form:with V and K – the linear vibronic coupling and the harmonic force constants, respectively (see Section 7 for the derivation of their values), and B – the electron transfer parameter (see eqn (1.1) and (1.2) and following text). Before computing and analyzing the effect of vibronic coupling (K and V) and electron transfer (B) on the vibronic levels it is quite instructive to take the resonance form Fe²⁺(μTe₂)Fe³⁺ (left) first and neglect electron transfer. This is the generic case of a hypothetical completely localized mixed valence (Class I) complex. It may apply for 1 and 2 but not for 3. Localization may be induced by the presence of two different coordination environments around the two centers (see examples for such Fe²⁺–Fe³⁺ dimers see ref. 17) and/or by static geometrical distortions caused by crystal packing effects, of from both. The Hamiltonian (eqn (8.2)) is set up in the basis of simple BO products |SM_S〉|n〉, with |SM_S〉 the functions M_S = 3/2,1/2, −1/2, −3/2 to the S = 3/2 spin and |n〉 -the wave functions of the quantum Harmonic oscillator up to a number N_vib(n = 0, 1, …, N_vib).

Ĥ_total(Class I) = Ĥ_vibr + Ĥ_ZFS + Ĥ_Zeeman (8.2)

In eqn (8.2)Ĥ_vibr (eqn (8.3)), Ĥ_ZFS (eqn (8.4)) and Ĥ_Zeeman (eqn (8.5)) are the vibrational, zero-field splitting and the Zeeman terms, respectively:

Ĥ_vibr = Eⁿ_vibr = ħω (n + 1/2) (8.3)

A Herzberg-Teller expansion of Ĥ_ZFS (eqn (8.6)) and Ĥ_Zeeman (eqn (8.7)) up to quadratic terms with respect to the symmetry breaking local distortion mode q has been used to compute the ZFS and Zeeman matrices, and the matrix Ĥ_total(Class I) using the master formula for its matrix elements (eqn (8.8)).

Ĥ_ZFS = Ĥ^o_ZFS + (∂Ĥ^o_ZFS/∂q)_oq + (1/2)(∂²Ĥ^o_ZFS/∂q²)_oq² (8.6)

Ĥ_Zeeman = Ĥ^o_Zeeman + (∂Ĥ^o_Zeeman/∂q)_oq + (1/2)(∂²Ĥ^o_Zeeman/∂q²)_oq² (8.7)

Eqn (8.6) and (8.7) account for the dependence of the six components of the D tensor and nine independent matrix elements of the g-matrix, respectively on q as illustrated in Fig. 10 (top for D, middle and bottom for g). A Herzberg-Teller expansion of each matrix element up to quadratic terms in q shows a strong non-linearity both in the diagonal and in the off-diagonal matrix elements. In what follows only linear terms of this expansion will be considered; quadratic terms have been used in order to more precisely approximate linear terms. The solution of the vibronic problem (eqn (8.8)) can be greatly simplified by noticing, that terms in eqn (8.6) and (8.7) linear in q are of the order of 0.1 cm⁻¹, and therefore much smaller then ħω (141 cm⁻¹) and the vibronic stabilization ΔE_stab (235 cm⁻¹, see Table 7). Therefore, one can eliminate linear vibronic coupling term in the APES (the Vq in eqn (8.1)) by shifting the origin of q from 0 to its value at the minimum of energy q_min = V/K = 0.05 Å) after properly transforming the q dependent terms in eqn (8.6) and (8.7). Energies of the lowest four vibronic levels and g-factors are listed in top half of Table 8. While no changes are computed concerning energies and g-factors of the two Kramers doublets resulting from the scalar relativistic S = 3/2 ground state, the next two excited states, originally a purely vibrational excitation at 141.2 cm⁻¹, now become split into two sublevels (141.3 and 145.4 cm⁻¹) the difference exactly matching the ground state zero-field splitting (4.1 cm⁻¹). Vibronic mixing between the ground and excited Kramers doublets KD1-KD2, 0.069 and 0.204 cm⁻¹, are of importance for the magnetic relaxation. This is outside the scope of the work herein.

Fig. 10 Dependence of the spin-Hamiltonian parameters of the S = 3/2 ground state of 3′ – the D tensor, top (in cm⁻¹) and g-matrix: diagonal-(middle) and off-diagonal (bottom) matrix elements on the symmetry breaking local mode q (in Å); linear and quadratic derivatives are given in cm⁻¹/Å and cm⁻¹/Ų respectively (for D) and 1/Å and 1/Ų (for g) When used in the context of eqn (8.6)–(8.8) these have to be converted into dimensionless q dividing by the conversion factor ; M_eff = 2M_Te = 255.2, ħω = 141.2 cm⁻¹ – the vibrational frequency; q → fq; (dD/dq)_o → (dD/dq)_o/f; (d²D/dq²)_o → (d²D/dq²)_o/f²; (dg/dq)_o → (dg/dq)_o/f; (d²g/dq²)_o → (d²g/dq²)_o/f².

Table 8 Lowest thermally accessible SOC vibronic energy levels and Kramers doublet g-factors from CAS(11,10) CASSCF/NEVPT2 calculations^a

Class I	ħω = 141.2 cm^−1, Nvib = 32
E/cm^−1	g 1	g 2	g 3
a KD = Kramers doublet.
KD1 0	1.401	1.882	5.767
KD2 4.2	1.588	2.159	5.421
KD3 141.3	1.401	1.882	5.767
KD4 145.4	1.588	2.159	5.421

---

Class III	V = −8369 cm^−1; ħω = 141.2 cm^−1, Nvib = 32, B = 163 cm^−1
KD1 0	1.402	1.883	5.776
KD2 4.2	1.587	2.159	5.422
KD3 44.2	1.401	1.882	5.761
KD4 48.4	1.588	2.158	5.420

---

We now invoke electron transfer into our consideration. The Hamiltonian of eqn (8.1) can be transformed from a localized (|Ψ_left〉 = |Ψ_A(d⁶) Ψ_B(d⁵)〉, |Ψ_right〉 = |Ψ_A(d⁵) Ψ_B(d⁶)〉) into a canonical (|Ψ_±〉) form following eqn (8.9). In this presentation electron transfer (2B) and vibronic coupling (Vq) change their roles; the first appears now at the diagonal while the latter enters into the off-diagonal part of the APES (eqn (8.10)). A configuration energy diagram is presented schematically in Fig. 11.

Fig. 11 Configuration energy diagram (schematic) for 3′ in the limit of a complete delocalization of the extra electron over the two iron sites (Fe^2.5+_μTe₂Fe^2.5+, Class III).

The total Hamiltonian (8.11) has been extended by one additional term, Ĥ^±_vibronic and the basis set is now |SM_S〉|n〉|l〉 of total dimension 4N_vib2 extended with an additional quantum number l = +1, −1, labelling basis functions corresponding to the |Ψ₊〉 and |Ψ₋〉 branches, respectively.

Ĥ_total = Ĥ_vibr + Ĥ_ZFS + Ĥ_Zeeman + Ĥ^±_vibronic (8.11)

The matrix of Ĥ_total in this basis consists of two diagonal matrices of dimension 4N_vib identical to Ĥ_total(Class I) (eqn (8.2)), their diagonal elements are now shifted up (for branch|Ψ₊〉) and down (for branch |Ψ₋〉) in energy by 2B, respectively, and a sub-diagonal matrix appears whose matrix elements are computed according to the master equation eqn (8.12).

The structure of the vibronic Hamiltonian of a mixed valence Class III complex is presented in Fig. 12.

Fig. 12 Structure of the vibronic Hamiltonian of a mixed valence Class III binuclear transition metal dimer complex; I -the 4N_vib × 4N_vib unit matrix.

The parameter B is not uniquely defined and will be subjected to an analysis in Section 9. Yet, we decided here to render B freely variable and obtained a value of B = 163 cm⁻¹ from a simulation of the χT vs. T data (Fig. 13); a value about five times smaller than previously reported (B = 750 cm^−1 10) results. Using this value, we computed energies of the lowest four Kramers doublets along with their anisotropic g-values (see bottom half of Table 8). The comparison between the results from the Class III and Class I modelling of 3′ shows the combined effect of electron transfer and vibronic coupling on the vibronic levels: the energies of the third and the fourth excited Kramers states change from 141.3 and 145.4 cm⁻¹, energies stemming from a zero-field split vibrational excitation (Class I) (the value of the nominal first vibrational excitation) to as much as three times smaller energies, 44.2 and 48.4 cm⁻¹. Energies of the lowest three excited levels of 3 are not known. However, spin-relaxation times determined from the temperature dependence of the Lorentzian broadening of spin-relaxation times have been reported (Fig. 14). From these data the authors obtained an Orbach energy gap of U = 60 cm⁻¹ ≈ 4D. Inspection of the Arrhenius ln(1/τ) vs. 1/kT plots in Fig. 14 allows one to recognize two different slopes, a larger one for the data at higher temperatures (T = 20, 13.3 K), and a smaller slope reflected by the data at lower temperatures (T = 10, 8, 6.7 K). With some precaution (in view of the inherent approximations), we were able to simulate the data using two exponents to afford two Orbach energies U₁ = 5.6 cm⁻¹ and U₂ = 30.5 cm⁻¹. These are in a reasonable agreement with the computed values −4.2 cm⁻¹ and 45 cm⁻¹, respectively. A more rigorous direct access to the energies of the lowest excitations would be possible using dedicated experiments, such as inelastic neutron scattering. We hope that the theoretical results of this work will stimulate such further studies on 3.

Fig. 13 Comparison between experimental (“Exp”) and theoretical χT vs. T dependencies using the static geometry of 3′ and vibronic wave functions computed using Class I and Class III models.

Fig. 14 Spin-relaxation rates from Lorentzian broadening of the low-field sharp line transition in the X-Band EPR spectrum of 3 (red points);¹⁰ simulation of the data using a single exponential Orbach relaxation mechanism (black line), and simulation of the data using a bi-exponential Orbach mechanism, blue line (T = 20, 13.3 K data points), green line (T = 10, 8, 6.7 K data points).

9. Critical view of the double-exchange model and its application for the interpretation of magnetic and spectroscopic data of binuclear mixed valence complexes

In this study we employed the double exchange vibronic coupling (DEVC) model (eqn (1.2), (8.2) and (8.10) to the interpretation of the spectroscopic and magnetic properties of the Class III Fe^2.5+μTe₂Fe^2.5+ mixed valence dimer 3. This model has its scope of validity and limitations, which we now discuss in some detail. The goal of the DEVC model is to describe only part of the energy spectrum of the total Hamiltonian, that part connected with the spin-exchange within the Fe²⁺(S = 2)–Fe³⁺(S = 5/2) pair and the transfer of one electron from the Fe²⁺ to the Fe³⁺ center. Vibronic coupling to one symmetry breaking effective normal mode q was considered. The effects of these complex processes were captured by four parameters only: spin exchange (J), electron transfer (B) and vibronic coupling (V,K). The DEVC model presupposes one state of each spin-multiplicity from S = 1/2 to S = 9/2. This tacitly assumes, that other states of the same spin are far apart in energy and cannot mix (i.e. the “giant spin approximation”) and it neglects spin–orbit coupling. However, SOC and low-symmetry yield leading contributions to the magnetic anisotropy.³⁵ Therefore the DEVC model intrinsically fails to relate observables – EPR spectra and the magnetic susceptibility – with multiplet fine structures. First principles based multireference (MR) ab initio calculations do not suffer from this drawback but are rather complex and difficult to analyze. Some compromise between DEVC and MR approaches is possible: (1) CASSCF/NEVPT2-SOC calculations using an orbital localization procedure allow for determination of the parameters of the spin-Hamiltonian – the zero field splitting tensor D and g-factors extracted from first principles multireference calculations. (2) Linear vibronic coupling and linear spin–phonon coupling parameters – the derivatives of ZFS D-tensor and the g-matrix with respect to q – were then obtained numerically. (3) The most critical point in the chain of approximations is the determination of B which in the DEVC model defines the electron transfer operator T̂₁₂. The action of T̂₁₂ is defined by eqn (9.1); it considers the transfer of one electron between the two resonance forms, in our case:

T̂₁₂|Ψ_left〉 = |Ψ_A(d⁶)Ψ_B(d⁵) 〉 = |Ψ_right〉 = |Ψ_A(d⁵)Ψ_B(d⁶)〉 (9.1)

The first principles equivalent of T̂₁₂ acts on the complex SOC components of two S = 3/2 ground state Kramers doublets (eqn (9.2)), consisting of one- and two-electron electron transfer operators:

In eqn (9.2) the operator T̂₁₂ is written down in second quantization with a⁺_jβ∈B, a⁺_jβ∈A and a_iβ∈A, a_iβ∈B pairs of creation and annihilation operators describing the transfer of the surplus β spin from the Fe²⁺(d⁶) center (left) to the half-filled shell of Fe³⁺(d⁵) (right, or vice versa). Analysis of the local Fock operator in the basis of converged CASSCF wave function, evaluated using CAS(11,10) converged natural orbitals allows for extraction of the one-electron transfer integrals illustrated in Fig. 15 (see ESI,† for computational details). It follows from this numerical example, that out of the five possible one-electron transfer pathways, three, from the lowest occupied orbital of Fe²⁺ to the second, third and fourth orbital of Fe³⁺, are comparable in magnitude und should be simultaneously accounted for. In other words, three different parameters (B₁₂, B₁₃ and B₁₄), rather than just only one B parameter, should be considered. Since four excited states of the d⁶ configuration may couple the ground state one (Fig. 15) via SOC, a total of as much as 25 charge transfer integrals (Table 9) would have to be considered. The quantities in Fig. 15 have been obtained using CASSCF wave functions. One expects that dynamical correlation will modify the parameters B_ij to render the problem at the present stage not tractable. Efforts in this direction are the subject of ongoing projects in our group.

Fig. 15 One-electron transfer (hopping) integrals for spin-down electron in the S = 2 ground state of Fe²⁺ from 5 × 5 block diagonalization of the 10 × 10 Fock-matrix of the tellurium bridged iron-pair extracted from converged CASSCF calculations with a model space distributing 11 electrons over the ten active localized MOs [CAS(11,10)].

Table 9 One electron (hopping) integrals quantifying the transfer of the excess spin-down (β-spin) electron from the Fe²⁺(d⁶) site to the half-filled shell of the Fe³⁺(d⁵) site from the block-diagonalization the 10x10 Fock matrix build up from the localized natural orbitals of the converged CASSCF calculation of 3′

Diagonal matrix elements (eV) Fe^III(d^5) Fe^II(d^6)^a↓	−5.100	−4.977	−4.499	−4.167	−4.116
a Energies of the off-diagonal matrix elements are given in cm^−1.
−1.257	−39	1917	1210	616	105
−1.251	214	1419	−1705	−866	−88
−1.023	773	−77	−5512	186	−42
−0.829	−688	44	185	4724	632
−0.475	−15	36	−23	103	−1548

---

10. Conclusions

(1) Binuclear mixed valence transition metal complexes represent a challenge to multireference ab initio methods. In this work, we studied theoretically a mixed valence Fe^2.5+μTe₂Fe^2.5+ complex with the aim of explaining its unusual S = 3/2 ground state and interpreting its recently reported X-band EPR spectra and magnetic properties. We have found that using localized orbitals in setting up the many particle basis (the CFSs) of the Fe²⁺(d⁶)–Fe³⁺(d⁵) pair in 3′, one is able to correctly reproduce the antiferromagnetic (negative) sign of Heisenberg exchange integral J, while extension of the active space from CAS(11,10) to CAS(23,16), including six more doubly occupied 5p_x, 5p_y and 5p_z MOs on the bridging tellurium ligands, further increase its negative value. Without localization, using a guess of canonical MOs, a S = 9/2 ground state along with large violations from the Landé interval rule for the exchange coupled spin systems was predicted.

(2) Scalar relativistic localized molecular orbital optimizations of each spin-ladder state – S = 1/2, 3/2, 5/2, 7/2 and 9/2 separately, was found to be mandatory in order to achieve reliable and quantitative correct results; depending on the spin, each of these sets of orbitals relax differently. NEVPT2 corrections to the CASSCF reference energies decrease with increasing spin, thus contributing to larger negative values of J.

(3) Spin–orbit coupling calculations using the entire Hilbert space of a given spin multiplicity allowed to compute g-factor anisotropies, which for S = 3/2 exactly match g-factors extracted from simulations of the X-band EPR spectrum. This might be used for diagnostic purposes by comparing experimental and computed spin-Hamiltonian parameters. However, magnetic susceptibilities could not be reproduced in this way. To this end, we felt it necessary to account for vibronic coupling.

(4) Vibronic coupling has been considered in the approximation of a valence localized Fe²⁺μTe₂Fe³⁺ complex (Class I) excluding electron transfer. Using a scan of the ground state energy along the effective symmetry breaking normal mode q, both scalar relativistic and spin–orbit split S = 3/2 ground states were studied. We have found that spin–orbit coupling greatly reduces geometrical distortions and vibronic stabilizations over the scalar relativistic results. Numerical values for the linear vibronic coupling constant (V) and the Harmonic force constant (K) allow for quantification of the effect. From the SOC calculations the first derivatives of ZFS tensor D and the g-matrix of the spin 3/2 ground state spin Hamiltonian have been also determined. Extending the model with electron transfer we could simulate the temperature dependence of the magnetic susceptibility using a single electron transfer parameter B, whose value can hardly be estimated otherwise when using ab initio methods. Taking this parameter as freely variable, we can adjust its value from a best fit to the magnetic susceptibility data (B = 163 cm⁻¹). It is 4.5 time smaller than the one previously reported.¹⁰ Employing this value in the calculation of vibronic levels, energies of the lowest two excited states, 4.2 and 43 cm⁻¹ were computed in reasonable agreement with numbers extracted from the temperature dependence of the spin-relaxation times (5 and 32 cm⁻¹, respectively).

(5) The expediency and limitations of the DEVC model have been discussed in the light of concurrent electron transfer pathways with parameters extracted from correlated wave functions using local Fock matrices. It is pointed out, that values of B obtained from interpretation of experimental data should be regarded with caution, leading in some cases to some unexpected observations. Thus, the value of (B/J) derived from the S = 3/2 SOC ground state herein (B/J = 163/63) = 2.58) is outside the range of 3 < B/J < 5 predicted by the simple model (eqn (1.1)).

(6) The computational protocol developed in this work, particularly the vibronic analysis part of the zero-field splitting D-tensor and g-matrix are equally applicable to any binuclear transition metal dimer and is potentially applicable to all types of spin-relaxation mechanisms.

(7) Higher-nuclearity clusters can be potentially treated using the same computational recipe by setting up a composed model spin-Hamiltonian. A prerequisite for such applications are well defined metal based oxidation states and d-electron counts at each paramagnetic center (classes I and II mixed valence complexes). Taking as an example a triangular Mn₁Mn₂Mn₃ cluster with Mn₁, Mn₂ and Mn₃ in oxidation state III (high-spin d⁴, stabilized in a Jahn-Teller distorted geometry), diamagnetic substitution of two of the manganese atoms by Ga^III at the manganese sites, say Mn₂ and Mn₃, CASSCF/NEVPT2 calculations will give access to the zero-field splitting tensor D₁ (similar substitutions of Mn₁ and Mn₃ to D₂, and Mn₁, Mn₂ to D₃, respectively). Diamagnetic substitution of only one Mn^III center, and applying the computational CASSCF/NEVPT2 protocol described here will in turn give access to the Heisenberg exchange parameters J₁₂ (by substituting Mn₃ with Ga) and similar J₁₃ (by substituting Mn₂ with Ga) and J₂₃ (substituting Mn₁ by Ga). This allows one to set up the spin Hamiltonian of the triangular complex in the spin-uncoupled basis of the three S = 2, Mn^III spin centers (eqn (10.1)) which can be solved by any program dealing with high nuclearity spin clusters, such as MagnPack³⁶ and PHI.³⁷

H_sp-ham = − 2J₁₂Ŝ₁·Ŝ₂ − 2J₁₃Ŝ₁·Ŝ₃ − 2J₂₃Ŝ₂·Ŝ₃ + Ŝ₁D₁Ŝ₁ + Ŝ₂D₂Ŝ₂ + Ŝ₃D₃Ŝ₃ (10.1)

By processing in the same way, one can also treat four nuclear clusters related, say to the oxygen evolving model complex of Photosystem II of natural photosynthesis.

Open Access funding provided by the Max Planck Society.

Conflicts of interest

There are no conflicts to declare.

References

1. D. A. Pantazis, ACS Catal., 2018, 8, 9477.
2. C. Van Stappen, L. Decamps, G. E. Cutsail III, R. Bjornsson, J. T. Henthorn, J. A. Birrell and S. DeBeer, Chem. Rev., 2020, 120, 5005.
3. (a) R. Bjornsson, F. A. Lima, T. Spatzal, T. Weyhermuller, P. Glatzel, E. Bill, O. Einsle, F. Neese and S. DeBeer, Chem. Sci., 2014, 5, 3096; (b) J. K. Kowalska, J. T. Henthorn, C. Van Stappen, C. Trncik, O. Einsle, D. Keavney and S. DeBeer, Angew. Chem., Int. Ed., 2019, 58, 9373.
4. L. Noodleman and D. A. Case, Inorg. Chem., 1992, 38, 423.
5. L. Noodleman, J. Chem. Phys., 1981, 74, 5737.
6. B. Benediktsson and R. Bjornsson, Inorg. Chem., 2020, 59, 11514.
7. B. Benediktsson and R. Bjornsson, J. Chem. Theory Comput., 2022, 18, 1437.
8. (a) J.-P. Malrieu, R. Caballol, C. J. Calzado, C. de Graaf and N. Guihéry, Chem. Rev., 2014, 114, 429; (b) R. Maurice, C. de Graaf and N. Guihéry, J. Chem. Phys., 2010, 133, 084307; (c) R. Maurice, R. Bastardis, C. de Graaf, N. Suaud, T. Mallah and N. Guihéry, J. Chem. Theory Comput., 2009, 5, 2977; (d) R. Maurice, N. Guihéry and C. de Graaf, J. Chem. Theory Comput., 2010, 6, 55; (e) R. Maurice, A. M. Pradipto, N. Guihéry, R. Broer and C. de Graaf, J. Chem. Theory Comput., 2010, 6, 3092.
9. R. Maurice, R. Broer, N. Guihéry and C. de Graaf, Handbook of Relativistic Quantum Chemistry, ed. W. Liu, Springer-Verlag, Berlin Heidelberg, 2016, p. 765.
10. J. Henthorn, G. E. Cutsail, T. Weyhermüller and S. DeBeer, Nat. Chem., 2022, 14, 328.
11. M. B. Robin and P. Day, Mixed Valence Chemistry, Advances in Inorganic Chemistry and Radiochemistry, 1968, vol. 10, pp. 247–422 DOI:10.1016/S0065-2792(08)60179-X.
12. P. W. Anderson and H. Hasegawa, Phys. Rev., 1955, 100, 675.
13. J.-J. Girerd, J. Chem. Phys., 1983, 79, 1766.
14. G. Blondin, S. Borshch and J.-J. Girerd, Comments Inorg. Chem., 1992, 12, 315.
15. S. B. Piepho, E. R. Krausz and P. N. Schatz, J. Am. Chem. Soc., 1978, 100, 2996.
16. Trough out this work it was assumed, that during the process the Fe–Fe distance does not change. This is justified to the extent of the constraints imposed by the crystal packing.
17. A. Domingo, C. Angeli, C. de Graaf and V. Robert, J. Comput. Chem., 2015, 36, 861.
18. F. Neese, Wiley Interdiscip. Rev.:Comput. Mol. Sci., 2018, 8, e1327.
19. F. Neese, F. Wennmohs, U. Becker and C. Riplinger, J. Chem. Phys., 2020, 152, 224108(1–18).
20. N. Spiller, V.-G. Chilkuri, S. DeBeer and F. Neese, Eur. J. Inorg. Chem., 2020, 1525.
21. M. Atanasov, D. Ganyushin, K. Sivalingam and F. Neese, A Modern First-Principles View on Ligand Field Theory Through the Eyes of Correlated Multireference Wavefunctions, Molecular Electronic Structures of Transition Metal Complexes II, ed. D. M. P. Mingos, P. Day and J. P. Dahl, Structure and Bonding, 2012, vol. 143, p. 149.
22. M. Atanasov, D. Aravena, E. Suturina, E. Bill, D. Maganas and F. Neese, Coord. Chem. Rev., 2015, 289–290, 177.
23. S. K. Singh, J. Eng, M. Atanasov and F. Neese, Coord. Chem. Rev., 2017, 344, 2.
24. B. O. Roos, P. R. Taylor and P. E. M. Sigbahn, Chem. Phys., 1980, 48, 157.
25. C. Angeli, S. Evangelisti, R. Cimiraglia and D. A. Maynau, J. Chem. Phys., 2002, 117, 10525.
26. C. Angeli, R. Cimiraglia, S. Evangelisti, T. Leininger and J.-P. Malrieu, J. Chem. Phys., 2001, 114, 10252.
27. C. Angeli, R. Cimiraglia and J.-P. Malrieu, Chem. Phys. Lett., 2001, 350, 297.
28. C. E. Schäffer and C. K. Jørgensen, Mol. Phys., 1965, 9, 401.
29. C. E. Schäffer, Struct. Bonding, 1968, 5, 68.
30. F. Neese, J. Chem. Phys., 2005, 122, 034107.
31. D. Ganyushin and F. Neese, J. Chem. Phys., 2006, 125, 024103.
32. F. Weigend and R. Ahlrichs, Phys. Chem. Chem. Phys., 2005, 7, 3297 , recontracted for DKH2 by D. A. Pantazis.
33. D. A. Pantazis, X.-Y. Chen, C. R. Landis and F. Neese, J. Chem. Theory Comput., 2008, 4, 908.
34. Because the two Kramers doublets originating from S = 3/2 ground state share a common spatial wave function, the excited Kramers doublet at 4.2 cm⁻¹ yields identical results.
35. Yet another source of the magnetic anisotropy is spin-spin coupling (SSC) see E. R. Davidson and A. E. Clarc, Phys. Chem. Chem. Phys., 2007, 9, 1881 Compared to SOC, SSC usually yields much smaller (yet non-negligible) contributions to the anisotropy, see F. Neese, J. Am. Chem. Soc., 2006, 128, 10213.
36. J. J. Borrás-Almenar, J. M. Clemente-Juan, E. Coronado and B. S. Tsukerblat, J. Comput. Chem., 2001, 22, 985.
37. N. F. Chilton, R. P. Anderson, L. D. Turner, A. Sonchini and K. S. Murray, J. Comput. Chem., 2013, 34, 1164.

---

Footnote

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

---