Predictive Computational Analysis of Magnetic Exchange Interactions in Gd³⁺ Complexes with Biradical Ligands

1. Introduction

Addressing the requirements of evolving technological paradigms necessitates the exploration of material spaces beyond the established canon. Molecular magnetism is a field where chemistry, physics, and materials science converge to explore paramagnetic molecular systems. After a century as a testbed for d- and f-electron physics, magnetic molecules are emerging as key components for technologies like data storage, spintronics, and quantum computing. Their main "ace in the hole" of this field is the ability to control magnetism through chemical synthesis. This chemical tuning knob allows for the creation of countless variations in the "structure-property" relationship. Even with significant accumulated knowledge in the field, synthetic chemists increasingly require clear guidelines. These are essential for selecting suitable paramagnetic precursors and rationally assembling them into functional molecular architectures with predetermined magnetic characteristics. Selecting an appropriate computational methodology for prediction is intrinsically linked to the dimensionality of the molecular magnetic material. Presently, research in this area is advancing along several distinct fronts. The rational design of novel molecular magnetic materials requires robust computational tools to screen prospective building blocks and assembly strategies prior to synthesis. Such preliminary assessment is crucial and can be executed via complementary theoretical approaches [1,2,3,4,5,6,7].

One powerful strategy is data-driven prediction, which utilizes spatial-electronic cooperativity via machine learning [8,9,10,11]. This inherently bottom-up approach mines complex, often non-intuitive relationships between a molecule's spatial structure and its emergent electronic properties to identify promising starting molecular objects. Conversely, the same computational paradigm facilitates a top-down design framework by predicting the properties of complex, multi-dimensional magnetic systems for rational selection of optimal molecular constituents [11].

At present, computational methods have achieved a level of accuracy sufficient for reliable modelling the magnetic behavior of synthesized single-molecule magnets (SMMs), which exhibit slow magnetization relaxation, a property stemming from magnetic anisotropy that creates an energy barrier (U_eff) between magnetic states with opposite sign of spin or total angular momentum projections [12]. As the key challenges in developing quantum magnets are to increase this energy barrier and prevent magnetization relaxation [13], the preliminary theoretical "screening" of the electronic and magnetic parameters of candidate molecular systems is of critical importance.

When designing SMMs, the main characteristics are the spin configuration of the ground state (ferro- or antiferromagnetic) and the energy difference, ΔE_J, between the antiferromagnetic and ferromagnetic spin states, which is determined by the exchange interaction between the paramagnetic centers. The greater the magnitude of this interaction, the lower the probability of magnetic relaxation. Quantum tunneling of magnetization (QTM) can be suppressed by sufficiently strong magnetic coupling. The latter acts as an internal magnetic field that lifts the degeneracy within a spin manifold. Notably, numerous complexes with a U_eff still exhibit a low blocking temperature. One key reason for this is rapid QTM [14]. Therefore, from the perspective of potential for designing molecular magnets, coordination compounds of paramagnetic metal ions with paramagnetic ligands are particularly attractive. This is especially valuable for designing SMMs based on complexes of organic radicals and f-elements, as in the latter the unpaired electrons reside on shielded orbitals [15,16].

2. The Purpose and Objectives of the Study

2.1. A Short Presentation of Donor–Acceptor Biradicals

The prior studied donor−acceptor (DA) biradicals (BRs), (SQNN)⁻ and (SQPhNN)⁻, (semiquinone–nitronylnitroxides) were isolated exclusively as d-metal complexes. The complexes Tp^CumM(SQNN) and Tp^CumM(SQPhNN), where M is a divalent 3d transition metal, were synthesized by reacting the corresponding catechol (Cat) precursors, CatNN or CatPhNN, with the starting complexes Tp^CumM(OH) under an inert atmosphere, followed by oxidation with atmospheric oxygen (see Scheme 1) [17,18,19]. Here, Tp^Cum is a tripodal ligand specifically designed [20] to protect the semiquinone coordination site from air oxidation. The DA BRs feature high-spin ground states (S = 1) arising from strong ferromagnetic exchange coupling between radical spins. This magnetic configuration enables high-resolution investigation of their electronic structures using techniques such as EPR, magnetic susceptibility, magnetic circular dichroism, and variable-temperature electronic absorption spectroscopy. Such organic heterospin systems can exhibit inherently stronger ferromagnetic coupling than their homospin analogues. This principle establishes a key structure–property relationship: designing molecules from a donor–acceptor perspective is a powerful strategy for achieving high-spin ground states [21]. These pioneering studies established an entire research direction in the coordination chemistry of d-metal complexes with DA BRs [22,23] and revealed their potential application in quantum information science [24] and molecular electronics [25].

