Multiconfiguration Pair-Density Functional Theory for Chromium(IV) Molecular Qubits.

Pseudotetrahedral organometallic complexes containing chromium(IV) and aryl ligands have been experimentally identified as promising molecular qubit candidates. Here we present a computational protocol based on multiconfiguration pair-density functional theory for computing singlet-triplet gaps and zero-field splitting (ZFS) parameters in Cr(IV) aryl complexes. We find that two multireference methods, multistate complete active space second-order perturbation theory (MS-CASPT2) and hybrid multistate pair-density functional theory (HMS-PDFT), perform better than Kohn-Sham density functional theory for singlet-triplet gaps. Despite the very small values of the ZFS parameters, both multireference methods performed qualitatively well. MS-CASPT2 and HMS-PDFT performed particularly well for predicting the trend in the ratio of the rhombic and axial ZFS parameters, |E/D|. We have also investigated the dependence and sensitivity of the calculated ZFS parameters on the active space and the molecular geometry. The methodologies outlined here can guide future prediction of ZFS parameters in molecular qubit candidates.

Introduction

Advances in quantum computing, quantum sensing, and high-density data storage require fine control of the superposition states in qubits—the fundamental building units of quantum information.[1−4] A variety of different materials have been proposed to actualize these qubits, including superconducting circuits,[5] trapped ions,[6] silicon-based qubits,[7] rare-earth ions in crystals,[8] and nitrogen-vacancy (NV–) defects in diamond.[9,10] Each of these materials have physical properties and features that can be synthetically controlled, making them potentially promising qubit candidates.

Some paramagnetic transition metal complexes have also been proven to be reliable spin qubits. They are optically addressable for photoluminescence readout, and they can be coherently manipulated with microwave radiation. Also, they can display coherence times longer than 100 μs; it is possible to increase their coherence times through the modification of the molecular environment.[9,11,12] These features make them promising spin qubits with properties tunable by chemical synthesis.[11,13−16]

Among these paramagnetic transition metal complexes, molecules with a triplet ground-state and small zero-field splitting (ZFS) provide two-level quantum systems addressable via microwave frequencies.[11,14,17] Strikingly, Bayliss et al.[17] showed that the ZFS parameters of pseudotetrahedral chromium(IV) aryl complexes can be tuned by modifying their organic ligands. In order to design new, better qubits, a very large chemical space of possible ligands needs to be explored. It is not feasible to achieve this by purely synthetic approaches. Therefore, a reliable computational methodology to predict the properties of target compounds needs to be developed.

Several quantum-mechanical theoretical models based on multireference methods or Kohn–Sham density functional theory (KS-DFT) have been used to predict molecular properties.[18,19] It is possible to calculate the ZFS parameters with a coupled-perturbed response treatment using density functional theory.[20] When dealing with open-shell transition metal complexes, KS-DFT performance can be questionable and is functional-dependent.[18,21−23] Due to close-lying excited states and the importance of the proper account of electron correlation of the d-orbitals in transition metal complexes, multireference methods are preferred. Multiconfiguration perturbation theory methods perform better than KS-DFT, but they are computationally expensive when dealing with large molecular systems.[24] In the present work, we used multiconfigurational pair-density functional theory (MC-PDFT) to investigate the spin-state gaps and ZFS parameters of three Cr(L)4 compounds with different aryl ligands, L = 2-methylphenyl, 2,3-dimethylphenyl, or 2,4-dimethylphenyl (Figure ). We refer to these Cr(IV) complexes as 1, 2, and 3, respectively. Multireference methods have been shown to be useful for understanding the electronic structure of these molecular qubits.[25−30]

Figure 1

Potential molecular spin-qubits based on Cr(IV) aryl complexes. (a) Distorted tetrahedral environment with the Cr(IV) center (blue) and the ligating ligand atoms (gray). The organic ligands are based on (b) toluene, (c) o-xylene, and (d) m-xylene molecules. The dashed line indicates the carbon atom bonded to the metal center.

