Toward Accurate Spin–Orbit Splittings from Relativistic Multireference Electronic Structure Theory

Most nonrelativistic electron correlation methods can be adapted to account for relativistic effects, as long as the relativistic molecular spinor integrals are available, from either a four-, two-, or one-component mean-field calculation. However, relativistic multireference correlation methods remain a relatively unexplored area, with mixed evidence regarding the improvements brought by perturbative treatments. We report, for the first time, the implementation of state-averaged four-component relativistic multireference perturbation theories to second and third order based on the driven similarity renormalization group (DSRG). With our methods, named 4c-SA-DSRG-MRPT2 and 3, we find that the dynamical correlation included on top of 4c-CASSCF references can significantly improve the spin–orbit splittings in p-block elements and potential energy surfaces when compared to 4c-CASSCF and 4c-CASPT2 results. We further show that 4c-DSRG-MRPT2 and 3 are applicable to these systems over a wide range of the flow parameter, with systematic improvement from second to third order in terms of both improved error statistics and reduced sensitivity with respect to the flow parameter.

The accurate description of the electronic structure of systems containing heavy atoms requires a balanced treatment of both strong electron correlation and relativistic effects. 1,2e former arises from the commonplace occurrence of near-degenerate configurations in such systems, and the latter arises from the Z 2 dependence of relativistic contribution to the the valence electron energy, 3 where Z is the atomic number.1][12] Four-component (4c) mean-field and electron correlation methods employ Hamiltonians that explicitly couple the spin and orbital angular momentum degrees of freedom, and are the methods of choice to treat SOC effects in an ab initio manner.However, relativistic electron correlation methods are not always able to consistently improve the description of SOC effects from relativistic mean-field references. 135][16][17][18][19][20][21][22][23][24][25][26][27][28][29][30][31][32][33] Most nonrelativistic correlation methods can be readily adapted to account for relativistic effects within the no-pair approximation, as long as the complex-valued molecular spinor integrals are available, [34][35][36] after either a two-or four-component mean-field calculation.As such, considerable work has been devoted to developing mean-field methods.8][49][50][51] Multireference relativistic methods have received proportionately more attention than their nonrelativistic cousins, as many systems where a relativistic study is warranted, such as late-row transition metal, lanthanide, and actinide complexes, exhibit strong configuration mixing in the ground state, and as a result require the zeroth order wave function to be multiconfigurational.Examples of relativistic multireference theories include Fock-space multireference coupled cluster (FSMRCC), 1 4c internally contracted MRCI (ic-MRCI), CASPT2, 2 NEVPT2, 52 generalized van Vleck PT2 (GVVPT2), 53 multireference Møller-Plesset (MRMP). 1,54In several ways, these methods reflect one or more limitations of their nonrelativistic counterparts: 1) lack of proper scaling with system size, 2) difficulties with scaling to large active spaces, 3) numerical divergences arising from the intruder state problem, and 4) insufficient accuracy in the description of dynamical electron correlation.
To address the abovementioned insufficiencies of relativistic multireference theories, we have explored a 4c generalization of the multi-reference driven similarity renormalization group (MR-DSRG) approach. 55,56The MR-DSRG is a size-extensive, low-scaling, numerically robust, and systematically improvable many-body formalism for treating dynamical electron correlation starting from strongly correlated reference states. 55,56In this paper, we show that a combination of 4c strongly correlated reference states and the MR-DSRG is a very promising route to systematically account for relativistic and correlation effects.We demonstrate this point by reporting the first implementation of 4c MR-DSRG methods truncated to second-and third-order in perturbation theory (4c-DSRG-MRPT2/3) and benchmarking them on the spin-orbit splittings of second-to fourth-row atoms, which previous work has shown to be a weak point for four-component correlation methods. 13e starting point for formulating the 4c-MR-DSRG is first-quantized quasi-relativistic many-electron Hamiltonian: 57 where h D (i) is the one-electron Dirac Hamiltonian for electron i, and g(i, j) is the two-electron interaction operator for an electron pair i and j.The latter can contain either the bare-Coulomb operator, leading to the Dirac-Coulomb (DC) Hamiltonian, or can additionally include the Gaunt and gauge terms, leading to the Dirac-Coulomb-Gaunt (DCG) or Dirac-Coulomb-Breit (DCB) Hamiltonians.These operators are defined as where α i is a 3-vector of the Dirac matrices for particle i.The Coulomb operator contains charge-charge and spin-same-orbit interactions, the Gaunt and Breit operators contain additional current-current and spin-other-orbit interactions, and a gauge term resulting from working in the Coulomb gauge. 34,36Following common practice, we use the DCB Hamiltonian in combination with the no-pair approximation, whereby the many-body basis is constructed from positive-energy spinors only. 57The four-component molecular spinors {ψ i (r)} are expanded in a two-spinor basis set {ϕ S/L µ } for the small/large components.To avoid variational collapse, we impose the kinetic balance condition, 2cϕ S µ = σ • pϕ L µ , where σ is a 3-vector of Pauli matrices. 58The large component 2-spinor basis functions comprise standard scalar Gaussian basis functions multiplied by 2-spinor spherical basis functions. 57 consider the most general, state-averaged formulation of the MR-DSRG approach, which can target one or multiple electronic states in one computation.We approximate the k-th electronic state to zeroth-order with 4c-CASSCF states of the form where N CAS is the number of CAS determinants.In 4c-CASSCF, the spinor orbitals {ψ i } that enter into Slater determinants {Φ µ } and the complex CASCI coefficients C (k) µ are simultaneously optimized via the minimax principle. 59No positronic orbitals enter into the determinants, in what is known as the "no dressed pair" or "no virtual pair" approximation (NVPA). 40,60Furthermore, to ensure the correct degeneracies of the electronic states, we use the state-averaged CASSCF formalism (4c-SA-CASSCF), 61 where the average energy of a chosen set of states is optimized with respect to the variational parameters. 62To account for dynamical electron correlation, we combine 4c-(SA-)CASSCF wave functions with the multireference driven similarity renormalization group (MR-DSRG) formalism. 55The starting point of the 4c-MR-DSRG formalism is the second-quantized version of the quasi-relativistic no-pair Hamiltonian: The orbital indexing convention and definition of orbital spaces used in this work reflects the partitioning of the orbitals into core, active, and virtual sets, as shown in Fig. 2. The quantities f q p that enter in the Hamiltonian are elements of the generalized Fock matrix, where γ i j is the one-body density matrix defined for a general state Ψ as γ i j = ⟨Ψ|â † i âj |Ψ⟩.The MR-DSRG performs a continuous unitary transformation of the Hamiltonian: Ĥ → H = e − Â Ĥe Â. ( The operator Â is anti-Hermitian, and is expressed in terms of a cluster operator as Â = T − T † .T is parameterized as in the internally-contracted generalization of coupled cluster theory [63][64][65] as T = T1 + . . .Tn , where a general k-body term is expressed as where {•} indicates normal-ordering with respect to a correlated reference state Ψ or ensemble density. 66The effect of the MR-DSRG transformation is to bring the Hamiltonian into a more band-diagonal form by zeroing the couplings between the multi-configurational reference |Ψ 0 ⟩ and the excited configurations, i.e., the set of internally-contracted determinants |Ψ ab... ij... ⟩ = {â ab... ij... } |Ψ 0 ⟩. 66,67To avoid numerical divergences caused by intruder states, in the DSRG we solve a regularized equation that enforces the decoupling of the off-diagonal elements of H, denoted as [ H] od .This equation reads as [ H(s)] od = R(s), where R(s) is a regularization term derived by a second-order perturbative analysis of the flow similarity renormalization group. 55,68,69Several authors have pointed out the benefits of introducing a parameterized denominator shift or a regularizer of the first-order amplitudes in second-order perturbation theory-both in its single and multireference versions. 70,71In the DSRG, this role is played by the flow parameter s, which in typical applications is found to be optimal for values in the range s ∈ [0.5, 2], depending on the truncation level. 72To generate equations and corresponding tensor contractions for the 4c relativistic extension of the MR-DSRG truncated to second-and third-order in perturbation theory, we have developed the automated code implementation pipeline depicted in Fig. 3.This approach enables the rapid implementation of highly complex electronic structure theories, removing the potential for human error.This workflow uses our code generator Wick&d 73 coupled to relativistic integrals obtained from PySCF. 44,45,74Since 4c-CASSCF is currently not available in PySCF, we generated the corresponding spinor coefficients using the ChronusQ package. 43Computational details are found in Section 1 of the Supplementary Information.To examine the ability of the 4c-DSRG-MRPT2/3 to improve systematically upon 4c-CASSCF, we computed spin-orbit splittings of second-to fourth-row p-block elements, and the spectroscopic constants for the hydroxyl radical.The p-block splittings span three orders of magnitude, from a few reciprocal centimeters (cm −1 ) to a few thousand, as shown in Fig. 4.
Therefore, they form a good set of benchmarks for both the absolute and relative accuracies of relativistic methods.
In Fig. 5, we compare error statistics for the splittings for the 15 p-block elements calculated with 4c-CASSCF, 4c-SA-DSRG-MRPT2/3, 4c-CASPT2, and 4c-MR-CISD+Q (taken  from the work by Zhang et al. 13 ) to the experimental splittings for three values of the flow parameter (0.2, 0.35, and 0.5 E −2 h ), which cover the range of flow parameters used in non-relativistic computations and that produce results competitive with 4c CASSCF and CASPT2.Table S1 reports splittings for each state with optimal flow parameters that achieve the smallest mean absolute error in cm −1 for the whole set of elements for each of the methods respectively (0.24 E −2 h for MRPT2 and 0.35 E −2 h for MRPT3).Looking at the distribution of errors (top plot in Fig. 5), both MRPT2 and 3 achieve better average performance on the test set of atoms than CASSCF and the two other multireference electron correlation methods (except for s = 0.5 E −2 h in MRPT2).Improving the treatment of dynamical electron correlation from second to third order results in a systematic reduction of the error in the splittings and a decrease in flow parameter sensitivity.Across all methods, the main outliers in percentage error (bottom plot in Fig. 5) are elements with small (10-150 cm −1 ) splittings, with the exceptions of arsenic (322 cm −1 ).The complete dataset for the splittings can be found in Sections 2 and 3 in the Supplementary Information.When compared according to signed errors, irrespective of the value of s, the main outliers are fourth-row elements with large splittings, with selenium and gallium being the worst offenders.We also note that DSRG-MRPT3 achieves sub-wavenumber accuracy for boron and carbon, and rather impressively, less than 2 cm −1 errors for the fourth-row selenium and bromine.4c-CASPT2 and 4c-MR-CISD+Q can be seen to consistently and severely underestimate spin-orbit splittings.This behavior was also pointed out in Zhang et al. 13 to occur for smaller and larger basis sets than the ones used here.In Table S2, we report additional splittings represented by dashed arrows in Fig. 4. We again observe that MRPT2/3 outperform CASPT2, which in turn outperforms CASSCF in this set of transitions.This shows that the 4c-SA-DSRG formalism is capable of accurately capturing SOC effects in higher excited states.
After assessing the quality of DSRG-based methods, we turn to the question of the dependence of our results on the flow parameter and its optimal choice for spin-orbit coupling computations.To this end, we studied the variation of the mean absolute percentage error in the splittings as the flow parameter is varied in the range of s ∈ [0, 1.3] E −2 h .This quantity is plotted in Fig. 6 for both the DSRG-MRPT2/3.In the middle inset, we first observe the correct limiting behavior for s = 0, where MRPT2 and MRPT3 reduce to CASSCF.For  Sections 5 and 6 in the Supplementary Information report summary error statistics broken down into groups and periods.An analysis of this data shows that second-row elements display significantly larger absolute and percentage errors, which cannot be simply attributed to the elements having smaller splittings.This is most likely due to the poor description of the core region by the relatively small uncontracted cc-pVTZ basis set for boron through fluorine.As SOC effects decay as 1/r 3 , 78 having variational flexibility in the core region is crucial: Visscher and Dyall 79 found that the cc-pVTZ basis augmented with an additional tight p function reduced the error in the spin-orbit splitting of F 2 tenfold to 1%.We have tested this hypothesis by following Visscher and Dyall and added a tight p function with the exponent of 250.83491 to the cc-pVTZ basis of fluorine.The resulting 2 P3 /2 → 2 P1 /2 splitting error of fluorine for MRPT3 using s = 0.35 E −2 h improved from −12.46 cm −1 (−3.1%) to −1.04 cm −1 (−0.3%).Analogous data for the remaining second-row elements is reported in Section 8 of the Supplementary Information.From these results, we can see that except for the N 2 D5 /2 → 2 D3 /2 splitting, the provision of more core flexibility significantly improves the description of second-row valence SOC effects with a negligible increase in basis size.The improvement in the underlying CASSCF splittings largely drives the improved description.
However, DSRG-MRPT2 and 3 still provide an improvement over the CASSCF results.
Another aspect we investigated is a reduction of the cost of the 4c-CASSCF procedure, which in our computations, especially on third-and fourth-row elements, is overwhelmingly the most time-consuming step.In Section 7 of the Supplementary Information, we report additional data for a mixed scheme aimed at reducing this cost that uses a combination of the DC and DCB Hamiltonians in the 4c-CASSCF and DSRG-MRPT2 computations.The data presented therein supports the use of such a mixed scheme, as the reduction in accuracy is negligible.
Lastly, we consider the accuracy of potential energy surfaces computed with the 4c-DSRG-MRPT methods.For example, in Table 1 we report the spectroscopic constants of the ground state of the hydroxyl radical computed with 4c-CASSCF and 4c-DSRG-MRPT3.
1][82] Our results show that although both 4c-CASSCF and 4c-DSRG-MRPT predict ZFSs in excellent agreement with the experimental value, 83 the spectroscopic constant from the 4c-DSRG-MRPT3 potential display significantly smaller errors.This point is particularly evident for the harmonic vibrational frequency of OH, which 4c-CASSCF underestimates by 210.9 cm −1 vs. 8.7 cm −1 for the 4c-DSRG-MRPT3.Finally, we computed the potential energy surfaces of the X 2 Π 3/2 and 2 Π 1/2 states.
Table 1: Spectroscopic constants of the X 2 Π ground state of the OH radical calculated with different methods.Shown are the differences with respect to experimental values, taken from the Diatomic Molecular Spectroscopy Database. 83,84he (adiabatic) zero-field splitting includes anharmonic zero-point vibrational energy contributions.