The prior studied donor−acceptor (DA) biradicals (BRs), (SQNN)⁻ and (SQPhNN)⁻, (semiquinone–nitronylnitroxides) were isolated exclusively as d-metal complexes. The complexes Tp^CumM(SQNN) and Tp^CumM(SQPhNN), where M is a divalent 3d transition metal, were synthesized by reacting the corresponding catechol (Cat) precursors, CatNN or CatPhNN, with the starting complexes Tp^CumM(OH) under an inert atmosphere, followed by oxidation with atmospheric oxygen (see Scheme 1) [17,18,19]. Here, Tp^Cum is a tripodal ligand specifically designed [20] to protect the semiquinone coordination site from air oxidation. The DA BRs feature high-spin ground states (S = 1) arising from strong ferromagnetic exchange coupling between radical spins. This magnetic configuration enables high-resolution investigation of their electronic structures using techniques such as EPR, magnetic susceptibility, magnetic circular dichroism, and variable-temperature electronic absorption spectroscopy. Such organic heterospin systems can exhibit inherently stronger ferromagnetic coupling than their homospin analogues. This principle establishes a key structure–property relationship: designing molecules from a donor–acceptor perspective is a powerful strategy for achieving high-spin ground states [21]. These pioneering studies established an entire research direction in the coordination chemistry of d-metal complexes with DA BRs [22,23] and revealed their potential application in quantum information science [24] and molecular electronics [25].

2.2. Theoretical Notes on Semiquinone — Nitronylnitroxide Donor-Acceptor Biradicals (BRs)

2.2.1. Complexes of the (SQNN)⁻ and (SQPhNN)⁻ BRs with 3d Metals

The idea that DA interactions could underpin molecular electronics began with a seminal assumption a donor–bridge–acceptor (DBA) molecular rectifier [26]. DA BRs of the Tp^CumM(SQNN) type (where M is Zn²⁺) are considered to be ground-state analogues for photo-induced charge-separated states. Spectroscopic and computational studies have extensively evaluated the pairwise exchange in the complexes Tp^CumM(SQNN) and Tp^CumM(SQBNN) to understand excited-state contributions to the ground state Explaining the ground-state exchange requires models that incorporate excited-state configurations, going beyond the active-electron approximation [27,31,32,33].

The complete set of SCRXD-derived bond lengths for the (SQNN)⁻ metal complexes is fully consistent with a +2 oxidation state of the Mn, Co, and Ni ions, supporting an M^IISQ electronic configuration rather than M^IIICat [18,19]. Moreover, the bond lengths of NN part differ from those in the corresponding CatNN compounds, further indicating an interaction between the SQ and NN moieties. This leads to a description with a very strong intraligand SQ–NN exchange (J_SQ–NN) and an effective spin Hamiltonian governed by a single metal–ligand exchange parameter, J_M–SQNN [19,34]. The exchange interactions between the SQNN⁻ organic BR ligand and the divalent transition metal ion are antiferromagnetic for Tp^CumNi^II(SQNN) (J = –87.8 cm⁻¹) and Tp^CumMn^II(SQNN) (J = –41.3 cm⁻¹), but ferromagnetic for Tp^CumCu^II(SQNN) (J = +75.6 cm⁻¹). The exchange parameter J can be expressed as the sum of ferromagnetic and antiferromagnetic contributions: J = J_F + J_AF [35,36]. In general, the ferromagnetic term J_F is small, and J becomes ferromagnetic only when J_AF is negligible. Accordingly, trends in metal–ligand antiferromagnetic exchange can be discussed solely in terms of J_AF. Hence, J values indicate that the magnitude of the antiferromagnetic exchange between (SQNN)⁻ and M^II decreases in the order: Ni > Mn > Cu, with the antiferromagnetic contribution being very weak in the copper complex. To assess qualitative magnetic trends in related systems, relative exchange couplings can be estimated from extended Hückel calculations [37]. For a dimer with one unpaired electron per site, magnetic orbitals φₐ and φ_b give J_AF = –2(t_ab)²/U_eff, where (t_ab)² is the hopping integral between them [31,35,38].

2.2.2. Frontier Molecular Orbital (MO) Description for (SQNN)⁻ and (SQPhNN)⁻ DA BRs

Density functional theory (DFT) calculations on the Tp^CumM(SQNN) BR and its fragments (Figure 1) reveal key orbital characteristics [27,28,34]. The singly occupied (SOMO) of the SQ fragment is delocalized across its π-system, whereas the orbital, SOMO_NN, has a node at the bridging carbon atom. This nodal structure prevents direct communication between the SQ and NN constituents of SOMOs active electron wavefunctions leading to that overlap integral ⟨SOMO_NN|SOMO_SQ⟩ ≈ 0, and thus direct exchange between the two radical spins. The MO picture shows that the SOMO2 orbital of the combined SQNN system arises from admixture the SOMO_SQ with the lowest unoccupied LUMO_NN. This mixing creates a net positive overlap density via the NN fragment's contribution to both spin-bearing SOMO1_SQNN and SOMO2_SQNN orbitals, resulting in a ferromagnetic contribution to the ground-state exchange. Owing to the intense ferromagnetic exchange in the DA BR, magnetic susceptibility measurements yielded solely a lower-limit estimate for the intraligand exchange parameter, yielded only a lower limit for the intraligand exchange of J_SQ–NN > +300 cm⁻¹ (assuming Ĥ =−2J(Ŝ_SQ∙Ŝ_NN)). Higher value of +550 cm⁻¹ was found for J_SQ–NN from the EPR spin densities [29].

