Designing molecular qubits: computational insights into first-row and group 6 transition metal complexes

Abstract

In the realm of optically addressable qubits, a previously synthesized and characterized Cr(IV) pseudo-tetrahedral complex, featuring four strongly donating ligands surrounding the chromium center, has demonstrated potential as a qubit candidate. This study proposes analogs of this complex through a metal substitution strategy, extending the investigation to different complexes based on metal centers selected from first-row and Group 6 transition metals. Computational modeling based on multiconfigurational methods CASPT2 and MC-PDFT was utilized to calculate energy gaps between ground and excited electronic spin states, and zero-field splitting parameters. Simulations were applied to each equilibrium geometry and related deformations based on vibrational modes. All results align with previous experimental findings, but also show that qubits based on V and Ti centers could be more electronically stable than the Cr one, suggesting a lower electronic features dependency from their related geometry. In some cases geometrical deformations provide changes in relative energy gaps between triplet and singlet excited state, that could potentially swap, offering a different initialization process, and some inspiration for ligand design based on such deformations. Additionally, this study identifies an unsynthesized Ti(II) compound as a promising candidate for molecular qubits. This finding highlights the role of computational multireference methods in the rational design of qubit systems.

---

1 Introduction

Quantum bits (qubits) are the fundamental building blocks of quantum information science and technology. Their development is expected to significantly enhance transformative applications, including secure communication,^1,2 complex financial modeling,^3,4 and efficient drug design.^5,6 Spin-based qubits have garnered significant attention due to their potential for long coherence times and compatibility with established experimental techniques. Among these, nitrogen-vacancy centers (NVs) in diamond have been extensively studied and are well-recognized for their applications in molecular sensing and quantum sensing technologies.^7–12 While NV centers remain the most advanced systems demonstrating these properties, transition metal complexes with metal centers like chromium and vanadium also show significant promise as alternatives. This is because of the difficulty in selectively addressing disordered defect centers. In contrast, molecular qubits offer structural reproducibility and modularity. Molecular systems with a triplet ground state and a singlet excited state have demonstrated the ability of mimicking the mechanisms of the NV centers for optical initialization and readout mechanisms through spin polarization.¹³ The optical transitions in these materials typically fall within the visible to near-infrared spectrum, corresponding to energy gaps suitable for efficient optical manipulation. Furthermore, they present unique advantages such as tunable physical and chemical properties enabled by chemical synthesis, as well as scalability in production.^14–16 Various platforms have been explored for molecular qubits, including organic radicals,¹⁷ rare-earth complexes,^18,19 transition metal complexes,^14,20–24 and hybrid systems featuring metal ions coordinated with organic radicals.^25,26

Designing optically addressable molecular qubits requires specific spin properties and controlled zero-field splitting (ZFS). A non-singlet spin state, often a triplet analogous to NV centers, can be achieved by incorporating an appropriate metal ion into symmetric environments such as octahedral or tetrahedral coordination. For instance, d² tetrahedral complexes^14,15 have demonstrated favorable magnetic anisotropy, crucial for stabilizing specific spin states. Additionally, higher spin multiplicities in organic systems, such as quartet states, can also be viable qubits through intersystem crossing with a chromophore.¹⁷

Another key requirement for coherent microwave control is a small ZFS parameter (|D|), enabling precise spin-state manipulation. For optically addressable molecular qubits, a |D| value below 20 GHz is preferred,^15,27 as it ensures that the spin polarization values fall within the X-band frequency range commonly utilized in electron paramagnetic resonance (EPR) spectroscopy. Triplet states have been widely studied due to their simplicity featuring only two unpaired electrons and three M_S levels (0 and ±1). Finally, an essential criterion is that the spin-lattice relaxation time (T₁) must exceed the optical emission lifetime (τ_opt) to maintain coherence.

In this study, we explore transition metal systems with a pseudo-tetrahedral (S₄) symmetry and a strong ligand field to achieve the desired electronic structure.^15,28–30 This structural choice optimizes microwave addressability by tuning the orbital configuration to yield a small |D| value, ensuring efficient spin control.^31–35

When starting from existing qubit structures, these features are addressed after structural modifications. These modifications involve metal or ligand substitution, and solvent or host matrix modification.¹⁵ ZFS changes upon ligand substitution were previously experimentally studied by Freedman and collaborators,³² and computationally by^36,37 using Cr(IV) as the metal center.³⁰

One challenge in performing electronic structure calculations on such systems is selecting methodologies that can accurately capture their complex electronic structure. Multireference methods and Kohn–Sham density functional theory (KS-DFT) are usually the methods of choice.^38–40 The zero-field splitting parameters can be computed using a coupled-perturbed response formalism within the framework of density functional theory.^41,42 In systems containing transition metals, excited states lie close in energy, and multireference methods are generally preferred.⁴³ While multiconfiguration perturbation theory methods like CASPT2 and NEVPT2 (ref. 44–46) provide improved accuracy over KS-DFT, their high computational cost limits their applicability to large molecular systems. Multiconfiguration pair-density functional theory^47–51 (MC-PDFT) provides results of PT2 quality at a significantly reduced cost because it does not require the higher order density matrices used in PT2-based methods.

In this study, we computationally investigate the effect of metal substitution, starting from the experimentally made pseudo-tetrahedral Cr(IV) complex.^32,36 We change the metal center while maintaining the same strong-field ligand environment formed by four o-tolyl ligands (Fig. 1). Seven complexes in addition to the Cr complex were considered; the metal centers were selected based on their similarity to Cr as a first-row transition metal and Group 6 elements. Specifically, we examined Ti, V, Fe, Co, Ni (first-row), and Mo and W (Group 6). The resulting complexes include titanium(II) (Ti(o-tol)₄²⁻), vanadium(III) (V(o-tol)₄⁻), chromium(IV) (Cr(o-tol)₄), iron(II) (Fe(o-tol)₄²⁻), cobalt(II) (Co(o-tol)₄²⁻), nickel(II) (Ni(o-tol)₄²⁻), molybdenum(IV) (Mo(o-tol)₄), and tungsten(IV) (W(o-tol)₄). This selection was intended to explore the chemical space for metals that might offer ideal properties for molecular qubits, such as a small triplet–singlet gap, within the visible to near-infrared spectrum, and without any other triplet state in between, and axial ZFS parameters |D| smaller than 20 GHz (EPR X-band limit).²⁷ Furthermore, given the potential toxicity of Cr, substituting it with less toxic metals such as Ti or Fe could present an advantage for biological or medical applications of molecular qubits. Multireference computational methods, namely CASPT2 and MC-PDFT, were used to assess the electronic and magnetic properties of these complexes, facilitating the identification of promising candidates for experimental validation. By systematically analyzing the impact of metal substitution on the electronic and magnetic properties of these compounds, we aim to provide valuable insights into the design principles of molecular qubits. Our major finding is that the Ti-based species investigated here for the first time has the requirements needed for a potential effective qubit (if it can be synthesized and stabilized). Furthermore, we demonstrate that MC-PDFT provides a quantitatively accurate description of these systems, making it a promising method for future studies due to its computational advantages over CASPT2.