Method
∆r e (pm) ∆ω e (cm −1 ) ∆ω e x e (cm −1 ) ∆D 0 (eV) ∆ZFS (cm In summary, we have presented second-and third-order four-component multi-reference perturbation theories based on the driven similarity renormalization group formalism, 4c-DSRG-MRPT2/3.We benchmarked these methods on the spin-orbit splittings of second-to fourth-row p-block atoms, showing that they outperform both the underlying 4c-CASSCF and other four-component electron correlation methods, namely 4c-MR-CISD+Q and 4c-CASPT2.We have further shown that 4c-DSRG-MRPT2 and 3 are applicable to these systems over a wide range of the flow parameter, with systematic improvements in error metrics and sensitivity with respect to s from second to third order.In our calculations, most of the wall time is spent in the integral transformation step, which is known to be the main drawback of four-component theories.Exact two-component (X2C) theories with atomic mean-field spin-orbit effects (amfX2C) [85][86][87][88]  For the p-block elements, whenever the neutral atom cannot be converged by the 4c-HF procedure, the corresponding ion with the configuration of the closest noble gas is used instead (e.g., Al 3+ , Se For the calculations involving the hydroxyl radical, a complete active space of 5 electrons in 8 spinor orbitals is used, which corresponds to the 2p σ , the doubly degenerate 2p π , and the 2p * σ MOs.The state-averaged calculations averages over the lowest 6 states with equal weights, corresponding to the doubly degenerate X 2 Π 3/2 , 2 Π 1/2 , and A 2 Σ + states.The default convergence criteria from Chronus Quantum are used.For all DSRG-MRPT2/3 computations, a flow parameter of 0.5 E −2 h is used.The spectroscopic constants are obtained from the psi4.diatomic.anharmonicityfunction provided by the open-source package Psi4, 5 which is an implementation of the algorithm described by Bender et al. 6 .A grid size of 0.005 Å is used in the range of [0.920, 1.020] Å, with a finer grid of 0.001 Å in the range of [0.960, 0.990], for a total of 45 points, to determine the constants.The dissociation energies, D 0 , are calculated using state-specific CASSCF references, as state-averaged CASSCF calculations using the same settings are not feasible due to multiple degenerate states coming in from higher (untracked) excited states. 7
2 The 4c-MR-CISD+Q computations were intractable for these atoms. 3The uncontracted cc-pVDZ basis set was used for these atoms.
4 Unavailable data points have been omitted from averaging.
TABLE S1.Comparison between the spin-orbit splittings of the 15 second-to fourth-row p-block elements calculated with 4c-SA-CASSCF, 4c-CASPT2, 4c-MR-CISD+Q, and MRPT2 and 3 to the experimental splittings.Flow parameters, s, are in units of E −2 h , and all results are reported in units of cm −1 , unless otherwise noted."Exp."stands for experimental splitting, "MAE" stands for mean absolute error, and "MAPE" stands for mean absolute percentage error. 1Two core spinors were frozen for nitrogen.
2 The 4c-MR-CISD+Q computations were intractable for these atoms. 3The uncontracted cc-pVDZ basis set was used for these atoms. 4As the most error-prone data points were intractable for MRCI, the average is not appropriate for comparison with other methods.

VII. RESULTS FOR CALCULATIONS USING DC COEFFICIENTS
We have experimented with using only the DC Hamiltonian for the 4c-CASSCF stage, using the resulting coefficients to transform the DCB atomic spinor integrals, which are then fed into 4c-DSRG-MRPT2.We can see that the performance is very similar to those using full DCB molecular spinor coefficients.This demonstrates that iterative reference relaxation can largely correct for inexact starting CI and MO coefficients, and can therefore be used to bypass expensive iterative DCB integral transformations in the SCF stages.This is further supported by the fact that, if the same MO coefficients used in the AO to MO transformation using DC atomic spinor integrals ("0.

VIII. RESULTS FOR USING THE AUGMENTED CC-PVTZ BASIS
We also augmented other second row elements with a p function each with exponents determined by scaling the fluorine exponent by the ratio between the exponents of the tightest p functions in the original basis set of a given element and fluorine, resulting in exponents of 68.8243, 106.9535, 152.2273, and 196.9866S8.Percentage error improvements in all 8 calculated spin-orbit splittings for second-row elements after augmenting the cc-pVTZ basis with a tight p function.A positive percentage signifies that the result after the augmentation is closer to the experimental splitting, and vice versa for a negative percentage.Flow parameters are in units of E −2 h , "Avg." stands for "average", "MAPE" stands for mean absolute percentage error, "orig."stands for "original basis", and "aug."stands for "augmented basis".All results are in percentage points.

Figure 1 :
Figure 1: Summary of the 4c-SA-DSRG-MRPT2/3 algorithm.The iterative reference relaxation process is terminated once a threshold is reached.

Figure 2 :
Figure 2: The partitioning of the orbital basis, and the orbital indices used in the MR-DSRG formalism.

Figure 3 :
Figure 3: Overview of the automated strategy used in this paper to implement the 4c DSRG methods.Starting from the left, the DSRG ansatz is perturbatively expanded to give working equations for the energy corrections.These are then turned into tensor contractions using the open-source Wick&d package. 73Contractions generated by Wick&d are then automatically compiled into executable Python code.

Figure 4 :
Figure4: A scheme of the spin-orbit splittings of the p-block elements discussed.Stateaveraging was performed over all states shown.The greyed-out splittings are excluded from analysis and are only included in the state-averaging procedure to ensure variational stability.A switch between the 2 D 3/2 and 2 D 5/2 states in the group 15 elements is indicated by negative splittings.Energy levels accessed from the NIST Atomic Spectra Database.75

Figure 5 :
Figure5: Distributions of the errors in the computed spin-orbit splittings.The dashed vertical lines at 0 represent experimental splittings.The plots were generated with the violinplot function of seaborn,76,77 where the data points and the kernel density estimates of the splitting errors of each method are shown.The outlier elements of each plot are labeled.
both PT2 and PT3, the mean absolute percentage errors are smaller than the CASSCF and CASPT2 values for a large range of s, up to s ≈ 0.75 E −2 h for PT3.For PT2, and PT3 especially, these values overlap with the commonly used range s ∈ [0.35, 1.0].The kinks in the main curves are artifacts of the absolute values, as the mean signed error curves (bottom insets) are smooth.These results show that the error profiles of MRPT2 and MRPT3 do not simply interpolate between those of CASSCF and CASPT2, and bring systematic improvements from a balanced treatment of dynamical correlation.

Figure 6 :
Figure 6: Variation of the mean average error of 4c-DSRG-MRPT2 and MRPT3 as a function of the flow parameter, s, for all of the splittings computed using the DCB Hamiltonian with DCB CASSCF molecular spinor coefficients.The main curves show absolute errors, and the bottom insets show signed errors.The middle insets zoom into the region close to s = 0.
are an attractive way to reduce the computational effort not only of integral transformations but also in the SCF iterations themselves, depending on the flavor of X2C used.In the future, we plan to explore the combination of amfX2C with two-component MR-DSRG methods to extend multireference relativistic computations to larger systems.Overall, these preliminary results show that a combination of relativistic Hamiltonians and multireference DSRG methods could open many new avenues in modeling the chemistry of open-shell species containing transition metals and heavy elements.Supplementary information: "Towards accurate spin-orbit splittings from relativistic multireference electronic structure theory" Zijun Zhao 1, a) and Francesco A. Evangelista 1, b) Department of Chemistry and Cherry Emerson Center for Scientific Computation, Emory University, Atlanta, Georgia 30322, USA (Dated: 13 May 2024) a) Electronic mail: zijun.zhao@emory.edub) Electronic mail: francesco.evangelista@emory.edu 1 arXiv:2405.06077v1[cond-mat.str-el]9 May 2024 I. COMPUTATIONAL DETAILS We used the Dirac-Coulomb-Breit (DCB) Hamiltonian with the uncontracted Dunning's correlation-consistent polarized valence triple-zeta (uc-cc-pVTZ) basis set, 1-3 finite Gaussian nuclei, 4 freezing non-valence orbitals in all third-and fourth-row atoms from the correlated calculations.