The MO diagram (Figure 2) based on the results of bonding calculations displays the main SQ, Ph (phenilen), and NN in frontier MO interactions in Tp^CumM(SQPhNN). The phenylene e₁ orbital have the suitable symmetry to mix with the SOMO_SQ and the LUMO_NN. Strong orbital mixing between the LUMO_NN and Ph(e₁) orbitals gives a stabilized LUMO_SQPhNN orbital with substantial LUMO_NN–Ph(e₁) bonding character. The energetic stabilization of the LUMO_SQPhNN orbital allows for stronger interaction with SOMO_SQ orbital. This stabilization enhances its interaction with the SOMO_SQ. Consequently, the resulting SQ–NN–Ph mixing clearly facilitates delocalization of positive spin density from the SQ donor onto both the Ph–ring and the NN nitrogen atoms. For the phenylene-bridged DA BR in its complex Tp^CumZn(SQPhNN), magnetic measurements yielded an exchange constant of J_SQ–PhNN = +100 cm⁻¹ from variable-temperature data analysis [19]. Note that the admixture of the LUMO_NN into the SOMO2_SQNN orbital points to a significant contribution of excited states to the ground state magnetic exchange, a concept that forms the basis of the valence bond configurational interaction approach.

2.2.3. Valence Bond Configuration Interaction Model (VBCIM)

For organic nondisjoint type BRs [39], intramolecular exchange coupling is typically described within the active electron approximation, involving only two magnetic orbitals. However, an alternative treatment for DA systems is provided by the valence bond configuration interaction model (VBCIM) [28,40], opposed to orbital description of the system. Thus, the VBCI approach incorporating parameters accessible from spectroscopic data, enables assessment of excited-state contributions to ground-state exchange. In this model, covalency is manifested by configurational mixing of charge-transfer excited states into the ground state via the electronic coupling matrix element, H_DA. Crucially, the VBCIM primarily highlights the key role of an SOMO_SQ → LUMO_NN charge transfer configuration (CTC). Its mixing into the electronic ground state is necessary to facilitate stabilization of the high-spin triplet, S = 1 ground state in Tp^CumM(SQNN) and Tp^CumM(SQPhNN) BRs. The intraligand SQ→NN CT-band in Tp^CumM(SQNN) is the organic analogue of metal-to-metal charge transfer (MMCT) bands observed in transition metal dimers [41,42]. One advantageous feature of these DA BRs is that there is only a single dominant superexchange pathway leading to a straightforward determination of H_DA for a DA BR within the VBCI formalism [41].

The electronic origin of the strong ferromagnetic exchange in Tp^CumZn(SQNN) [27,28] was interpreted using VBCIM shown in Figure 3. The active space comprises three orbitals: SOMO_SQ, SOMO_NN and LUMO_NN orbitals (SOMO¹_NNSOMO¹_SQLUMO⁰_NN; for 110 ground-state configuration). The coupling constant J_SQ-NN is related to spectroscopic parameters via: ${2J}_{SQ-NN}={\displaystyle \frac{{2H}_{SQNN}^2{K}_0}{\left(U^2-K_0^2\right)}}$ where H_SQNN is the electronic coupling matrix element connecting SOMO_SQ with LUMO_NN, K₀ is the single-site exchange integral that splits the singlet and triplet configurations arising from a one-electron SOMO_SQ→LUMO_NN promotion, and U is the mean charge transfer energy for this process[34]. The ferromagnetic exchange interaction in DA BR is maximized when H_SQNN and K₀ are both large, and the average energy of the intraligand DA CT energy (U–K₀) is low.

The electronic coupling parameter H_SQ-NN (or H_DA) can be determined experimentally for DA BRs by combining magnetic susceptibility and EPR measurements to obtain J, and optical spectroscopies to determine both U and K₀. However, as noted in the original work [27], the experimental determination of the ³CTC–¹CTC splitting (2K₀) is particularly challenging, as it requires the observation of the spin-forbidden ³CTC–¹CTC intraligand charge-transfer transition. Thus, evaluating H_DA in a simple DA biradical requires three parameters:J, obtainable from magnetic susceptibility/EPR), and U and K₀ (derivable from optical spectroscopy).

2.2.4. Problem Statement for a Predictive Computational Study

