{'fraction_non_alphanumeric': 0.029346012154035078, 'fraction_numerical': 0.01275229068925344, 'mean_word_length': 4.995846705682462, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Zero-Field Splitting Parameters from Four-Component Relativistic Methods

We report an approach for determination of zero-field splitting parameters from four-component relativistic calculations. Our approach involves neither perturbative treatment of spin-orbit interaction nor truncation of the spin-orbit coupled states. We make use of a multi-state implementation of relativistic complete active space perturbation theory (CASPT2), partially contracted N-electron valence perturbation theory (NEVPT2), and multi-reference configuration interaction theory (MRCI), all with the fully internally contracted ansatz. A mapping is performed from the Dirac Hamiltonian to the pseudospin Hamiltonian, using correlated energies and the magnetic moment matrix elements of the reference wavefunctions. Direct spin-spin coupling is naturally included through the full 2-electron Breit interaction. Benchmark calculations on chalcogen diatomics and pseudotetrahedral cobalt(II) complexes show accuracy comparable to the commonly used state-interaction with spin-orbit (SI-SO) approach, while tests on a uranium(III) single-ion magnet suggest that for actinide complexes the strengths of our approach through the more robust treatment of spin-orbit effects and the avoidence of state truncation are of greater importance.


I. INTRODUCTION
The past two decades have witnessed a continuous growth of interest in molecular magnets due to their potential for use in quantum information processing, high-density information storage, and spintronics. 1-3 The magnetic properties of open-shell complexes are conventionally modeled using a phenomenological spin Hamiltonian, whose parameters are determined experimentally via electron paramagnetic resonance spectroscopy, magnetometry, inelastic neutron scattering, and other methods. 4-6 To second order, the spin Hamiltonian is expressed aŝ with axial and rhombic zero-field splitting (ZFS) parameters D and E, the Bohr magneton µ B , external magnetic field B, the g-matrix g expressing the strength and anisotropy of the electronic Zeeman effect, and spin operatorŜ. Additional terms can be added for smaller effects such as hyperfine coupling, interaction with the nuclear quadrupole moment, higher-order ZFS contributions, and so on. One of the key parameters necessary for single-ion magnetism is easy axis ZFS, which contributes to the energy barrier to reversal of spin magnetization. The ZFS is a breaking of degeneracy of the ground spin multiplet of an open-shell complex with total spin of 1 or greater, which results from spin-orbit and spin-spin coupling, and its accurate determination is essential for an understanding of the behavior of molecular magnets. Several methods are already available for the prediction of spin Hamiltonian parameters. The anisotropic g-matrix has been predicted using both relativistic density functional theory and wavefunction-based methods. [7][8][9][10][11][12][13][14] For the prediction of zero-field energy splittings, the state-interaction with spin-orbit (SI-SO) approach has been widely used. 15 Van den Heuvel et al. 16 have introduced a simplified approach where configuration-averaged Hartree-Fock orbitals are used in a complete active space configuration interaction (CASCI) diagonalization with spin-orbit effects included. This gives comparable results to complete active space self-consistent field (CASSCF) followed by SI-SO for lanthanide systems with weak crystal field effects, although neither method could accurately reproduce the experimental low-energy excitation spectra. 17 Contributions from direct spin-spin coupling have been introduced in the context of density functional theory (DFT). 18 In this work, we demonstrate a robust method for the ab initio determination of ZFS parameters of mononuclear complexes with strong spin-orbit and spin-spin coupling. Our approach makes use of four-component relativistic CASSCF, [19][20][21] which can be extended with perturbation theory for dynamical correlation 22,23 or the Gaunt or full Breit interaction 24 to include direct spin-spin coupling. To extract spin Hamiltonian parameters from the ab initio energies and wavefunctions, we have adapted the pseudospin mapping approach of Chibotaru and Ungur 12 for the four-component wavefunctions used here.
We begin by reviewing the relevant theory for density fitted, 4-component relativistic CASSCF, fully internally contracted complete active space perturbation theory (CASPT2), and partially contracted N-electron valence perturbation theory (NEVPT2). We also review the procedure for mapping to the pseudospin Hamiltonian before going on to benchmark the methods for several series of molecules. We start with a series of chalcogen diatomics, comparing to existing methods and examining the importance of multi-state approaches in each case. We then present example calculations of ZFS parameters for single-ion magnets with moderate and very strong spin-orbit effects. We begin with the pseudotetrahedral cobalt complexes [Co(XPh) 2− 4 for X = O, S, and Se, [25][26][27][28] , then move on to the actinide complex U(H 2 BPz 2 ) 3 , [29][30][31][32] where spin-orbit effects are much stronger and should ideally be included concurrently with orbital optimization.

A. Relativistic CASSCF
Throughout this section, we use I, J, and K to label Slater determinants, while L, M, and N refer to CASSCF states. The molecular orbital indices i and j refer to closed orbitals, r, s, and t to active orbitals, a and b to virtual orbitals; w, x, y, and z are used as generic orbital indices, and µ is used for atomic basis functions. Our four-component molecular spinors are constructed from a basis set of 2-spinor atomic basis functions generated using restricted kinetic balance (RKB). 33,34 Starting orbitals are generated from a converged Dirac-Hartree-Fock wavefunction, 24 with active electrons removed when appropriate.
At the level of CASSCF, the multi-configuration wavefunction is parametrized as where we seek to find eigenfunctions of the Dirac Hamiltonian within the CASCI space: The relativistic Dirac-Breit Hamiltonian can be expressed in second quantized form aŝ with the excitation operators defined and the nuclear potential modeled with a Gaussian charge distribution with exponent ζ A and total charge Z A . 35 We can also obtain the Dirac-Coulomb-Gaunt or Dirac-Coulomb Hamiltonian by neglecting the last one or two terms (respectively) of Eq. (6). 24 We perform a two-step optimization procedure to obtain the CASSCF wave functions, where for each macroiteration the state-specific energies are made stationary with respect to c I,M coefficients, i.e., The c I,M coefficients are obtained using the relativistic analogue of full configuration interaction, described previously by one of the authors. 19,20 The molecular coefficients are optimized through a series of macroiterations, which makes the state-averaged energy stationary with respect to orbital coefficients C µp with appropriate normalizations, The orbital optimization is performed using a second-order minimax algorithm, in which explicit construction of the orbital Hessian matrix is avoided by direct multiplication of density fittted integrals to trial vectors; further details are given in Ref. 21.

B. Relativistic XMS-CASPT2 and NEVPT2
For the relativistic dynamical correlation theories, we employ a fully internally contracted (FIC) ansatz for the first-order wavefunction, using either the multi-state, multireference (MS-MR) parametrization: or the single-state, single-reference (SS-SR) parametrization: The former choice comes at higher cost but ensures the invariance of the correlated wavefunction with respect to rotation among reference states. For post-CASSCF inclusion of dynamical correlation, we make use of the extended multi-state (XMS) algorithm for quasi-degenerate perturbation theory proposed by Granovsky. 36,37 This is necessary to ensure invariance among rotations between reference states, and can be achieved by diagonalizing the zeroth-order Hamiltonian in the basis of reference states: Using the rotated reference functions |M , the implementation of XMS-CASPT2 is the same as conventional MS-CASPT2. There is no need for an extended multi-state algorithm for NEVPT2, because in this case the zeroth-order Hamiltonian is already diagonal. The zeroth-order Hamiltonian for (X)MS-CASPT2 is given by the state-averaged generalized Fock operator: using with a projection to the reference spacê and state-averaged first-order density matrix where N st is the number of states included. E s is a vertical shift parameter. The working equations for single-state relativistic fully internally contracted CASPT2 and multi-reference configuration interaction (MRCI) were reported by one of us. 22 The working equations for the residual, source, and norm terms were implemented using the code-generator SMITH3 38 and interfaced to the XMS-CASPT2 architecture using the procedure outlined in Ref. 37.
The effective Hamiltonian is constructed including shift corrections for all elements: the eigenvalues of which correspond to our XMS-CASPT2 second-order energies.
Dyall proposed an alternative formulation of multireference perturbation theory, where 2-electron interactions among active orbitals are fully included in the zeroth-order Hamiltonian, which he called the CAS/A operator. 39 Using only the diagonal blocks of Dyall's zeroth-order Hamiltonian allows us to run partially contracted NEVPT2. 40,41

C. Mapping to Pseudospin Hamiltonian
For complexes with spin larger than 1, the energies alone are insufficient to define the ZFS parameters, so the ab initio states must be mapped to a model Hamiltonian. 12,42 In order to determine the spin Hamiltonian parameters, we require a mapping from the eigenfunctions of the Hamiltonian to a suitable basis of effective or fictitious spin states. We first compute the matrix elements of the magnetic moment operator using the four-component analogue of the Pauli spin operator, where σ contains the three standard 2-component Pauli spin matrices. 43 We build the matrix representation ofμ z within the selected manifold of states and diagonalize it. By convention, z is chosen to be one of the main magnetic axes of the complex, which can be determined as described by Chibotaru and Ungur. 12,44 The column eigenvectors of µ z give us the appropriate combinations of Hamiltonian eigenfunctions that should be put in correspondence of pseudospin eigenfunctions, with one complication. The phases of the column eigenvectors are arbitrary, while certain phase restrictions on the pseudospin eigenfunctions must be enforced. The first requirement is time-reversal symmetry, the restoration of which is considered in the Supporting Information. Apart from time-reversal symmetry, the relative phase between pseudospin eigenfunctions is in principle arbitrary; however, we must ensure that the relative phase of true angular momentum matrix elements matches that of the pseudospin operators. The conventional choice is that which makes the Clebsch-Gordon coefficients real; 45 it can be enforced in our context by choosing a phase transformation such that all the m s |μ x |m s +1 matrix elements are real and positive. This requirement ensures that the pseudospin reduces to the true spinŜ at the limit of zero spin-orbit and spin-spin coupling and leads to consistent behavior under rotation of the coordinate system.

D. Extracting the Zero-Field Splitting Parameters
The numerical pseudospin Hamiltonian is constructed from eigenstates of the ab initio Hamiltonian.
where E i is the energy of state i, and U contains the phaseadjusted eigenvectors of µ z as discussed above and further in the Supporting Information. When the pseudospin Hamiltonian is built up in terms of extended Stevens operators, we have in matrix form The extended Stevens operators are built from pseudospin operators; for ZFS, the nonzero contributions will have even k between 2 and 2S , while q runs from −k to k. The general algorithm for generation of the extended Stevens operators has been published by Ryabov 46,47 and is concisely recapitulated in the Supporting Information. Setting Eqs. (32) and (33) equal leads to a set of linear equations for the coefficients B q k . The second-order contributions are more often expressed in terms of the D-tensor with the parameters and the conventional D and E parameters derived from the diagonalized form of D.
For dynamically correlated calculations, we approximate the mapping to pseudospin states by using rotated reference CASSCF reference functions to determine the magnetic moment matrices, Using these matrices, the transformation matrix U is determined as normal, but the dynamically correlated energies are used when building the Hamiltonian in Eq. (31). In this way, we capture the effects of correlation on energy splittings as well as any mixing that occurs in (extended) multi-state methods, while avoiding the computation of relaxed density matrices associated with the dynamically correlated methods. We consider this approach similar in spirit to the common practice with state interaction methods of using correlated energies along only the diagonal elements of the spin-orbit Hamiltonian. 48

III. COMPUTATIONAL DETAILS
Four-component relativistic calculations were performed with the BAGEL electronic structure package 49,50 developed by our group. We use spinor basis functions and make use of density fitting with a customized fitting basis set described below. Except where otherwise noted, the Dirac-Coulomb relativistic Hamiltonian was used. Reference data using the SI-SO approach were computed in MOLCAS 51 and made use of density fitting with an atomic compact fitting basis set 52 generated on the fly.
There is some difficulty in the use of contracted orbital basis sets with four-component relativistic calculations. Part of the basis set incompleteness error stems from the fact that convetional basis sets are optimized using a different Hamiltonian operator, but there are additional complications due to restricted kinetic balance. Because the RKB condition does not exactly reproduce the relationship between large and small components in the presence of relativistic effects, a strongly contracted basis function optimized for the large component will generate a small component basis function that is not optimal. 53 Dyall has made some effort to develop basis sets for fully relativistic calculations with different contractions for each of the four components. [54][55][56][57] Currently, our program can only accommodate one set of contraction coefficients in accord with the RKB condition, and in our previous work we found Dyall's spin-free basis sets to be inadequate for our needs. 58 The safest approach is to fully decontract the basis set, but many calculations would be intractable with the large primitive basis sets that result.
In this work, the Atomic Natural Orbitals -Relativistically Core Correlated (ANO-RCC) basis sets 59 are used for all atoms. In Section IV A, results for a series of small molecules are given for both the uncontracted form and the triple-ζ plus polarization (TZP) contraction of this basis set. For fourcomponent calculations on the larger molecules studied in this work, we compromise by using an uncontracted basis on the central metal atom and contracted basis sets on the ligands, with the set of contractions specified in the appropriate section. For the SI-SO comparison data, the TZP contraction is used for central metal atoms, because the atomic mean field integrals used in the SI-SO cannot be accurately computed with a primitive basis set in MOLCAS. 60 We require a customized auxiliary basis set for density fitting because standard fitting basis sets do not include heavy atoms or contain the flexibility needed to capture the Gaunt for full Breit interaction. 24 To meet these demands, we have augmented standard basis sets in an even-tempered manner by adding primitive functions with Gaussian exponents increasing by a factor of 2.5. For H atoms we have simply decontracted the standard TZVPP-JKFIT auxiliary basis set. For B, C, N, O, S, Co, and Se, we have decontracted the TZVPP-JKFIT auxiliary basis set, and we added one tight function to each of the p, d, and f shells. For Te and U, we begin with the decontracted ANO-RCC basis set. We add polarization functions by applying the exponents of the highest angular momentum functions to higher shells up to i-type (i.e., for Te the exponents of 3 g functions are copied to the h and i shells, and for U the 3 h exponents are copied to i). We then add several tight functions to each shell, with each one's exponent greater than the last by a factor of 2.5. The numbers of functions added are as follows: 3 s, 3 p, 6 d, 7 f , 5 g, and 1 h for Te, and 4 s, 4 p, 7 d, 8 f , 6 g, 6 h, and 5 i for U.
For dynamically correlated calculations, core orbitals were frozen to reduce the cost of the calculation. In some cases, the higher-energy virtual orbitals resulting from basis set decontraction were deleted as well; the energy cutoffs are specified in the relevant results section. A vertical shift of 0.2 was used in four-component XMS-CASPT2 and NEVPT2 calculations, with correction terms included in all elements of the effective Hamiltonian. 37,61 For CASPT2 + SI-SO, an imaginary shift of 0.2 was used along with the standard 0.25 IPEA shift. 62,63 All computations were performed at the experimental structures. 26,28,29,64 The bond lengths (for diatomics) and molecular coordinates (for larger molecules) are provided in the Supporting Information. For the larger molecules, heavy atom positions were obtained from the crystal structure, and hydrogen positions were optimized using DFT. The DFT optimizations were performed using the Turbomole package. 65,66 We used the PBE functional 67 along with the def2-TZVPP basis set for all atoms except U, for which we used def-TZVPP with effective core potential. [68][69][70][71][72] We also used the resolution of the identity approximation along with the corresponding RI-J auxiliary basis sets. 73,74 IV. RESULTS AND DISCUSSION

A. Benchmarking Dynamical Correlation Methods with
Chalcogen Diatomics Table I presents ZFS results for a series of chalcogen diatomics using the methods reported in the previous sections. The diatomics in this series have been previously used to benchmark other methods by Rota et al. 78 For dynamically correlated calculations, deep core orbitals were frozen (2 for S 2 , 10 for Se 2 , and 18 for Te 2 ). As indicated in the table, an additional 18 semi-core orbitals were frozen in some Te 2 calculations to reduce the cost. These benchmark tests were performed with the MS-MR parametrization.
We have also tested the influence of basis set contraction upon ZFS results, because ANO-RCC basis contractions are not optimized for the four-component treatment of relativistic effects. Due to cost constraints, MRCI+Q results are included only for the contracted basis set, and in some cases high-energy virtual orbitals were deleted. In most cases, the results with the contracted basis set are further from the experimental values, with errors of 10-20% in the cases of S 2 and Se 2 , and smaller errors for O 2 and Te 2 . This suggests that in the pursuit of high-quality results without the large costs associated with decontraction, there is use for further development of basis sets suitable for four-component relativistic calculations. For the time being, we prefer using uncontracted basis sets and compromise when necessary by contracting basis sets on ligands but not spin centers.
The XMS-CASPT2 results in Table I show an extremely strong dependence upon the choice of states to include in the calculation, with large errors introduced when excluding some low-lying excited singlets. The choice of which excited states to include affects the calculation primarily in two ways: through the orbital optimization in state-averaged CASSCF, and through the zeroth-order Hamiltonian which uses a stateaveraged density matrix. Table S2 in the Supporting Information separates these effects by correlating the first three or six states using both sets of orbitals, and clearly shows that the zeroth-order Hamiltonian is the main source of error; the change to the state-averaged orbitals alters the XMS-CASPT2 splitting by a few percent only. It is unsurprising, then, that the strong dependence on state averaging is not observed in the other methods; Dyall's Hamiltonian used by NEVPT2 replaces the active-active part of the state-averaged Fock matrix with the full Hamiltonian, and MRCI+Q does not use it at all. MRCI+Q is unsurprisingly the most reliable method tested here, although for large problems it will be intractable. When a more approximate dynamically correlated method is needed, NEVPT2 is preferred due to its relatively low sensitivity to the set of states included, although for this series the accuracy of CASSCF is competitive with the perturbative methods, in comparison to experimental values.
We briefly mention one other complication which has been observed. For calculations with a very minimal active space (e.g., 2 electrons in 2 orbitals with all six states included), it is possible for states which are nondegenerate underĤ to become degenerate underĤ 0 and, therefore, to freely mix in the XMS rotation. In such cases, if a vertical shift is used with only the diagonal correction terms, unphysical breaking of Kramers degeneracy will result, but if no shift is used or we correct the entire effective Hamiltonian, degeneracy is restored in the final multi-state mixing step.

B. Zero-Field Splitting of Pseudotetrahedral Cobalt(II) (XPh) 4 Complexes
To gauge the performance of relativistic correlation methods for mononuclear transition metal complexes, we look to the series of pseudotetrahedral cobalt(II) (XPh) 4 complexes with X = O, S, Se whose structures are shown in Figure 1. Molecule (c), (PPh 4 ) 2 Co(SPh) 4 was one of the first field-free transition metal single-ion magnets to be discovered; 25 the analogous complexes here have been characterized both experimentally and theoretically to probe the magneto-structural correlations present. [26][27][28] Structures (a) and (b) as well as (c) and (d) differ only in the counter ion and solvent present in the crystal, which dramatically alters both the conformation and ZFS of the anion.
CASSCF-level ZFS parameters computed for the whole series are compared along with experimental values in Table II. A mixed basis set is used, with the central cobalt atom decontracted in four-component results, the TZP contraction for chalcogen atoms, the double-ζ plus polarization (DZP) contraction for carbon atoms, and double-ζ (DZ) for hydrogen atoms. Agreement with experiment is reasonably good for this series, with qualitative agreement throughout. The major exception is the Dirac-CASSCF results with only the ground quartet included in the state-averaging procedure; in most cases more accurate results are obtained when the orbitals are balanced by including all d 5 excitations (10 quartets and 40 doublets). The choice of whether to include a second, correlating set of d-orbitals is of lesser importance, typically altering the observed splitting by a few wavenumbers or less.
Four-component NEVPT2 calculations have also been performed on structures (a), (c), and (e) of this series; the results are summarized in Table III. Due to cost constraints, it was necessary to reduce the basis set combination used: for these data, the basis sets for cobalt and hydrogen are the same as above, but the DZP contraction was used for chalcogen atoms and the DZ contraction for carbon atoms. All results in this table use an active space of 7 electrons in 5 active orbitals. For Dirac-NEVPT2, the MS-MR contraction is used while cor- relating the four lowest states, but using orbitals that were obtained by finding a minimax state-averaged energy for all excitations of d-electrons. The number of frozen core orbitals is 37, 53, and 69 for structures (a), (c), and (e), respectively. 89 high-energy virtual orbitals were deleted, corresponding to an energy cutoff threshold of 11 E h to avoid splitting up neardegenerate orbitals at 10 E h . We used the MS-MR internal contraction scheme in this and the previous section. It is tempting to switch to the SS-SR parametrization, which would reduce the memory cost by roughly a factor of N st and the timing cost by a greater amount. Some tests have been performed using truncated models for this molecular series; the details are presented in the Supporting Information. We found that reducing the expansion set to the SS-SR scheme introduced errors in the splitting energies of several percent, and that the ratio |E/D| was much more sensitive. Therefore, we prefer the MS-MR parametrization when possible, while noting that the errors introduced by SS-SR in the energy spacings were moderate. We see corrections due to dynamical correlation of 2 to 20 cm −1 , slightly larger than are observed in the CASPT2 + SI-SO reference data. The basis set reduction necessary to make the NEVPT2 calculations tractable introduced errors of up to 3.6 cm −1 at the Dirac-CASSCF level, showing somewhat higher basis set sensitivity than the SI-SO results; one might expect a somewhat greater basis set incompleteness error in the correlated data, which would be consistent with the results in Table I for smaller test systems. For structure (a), the Dirac-NEVPT2 best reproduces the experimental findings, but for the other two complexes it appears to overcorrect.

C. Zero-Field Splitting of a Uranium(III) Pyrazole Complex
To test the performance of these methods for magnetic complexes with very heavy elements, we have computed the ZFS of the ground J = 9 2 multiplet in the uranium(III) dihydrobispyrazolylborate complex U(H 2 BPz 2 ) 3 , shown by Long and co-workers to act as a single-ion magnet in magnetically dilute samples. 29,30 The structure is depicted in Figure 2, and results are presented in Table IV. Previous theoretical calculations have been carried out with the SI-SO approach by Spi-    vak et al. 32 and with a semi-empirical crystal-field approach by Baldovi et al. 31 The results reported by Spivak et al. show a somewhat smaller splitting than that in Table IV, which might result from a different choice of molecular geometry; further details are discussed in the Supporting Information.
The results in Table IV were computed with either the decontracted basis or TZP contraction for the uranium atom as appropriate, the DZP basis set for boron, carbon, and nitrogen atoms, and the DZ contraction for hydrogen atoms. An exception is the Dirac-NEVPT2 result, for which the polarization functions were not used for main group elements; one column also shows CASSCF results with the same basis set. The active space in all cases is 3 active electrons in the 7 forbitals. Orbital optimizations were perfomed for the minimax state-averaged energy for all excitations within the forbitals. NEVPT2 calculations included the lowest 10 states and used the SS-SR parametrization; 72 core orbitals were frozen, and 160 virtual orbitals were deleted, corresponding to an energy cutoff of 10 E h .
We computed an excitation energy of 161.2 cm −1 with Dirac-CASSCF and 103.3 cm −1 with NEVPT2, somewhat lower than the value of 230 cm −1 predicted by a semiempirical crystal-field approach. 31 We obtain qualitatively similar results with Dirac-CASSCF vs. CASSCF + SI-SO. The Breit correction at the CASSCF level was determined to be only a few wavenumbers, smaller than other errors relating to electron correlation and the molecular geometry dependence.
The excitation spectrum of U(H 2 BPz 2 ) 3 is not available for comparison. The closest analogue we have is a measurement of the magnetic relaxation barrier, which was found to be 16 cm −1 in a magnetically dilute sample. 30 All theoretical meth-FIG. 2. Geometry of U(H 2 BPz 2 ) 3 . 29 ods examined here give a first excited state at least several times higher. There are several competing relaxation mechanisms in this molecule, 30 so the observed barrier is not directly related to the energy splitting computed in this work.
The CASPT2 + SI-SO approach requires computing the second-order correlation energy for a large number of excited states, and there are practical limits to how many states can be included. We have truncated our CASPT2 space at 21 quartets and 11 doublets, out of 35 quartets and 112 doublets accessible by exciting within the f -orbital manifold, which is the same cutoff used in the earlier study. 32 As shown in Table IV, truncating the coupling of CASSCF states at the same level shifts the energy splittings by 70-160 cm −1 , so the aggressive truncation needed for large active spaces comes at a substantial accuracy cost for ZFS calculations. One of the advantages of the four-component approach is that spin-orbit effects are included from the start and can be captured without computing many excited states, so it avoids this particular problem. Of course, there is a tradeoff; the large memory costs associated with the four-component NEVPT2 required us to make compromises to limit the number of orbitals included in the correlation, although in this case the errors associated with reducing the basis set are manageable, at least at the CASSCFlevel. For this molecule, the Dirac-NEVPT2 correction to the enery splittings is of comparable magnitude but opposite sign relative to the CASPT2 + SI-SO reference.

V. CONCLUSIONS
We have developed a computational tool for the determination of ZFS parameters in relatively large molecules containing strong relativistic effects. Specifically, we report the development of (extended) multistate implementations of four-component relativistic CASPT2, NEVPT2, and MRCI as well as a mapping from the four-component wavefunctions to a pseudospin Hamiltonian for the extraction of ZFS parameters. We observed a surprisingly large dependence upon the state-averaging procedure used to obtain the orbitals, which is troublesome for XMS-CASPT2 but reduced in CASSCF, NEVPT2, and MRCI+Q. Benchmark calculations performed on the series of chalcogen diatomics and pseudotetrahedral Co(XPh) 4 complexes showed similar results to the less costly scalar relativistic + SI-SO approach, suggesting that for light elements the fully relativistic methods can serve as a useful benchmark but are not always needed. Tests on the actinide single-molecule magnet U(H 2 BPz 2 ) 3 showed larger effects, indicating that for very heavy metals the simultaneous treatment of relativity and electron correlation, as well as avoidance of state truncation sometimes required for SI-SO, can be helpful when attempting high-accuracy calculations. The methods reported in this manuscript are publicly available in the BAGEL program package 49 and distributed under the GNU General Public License.

A. Supporting Information
The following content is available in the Supporting Information: Discussion of the theory and working equations used in the implementation of the mapping to the pseudospin Hamiltonian; ZFS results for chalcogen diatomics with a reduced active space, including an investigation to separate the effects of state-averaged orbitals and state-averaged Fock operator on dynamically correlated results; comparison of ZFS results using the SS-SR and MS-MR parametrizations for truncated models of the Co(XPh) 2− 4 series; discussion of the molecular geometry of U(H 2 BPz 2 ) 3 and the effects of optimization of hydrogen positions; xyz coordinates of the molecules used for testing.