The ZFS parameters and the triplet-singlet gap are crucial descriptors for determining whether molecular qubit candidates can be initialized and addressed. In this paper, we use multireference calculations to study the ZFS parameters and vertical excitation energies of the chromium(IV) aryl complexes shown in Figure . Because the ground triplet state is orbitally nondegenerate, its zero-field spin–orbit splitting only becomes nonzero at second order, and therefore, it is small and depends strongly on the spin–orbit couplings to the excited states. Janicka et al.[31] showed that the spin–spin contribution to the ZFS is negligibly small, and we shall not consider it here. Our paper differs from their in that they used perturbation theory, whereas we demonstrate the applicability of the MC-PDFT method, which, due to its lower cost, may be used efficiently in future work to study large chemical spaces of ligands. Another difference from previous work is that, in order to better understand the multiconfigurational character of the ZFS parameters, we examine the dependence of the calculated ZFS parameters on the choice on active space; we consider a wide range of choices.

We investigate changes in the ZFS parameters with respect to (i) molecular geometry, (ii) the choice of active space, and (iii) the choice of multireference method. We also study the origin of the small ZFS values in the examined complexes. Overall, the analysis reported here demonstrates the remarkable performance of hybrid multistate pair-density functional theory (HMS-PDFT) for the computation of triplet-singlet gaps and ZFS parameters with accuracy comparable to the well-established multistate multiconfiguration complete active space second-order perturbation theory (MS-CASPT2).

Theoretical Framework

The quantum mechanical treatment of the ZFS relies on the phenomenological Hamiltonian:where Ŝ is the pseudospin operator, and D is a real symmetric tensor that characterizes the ZFS.[32,33] After diagonalization, the Hamiltonian in eq can be rewritten as[34]By using the relation = + + and using the standard convention that the trace of D is zero, the Hamiltonian can be expressed as,[34]where the axial (D) and rhombic (E) ZFS parameters are defined to be

Different quantum technology applications place different demands on the value of D. For example, D needs to be large and negative in single-molecule magnets used for data storage applications.[1,35−37] However, for optically addressable molecular qubits, D should be small, allowing for control of the electron spin by microwave radiation.[11,17] One can try to design molecules with small D by changing the ligand field; the magnitude of D can then be measured in EPR experiments.[11]

We use multiconfiguration wave functions based on active spaces, and choosing an active space requires the selection of a set of n active electrons in m active orbitals; the resulting active space is denoted (n, m). Active space selection is not a black box operation; there are several procedures for active space selection,[38−43] and specialized considerations may be advisible when dealing with magnetic properties.[44−46] Following the protocol of ref (47), a basis of states computed without including spin–orbit coupling (called spin–orbit-free states or spin–orbit-free roots) are computed first. The spin–orbit-free states include the ground-state and various excited states with different spin multiplicities. For transition metal systems, the valence 3d orbitals are an obvious choice for inclusion in the active space; however, the resulting minimal active space is not large enough to capture a sufficient portion of the multiconfigurational character when the covalency of the metal–ligand bonds become significant. Alternately, inclusion of too many metal-to-ligand or ligand-to-metal charge transfer states may result in an unaffordable state-averaged CASSCF calculation. Thus, it is necessary to fine-tune the number of excited states used to compute the ZFS parameters.[44] The spin–orbit matrix elements are computed under the restricted active space state interaction spin–orbit (RASSI-SO) formalism.[47−49] The resulting matrix elements allow the calculation of the magnetic properties.

We calculate the spin–orbit-free state energies with the state-averaged complete active space self-consistent field method (SA-CASSCF).[50] This method provides a good qualitative description of systems with static correlation. To quantitatively describe dynamic electron correlation in chemical systems, a post-SA-CASSCF calculation is performed, using multistate complete-active-space second-order perturbation theory (MS-CASPT2),[51] multiconfiguration pair-density functional theory (MC-PDFT), or multistate pair-density functional theory (MS-PDFT). MC-PDFT and MS-PDFT are methods that provide CASPT2-quality and MS-CASPT2-quality accuracy at reduced computational cost.[52−57] The MC-PDFT or MS-PDFT computational cost is similar to that of computing the reference wave function.[58] MS-PDFT has recently been successfully used to study the g tensor[59] and the ZFS[46] of several transition metal complexes.

The MC-PDFT and MS-PDFT performance can be further improved by using a hybrid MC-PDFT (HMC-PDFT) and hybrid MS-PDFT (HMS-PDFT).[60] Similarly to Kohn–Sham DFT, in which local functionals are combined with Hartree–Fock exchange, the HMC-PDFT energy can be expressed aswhere the parameter λ indicates the amount of CASSCF energy included in the hybridization. In the present work, we consider the hybrid translated functional tPBE0 defined previously,[60] for which the λ parameter is set equal to 0.25.