The main motivation for this predictive investigation is to pursue a synthetic quest for the preparation of coordination compounds comprising DA BRs and rare earth ions. The lanthanide ions are anticipated to impart high magnetic anisotropy to the molecules, owing to their strong intrinsic spin-orbit coupling. Moreover, we previously synthesized neutral complexes of the Tp₂LnSQ^DB type, where SQ^DB denotes tripodal ligand, 3,5-di-tertbutylsemiquinonate [43,44,45], whose coordination site and charge are identical to those of the DA BRs. Moreover, the class of Tp₂LnSQ^DB compounds still remains a subject of active research interest [46,47,48].

Regarding the magnetic properties of the Tp₂LnSQ^DB complexes, the coupling constant between Gd and SQ^DB, determined from magnetic data [43], reaches −5.7 cm⁻¹ (2J formalism, which will be used for all coupling constants throughout the text). In contrast, a recent estimate of the J value, obtained from both spectroscopic (inelastic neutron scattering) and magnetic studies, gave a value of −1.8 cm⁻¹ for Tp₂ErSQ^DB [47].

We performed computational studies on the hypothetical structures of gadolinium complexes with DA radicals, Tp₂Gd(SQNN) (1) and Tp₂Gd(SQPhNN) (2). This was done to evaluate how the magnitude of the metal–radical magnetic exchange interaction would change in the case of coordination of a semiquinone–nitronylnitroxide BR, (SQNN)⁻ or (SQPhNN)⁻ to an Ln–center.

3. Computational Methods

3.1. Generation of Molecular Structures for Computational Calculations

The starting molecular structures were obtained from the SCXRD geometry of Tp₂Gd(SQ) [43] by replacing the t-Bu group at the 5-position with the –NN or –PhNN moiety. Partial geometry optimizations of Tp₂Gd(SQNN) (1) and Tp₂Gd(SQPhNN) (2) were performed for the triplet state at the (U)BP86 [49,50]/def2-TZVP [50] level. To reduce the computational costs, the Gd³⁺ was replaced by a diamagnetic Y³⁺ while retaining the geometric constraints of the Tp₂Y(SQ) fragment; the def2-ECP [51] was used for Y³⁺. The identity resolution (RI) method with the auxiliary def2/J basis set was applied [52].

3.2. Methods and Basis Sets Employed in the Calculations

The SA-CASSCF/QD-NEVPT2 [53] calculations with the different number of excited states were performed for non-coordinated forms of biradicals (BRs) and their complexes under study. The electronic structure calculations were accomplished at the state averaged (SA) CASSCF level [54,55] taking into account dynamic correlation at NEVPT2 [56] level for one dectet, two octet and one sextet states for 1 and 2. The active space of SA calculations consists of 11 or 13 molecular orbitals (MOs) with 11 or 13 electrons on them correspondingly; the orbital set is composed of seven Gd³⁺ 4f atomic orbitals (AOs), HOMO_SQ, SOMO_SQ, LUMO_SQ and SOMO_NN for 1 and 2. The complete sets of MOs obtained from the calculations are provided in Table S1 ÷ S4.

Additionally, the SA-CASSCF/NEVPT2 calculations were performed for one triplet and one singlet states of noncoordinated forms of (SQNN)⁻ and (SQPhNN)⁻ BRs; as well as for the optimized geometries of complexes Tp₂Y(SQNN) (1Y) and Tp₂Y(SQPhNN) (2Y), constructed from 1 and 2 by substitution the Gd center by Y³⁺ ion their optimized geometries. The active space was constructed by HOMO_SQ, SOMO_SQ, LUMO_SQ, and SOMO_NN with six electrons on them. To determine the electronic coupling matrix element in VBCI formalism [41], the SA-CASSCF(6,5)/QD-NEVPT2 [53] calculations were performed with 10 triplet and 13 singlet states for two BR-anions, their yttrium complexes and with 10 triplet and 14 singlet states for Tp^CumZn(SQNN) (3) and the same active space. The XRD geometry of 3 was used for all calculations [44]. As well, the SA-CASSCF(10,10)/NEVPT2 calculations were performed for nonet and septet states of complex Tp^CumGd(SQ) where seven 4f-AOs of Gd³⁺ and HOMO_SQ, SOMO_SQ, LUMO_SQ made up the active space.

The scalar relativistic effects were taken into account using DKH2 Hamiltonian [57]. The SARC2-DKH-QZVP and auxiliary SARC2-DKH-QZVP/JK basis sets [58] for Gd³⁺ and the DKH-def2-TZVP(-f) [59] together with the def2/JK [60] and def2-TZVP/C [61] basis sets for other atoms were used. The RI-JK approximation [60,62] was applied with DefGrid2 [63]. ORCA 6.0 and 6.1.0 program package [64] were used for all computations, including the calculation of the temperature dependence of magnetic susceptibility. The approximation of the state energies was carried out in PHI program [65] using model spin-Hamiltonian:

$$\mathcal{H}={\mu }_{B}g\mathbf{B}·\left({\mathit{S}}_{Gd}+{\mathit{S}}_{SQ}+{\mathit{S}}_{NN}\right)-2J_{SQNN}{\mathit{S}}_{SQ}·{\mathit{S}}_{NN}-2J_{GdSQ}{\mathit{S}}_{Gd}·{\mathit{S}}_{SQ}-2J_{GdNN}{\mathit{S}}_{Gd}·{\mathit{S}}_{NN}$$

(1)

4. Discussion

The electronic structures of complexes Tp₂Gd(SQNN) (1) and Tp₂Gd(SQPhNN) (2) were computed at the SA-CASSCF/NEVPT2 level to investigate the exchange interaction between the Gd³⁺ ion (S_Gd = 7/2) and the ferromagnetically coupled heterospin BR (S_BR = 1), which affords one decet, two octet, and one sextet state. This number of states is sufficient to obtain information about the exchange coupling between Gd³⁺ ion and the biradical as well as to assess the value of exchange interaction between the SQ and NN constituents of the DA biradicals. The energies of the states for 1 and 2 calculated at SA-different levels are presented in Table 1. The active space orbitals for these calculations are presented in Table 2 for both complex 1 and complex 2. The order of the lowest three states at the SA-CASSCF and SA-CASSCF/NEVPT2 differs for both complexes. This indicates that SA-CASSCF gives an incorrect sign and underestimates the strength of the exchange interactions between Gd³⁺ and SQNN^– or SQPhNN^– entities. SA-CASSCF thus predicts a ferromagnetic nature of the interaction whereas the inclusion of dynamic of dynamic correlation leads to the correct antiferromagnetic coupling, which is consistent with the experimental magnetic data. The position of the second excited octet state correlates with the strength of the exchange interaction between SQ and NN moieties of the biradicals and corresponds to their singlet forms.

The calculated temperature dependencies of magnetic susceptibility for complexes Tp₂(GdSQNN) (1) and Tp₂Gd(SQPhNN) (2) are shown in Figure 4. Their behavior correlates with the energy levels (Table 1), showing a gradual increase with temperature for both complexes. The overall shape of the χT curves for complexes 1 and 2 in the temperature range 2–300 K closely resembles the magnetic behavior of complex Tp₂Gd(SQ) (3): a gradual decrease in χT values with decreasing temperature, with a step-like feature in the low-temperature region [43]. Measured at 245 K, experimental χT value for 3 is 7.9 emu·K·mol⁻¹. This is slightly lower than the theoretical value of 8.250 emu·K·mol⁻¹ expected for uncoupled spin carriers: S_Gd = 7/2, S_SQ = 1/2 and g = 2 (7.875 emu·K·mol⁻¹ for the free Gd³⁺ ion [66] and 0.375 emu·K·mol⁻¹ for one S = 1/2 spin from SQ^–). In the temperature range of 0 ÷ 20 K, Tp₂Gd(SQ) shows a decrease in χT to 5.98 emu·K·mol⁻¹ at about 10–15 K, corresponding to a state with S = 3. The reduced χT values were also observing for computed in our study complexes 1 and 2. At 300 K, they reach 8.627 and 8.545 emu·K·mol⁻¹, respectively. This is somewhat lower than the value of 8.875 emu·K·mol⁻¹ expected for a sum of uncoupled paramagnetic contributions (S_Gd = 7/2, g_Gd = 2.00 and S_BR = 1), but it is very close to 8.625 emu·K·mol⁻¹, which is in better agreement with the uncorrelated S_Gd = 7/2, S_SQ = 1/2, and S_NN = 1/2 situation. This is most probably caused by the admixture of excited states to the ground spin state [34]. Furthermore the reduced χT (T = 300 K) values (compared to expected theoretical values) have previously been observed for Gd³⁺ radical complexes [67,68,69,70,71,72,73,74,75,76]. However, for both complexes 1 and 2 the χT value at 5 K is 4.399 emu·K·mol⁻¹, which corresponds well to the theoretical value of 4.375 emu·K·mol⁻¹ for an S = 5/2 state. The almost one and a half times discrepancy in the characteristic temperatures in magnetic behavior for complex 3 compared to that for complexes 1 and 2 arises from the unequal exchange coupling strengths of the Gd–SQ and Gd–SQNN exchange interactions. The further decrease down to 0 K can be attributed to the influence of zero-field splitting (ZFS) for all complexes.