Fig. 1 Molecular qubit structure with pseudo-tetrahedral ligand field. The M label indicates the considered metal centers, while the labels C₁, C₂, C₃, and C₄ the carbon atoms of each o-tolyl ligand bonded to the metal.

2 Results and discussion

The electronic configuration, formal charge of the metal ions, and spin multiplicity of the complexes are reported in Table 1.

Table 1 Electronic configuration (El. Conf.) and formal charge of the metal ion, and spin considered for each compound

Complex	El. Conf.	Charge	Spin
a Notice that the ground state of Ni(o-tol)4^2− is a singlet state, but for comparative purposes the triplet has been considered in this work.
Ti(o-tol)4^2−	3d^2	II	Triplet
V(o-tol)4^−	3d^2	III	Triplet
Cr(o-tol)4	3d^2	IV	Triplet
Fe(o-tol)4^2−	3d^6	II	Quintet
Co(o-tol)4^2−	3d^7	II	Quartet
Ni(o-tol)4^2−	3d^8	II	Triplet^a
Mo(o-tol)4	4d^2	IV	Triplet
W(o-tol)4	5d^2	IV	Triplet

---

The molecular geometries of all species at any accessible spin state were optimized using density functional theory (DFT) with four functionals: PBE-D3BJ, B3LYP-D3BJ, TPSSh-D3BJ, and M06 (see Computational details). Upon optimization, these complexes adopt a pseudo-tetrahedral structure in all accessible spin states (beside the Ni singlet state, see below), which is an optimal geometry for optically addressable molecular qubits.^30,32 The adiabatic energy gaps between spin states were evaluated with all the listed functionals. The geometric deviations among the selected functionals were minor, therefore the TPSSh-D3BJ optimized structures were selected for further multireference calculations due to its accuracy with most transition metal complexes.^52–54 Vertical triplet–singlet gaps were calculated with ab initio multireference methods (MR) state-specific complete active space followed by second-order perturbation theory (CASPT2),^55,56 and multiconfigurational pair-density functional theory (MC-PDFT)^47,57 with the tPBE and tPBE0 functionals.

In the following, such complexes are discussed in two separate sections: first the complexes with electronic configuration d², which have a triplet ground state (Cr(o-tol)₄, V(o-tol)₄⁻, Mo(o-tol)₄, W(o-tol)₄, Ti(o-tol)₄²⁻); and subsequently those with a different ground-state spin multiplicity (Fe(o-tol)₄²⁻, Co(o-tol)₄²⁻, Ni(o-tol)₄²⁻).

For the MR calculations, we explored active spaces of increasing size (see Computational details). In the following, we present results with the largest active space for each system, while we discuss the effect of increasing the active space in the Active Space Dependency section in the ESI.†

2.1 Triplet ground state: Ti, Cr, V, Mo, and W systems

In this section, we discuss the Cr(o-tol)₄, Mo(o-tol)₄, W(o-tol)₄, Ti(o-tol)₄²⁻, and V(o-tol)₄⁻ complexes, which exhibit a d² electronic configuration and a triplet ground state. The Cr, Mo, and V complexes have been previously investigated.³¹ The triplet–singlet energy gap of Cr(o-tol)₄ has been experimentally determined as 1.209 eV (1025 nm).³⁰ For Mo(o-tol)₄ and V(o-tol)₄⁻, the triplet–singlet gaps have been computationally estimated using periodic boundary condition density functional theory (PBC-DFT) with the PBE functional, yielding values of 0.69 eV (1800 nm) and 0.80 eV (1550 nm), respectively.^31,58 In Table 2 we report the triplet–singlet gaps for these species.

Table 2 Triplet–singlet gaps (ΔE_TS) for the Cr(o-tol)₄, V(o-tol)₄⁻, Mo(o-tol)₄, W(o-tol)₄, and Ti(o-tol)₄²⁻ complexes. The DFT values correspond to adiabatic energy differences. The multireference (MR) values correspond to vertical energy differences, computed at the TPSSh triplet optimized geometries. All values are in eV. All the species have a d² electronic configuration and a triplet ground state

Complex^a	MR vertical ΔETS			DFT^b adiabatic ΔETS				Ref.
	CASPT2	tPBE	tPBE0	PBE	B3LYP	TPSSh	M06
a For Cr(o-tol)4, Mo(o-tol)4, W(o-tol)4, and Ti(o-tol)4^2−, the data corresponds to the (10,15) active space, and for V(o-tol)4^−, the (8,13) active space.b All DFT functionals except M06 include D3BJ empirical correction. The RMSD values between optimized geometries is lower than 0.4 Å, thus DFT adiabatic triplet–singlet differences can be compared with MR vertical differences.c Experimental value from ref. 30.d These values are estimations based on periodic PBE diabatic calculations from ref. 31 corrected using the experimental values from ref. 30.
Cr(o-tol)4	1.44	1.04	1.20	0.51	1.52	1.55	1.67	1.209^c
V(o-tol)4^−	1.18	0.65	0.85	0.92	1.42	1.36	1.37	0.80^d
Mo(o-tol)4	0.82	0.54	0.67	0.87	1.03	1.07	1.01	0.69^d
W(o-tol)4	0.65	0.43	0.55	0.74	0.89	0.91	0.87	—
Ti(o-tol)4^2−	0.56	0.49	0.55	0.20	0.43	0.41	0.44	—

---

The multireference values are vertical differences computed at the ground-state triplet geometry, while the DFT ones are adiabatic. The geometrical differences between the optimized singlet and the triplet structures are minimal, (root mean square deviation values are lower than 0.4 Å see ESI†) therefore comparing DFT adiabatic values with vertical MR values is justified. Geometric parameters values are reported in the ESI in Tables S1 and S2.†

Our calculations predict a triplet ground state for all species. Table 2 shows that, among the multireference methods, the tPBE0 relative energies closely align with reference values, namely the experimental value for the Cr(o-tol)₄ complex, and computational values for Mo(o-tol)₄, and V(o-tol)₄⁻,³¹ the values obtained are 1.20, 0.85, and 0.67 eV (27.7, 19.6, 15.5 kcal mol⁻¹) respectively. In contrast, CASPT2 and tPBE exhibit a slightly larger discrepancy of approximately 0.2 eV in comparison to tPBE0 and the reference values (Table 2). The triplet–singlet energy gaps calculated, for the same list of complexes, are 1.44, 1.18, 0.82 eV (33.2, 27.2, 18.9 kcal mol⁻¹) with CASPT2, and 1.04, 0.65, 0.54 eV with tPBE. We also considered the W(o-tol)₄ and Ti(o-tol)₄²⁻ complexes, for which there are no experimental reference values; the triplet–singlet energy gap is 0.55 eV (12.7 kcal mol⁻¹) for both complexes with tPBE0, 0.65, and 0.56 eV (15.0, 12.9 kcal mol⁻¹) respectively with CASPT2, and 0.43, and 0.49 eV respectively with tPBE (Table 2).