Computational Methods

We performed various kinds of calculations, which we now describe in turn. All calculations were performed without imposing any symmetry.

Crystal Calculations

We optimized the geometries of the chromium(IV) aryl complexes (Figure ) in periodic calculations on the molecular crystals to consider the experimental environment. These calculations were carried out using the VASP package[61−64] with the PBE exchange-correlation functional[65] and the D3 dispersion term damped by the Becke–Johnson damping factor.[66] A plane wave kinetic energy cutoff of 600 eV was used along with the PAW pseudopotential.[67,68] A triplet spin multiplicity was considered since all the complexes studied herein have triplet ground states. The Brillouin zone was sampled only at the Γ point. The experimental crystal structure[17] was used as an initial guess, and the ionic positions and the lattice coordinates were allowed to relax without any space group symmetry imposed.

Subsequent calculations (described below) were out for a single molecule using both the molecular structure extracted from this optimized crystal and the molecular structure extracted from the experimental crystal structure.

Gas-Phase Calculations

Using the crystal structure as input, we optimized the molecular geometries for the triplet state of the complexes studied herein. The calculations were done using Orca 5.0.[69] To be consistent with the crystal calculations, we used the PBE Kohn–Sham functional and the D3 dispersion model with Becke–Johnson damping factor,[66] and Def2-TZVPP basis set.[70]

The optimized geometries were used as inputs for the following multiconfigurational calculations.

SA-CASSCF Calculations

The optimized and crystal structures were used as inputs for SA-CASSCF calculations of the triplet ground state and first singlet excited state energy. Calculations for both spin states were performed at the ground-state (triplet) geometries without spin–orbit coupling. These calculations were performed using OpenMolcas v21.10.[71] We used the ANO-RCC-VTZP basis set for the chromium atom and the ANO-RCC-VDZ basis set[72] for the carbon and hydrogen atoms. The second-order Douglas–Kroll–Hess Hamiltonian to was used to account for scalar relativistic effects.[73,74] We used the resolution of the identity and the Cholesky decomposition to speed up the two-electron integral calculations.[75]

We considered several active spaces (Figure ). The Cr(IV) ion has a 3d2 electron configuration; the minimal active space for the ion in tetrahedral geometry would be (2,2), consisting of a doubly degenerate e orbital. A (2,5) active space is constructed by adding the empty t2* orbitals. Adding a correlating subshell of d orbitals yields a (2,10) active space. Alternatively, one can add Cr–ligand occupied molecular orbitals, t2, to the to the (2,5) active space to construct the (8,8) active space that accounts better for metal–ligand covalency. Finally, for the largest active space, we added the correlating d subshell and the chromium s bonding and antibonding orbitals to the (8,8) active space to obtain the (10,15) active space. (See Figures and S7–S11. Figures and tables with a prefix S are given in the Supporting Information.)

Figure 2

Molecular orbitals of molecule 1 are included in the different active spaces considered for the multireference calculations. Left: In the diagram, the red, blue, purple, orange, and green boxes correspond to the (2,2), (2,5), (2,10), (8,8), and (10,15) active spaces. Right: Drawings of all the molecular orbitals considered for the aforementioned active spaces. The orbitals are displayed in the same order as in the diagram on the left.

The number of states averaged was 10 triplets and 15 singlets for the (2,5) and (2,10) active spaces and 7 triplets and 9 singlets for the (8,8) and (10,15) active spaces. Equal weights were used for all averaged energies.

Multiconfiguration Calculations of the Singlet–Triplet Splitting

The first excited singlet state is nearly doubly degenerate due to the near-tetrahedral geometry. Using the optimized spin-triplet structures, we calculated the vertical excitation energies to the lower of the two states. These excitation energies were calculated without spin–orbit coupling by SA-CASSCF (described above), MS-CASPT2, single-state MC-PDFT, and single-state hybrid MC-PDFT (HMC-PDFT). The reference functions for the latter three kinds of calculation were the SA-CASSCF calculations described above. The MS-CASPT2 calculations employed an imaginary shift of 0.3 au and an ionization-potential electron-affinity shift of 0.25 au;[59] the model space sizes were the same as the numbers of states averaged in the SA-CASSCF steps. The MC-PDFT calculations employed the tPBE on-top functional.[52] The HMC-PDFT calculations employ the tPBE0 functional,[60] which combines 25% CASSCF nonclassical energy with 75% tPBE nonclassical energy. A fine grid was used for the MC-PDFT and HMC-PDFT calculations.