TABLE S2 .
Comparison between additional spin-orbit splittings (those indicated with dashed arrows in ??) of the 9 second-to fourth-row p-block elements calculated with 4c-SA-CASSCF, 4c-CASPT2, 4c-MR-CISD+Q, and MRPT2 and 3 to the experimental splittings.Flow parameters, s, are in units of E −2 h , and all results are reported in units of cm −1 , unless otherwise noted."Exp."stands for experimental splitting, and "MAPE" stands for mean absolute percentage error.

TABLE S4 .
Summary error statistics for selected four-component methods across groups 13 through 17.The mean absolute errors (abs.) are in cm −1 , and the mean absolute percentage errors (pct.) are in percentage points.
DCDC') to the experimental splittings.All calculations employed a flow parameter, s, of 0.2 E −2 h , and all results are reported in units of cm −1 , unless otherwise noted.'Exp.' stands for experimental splitting, 'MAE' stands for mean absolute error, and 'MAPE' stands for mean absolute percentage error.
for boron, carbon, nitrogen, and oxygen respectively.

TABLE S7 .
Splittings in all 8 calculated spin-orbit splittings for second-row elements before and after augmenting the cc-pVTZ basis with a tight p function.Flow parameters are in units of E −2 h , 'orig.' stands for 'original basis', and 'aug.' stands for 'augmented basis'.All results are in cm −1 .