Multireference calculations can offer insights that can help rationalize the experimental results on the V(o-tol)₄⁻ complex, which exhibits no measurable emission in the 900–1700 nm range.³¹ It has been suggested³¹ that this lack of emission may be due to a low-lying S₁ excited state (>1200 nm, ∼1 eV) whose relaxation overlaps spectrally with high-energy C–H stretching overtones. This overlap facilitates a non-radiative decay via multi-phonon mediated relaxation pathways. The CASPT2 and tPBE0 singlet–triplet energy gaps of 1.18 eV and 0.85 eV respectively, fall within this energy range, supporting the proposed explanation.

All the DFT functionals predict a triplet–singlet gap larger than the previously reported values³¹ by an average of 0.3 eV. All values are reported in Table 2 and follow the same trend as the MR methods, with the exception of PBE-D3BJ. Indeed, PBE for the Cr(o-tol)₄ complex provides an energy gap of 0.51 eV, which is 0.7 eV lower than the experimental reference. The values obtained for the list of complexes Cr(o-tol)₄, V(o-tol)₄⁻, Mo(o-tol)₄, W(o-tol)₄, and Ti(o-tol)₄²⁻ are, with PBE-D3BJ 0.51, 0.92, 0.87, 0.74, 0.20 eV, while with TPSSh-D3BJ 1.55, 1.36, 1.07 0.91, and 0.41 eV (35.7, 31.4, 24.7, 21.0, 9.5 kcal mol⁻¹) respectively (B3LYP-D3BJ, and M06 values are very similar to the ones obtained with TPSSh-D3BJ therefore they are only reported in Table 2). In conclusion, CASPT2, tPBE0, and TPSSh-D3BJ agree better with the experimental value in the Cr(o-tol)₄ case, these methods will be used for the rest of the study. The triplet–singlet gaps follow the trend Cr(o-tol)₄ > V(o-tol)₄⁻ > Mo(o-tol)₄ > W(o-tol)₄ > Ti(o-tol)₄²⁻ with all these three methods (Table 2, and Fig. S2†).

2.2 Non-triplet ground state: Fe, Co, and Ni systems

The Fe(o-tol)₄²⁻ and Co(o-tol)₄²⁻ complexes were optimized in their quintet and quartet ground states, respectively, resulting in pseudo-tetrahedral structures. Instead, Ni(o-tol)₄²⁻ has a singlet ground state and a planar structure, however its triplet-optimized geometry is pseudo-tetrahedral; thus for the sake of comparisons with the other complexes, this non-ground state geometry is used for further calculations (even though it is not the lowest energy state for the current ligand environment). The complexes discussed in this section do not have a triplet ground state, and moreover a direct spin-flip first optical transition is missing due to the presence of other states between the ground state and the first excited state (see ESI†). Therefore they are not ideal qubit candidates. However, a brief description is reported next.

For the Ni(o-tol)₄²⁻ species, at the triplet optimized geometry, (Fig. S2†) the CASPT2 and tPBE0 triplet–singlet gaps are 1.10 eV and 0.91 eV (25.4, 21.0 kcal mol⁻¹) respectively. Different ligand environments for Ni(II) could eventually be explored to have a triplet ground state, for qubit design purposes.

The Fe(o-tol)₄²⁻ complex has a quasi doubly-degenerate quintet ground state, and the first excited state is also a quintet state, positioned at about 0.7 eV (16.1 kcal mol⁻¹), according to all multireference methods (Fig. S9 and S17†). The lowest triplet excited state lies 1.53 eV, and 1.02 eV higher in energy using the CASPT2 and tPBE0 methods (35.5, 23.5 kcal mol⁻¹), respectively. Moreover, this triplet is higher than several quintet states; therefore the qubit initialization cannot take place since there are no close-lying states of other spin multiplicities (see ESI†).

For the Co(o-tol)₄²⁻ complex, similarly to the Fe(o-tol)₄²⁻ complex, both the ground state and the first excited state are quartet states (Fig. S10 and S18†). The first doublet state is located after several quartet excited states and similarly the qubit initialization cannot take place (see ESI† for full details).

To summarize, due to their electronic structures, the Fe(o-tol)₄²⁻, Co(o-tol)₄²⁻, and Ni(o-tol)₄²⁻ complexes are not relevant as potential molecular qubits. In the following we will consider only the Cr(o-tol)₄, V(o-tol)₄⁻, Mo(o-tol)₄, Ti(o-tol)₄²⁻, and W(o-tol)₄ complexes for the zero-field splitting discussion.

2.3 Zero-field splitting parameters of the Cr, V, Mo, and Ti complexes

Optically addressable molecular qubits are considered for practical use when the axial ZFS parameter |D| falls into the X-band range, thus below 20 GHz (0.67 cm⁻¹).^15,27 The parameter |D| reflects the splitting of the triplet ground state M_S levels, which directly impacts the qubit's initialization and control, and the |E| parameter indicates whether a molecular system is uniaxial⁵⁹ (see ESI† for further details). The ZFS depends on relativistic effects and electronic state interactions, both of which are incorporated into our calculations (see Computational details section).

The experimentally reported |D| values for the Cr(o-tol)₄, V(o-tol)₄⁻, and Mo(o-tol)₄ complexes are 3.53, 5.62, and 7.30 GHz, respectively.^30,31 We computed ZFS parameters for these complexes as well as for the two unsynthesized species, Ti(o-tol)₄²⁻ and W(o-tol)₄, using multireference methods (see Computational details). The computed CASPT2 and tPBE0 |D| values are presented in Fig. 2.

Fig. 2 Calculated axial ZFS parameter |D| using CASPT2 (solid bars) and tPBE0 (striped bars) for the complexes: Cr(o-tol)₄ red, V(o-tol)₄⁻ orange, Mo(o-tol)₄ teal, Ti(o-tol)₄²⁻ gray, and W(o-tol)₄ cyan. Active spaces of (8,13) for V(o-tol)₄⁻ and (10,15) for Cr(o-tol)₄, Mo(o-tol)₄, Ti(o-tol)₄²⁻, and W(o-tol)₄ were used. Values are in GHz. Dashed lines indicate experimental data from ref. 30 and 31.

The CASPT2 |D| values for the Cr(o-tol)₄, V(o-tol)₄⁻, Mo(o-tol)₄, Ti(o-tol)₄²⁻, and W(o-tol)₄ complexes are 2.72, 6.87, 6.82, 15.11, and 145.50 GHz, respectively and the corresponding tPBE0 values are 2.39, 2.95, 4.09, 13.82, and 124.04 GHz. The W(o-tol)₄ complex exceeds the optimal threshold for |D|, disqualifying it as a viable molecular qubit candidate. The tPBE0 method aligns closely with CASPT2 for Cr(o-tol)₄, which exhibits the smallest |D| among the studied complexes. For the V(o-tol)₄⁻ and Mo(o-tol)₄ complexes, CASPT2 provides better agreement with experimental data than tPBE0, but both methods reproduce the experimental trend. Finally, we also calculated the rhombic parameter |E| which turned out to be 0 GHz in all cases (see ESI in Table S10†).