Inclusion of Spin–Orbit Coupling

The ZFS parameters were calculated for all mentioned active spaces except for the minimal one. The ZFS calculations include spin–orbit coupling by methods describe previously;[46,49,59,76] in particular, spin–orbit coupling was included in SA-CASSCF, MS-CASPT2, MS-PDFT,[77] and HMS-PDFT calculations by using the atomic-mean field approximation to the Breit–Pauli Hamiltonian, where the spin–orbit coupling operator is diagonalized using the eigenstates obtained with the multiconfiguration methods aforementioned, to obtain the spin–orbit coupled eigenstates and energies. This procedure was widely discussed in ref (59), and the available implementation on the RASSI module[48,49] of OpenMolcas. This allows us to use molecular orbitals optimized separately for the singlet and triplet manifolds. For the latter three kinds of calculation, the states included in the model space were the same as the states averaged in the SA-CASSCF calculations. In the MS-PDFT and HMS-PDFT calculations, the basis for the model space calculations was the set of SA-CASSCF eigenvectors (i.e., we did not transform to an intermediate basis). Additional data for the on-top functionals tested is available in the Supporting Information (Figures S15–S17).

NEVPT2 Calculations

For comparison, we also performed SA-CASSCF calculations and N-electron valence-state second-order perturbation theory[78−80] (NEVPT2) calculations (with these SA-CASSCF calculations providing the reference states) using the Orca 5.0 code.[69] In these calculations, the orbitals are optimized simultaneously for the singlets and triplets, resulting in a single set of molecular orbitals. Ten triplets and 15 singlets were included in the state averaging with equal weights. We used an active space of 2 active electrons in 5 active orbitals (the 3d orbitals). Representative active space orbitals are shown in Figure S2. The Douglas–Kroll–Hess[73,74,81−83] scalar relativistic Hamiltonian and the DKH-def2-TZVP basis set were used. SCF convergence criteria with an energy tolerance of 10–7 hartrees were applied. The NEVPT2 calculations were performed on the experimental and periodic optimized Cr(IV) structures (see Tables S7 and S8).

Kohn–Sham Density Functional Calculations

The Orca 5.0(69) package was used to determine the ZFS parameters with various Kohn–Sham functionals. Details are in the Supporting Information.

Spin–Spin Coupling

We also considered the effect of spin–spin coupling (SSC) in the SA-CASSCF/NEVPT2 calculations and Kohn–Sham density functional calculations. The SA-CASSCF/NEVPT2 calculations were done without and with SSC, and this showed that SSC has only a small effect on the results (see Tables S7 and S8). SSC was not included in the OpenMolcas calculations.

Results and Discussion

Structures

It is well-known that the spin properties of the metal ion depend significantly on its coordination environment.[16,17,28,84−88] In order to study the sensitivity of the results to the geometry, we compared calculations with (i) the experimental crystal structures (measured at 100 K) from ref (17), (ii) the structures optimized with periodic boundary conditions on the unit cell, and (iii) the gas-phase optimized structures. The optimized structures differ from the former by lack of the zero-point and temperature effects.

The experimental and optimized structures showed slight differences in the structural parameters. The four Cr–ligand bond lengths in the periodic optimized complex 1 structure are 1.97 Å, in the gas-phase structure are 1.98 Å, while the bond lengths range from 1.98 to 2.00 Å in the structure obtained from the crystallographic data. The Cr–ligand bond distances are 1.98 Å in the complex 2 optimized geometries, whereas the two bond distances are 2.02 and 2.00 Å in the crystal structure. For 3, there are two different Cr–ligand bond lengths, 1.98 and 1.97 Å in the optimized geometries as compared to 2.00 and 1.99 Å in the experimental structure. The geometric parameters differ the most in the compound 1, where the largest difference in the Cr–ligand bond distances is 0.03 Å, and the largest difference in the R–Cr–R bond angles is 2.8°. For the complexes 2 and 3, the corresponding values between the experimental and periodic optimized structures are 0.04 Å and 0.1°, and 0.03 Å and 1.4°, respectively. However, between the gas-phase and experimental structures, the differences between the bond angles increase up to 5°. The experimental structures and those optimized with periodic calculations show more similarity between them (see Figures S4–S7).