However, the χT product for complex 1 continues to rise with temperature above 300 K, whereas for complex 2 it approaches saturation. Thus, the magnetic susceptibility tends to higher values for complex 1, which is consistent with the interaction between S_Gd = 7/2 and S_BR = 1. The deviation between the calculated and anticipated χT values at 300 K may arise from several factors. As reported in reference [17], the χT value for Tp^CumZn(SQNN) at 300 K is approximately 0.98 emu·K·mol⁻¹. Hence, in complex 1, the contribution of Gd³⁺ is equal 7.647 emu·K·mol⁻¹. This corresponds to g_Gd = 1.97, which differs from the isotropic value of 2.0. This difference could be due to weak magnetic anisotropy or to uncertainties in the experimental data processing that affect the g-value determination. Alternatively, antisymmetric exchange interactions between Gd³⁺ and the SQNN^– biradical may also be operative [19] due to the admixture of excited states to the ground spin state [34]. For complex 2, the χT value at 300 K is closer to that expected for a paramagnetic BR and corresponds to g_Gd = 1.99, which is consistent with the weaker exchange interaction in the (SQPhNN)^– moiety compared to the BR fragment in complex 1.

To retrieve the exchange interaction values, the approximation of the calculated energies was performed using the model spin-Hamiltonian (1). As the exchange coupling is strong in the DA biradicals, the parameters J_Gd–SQ and J_Gd–NN are degenerate, and a more appropriate parameter, J_Gd–SQNN, is considered. The calculated exchange interaction parameters for the metal–biradical and intra-biradical interactions obtained in this study, together with those reported for related systems in the literature, are summarized in Table 3. Complexes 1 and 2 exhibit similar strength of Gd³⁺–biradical exchange interaction, whereas the parameter J_SQ–NN is an order of magnitude higher than J_SQ–PhNN. Additional SA-CASSCF(10,10)/NEVPT2 calculations performed for complex Tp₂Gd(SQ) yielded an exchange interaction value J_Gd–SQ = −5.4 cm⁻¹ using Ĥ = –2JŜ₁Ŝ₂ formalism, which correlates well with the value determined from experiment [43]. Thus, the interaction of the Gd³⁺ ion with the semiquinone site of strongly coupled biradicals leads to a twofold decrease in the exchange parameter compared with that in the monoradical complex Tp₂Gd(SQ).

Compared with the DFT broken-symmetry approach, the CASSCF method offers a more detailed description of the exchange coupling in terms of configurational interaction (CI) contributions. To elucidate the contributions of the NN and PhNN moieties to the intraligand exchange interactions and to understand the effect of Gd³⁺ coordination on the coupling constants in the biradicals, SA-CASSCF(6,5)/NEVPT2 calculations were performed on the (SQNN)⁻ and (SQPhNN)⁻ species, as well as on the diamagnetic metal complexes Tp₂Y(SQNN) (1Y) and Tp₂Y(SQPhNN) (2Y). The calculated with different number of triplet and singlet states exchange parameters for the BR^– and their Y³⁺ complexes are listed in Table 4. The computations performed for the DA^– biradicals gave the coupling constants J_SQ-NN = 633.7 cm⁻¹ and J_SQ-PhNN = 109.9 cm⁻¹, whereas the same parameters for the yttrium complexes 1Y and 2Y are 618.4 cm⁻¹ and 63.9 cm⁻¹, respectively, while the parameter J_SQ–NN for complex 3 achieving 930.9 cm⁻¹. The latter value having been recently estimated by ab initio computation [34], which gave a value of 743.3 cm⁻¹. This difference between the J-values is related to the geometry of the biradical. The torsion angle is 17–19° in BR and complex 1Y, 0.5–1.5° for complex 3, compared with 2.0° reported in reference [34]. This discrepancy arises from the different geometry optimization approaches used: full optimization of the previously studied complex, partial optimization of the nitronylnitroxide fragment for complex 1, and the use of the experimental XRD geometry for complex 3. Complex 1Y displays a slight deviation in the exchange interaction J_SQ-NN value from corresponding uncomplexed (SQNN)⁻ form, whereas for complex 2Y, the parameter J_SQ–PhNN is considerably lower than that for free (SQPhNN)⁻ (Table 3). This discrepancy likely arises from a stronger redistribution of spin density from the SQ fragment toward the PhNN moiety, combined with the coordination effect of the Y³⁺ ion. At the same time, complex Tp^CumZn(SQNN) (3) exhibit a 1.5-fold increase in J. This change is attributed to the difference in dihedral (or torsion) angle, which is about 0.5–1.5° for 3 and lies in the range of 17–19° for complex 1Y and (SQNN)^–.

Comparative spin density maps for complexes 1 and 2 and their Y³⁺ analogues are displayed in Figure 5. Additionally, Table S5 summarizes the natural charge and Mulliken spin density distributions on the constituent fragments (Gd³⁺, SQ– and –NN) of complexes 1 and 2. The natural charge distributions are virtually identical across all complexes in all lowest states. The most pronounced variations in spin density are observed on the SQ and NN moieties of the biradical upon transition from the high-spin to the low-spin states, while the spin density on Gd³⁺ undergoes only a minor decrease. This behavior confirms that the spin inversion occurs predominantly on the biradical fragments, whereas the spin density on the metal ion remains essentially conserved. The small changes in spin density on Gd³⁺ are attributable to spin back-donation from the SQNN fragment in the sextet state.