Although we reported the absolute |D| values thus far, the sign of D is also critical, but often challenging to determine experimentally. Computationally, we can resolve the sign: the computed axial parameters for Cr(o-tol)₄, V(o-tol)₄⁻, Mo(o-tol)₄, and Ti(o-tol)₄²⁻ are all negative. This agrees with recent literature,⁶⁰ which reports D = −3.53 GHz for the Cr(o-tol)₄ complex.

The observed increase in |D| values within the Group 6 complexes (Cr(o-tol)₄, Mo(o-tol)₄, W(o-tol)₄) can be rationalized using eqn (S5) in the ESI,† which predicts a trend of increasing |D| with the spin–orbit coupling constant (ζ) and the atomic number (Z). Specifically, |D| scales proportionally to ζ² and Z⁴. The spin–orbit coupling constants for Cr(o-tol)₄, Mo(o-tol)₄, and W(o-tol)₄ are 325, 950, and 2300 cm⁻¹, respectively,⁶¹ illustrating the substantial increase with heavier elements.

This relationship is visualized in Fig. S3,† where panel (a) shows the |D| parameters in cm⁻¹ for the Cr(o-tol)₄, Mo(o-tol)₄, and W(o-tol)₄ complexes, and panel (b) displays their corresponding squared spin–orbit coupling constants. These results confirm that relativistic effects significantly influence the |D| values in heavier metal centers.

The |D| values increase from Cr(o-tol)₄ to V(o-tol)₄⁻ and Mo(o-tol)₄ to Ti(o-tol)₄²⁻ (Fig. 2). This trend is explained considering the relative energies shown in Fig. S4–S8.† According to eqn (S5) in the ESI,† when the energy difference between two states is small, the contribution to D is large.

The axial parameters for the Cr(o-tol)₄, V(o-tol)₄⁻, Mo(o-tol)₄, and Ti(o-tol)₄²⁻ complexes are all below the 20 GHz threshold, supporting their potential as molecular qubits. Although this is consistent with previous findings for Cr(o-tol)₄, V(o-tol)₄⁻, and Mo(o-tol)₄, the identification of Ti(o-tol)₄²⁻ as a promising candidate is noteworthy. To the best of our knowledge, titanium has not previously been reported as a potential metal center for molecular qubit applications. However, Ti in the oxidation state II occurs in reactive species and catalytic processes. Notably, there are experimentally stable Ti(II) complexes,⁶² which could serve as valuable inspiration for future qubit materials. Building on these insights, our future work will focus on developing and exploring molecular qubit candidates centered on Ti(II).

2.4 Effect of geometric distortion on the zero-field splitting parameters

The computations discussed so far were performed at equilibrium geometries. The systems, however, are in a host matrix and they are subject to distortions. In previous studies^63,64 electrostatic embedding approaches have been employed to incorporate the electrostatic influence of the host matrix in multireference calculations. In contrast, our approach focuses on deforming the molecules by following two vibrational modes, the symmetric stretch and scissoring (Fig. S23†), from which we generated eight distorted structures. We applied these distortions to the Cr(o-tol)₄, V(o-tol)₄⁻, and Ti(o-tol)₄²⁻ complexes and analyzed how they affect the energy gaps between the ground state (T₀) and excited states (S₁ or T₁), as well as the zero-field splitting parameter |D|. In this section, we discuss only the CASPT2 results, which have a similar trend to the tPBE0 ones.

For the Cr(o-tol)₄ complex, the energy gaps are only slightly affected by stretching/compressing the Cr–C bonds by 0.1 Å and the ZFS parameter |D| varies only by approximately 1 GHz (Fig. 3 and Table S77†). However, bond angle variations of ±7° along the scissoring vibration have a more pronounced effect on |D| up to about 12 GHz (Fig. 3 and Table S78†). This highlights the higher sensitivity to angular deformations than bond distortion.

Fig. 3 CASPT2 and tPBE0 axial parameter |D| of Cr(o-tol)₄ complex for (a) the Cr–C symmetric stretching, and (b) the C–Cr–C scissoring normal modes. The equilibrium bond angle is θ = 104.9°, and the equilibrium bond distance is d = 1.98 Å.

For the V(o-tol)₄⁻ complex, variations in the V–C bond length have little effect on the energy of the S₁ state. The T₁ state instead lowers in energy with bond elongation and eventually crosses below S₁, becoming the lowest-energy excited state (Fig. S26†). The T₁–S₁ inversion could influence the initialization process in the distorted geometry, potentially leading to fluorescence instead of non-radiative decay. As a result, the system may fail to achieve the required spin polarization, compromising initialization. Such a reversal in singlet–triplet stability warrants experimental validation, perhaps in a ligand environment that allows for a slightly larger V—ligand distance. In contrast, bond angle distortions have little effect on the energy gaps between the excited states and the ground state (Fig. S27†). The ZFS parameter |D| varies by approximately ±2 GHz during bond stretching and compression, while bond-angle variations affect |D| in a trend similar to the Cr(o-tol)₄ complex: slight deformations up to ±3° reduce |D| from 6 to 2 GHz, whereas larger distortions, up to ±6° increase |D| to 12 GHz (Fig. S38†).

For the Ti(o-tol)₄²⁻ complex, the change in relative energies of the T₁ and S₁ states is reversed to the V(o-tol)₄⁻ complex: variations in the Ti–C bond length have effect on the S₁ state which increases in energy with bond elongation and eventually crosses above T₁ which is constant in energy despite elongation (Fig. S28†). Variations in bond angles along the scissoring vibration do not have effects on the relative energies of the excited states (Fig. S29†). For the ZFS parameter |D|, when the Ti–C bond is compressed |D| remains unchanged, while upon stretching by 0.08 Å, it decreases from its equilibrium value of 14 GHz to 7 GHz. Bond angle modifications affect |D| by at most 4 GHz (Fig. S40 and S41†).

In conclusion, the relative energies of the lower electronic states of the Cr(o-tol)₄ complex remain largely unaffected by bond length variations, while bond angle variations significantly influence the ZFS parameter |D|, which increases with larger bond angles, consistent with experimental observations.³¹ V(o-tol)₄⁻ shows more pronounced effects on the relative energies of T₁ and S₁ compared to the Cr(o-tol)₄ complex, with a potential swap between the two states at larger bond elongation. Finally, the Ti(o-tol)₄²⁻ complex exhibits changes in both excited-state energetics and |D| with bond length variations, whereas bond angle modifications have little to no impact. Overall, since molecular deformations in a host matrix are primarily associated with bond angle distortions, our findings suggest that the Ti(o-tol)₄²⁻ complex is less prone to these deformations, and thus small |D| variations than the Cr(o-tol)₄, and V(o-tol)₄⁻ complexes.

3 Conclusion

We used multiconfigurational CASPT2 and tPBE0 methods to investigate the properties of a series of organometallic complexes as potential candidates for molecular qubits. The features under study were the energy gaps between states with different spin and zero-field splitting parameters |D| and |E|. The metals considered were first-row and Group 6 transition metals. Our results confirm prior experimental findings for Cr, V, and Mo complexes while extending the scope to new candidates. Notably, the Ti(II) complex is identified as a promising, novel candidate for molecular qubits.