Triplet–Singlet Gaps

Bayliss et al. reported for complexes 1, 2, and 3 the experimental zero-phonon lines ranging from 1.20 to 1.22 eV (1009 to 1025 nm).[17] We evaluated the capability of the methods explored herein to reproduce these reported experimental singlet–triplet gaps.

According to our calculations, the first singlet excited state is double degenerate, and the optimization of such structure is complicated with Kohn–Sham DFT based methods. Thus, we approximate the triplet-singlet gap as a vertical excitation. Table reports the triplet–singlet energy differences (ΔET-S) for the molecule 1 using the experimental, and periodic optimized geometries. We computed the vertical excitation energy at the ground (triplet) state geometry. The experimental gap from the zero-phonon line would correspond to a difference between the lowest vibrational levels of the two different spin states, which have different geometries. In Table , the MC-PDFT and HMC-PDFT results are labeled tPBE and tPBE0, respectively. The close agreement between the (2,2) MS-CASPT2 gap and the experimental value may be attributed to some cancellation of errors. When virtual orbitals are added to form the (2,5) and (2,10) active spaces, the triplet-singlet MS-CASPT2 energy gap is reduced by about 0.2–0.3 eV, with respect to the minimal active space. In contrast, tPBE0 is less affected by adding empty orbitals. Further inclusion of Cr–ligand bonding and antibonding orbitals causes both the MS-CASPT2 and tPBE0 gaps to increase. However, the tPBE0 value displays a better agreement with the experimental findings when using the (10,15) active space. With the largest active space, SA-CASSCF overestimates the gap by 0.29–0.46 eV, depending on the geometry, MS-CASPT2 overestimates it by 0.21–0.31 eV, MC-PDFT underestimates it by 0.18 or 0.04 eV, and HMC-PDFT is accurate within 0.04 eV at either geometry.

Table 1

Calculated Triplet–Singlet Gaps (ΔET-S, in eV) Using the Periodic Optimized Structure of Molecule 1 with the SA-CASSCF, MS-CASPT2, tPBE, and tPBE0 Methods

ΔET-S for optimized geometry
active space	SA-CASSCF	MS-CASPT2	tPBE	tPBE0	expt.
(2,2)	1.87	1.23	0.55	0.88
(2,5)	1.89	0.94	0.48	0.83
(2,10)	1.81	0.98	0.51	0.83
(8,8)	1.83	1.41	0.94	1.16
(10,15)	1.66	1.41	1.02	1.18	1.20

The ΔET-S values for 1, 2, and 3 are very similar, 1.18, 1.22, and 1.21 eV, respectively, as indicated in Tables , 2, and S9–S11. In agreement with this trend, our HMC-PDFT calculations also show only slight differences among the three complexes, as shown in Table . This table also shows Kohn–Sham density functional calculations from the literature[88] that are consistently larger than experiment and show much larger variation with changing ligands. We performed additional Kohn–Sham theory with several functionals; the results are given in Tables S1–S3, and they show systematic overestimation of the energy gap for all functionals tested. We conclude that it is better to use multireference methods to study this kind of problem. Table also shows the MS-CASPT2 for all three complexes. MS-CASPT2 agrees that the variations in the gap with the changes of ligands are small, but it systematically overestimates the magnitudes of the gaps. The finding that HMC-PDFT agrees better with experiment than does MS-CASPT2 is very encouraging, since HMC-PDFT has a lower computational cost with respect to MS-CASPT2.

Table 2

Calculated Triplet–Singlet Gaps (ΔET-S) with the (10,15) Active Space for the Periodic Optimized Geometry of Molecules 1, 2, and 3 Using MS-CASPT2, tPBE0, and B3LYP+D3(BJ) Methodsa

complex	MS-CASPT2	tPBE0	B3LYP+D3(BJ)b	expt.
1	1.41	1.18	1.63	1.20
2	1.46	1.22	1.47	1.22
3	1.46	1.21	1.31	1.20

The energy gaps are reported in eV.

The DFT values were obtained from ref (88).

Zero-Field Splittings

Figure shows the |D| and |E| parameters of the complex 1, in blue for the experimental structure and in green for the structure optimized in the crystal, with both employing the (2,5) active space. The HMS-PDFT results are labeled tPBE0 in the figures. To put the size of the splitting into perspective, we note that the experimental |D| value of 3.53 GHz is only 0.118 cm–1 = 0.0146 meV.