To elucidate the influence of various factors on the strength of exchange interactions (Table 3), the energy gaps for the triplet and singlet spin states (Figure 7) were calculated at the SA-CASSCF(11,11)/NEVPT2 level for complexes 1 and 2, and at the SA-CASSCF(6,5)/NEVPT2 level for the anion-biradicals (SQNN)⁻, (SQPhNN)⁻, and their metal complexes with triplet and singlet states. For both biradicals, the anionic forms exhibit higher energy gaps compared to the coordinated forms 1Y and 2Y, reflecting a redistribution of spin density upon coordination. However, in the case of the Gd³⁺ complexes, coordination of the biradical leads to completely different results. For complex 1, the second excited octet state, corresponding to the singlet form of the SQNN ligand, has a higher energy than the singlet state for 1Y (166 vs. 128 cm⁻¹). A possible reason for this could be spin-density donation from the metal to the ligand (M→L). At the same time, compound 2 shows an intermediate energy value for the corresponding excited spin states of (SQPhNN)⁻ and 2Y. It is noteworthy that M→L spin-density donation results in stronger SQ–NN exchange interactions in complex 2 relative to complex 2Y, where the metal ion is diamagnetic (71 cm⁻¹ vs. 64 cm⁻¹, respectively). Nevertheless, the phenylene bridge renders this interaction less efficient than that in compound 1. Furthermore, the largest energy gap for TpᶜᵘᵐZn(SQNN) (3) is due to the smallest torsional angle in the SQ–NN ligand among the compounds shown in Figure 6. For TpᶜᵘᵐZn(SQPhNN) (4), the torsional angle in SQ–Ph is slightly larger than that for compound 2 (see Table 3). At the same time, the Ph–NN angle in (4) is almost twice as small as in 2. This results in an energy-gap value that is slightly smaller than that for (SQPhNN)⁻ and larger than those for 2 and 2Y, due to stronger coupling between the SQ and NN monoradicals.

According to the previously described VBCI formalism [41], the electronic coupling matrix element (H_DA) can be determined using the spectroscopic observables, namely, the positions of the excited states determined by the intraligand charge transfer configuration (CTC) with SQ_SOMO → NN_LUMO one-electron promotion. To take into account the desirable charge transfer state, the SA-CASSCF/QD-NEVPT2 calculations were performed for the previously mentioned complexes with increased number of states.

Firstly, SA-CASSCF(6,5)/QD-NEVPT2 calculations were performed for the anion forms of the two biradicals, using 10 triplet and 13 singlet states. This number of states was chosen to target the desired charge-transfer configurations (CTC) corresponding to electron transfer from SOMO_SQ to LUMO_NN, with a dominant contribution to the wavefunctions of the excited states. When a different number of singlet and triplet states was used, the CTC became difficult to distinguish because of strong admixture with other configurations corresponding to other types of charge transfer or intraligand transitions. The energies of the desired charge transfer states are presented in Table 5 (³^,¹CTC).

For the (SQNN)⁻ specie, the target CT states correspond to the 3^rd triplet and 5^th singlet excited states at 18115 and 25028 cm⁻¹, with electron-promotion configurational contributions of 87% and 80%, respectively. For the (SQPhNN)⁻ biradical, the corresponding CT states are the 1st triplet and 5th singlet excited states at 12902 and 16639 cm⁻¹, with CT configurational contributions of 91% and 82% [77]. The calculated H_DA values are presented in the Table 6.

Further, the ab initio calculations at the same level of theory were performed for the complexes Tp₂YSQNN (1Y) and Tp₂YSQPhNN (2Y) with the same number of states. In the case of 1Y, the desirable states are the 4^th triplet and 7^th singlet excited states at 24538 and 28618 cm^-1 with 63% and 81% of CTC, correspondingly. The complex 2Y gives the 3^rd triplet and 6^th singlet excited states at 20344 and 23424 cm^-1 having the CTC contribution of 65% and 88%, correspondingly.

According to the obtained results, the excited-state energies are lower for the non-coordinated BRs than for complexes 1Y and 2Y. This correlates with the smaller energy gap between SOMO_SQ and LUMO_NN for the BRs. Specifically, the gaps are 4.54 eV for (SQNN)⁻ and 5.47 eV for complex 1Y, and 3.11 eV for (SQPhNN)⁻ and 4.59 eV for complex 2Y, respectively. Thus, a possible reason for this is the stabilization of SOMO_SQ upon coordination. For complex Tp^CumZn(SQNN) (3), calculations using 10 triplet and 14 singlet states gave improved results: the 4^th triplet and 7^th singlet excited states are located at 24556 and 28249 cm^–¹, with CTC contributions of 51% and 81%, respectively. The experimental CT transition energies are reported at about 24500 and 23000 cm⁻¹ for Tp^CumZn(SQNN) [28] and Tp^CumZn(SQPhNN) [77], respectively, which is in excellent agreement with the calculated results.