Through a systematic evaluation of vertical excitation energies and zero-field splitting parameters, we identified key trends across the transition-metal series. The vertical excitation energy of the lowest excited state (of different spin than the ground state) follows the sequence Ti(o-tol)₄²⁻ < W(o-tol)₄ < Mo(o-tol)₄ < V(o-tol)₄⁻ < Cr(o-tol)₄ reflecting the progressive filling of d orbitals and corresponding changes in electronic structure. In contrast, the electronic structures of Fe(o-tol)₄²⁻, Co(o-tol)₄²⁻, and Ni(o-tol)₄²⁻ suggest that they are unsuitable as molecular qubits, due to high-energy and unfavorable different-spin excitations. To explore these transition metals as viable candidates, alternative ligand environments would need to be considered.

Zero-field splitting calculations provided further refinement of qubit suitability. For Cr(o-tol)₄, V(o-tol)₄⁻, and Mo(o-tol)₄ complexes, both CASPT2 and tPBE0 reproduce the experimental trend for |D| values. The results also highlight that the W(o-tol)₄ complex does not meet the optimal threshold for |D|, because of its high atomic number, thereby eliminating it from consideration. Importantly, the Ti(o-tol)₄²⁻ complex has a predicted |D| value below 20 GHz, supporting its candidacy alongside Cr(o-tol)₄, V(o-tol)₄⁻, and Mo(o-tol)₄.

We investigated the effect of geometric distortions induced by the host matrix on the ZFS parameter |D| by varying the M–C bond length and C–M–C bond angle. Our analysis shows that the Ti(o-tol)₄²⁻ complex is less sensitive to geometric distortions compared to Cr(o-tol)₄ and V(o-tol)₄⁻, whose |D| parameter varies significantly with bond angle. These results highlight the importance of incorporating matrix effects when modeling molecular complexes as potential qubits.

The use of Ti(II) as a metal center for molecular qubit candidates is a promising avenue of future investigation. Stable Ti(II) complexes documented in the literature provide a basis for further exploration.^62,65 Future efforts from our side will focus on computational modeling and synthesizing Ti-based systems, as well as extending computational studies to explore other ligand environments and their influence on qubit properties.

4 Computational methods

4.1 Geometry optimization

Geometry optimization calculations were performed using Kohn–Sham density functional theory (KS-DFT) with Gaussian 16 rev. A03.⁶⁶ All calculations were carried out for isolated molecules, employing the following functionals: meta-GGA hybrid TPSSh,^67,68 hybrid B3LYP,^69–71 meta-GGA M06,⁷² and GGA PBE,⁷³ all including the empirical dispersion correction D3BJ based on Becke and Johnson damping, except for M06.⁷⁴ The Ahlrichs basis set def2-TZVP basis set⁷⁵ was used for all atoms, with pseudopotential def2-ECP for Mo, and W. Vibrational calculations were performed at the TPSSh-D3BJ/def2-TZVP level of theory.

4.2 SA-CASSCF calculations

The ground-state optimized structures served as inputs for the state-average complete active space self-consistent field (SA-CASSCF)⁷⁶ calculations, which were employed to calculate the electronic excitation energies. These calculations were carried out using OpenMolcas v23.10.⁷⁷

The ANO-RCC-VTZP basis set was employed for the metal centers, while ANO-RCC-VDZ basis sets were assigned to the carbon and hydrogen atoms.⁷⁸ Scalar relativistic effects were incorporated using the second-order Douglas–Kroll–Hess Hamiltonian.^79–83 To reduce computational costs, the resolution of the identity and Cholesky decomposition techniques were utilized.⁸⁴

For the Ti(o-tol)₄²⁻, Cr(o-tol)₄, Mo(o-tol)₄, and W(o-tol)₄ complexes, SA-CASSCF wave functions were computed using the active spaces shown in Fig. 4. The Ti, Cr, Mo and W metal centers have electronic configuration nd², where n = 3, 3, 4, 5 respectively. The final active spaces were obtained by progressively adding orbitals starting from the smallest active space, consisting of 2 electrons in the corresponding d-shell, (2,5), formed by the nd_z² and nd_x²−y² orbitals (we refer to them as e pair implying a tetrahedral symmetry, even if these species are only pseudo tetrahedral) and the nd_xy, nd_xz, and nd_yz orbitals (collectively referred to as ) shown in Fig. 4. A correlating set of d orbitals ((n + 1)d) was included to expand the active space, resulting in (2,10). By incorporating three metal–ligand molecular orbitals (t₂) into the active space starting from (2,5), the configuration was extended to (8,8). Finally, the addition of the correlating d-shell along with the metal s bonding (σ) and antibonding (σ*) orbitals expanded the active space to (10,15).

Fig. 4 Active spaces used for the multiconfigurational calculations of molecules studied herein. On the left energy scale, it is shown the active spaces used for molecules Ti(o-tol)₄²⁻, Cr(o-tol)₄, Mo(o-tol)₄, and W(o-tol)₄. The right energy scale shows only the largest active space used for complexes V(o-tol)₄⁻, Fe(o-tol)₄²⁻, Co(o-tol)₄²⁻, and Ni(o-tol)₄²⁻.

The relative energies with increasing active spaces are reported in Tables S44–S70.† There is a minor active space dependence in the results, a reason for which we presented the results with the largest active space, which is more balanced. In general the composition of the wave function does not vary significantly. For example, in the systems with a triplet ground state, the triplet is dominated by the single configuration (σ)²(t₂)⁶(nd_z²)^↑(nd_x²−y²)^↑((n + 1)d)⁰(σ*)⁰ while the lowest singlet has two dominant configurations, (i) (σ)²(t₂)⁶(nd_z²)⁰(nd_x²−y²)²((n + 1)d)⁰(σ*)⁰ and (ii) (σ)²(t₂)⁶(nd_z²)²(nd_x²−y²)⁰((n + 1)d)⁰(σ*)⁰.

Among the first-row transition metals ions, the electronic configurations are 3d^m, where m = 2 for V³⁺, m = 6 for Fe²⁺, m = 7 for Co²⁺, and m = 8 for Ni²⁺. Using SA-CASSCF, the (m, 5) and (m, 10) active spaces were obtained. On the other hand, the inclusion of the σ and σ* orbitals in the active spaces of the V(o-tol)₄⁻ complex faced some complications related to the inclusion of such orbitals in the active space. Also, due to the unsuitability of the Fe(o-tol)₄²⁻, Co(o-tol)₄²⁻, and Ni(o-tol)₄²⁻ complexes as molecular qubits a larger active space was not explored. Therefore, the SA-CASSCF calculations considered 13 orbitals rather than 15, as depicted in the right part of Fig. 4. The active spaces considered consist of the 3d orbitals (m, 5), the inclusion of a second d–shell (m,10), the alternative active space where we add three metal-bonding molecular orbitals to the 3d orbitals (m + 6,8), and the inclusion of the second d–shell to the latter active space to obtain (m + 6,13).