Figure 3

Computed ZFS parameters |D| and |E| of the complex 1 using the (2,5) active space and the SA-CASSCF, MS-CASPT2, and HMS-PDFT(tPBE0) methods. The vertical bars represents the absolute magnitude of the parameters for the experimental structure (blue) and the optimized one (green). The horizontal line represents the experimental values.[17]

The rhombic parameter (E) is zero by tensorial symmetry for a system with a 3-fold or higher rotation axis. The present system is slightly distorted, and we calculate a small but nonzero E parameter with the experimental structure, although the result reported experimentally[17] is zero. However, when we use the optimized geometry, our calculated E value is indeed zero. This shows the importance of the geometry.

The results Figure also indicate that |D| is also closer to experiment when calculated at the optimized geometry. The SA-CASSCF calculations overestimate |D| by more than a factor of 2, but MS-CASPT2 and HMS-PDFT are accurate within 5%. In ref (17), the authors were not able to obtain the sign of D. However, our calculations consistently predict that D is negative, meaning that for the triplet ground-state, the M = 0 level is above the M = ± 1 levels).

Figure shows the computed |D| parameter for the system 1 at the optimized geometry with different active spaces. We see that the value of |D| is not sensitive to the change of active space for MS-CASPT2 and for three of the four active spaces for HMS-PDFT, even though the SA-CASSCF results differ from one another by as much a factor of 2.5. Although the crystallographic data suggest the complex 1 has a S4 symmetry, we obtained nonzero |E| values when using the experimental structure (Figure ). This can be attributed to the disorder in the crystal structure. When using the optimized structure, all active spaces in Figure , yield an |E| value of zero–matching the experimental observations.

Figure 4

Calculation of |D| parameter of optimized complex 1 geometry with the SA-CASSCF, MS-CASPT2, and HMS-PDFT(tPBE0) methods. The active spaces used are (2,5) in blue, (2,10) in gray, (8,8) in green, and (10,15) in purple. The dashed line corresponds to the experimental data.[17]

The inclusion of more triplets and singlets in the RASSI calculation does not improve the computed value of the axial parameter |D| of molecule 1, as shown in Tables S12–S15. The (8,8) MS-CASPT2 value for |D| is 3.33 GHz when using an energy cutoff of 4.3 eV and including 7 triplets and 9 singlets. When the cutoff is increased to 5.0 eV, thereby including 18 triplets and 16 singlets, there is a decrease in the magnitude of D, resulting in worse agreement with experiment. This behavior is also found for HMS-PDFT. This is also true for the (10,15) active space; when using more triplet and singlet excited states there is a decrease of the magnitude of |D|. For both active spaces, the computed |E| parameter is always zero regardless of the number of spin–orbit-free states included.

For the compounds 2 and 3, there is also a decrease in the quality of the calculated ZFS parameters when going beyond 7 triplets and 9 singlets in the RASSI calculation (Tables S16–S23). Because these 7 triplets and 9 singlets are mostly metal-to-ligand charge transfer (MLCT) states, we conclude that the MLCT states rather than ligand-to-ligand or ligand-to-metal charge transfer states are most useful to obtain a qualitative and semiquantitative description of the axial parameter.

According to the experimental data, the axial ZFS parameter values are in the following order: |D| > |D| > |D|. However, when using the optimized geometry in the crystal, all of the theoretical methods showed a trend of |D| > |D| > |D|, as shown in Figure a. The figure shows that all theoretical methods correctly predict that D is smallest for 2, but they incorrectly predict the order for 1 and 3. Since D is inversely related to the energy gaps between the ground and excited states, a possible rationalization of the trend in the computed D values of the three complexes might be based on the calculated values of the spin–orbit-free energies. The calculated spin–orbit-free triplet spectrum is given in Tables S24–S27, and we observe that the order of the |D| for the three complexes is opposite to the order of the calculated energy difference between the ground state and the first triplet excited state.

Figure 5

Comparison of experimental and calculated values using the (10,15) active space for (a) the axial ZFS parameter |D|, (b) the rhombic ZFS parameter |E|, and (c) the ratio |E/D|. The blue bars correspond to the experimental measurements, and the other bars are for the following electronic structure methods: SA-CASSCF (brown), MS-CASPT2 (green), and HMS-PDFT(tPBE0) (red).