Subsequently, the electronic coupling matrix elements H_AD was calculated according to the scheme shown in Figure 3, using the first excited singlet state energies, ¹GC, from Table 5. The obtained values are presented in Table 6. These values indicate that not only the structure of the biradicals (torsion angles, presence and nature of the bridge) but also coordination influences the coupling strength. Moreover, the calculated values are in good agreement with those determined from experimental data for Tp^CumZn(SQNN) 13460 cm⁻¹ [28] and 5740 cm⁻¹ for Tp^CumZn(SQPhNN) [77].

5. Conclusions

The electronic structure and magnetic exchange interactions of the hypothetical Gd³⁺ complexes Tp₂Gd(SQNN) (1) and Tp₂Gd(SQPhNN) (2) with DA biradicals were systematically investigated using high-level ab initio calculations and compared with the previously experimentally characterized monoradical complex Tp₂Gd(SQ), which features an identical coordination environment. The common organization of the coordination sphere across all three complexes allows for a direct assessment of how the replacement of a simple semiquinone radical by a strongly coupled donor–acceptor biradical affects the magnetic properties.

The comparison between SA-CASSCF and SA-CASSCF/NEVPT2 calculations demonstrated that the inclusion of dynamic correlation is essential for correctly reproducing the antiferromagnetic exchange between Gd³⁺ and the BR-ligand. While SA-CASSCF alone yielded a ferromagnetic sign and underestimated the coupling strength, the NEVPT2-corrected energies were consistent with the expected antiferromagnetic behavior. The exchange interaction between Gd³⁺ and the SQ moiety in complexes 1 and 2 (J_SQ-NN ≈ –2.6 cm⁻¹) was found to be approximately two times weaker than that in the monoradical complex Tp₂Gd(SQ) (J_Gd–SQ= –5.4 cm⁻¹, as confirmed by additional SA-CASSCF(10,10)/NEVPT2 calculations); suggesting that the presence of the strongly coupled –NN or –PhNN fragment partially quenches the metal–ligand exchange. This reduction is attributed to the redistribution of spin density within the biradical moiety upon coordination to the lanthanide center.

The intraligand exchange parameters J_SQ–NN = 679.1 cm⁻¹ and J_SQ–PhNN = 71.0 cm⁻¹ for complexes 1 and 2, respectively, reflect a strong ferromagnetic interaction within the biradical moieties, with the phenylene bridge significantly reducing the coupling strength. Comparison of the yttrium analogues 1Y and 2Y with the corresponding free biradicals revealed that coordination to the metal center similarly leads to a redistribution of spin density and alters the intraligand exchange parameters. For the phenylene-bridged system, a more pronounced reduction in J_Gd–PhSQ was observed, attributed to stronger spin-density delocalization from the SQ– fragment toward the –PhNN moiety upon coordination. Application of the VBCIM formalism to the computed charge-transfer excited states yielded electronic coupling matrix elements H_DA~ of 15283 cm⁻¹ for 1Y and 6546 cm⁻¹ for 2Y, which correlated with the computed intraligand exchange parameters and demonstrated that both the biradical structure (torsion angle, bridge presence) and the coordination environment are important factors governing the coupling strength.

The calculated χT versus T curves for complexes 1 and 2 showed a gradual decrease with decreasing temperature, closely resembling the magnetic behavior of Tp₂Gd(SQ). The reduced χT values at 300 K compared to the uncoupled spin-only expectations were attributed to the admixture of excited states into the ground spin state, a phenomenon previously observed for other Gd³⁺–radical systems. Notably, the characteristic temperatures of the magnetic behavior for complexes 1 and 2 differ significantly from those of Tp₂Gd(SQ), reflecting the unequal exchange coupling strengths of the Gd–SQ and Gd–SQ–NN exchange interactions.

This predictive study demonstrates that the replacement of a monoradical semiquinone ligand by a strongly coupled donor–acceptor biradical in the Tp₂Ln(SQ) coordination platform leads to a significant reduction in the Gd³⁺–ligand antiferromagnetic exchange, while preserving the overall magnetic behavior pattern. The results provide valuable guidelines for the rational design of lanthanide-based molecular magnetic materials and suggest that the Tp₂Ln coordination environment, combined with properly designed DA biradicals, offers a promising route for tuning magnetic properties through chemical modification of the ligand framework while maintaining a constant coordination geometry around the metal center. These findings are expected to facilitate the development of new coordination compounds with tailored magnetic characteristics for potential applications in molecular spintronics and quantum information science.