For the Cr(o-tol)₄, Mo(o-tol)₄, and W(o-tol)₄ complexes, 7 triplets and 9 singlets were included in the state-averaged calculations for the (10,15) active spaces. For the first-row transition metal complexes, state-averaged calculations of Ti(o-tol)₄²⁻ and Ni(o-tol)₄²⁻ were performed using 10 triplets and 15 singlets for the corresponding active space. For V(o-tol)₄⁻ complex, 7 triplets and 9 singlets were calculated using the (8,13) active space. In the case of the Fe(o-tol)₄²⁻ complex, 5 quintets, 15 triplets, and 5 singlets were considered for the (12,13) active space. For the Co(o-tol)₄²⁻ complex, state-averaged calculations included 10 quartets and 30 doublets. The full discussion about these choices is in the ESI in the Zero-field splitting parameters section.†

4.3 CASPT2 and MC-PDFT calculations of energy gaps

The triplet–singlet energy gaps for the molecules Cr(o-tol)₄, Mo(o-tol)₄, W(o-tol)₄, Ti(o-tol)₄²⁻, V(o-tol)₄⁻, and Ni(o-tol)₄²⁻ were computed as vertical excitations from the triplet ground state geometry. On the other hand, the Fe(o-tol)₄²⁻ and Co(o-tol)₄²⁻ complexes do not have a triplet ground state but instead a quintet and quartet respectively. Therefore in these latter cases the quartet-doublet energy gap was calculated for the Co(o-tol)₄²⁻ complex, and the quintet–triplet gap for the Fe(o-tol)₄²⁻ complex, in both cases as vertical excitation using the ground state geometry. The reported energy gaps represent spin–orbit-free energies obtained through SA-CASSCF and post-CASSCF methods. Additionally, single-state complete active space second-order perturbation theory (CASPT2),⁵⁵ multiconfiguration pair-density functional theory (MC-PDFT) with the tPBE functional,⁴⁷ and hybrid MC-PDFT (HMC-PDFT) with the tPBE0 functional⁵⁷ were employed (see ESI†). The tPBE0 functional incorporates a mix of 25% CASSCF energy with 75% tPBE energy.⁵⁷

For the CASPT2 calculations, an ionization-potential electron affinity (IPEA) value of 0.25 a.u. and an imaginary shift of 0.30 a.u. were applied.

4.4 Zero-field splitting calculations

The zero-field splitting (ZFS) parameters were computed by incorporating spin–orbit coupling into the SA-CASSCF, CASPT2, tPBE and tPBE0 methods, utilizing the restricted active space state interaction (RASSI) method as implemented in OpenMolcas. The ZFS parameters were calculated for all the aforementioned active spaces using an effective Hamiltonian in a pseudospin basis approximation,⁸⁵ as implemented in the SINGLE_ANISO module of OpenMolcas. The procedure for computing the ZFS parameters is extensively described in ref. 36, 86 and 87.

Data availability

The data supporting this article have been included as part of the ESI.† OpenMolcas input files, and Cartesian coordinates of all structures are available in the all_structures.zip file.

Author contributions

A. S. and A. D. contributed equally to this work.

Conflicts of interest

Authors declare no conflict of interest.

Acknowledgements

This material is based upon work supported by the U.S. Department of Energy, Office of Science, National Quantum Information Science Research Centers. The authors thank the Research Computing Center (RCC) from the University of Chicago for computational resources. The authors especially appreciate the insightful discussions with Prof. John S. Anderson, and Dr María J. Beltrán-Leiva; and dedicate this work to the memory of José Luis Sauza Vega.

Notes and references