Figure b shows that none of theoretical calculations reproduces the experimentally observed trend that 2 has a larger rhombic parameter than 3, although they all correctly predict that |E| = 0 for compound 1.

Figure c compares the ratios |E/D|. The experimental trend |E/D| > |E/D| > |E/D| is reproduced by all methods when using (10,15) active space. The figure shows that the ratio |E/D|/|E/D| is 2.4 experimentally, 1.6 by SA-CASSCF, 2.4 by MS-CASPT2, and 1.2 by HMS-PDFT(tPBE0); thus, the HMS-PDFT value is semiquantitative. For the (2,5) and (2,10) active spaces, the MS-CASPT2 method shows an opposite trend (ratio less than 1). Figures S21–S24 show that there is no major change upon varying the active space; in particular, we obtain similar qualitative results when using either the (2,5) or (10,15) active spaces with HMS-PDFT(tPBE0) and MS-CASPT2. This observation may indicate some error cancellation while computing the quotient |E/D|, which allows HMS-PDFT to obtain a good estimate of this ratio. This result is important because knowledge of |E/D| and the g tensor is useful to recognize important spin systems[37,89−91] Thus, using the (2,5) active space and MS-CASPT2 or the HMS-PDFT(tPBE0) method can provide a good initial exploration of molecular qubit candidates. It is encouraging that HMS-PDFT(tPBE0) provides the same qualitative information as MS-CASPT2, but at a lower computational cost as observed in several studies for a wide range of molecules and properties (more information available in ref (58)). Hence, it may be affordable to explore larger molecular systems with this method.

In addition, when the gas-phase optimized geometries are used to compute the ZFS parameters, we find the experimental trend for |D| is satisfied for all methods when using (8,8) and (10,15) active spaces. However, the values for the rhombic parameter, |E|, are approximately zero (Figures S25–S28). Therefore, we obtain slightly different structures depending on the approximation used to optimize the molecular structures. Still, these differences significantly affected the ZFS parameters due to their small magnitude.

Tables S4–S6 shows that Kohn–Sham density functional calculations of the ZFS parameters using the coupled-perturbed method[20] do not improve over the results we obtained with multiconfigurational methods. Moreover, Figure S1 shows that the Kohn–Sham density functional calculations show strong dependence on the choice of functional, which would lead to uncertainty in using calculations to guide the synthesis of molecular qubit candidates. In contrast, the tested multiconfigurational methods do not show such a strong strong method dependence for the computation of the axial and rhombic parameters (Figures S18–S20).

Concluding Remarks

We have presented the results of multireference calculations for the axial and rhombic ZFS parameters, |D| and |E|, of three Cr(IV) aryl complexes.

The accurate prediction of ZFS parameters is a challenge for electronic structure theory. Our calculations show that the molecular geometry plays a vital role in obtaining values for |D| and |E|; small changes in the molecular geometry generate significant changes in these ZFS parameters. Gas-phase methods are useful to obtain molecular geometries. Still, periodic methods may be needed because of the geometric distortions resulting from crystal packing that may play a major role in tuning |D| and |E|. Therefore, we recommend the optimization of molecular geometries using periodic calculations followed by multireference calculations of the spin properties of molecular qubit candidates.

The HMS-PDFT(tPBE0) and MS-CASPT2 zero-field splitting values do not show a significant dependence on the active space. Furthermore, since it is impossible to consider all the spin–orbit-free roots for active spaces larger than (2,5), one has to choose a cutoff for the included excited states. The choice of cutoff for the studied complexes should consider the MLCT character of the excited states. For the compounds studied here, the energy cutoff was taken as 4.3 eV; the inclusion of higher-energy roots in the RASSI calculations did not improve in the results.

An appropriate singlet–triplet gap is essential for achieving photoluminescence readout of molecular spin qubits. The tPBE0 hybrid on-top functional with the (10,15) active space accurately reproduces the energy difference between the triplet ground state and first excited singlet state. Furthermore, all the examined multireference methodologies qualitatively reproduced the general magnitudes of |D| and |E|, although not all the trends between complexes; however, the |E/D| trend is in agreement with experiment. Hence, we expect that hybrid multistate pair-density functional theory will be useful to predict the optical and microwave addressability of unsynthesized molecular qubits via the accurate prediction of ZFS parameters at an affordable computational cost.