1. M. A. Nielsen and I. L. Chuang, Quantum Computation and Quantum Information: 10th Anniversary Edition, Cambridge University Press, 2010.
2. J. Biamonte, P. Wittek, N. Pancotti, P. Rebentrost, N. Wiebe and S. Lloyd, Nature, 2017, 549, 195–202.
3. A. Harrow and A. Montanaro, Nature, 2017, 549, 203–209.
4. P. Rebentrost, B. Gupt and T. R. Bromley, Phys. Rev. A, 2018, 98, 022321.
5. M. S. Fataftah and D. E. Freedman, Chem. Commun., 2018, 54, 13773–13781.
6. R. Santagati, A. Aspuru-Guzik, R. Babbush, M. Degroote, L. González, E. Kyoseva, N. Moll, M. Oppel, R. M. Parrish, N. C. Rubin, M. Streif, C. S. Tautermann, H. Weiss, N. Wiebe and C. Utschig-Utschig, Nat. Phys., 2024, 20, 549–557.
7. A. Gruber, A. Driäbenstedt, C. Tietz, L. Fleury, J. Wrachtrup and C. von Borczyskowski, Science, 1997, 276, 2012–2014.
8. N. Zhao, S. W. Ho and R. B. Liu, Phys. Rev. B: Condens. Matter Mater. Phys., 2012, 85, 115303.
9. J. Achard, J. Beck, J. Tissler, V. Jacques, P. R. Hemmer, F. Jelezko and J. Wrachtrup, Nat. Mater., 2009, 8, 383–387.
10. L. T. Hall, J. H. Cole and L. C. L. Hollenberg, Phys. Rev. B: Condens. Matter Mater. Phys., 2014, 90, 075201.
11. T. Zhang, G. Pramanik, K. Zhang, M. Gulka, L. Wang, J. Jing, F. Xu, Z. Li, Q. Wei, P. Cígler and Z. Chu, ACS Sens., 2021, 6, 2077–2107.
12. E. Janitz, K. Herb, L. Völker, W. S. Huxter, C. Degen and J. Abendroth, J. Mater. Chem. C, 2022, 10, 13533–13569.
13. V. Dobrovitski, G. Fuchs, A. Falk, C. Santori and D. Awschalom, Annu. Rev. Condens. Matter Phys., 2013, 4, 23–50.
14. J. M. Zadrozny, J. Niklas, O. G. Poluektov and D. E. Freedman, ACS Cent. Sci., 2015, 1, 488–492.
15. M. S. Fataftah and D. E. Freedman, Chem. Commun., 2018, 54, 13773–13781.
16. T. Mani, Chem. Phys. Rev., 2022, 3, 021301.
17. M. Mayländer, S. Chen, E. R. Lorenzo, M. R. Wasielewski and S. Richert, J. Am. Chem. Soc., 2021, 143, 7050–7058.
18. R. Hussain, G. Allodi, A. Chiesa, E. Garlatti, D. Mitcov, A. Konstantatos, K. S. Pedersen, R. De Renzi, S. Piligkos and S. Carretta, J. Am. Chem. Soc., 2018, 140, 9814–9818.
19. M. Atzori and R. Sessoli, J. Am. Chem. Soc., 2019, 141, 11339–11352.
20. M. F. Islam, J. F. Nossa, C. M. Canali and M. Pederson, Phys. Rev. B: Condens. Matter Mater. Phys., 2010, 82, 155446.
21. A. I. Johnson, F. Islam, C. M. Canali and M. R. Pederson, J. Chem. Phys., 2019, 151, 174105.
22. K. Bader, D. Dengler, S. Lenz, B. Endeward, S.-D. Jiang, P. Neugebauer and J. van Slageren, Nat. Commun., 2014, 5, 5304.
23. M. Atzori, L. Tesi, E. Morra, M. Chiesa, L. Sorace and R. Sessoli, J. Am. Chem. Soc., 2016, 138, 2154–2157.
24. M. J. Graham, J. M. Zadrozny, M. Shiddiq, J. S. Anderson, M. S. Fataftah, S. Hill and D. E. Freedman, J. Am. Chem. Soc., 2014, 136, 7623–7626.
25. P. Roseiro, S. Yalouz, D. J. R. Brook, N. Ben Amor and V. Robert, Inorg. Chem., 2023, 62, 5737–5743.
26. P. Roseiro, L. Petit, V. Robert and S. Yalouz, ChemPhysChem, 2023, 24, e202200478.
27. Coherent Control of NV Centers, https://www.zhinst.com/ch/en/applications/quantum-technologies/coherent-control-of-nv-centers.
28. J. Krzystek, J.-H. Park, M. W. Meisel, M. A. Hitchman, H. Stratemeier, L.-C. Brunel and J. Telser, Inorg. Chem., 2002, 41, 4478–4487.
29. P. J. Alonso, J. Forniés, M. A. García-Monforte, A. Martín and B. Menjón, Chem. Commun., 2001, 2138–2139.
30. S. L. Bayliss, D. W. Laorenza, P. J. Mintum, B. D. Kovos, D. E. Freedman and D. D. Awschalom, Science, 2020, 370, 1309–1312.
31. D. W. Laorenza, K. R. Mullin, L. R. Weiss, S. L. Bayliss, P. Deb, D. D. Awschalom, J. M. Rondinelli and D. E. Freedman, Chem. Sci., 2024, 15, 14016–14026.
32. D. W. Laorenza, A. Kairalapova, S. L. Bayliss, T. Goldzak, S. M. Greene, L. R. Weiss, P. Deb, P. J. Mintun, K. A. Collins, D. D. Awschalom, T. C. Berkelbach and D. E. Freedman, J. Am. Chem. Soc., 2021, 143, 21350–21363.
33. W. T. Morrillo, H. I. J. Cumming, A. Mattioni, J. K. Staab and N. F. Chilton, J. Am. Chem. Soc., 2024, 146, 25841–25851.
34. S. L. Bayliss, P. Deb, D. W. Laorenza, M. Onizhuk, G. Galli, D. E. Freedman and D. D. Awschalom, Phys. Rev. X, 2022, 12, 031028.
35. M. J. Amdur, K. R. Mullin, M. J. Waters, D. Puggioni, M. K. Wojnar, M. Gu, L. Sun, P. H. Oyala, J. M. Rondinelli and D. E. Freedman, Chem. Sci., 2022, 13, 7034–7045.
36. A. Sauza-de la Vega, R. Pandharkar, G. Stroscio, A. Sarkar, D. Truhlar and L. Gagliardi, JACS Au, 2022, 2, 2029–2037.
37. K. Janicka, A. L. Wysocki and K. Park, J. Phys. Chem. A, 2022, 126, 8007–8020.
38. F. Neese, Coord. Chem. Rev., 2009, 253, 526–563.
39. M. Atanasov, D. Aravena, E. Suturina, E. Bill, D. Maganas and F. Neese, Coord. Chem. Rev., 2015, 289–290, 177–214.
40. M. R. Pederson and S. N. Khanna, Phys. Rev. B: Condens. Matter Mater. Phys., 1999, 60, 9566–9572.
41. F. Neese, J. Am. Chem. Soc., 2006, 128, 10213–10222.
42. F. Neese, J. Chem. Phys., 2007, 127, 164112.
43. R. Maurice, R. Bastardis, C. de Graaf, N. Suaud, T. Mallah and N. Guihéry, J. Chem. Theory Comput., 2009, 5, 2977–2984.
44. M. Radoń and G. Drabik, J. Chem. Theory Comput., 2018, 14, 4010–4027.
45. D. Jacquemin, F. Kossoski, F. Gam, M. Boggio-Pasqua and P.-F. Loos, J. Chem. Theory Comput., 2023, 19, 8782–8800.
46. A. Heil, M. Kleinschmidt and C. M. Marian, J. Chem. Phys., 2018, 149, 164106.
47. G. Li Manni, R. K. Carlson, S. Luo, D. Ma, J. Olsen, D. G. Truhlar and L. Gagliardi, J. Chem. Theory Comput., 2014, 10, 3669–3680.
48. S. J. Stoneburner, D. G. Truhlar and L. Gagliardi, J. Chem. Phys., 2018, 148, 064108.
49. S. J. Stoneburner, D. G. Truhlar and L. Gagliardi, J. Phys. Chem. A, 2020, 124, 1187–1195.
50. P. Sharma, V. Bernales, D. G. Truhlar and L. Gagliardi, J. Phys. Chem. Lett., 2019, 10, 75–81.
51. P. R. Hertler, A. Sauza-de la Vega, A. Darù, A. Sarkar, R. A. Lewis, G. Wu, L. Gagliardi and T. W. Hayton, Chem. Sci., 2024, 15, 16559–16566.
52. A. Darù, C. Martín-Fernández and J. N. Harvey, ACS Catal., 2022, 12, 12678–12688.
53. N. Jiang, A. Darù, v. Kunstelj, J. G. Vitillo, M. E. Czaikowski, A. S. Filatov, A. Wuttig, L. Gagliardi and J. S. Anderson, J. Am. Chem. Soc., 2024, 146, 12243–12252.
54. M. L. Ploeger, A. Darù, J. N. Harvey and X. Hu, ACS Catal., 2020, 10, 2845–2854.
55. K. Andersson, P. Malmqvist and B. O. Roos, J. Chem. Phys., 1992, 96, 1218–1226.
56. P. Pulay, Int. J. Quantum Chem., 2011, 111, 3273–3279.
57. R. Pandharkar, M. R. Hermes, D. G. Truhlar and L. Gagliardi, J. Phys. Chem. Lett., 2020, 11, 10158–10163.
58. K. R. Mullin, D. W. Laorenza, D. E. Freedman and J. M. Rondinelli, Phys. Rev. Res., 2023, 5, L042023.
59. J. Telser, EPR Interactions – Zero-field Splittings, Wiley, 2017.
60. N. P. Kazmierczak, K. M. Luedecke, E. T. Gallmeier and R. G. Hadt, J. Phys. Chem. Lett., 2023, 14, 7658–7664.
61. S. Koseki, N. Matsunaga, T. Asada, M. W. Schmidt and M. S. Gordon, J. Phys. Chem. A, 2019, 123, 2325–2339.
62. E. M. Zolnhofer, G. B. Wijeratne, T. A. Jackson, S. Fortier, F. W. Heinemann, K. Meyer, J. Krzystek, A. Ozarowski, D. J. Mindiola and J. Telser, Inorg. Chem., 2020, 59, 6187–6201.
63. L. Baldinelli, D. Sorbelli, M. Toriyama, G. Bistoni, F. De Angelis and G. Galli, ChemRxiv, 2025, preprint,  DOI:10.26434/chemrxiv-2025-79sw7-v2.
64. C. M. Breneman and K. B. Wiberg, J. Comput. Chem., 1990, 11, 361–373.
65. A. Reinholdt, D. Pividori, A. L. Laughlin, I. M. DiMucci, S. N. MacMillan, M. G. Jafari, M. R. Gau, P. J. Carroll, J. Krzystek, A. Ozarowski, J. Telser, K. M. Lancaster, K. Meyer and D. J. Mindiola, Inorg. Chem., 2020, 59, 17834–17850.
66. M. J. Frisch, G. W. Trucks, H. B. Schlegel, G. E. Scuseria, M. A. Robb, J. R. Cheeseman, G. Scalmani, V. Barone, G. A. Petersson, H. Nakatsuji, X. Li, M. Caricato, A. V. Marenich, J. Bloino, B. G. Janesko, R. Gomperts, B. Mennucci, H. P. Hratchian, J. V. Ortiz, A. F. Izmaylov, J. L. Sonnenberg, D. Williams-Young, F. Ding, F. Lipparini, F. Egidi, J. Goings, B. Peng, A. Petrone, T. Henderson, D. Ranasinghe, V. G. Zakrzewski, J. Gao, N. Rega, G. Zheng, W. Liang, M. Hada, M. Ehara, K. Toyota, R. Fukuda, J. Hasegawa, M. Ishida, T. Nakajima, Y. Honda, O. Kitao, H. Nakai, T. Vreven, K. Throssell, J. A. Montgomery Jr, J. E. Peralta, F. Ogliaro, M. J. Bearpark, J. J. Heyd, E. N. Brothers, K. N. Kudin, V. N. Staroverov, T. A. Keith, R. Kobayashi, J. Normand, K. Raghavachari, A. P. Rendell, J. C. Burant, S. S. Iyengar, J. Tomasi, M. Cossi, J. M. Millam, M. Klene, C. Adamo, R. Cammi, J. W. Ochterski, R. L. Martin, K. Morokuma, O. Farkas, J. B. Foresman and D. J. Fox, Gaussian 16 Revision A.03, Gaussian Inc. Wallingford CT, 2016.
67. J. Tao, J. P. Perdew, V. N. Staroverov and G. E. Scuseria, Phys. Rev. Lett., 2003, 91, 146401.
68. V. N. Staroverov, G. E. Scuseria, J. Tao and J. P. Perdew, J. Chem. Phys., 2003, 119, 12129–12137.
69. A. D. Becke, J. Chem. Phys., 1993, 98, 5648–5652.
70. C. Lee, W. Yang and R. G. Parr, Phys. Rev. B: Condens. Matter Mater. Phys., 1988, 37, 785–789.
71. S. H. Vosko, L. Wilk and M. Nusair, Can. J. Phys., 1980, 58, 1200–1211.
72. Y. Zhao and D. G. Truhlar, Theor. Chem. Acc., 2008, 120, 215–241.
73. J. P. Perdew, K. Burke and M. Ernzerhof, Phys. Rev. Lett., 1996, 77, 3865–3868.
74. H. Schröder, A. Creon and T. Schwabe, J. Chem. Theory Comput., 2015, 11, 3163–3170.
75. F. Weigend and R. Ahlrichs, Phys. Chem. Chem. Phys., 2005, 7, 3297–3305.
76. B. O. Roos, Adv. Chem. Phys., 1987, 69, 399–445.
77. G. Li Manni, I. F. Galván, A. Alavi, F. Aleotti, F. Aquilante, J. Autschbach, D. Avagliano, A. Baiardi, J. J. Bao, S. Battaglia, L. Birnoschi, A. Blanco-González, S. I. Bokarev, R. Broer, R. Cacciari, P. B. Calio, R. K. Carlson, R. Carvalho Couto, L. Cerdán, L. F. Chibotaru, N. F. Chilton, J. R. Church, I. Conti, S. Coriani, J. Cuéllar-Zuquin, R. E. Daoud, N. Dattani, P. Decleva, C. de Graaf, M. G. Delcey, L. De Vico, W. Dobrautz, S. S. Dong, R. Feng, N. Ferré, M. Filatov(Gulak), L. Gagliardi, M. Garavelli, L. González, Y. Guan, M. Guo, M. R. Hennefarth, M. R. Hermes, C. E. Hoyer, M. Huix-Rotllant, V. K. Jaiswal, A. Kaiser, D. S. Kaliakin, M. Khamesian, D. S. King, V. Kochetov, M. Krośnicki, A. A. Kumaar, E. D. Larsson, S. Lehtola, M.-B. Lepetit, H. Lischka, P. López Ríos, M. Lundberg, D. Ma, S. Mai, P. Marquetand, I. C. D. Merritt, F. Montorsi, M. Mörchen, A. Nenov, V. H. A. Nguyen, Y. Nishimoto, M. S. Oakley, M. Olivucci, M. Oppel, D. Padula, R. Pandharkar, Q. M. Phung, F. Plasser, G. Raggi, E. Rebolini, M. Reiher, I. Rivalta, D. Roca-Sanjuán, T. Romig, A. A. Safari, A. Sánchez-Mansilla, A. M. Sand, I. Schapiro, T. R. Scott, J. Segarra-Martí, F. Segatta, D.-C. Sergentu, P. Sharma, R. Shepard, Y. Shu, J. K. Staab, T. P. Straatsma, L. K. Sørensen, B. N. C. Tenorio, D. G. Truhlar, L. Ungur, M. Vacher, V. Veryazov, T. A. Voß, O. Weser, D. Wu, X. Yang, D. Yarkony, C. Zhou, J. P. Zobel and R. Lindh, J. Chem. Theory Comput., 2023, 19, 6933–6991.
78. K. Pierloot, B. Dumez, P. O. Widmark and B. O. Roos, Theor. Chim. Acta, 1995, 90, 87–114.
79. D. Peng and M. Reiher, Theor. Chem. Acc., 2012, 131, 1081.
80. D. Peng and K. Hirao, J. Chem. Phys., 2009, 130, 044102.
81. M. Reiher and A. Wolf, J. Chem. Phys., 2004, 121, 2037–2047.
82. M. Reiher and A. Wolf, J. Chem. Phys., 2004, 121, 10945–10956.
83. M. Reiher, Theor. Chem. Acc., 2006, 116, 241–252.
84. F. Aquilante, T. B. Pedersen and R. Lindh, J. Chem. Phys., 2007, 126, 19416.
85. L. F. Chibotaru and L. Ungur, J. Chem. Phys., 2012, 137, 064112.
86. C. Zhou, D. Wu, L. Gagliardi and D. G. Truhlar, J. Chem. Theory Comput., 2021, 17, 5050–5063.
87. D. Wu, C. Zhou, J. J. Bao, L. Gagliardi and D. G. Truhlar, J. Chem. Theory Comput., 2022, 18, 2199–2207.

---

Footnotes

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

‡ These authors contributed equally to this work.

---