On the performance of DFT/MRCI-R and MR-MP2 in spin–orbit coupling calculations on diatomics and polyatomic organic molecules

ABSTRACT We have investigated the performance of different multi-reference quantum chemical methods with regard to electronic excitation energies and spin–orbit matrix elements (SOMES). Among these methods are two variants of the combined density functional theory and multi-reference configuration interaction method (DFT/MRCI and DFT/MRCI-R) and a multi-reference second-order Møller–Plesset perturbation theory (MR-MP2) approach. Two variants of MR-MP2 have been tested based on either Hartree–Fock orbitals or Kohn–Sham orbitals of the BH-LYP density functional. In connection with the MR-MP2 approaches, the first-order perturbed wave functions have been employed in the evaluation of spin–orbit coupling. To validate our results, we assembled experimental excitation energies and SOMES of eight diatomic and fifteen polyatomic molecules. For some of the smaller molecules, we carried out calculations at the complete active space self-consistent field (CASSCF) level to obtain SOMEs to compare with. Excitation energies of the experimentally unknown states were assessed with respect to second-order perturbation theory corrected (CASPT2) values where available. Overall, we find a very satisfactory agreement between the excitation energies and the SOMEs obtained with the four approaches. For a few states, outliers with regard to the excitation energies and/or SOMEs are observed. These outliers are carefully analysed and traced back to the wave function composition.


Introduction
The combined density functional theory and multireference configuration interaction (DFT/MRCI) method by Grimme and Waletzke [1] has proven itself as a very useful approach to an efficient treatment of static and dynamic correlation in larger molecules. In combination with Spock [2][3][4], which allows for the inclusion of electronic spin-orbit coupling, and Vibes [5][6][7], capable of treating the vibronic part of the problem, it has been used for the successful elucidation of many spin-dependent properties in molecules CONTACT Christel M. Marian Christel.Marian@hhu.de Supplemental data for this article can be accessed at http://dx.doi.org/./... [2,[8][9][10][11][12][13][14][15]. However, the DFT/MRCI method is not without shortcomings. The impetus for the present study came when we realised that the original DFT/MRCI method has some difficulties providing a balanced description of configurations with four open shells. Obviously, this is the case in electronically excited dimers and other weakly coupled bi-chromophores [16]. Less obvious are cases such as nitrobenzene where artificially low-lying excited singlet states with four open shells are encountered in the DFT/MRCI spectrum [17].
Recently, some of these shortcomings of the original DFT/MRCI were remedied by the parameterisation of a redesigned DFT/MRCI-R Hamiltonian [16], dubbed DFT/MRCI-R in the following. A completely different hybrid multi-configuration wave function and density functional theory (DFT) ansatz has been presented by Jensen and coworkers [18][19][20][21][22]. It combines a long-range multi-configuration self-consistent field (MCSCF) treatment with an adiabatic short-range DFT description of electron correlation. Time-dependent MCSCF shortrange DFT (TD-MC-srDFT) has been implemented in the Dalton program package [23] and might, therefore, be used for computing spin-dependent properties. To our knowledge, spin-orbit coupling calculations based on this method have not yet been reported, however.
Alternative ab initio methods that are able to treat static and dynamic electron correlation in sizeable molecules and allow for the computation of spinorbit matrix elements (SOMEs) are sparse. Spin-orbit coupling constants have been evaluated in a multiconfiguration linear response approach a long time ago [24,25], but the amount of dynamic electron correlation that can the covered by MCSCF is very limited. State-averaged complete active space self-consistent field (SA-CASSCF) calculations followed by second-order perturbation theory (CASPT2) have become something like a gold standard for computing excited-state energies when the active space can be made large enough to include all important excitations [26][27][28]. Spindependent properties are then evaluated by the so-called CASPT2/CASSI-SO procedure where SA-CASSCF wave functions are employed for computing SOMEs while adapting the energies of the states by shifting them to their CASPT2 energies [29][30][31]. The same applies to CASSCF/NEVPT2 calculations of spin-orbit interactions [32]. In that way, the perturbation corrections to the CASSCF wave functions are not accounted for when computing SOMEs. A slightly more general multireference second-order Møller-Plesset perturbation theory (MR-MP2) approach, albeit with configuration selection, was devised by Grimme and Waletzke [33] and is implemented in their mrci code. In this approach, all configurations, whose interactions with one of the reference vectors exceed a user-defined threshold, are included in the first-order space and contribute to the first-order perturbed wave function that can subsequently be used to determine SOMEs. The discarded configurations contribute solely to the MR-MP2 energy.
In this work, we investigate how the DFT/MRCI-R and MR-MP2 wave functions perform in spin-orbit coupling calculations on light-element compounds. To validate our results, we assembled data on energies and SOMEs of eight diatomic and fifteen polyatomic organic molecules. For the diatomic molecules, the excitation energies and spin-orbit coupling constants can be checked directly against experimental data. Energies of the polyatomic molecules can be compared with experimental data to a limited extent only. In particular, experimental data on optically dark states -that are often essential for the photophysical behaviour of a molecule and in particular for intersystem crossing -are missing for comparison. For some of the smaller molecules (formaldehyde, thioformaldehyde, nitromethane, furan, and thiophene), we carried out benchmark calculations at the CASSCF/CASPT2 level to obtain SOMEs and excitation energies to compare with. Excitation energies of the experimentally unknown states of the other molecules were assessed with respect to CASPT2 values from the literature [28,[34][35][36][37][38] where available.
Moreover, numerous examples exist in which the combination of DFT/MRCI and Spock yielded higher order spin-dependent properties such as phosphorescence and intersystem crossing rates in excellent agreement with experimental data [2,[8][9][10]12,13]. Hence, these SOMEs are considered reliable and may serve as a reference.
Furthermore, we shed light on an earlier assumption that SOMEs are more robust with regard to the wave function quality than energies [17]. In that work, the relaxation pathways of photoexcited nitrobenzene were studied. As mentioned already above, artificially low-lying excited singlet states with four open shells were encountered in the DFT/MRCI spectrum. For that reason, energies from the algebraic diagrammatic construction scheme ADC(3) were combined with SOMEs from DFT/MRCI wave functions. The ab initio MR-MP2 method treats the four open-shell configurations in a qualitatively correct manner which allows us now to quantitatively address this assumption.

Geometries and general remarks
The ground-state geometries of all polyatomic molecules were optimised at the DFT level, employing Becke's three-parameter exchange functional combined with the Lee-Yang-Parr correlation functional (B3LYP) [39][40][41]. Throughout, a valence triple zeta basis set with polarisation function (TZVP) [42] was used with additional polarisation functions on sulphur atoms (TZVPP) [43]. These calculations were performed with the Turbomole program package, version 6.5 [44]. The ground-state geometries were employed in all computations of excitedstate energies (which corresponds to vertical transitions) as well as the SOMEs. For some molecules with lowlying Rydberg states, we augmented the atomic orbital (AO) basis by a diffuse Gaussian s, p, and d function [45] after the geometry optimisation. These bases have been marked as augmented in contrast to valence for the pure TZVP/TZVPP basis sets.
For the diatomic molecules, experimental equilibrium bond distances of the ground and the first excited 3 (first 3 g for P 2 ) states were used [46] to enable comparison with experimental fine-structure splitting constants. In these cases, adiabatic energies were calculated.
Excited-state energies and wave functions were obtained by methods capable of treating multi-reference character in an efficient way (see below). The mrci program makes use of an interface to Turbomole that provides spin-independent one-electron integrals and three-index-two-electron integrals. For the resolutionof-the-identity (RI) approximation of the four-indextwo-electron integrals, optimised auxiliary basis sets (TZVP, TZVPP) from the Turbomole library [43] were utilised.

Electron correlation
Electron correlation calculations were carried out in four different manners as described in more detail below. All methods are implemented in the mrci code originally written by Grimme and Waletzke [1,33], later extended and parallelised in our laboratory [16,47].
In the electron correlation calculations, all electrons of the diatomic molecules and only the valence electrons of the polyatomic molecules were explicitly correlated. Unless noted otherwise, the initial reference space was generated in a restricted active space (RAS) manner with a maximum of two excitations out of ten electrons in ten orbitals. In total, 21 singlet and 20 triplet roots were obtained for all molecules. When molecular point-group symmetry is exploited, these roots do not necessarily relate to the lowest eigenvalues. In C s symmetric molecules, we solved for 11 1 A and 10 3 A states as well as for 10 roots each of 1, 3 A symmetry. In C 2v point-group symmetry, we computed six totally symmetric roots with singlet multiplicity and five roots each for other combinations of spin and spatial symmetry.

... DFT/MRCI calculations
The DFT/MRCI method designed by Grimme and Waletzke [1] is a combination of DFT (which gives information about dynamic correlation) and truncated MRCI expansions (to take the static correlation into account). It makes use of a semi-empirical Hamiltonian that has been parameterised in combination with the BH-LYP functional. In the original work, five parameters that depend on the multiplicity, the excitation class with respect to a closed-shell reference and the number of open shells of the configuration were fitted to experimental excitation energies. Only a small number of configuration state functions (CSFs), instead of the full MRCI space, is used in the actual calculation. Double counting of correlation is avoided by exponential scaling of the offdiagonal matrix elements. Configurations whose contributions would have been scaled away are discarded according to an energy gap criterion δE sel . δE sel was set to 0.8 E h for constructing the final reference space where all configurations with a squared coefficient exceeding 0.003 in one of the roots were kept. In the subsequent DFT/MRCI calculations that generate wave functions and excitation energies, the threshold was increased to 1.0 E h . RAS space between degenerate orbitals, in most of the diatomic molecules (AlCl, AlF, BCl, BF, P 2 and SiO), the 10-11-2 restricted active space was employed. Because of the same reason, AlBr needed the inclusion of yet another active orbital to obtain eventually a 10-12-2 RAS. In the MR-MP2 calculations, orbitals with energies larger than 10 E h were discarded. After every iteration, all configurations with a squared coefficient exceeding 0.003 in one of the roots were added to the reference space. The threshold for the truncation of the first-order space was set to 10 −7 E h , as recommended in the original publication [33]. This general approach had to be altered for some difficult cases when using the augmented basis. The calculations on thiophene were performed with the 10-12-2 RAS and the admixture of HF Rydberg and valence orbitals in thioformaldehyde and nitromethane led to the use of 10-11-2 as the initial reference space. The biggest starting RAS (12-12-2) was chosen for bithiophene in order to include both linear combinations of the lone-pair n(S) orbitals.
We employed the g0 operator as the unperturbed zeroth-order Hamiltonian, which is given by a sum of generalised one-particle Fock operators [33]. The MO basis and Fock matrix elements were generated in two different ways. First, Hartree-Fock (HF) orbitals were used, but this approach has possible difficulties because of the diffuse character of the virtual orbitals since these are generated in the field of N (instead of N−1) electrons. In their original publication, Grimme and Waletzke recommended to employ improved virtual orbitals (IVOs) [33]. Instead, we use Kohn-Sham (BH-LYP) orbitals as a second set of one-particle basis which do not suffer from the same difficulties as the HF orbitals. This additionally allows for a direct comparison of the performance between the methods lifting the orbital influence.

... CASSCF/CASPT calculations
State-averaged CASSCF and subsequent CASPT2 calculations were carried out with the MOLPRO package [48][49][50][51] using exactly the same AO basis sets as in the other calculations. We employed the minimum number of active electrons and active orbitals able to describe the desired molecular states appropriately. As the chosen active space varies from molecule to molecule, detailed specifications will be given below where the results for the individual molecules are presented.

Electronic spin-orbit coupling
Spin-orbit matrix elements for the DFT/MRCI and MR-MP2 wave functions were calculated using the spin-orbit coupling kit (Spock) developed in our laboratory [2][3][4]. It employs an effective one-electron mean-field approximation to the Breit-Pauli Hamiltonian [52]. A further simplification is introduced by neglecting all multicentre integrals, and atomic mean-field integrals are calculated with the AMFI program [53]. These approximations introduce errors which are usually lower than 5% [54,55]. To match the DFT/MRCI and MR-MP2 SOMEs as closely as possible, SOMEs of the CASSCF wave functions were computed with the one-centre approximation (als option) for the one-and two-electron Breit-Pauli spin-orbit integrals in MOLPRO [31].
The spin-orbit-free mrci code [1,33] computes only the M s = S sublevel of a multiplet state. This allows the calculation of triplet-triplet SOMEs in addition to singlet-triplet SOMEs which would not be possible if only the M s = 0 wave function were available. To enable the fast determination of SOMEs between the other multiplet sublevels by means of the Wigner-Eckart theorem [56], Spock expresses the spin part of the Hamiltonian in terms of components of a first-rank tensor operator, i.e.Ŝ +1 ,Ŝ 0 , andŜ −1 . In contrast, the spatial part of the spin-orbit Hamiltonian is adapted to the point-group symmetry of the molecule which corresponds, in most cases, to a Cartesian representation. Because of the relationŝ it is easy to transform the mixed representation of the Hamiltonian to a purely Cartesian one [57]. In the following, the x-and y-components of SOMEs between the M s =

Test set
The selected test molecules include eight diatomic molecules (AlBr, AlCl, AlF, BCl, BF, CS, P 2 , SiO) and fifteen polyatomic molecules (see Figure 1). The diatomic molecules primarily consist of first-and second-row elements, where first-order spin-orbit interactions still dominate the fine-structure splittings of spatially degenerate states. Their excitation energies and SOMEs will be compared directly to experimentally derived adiabatic transition energies and spin-orbit coupling constants. Unfortunately, such a comparison is not possible for spatially non-degenerate states of polyatomic molecules. In these cases, only off-diagonal SOMEs may adopt significant values. However, numerous examples have been presented in which the combination of DFT/MRCI and Spock yielded higher order spin-dependent properties such as phosphorescence and intersystem crossing rates in excellent agreement with experimental data [2,[8][9][10]12,13]. Hence, these results can serve as a reference. In addition, CASSCF/CASPT2 calculations on some of the smaller polyatomic molecules were carried out for comparison. The polyatomic molecules in the test set are chosen such that substantial spin-orbit coupling is expected which involves the presence of heteroatoms and states of different electronic characters. The states selected for comparison are mainly the ground and low-lying excited, photophysically and photochemically important states with ππ * and nπ * characters, which allows for an easier correspondence between the methods. In some cases, higher lying excited states with significant SOMEs are also included. Because their wave functions tend to be mixed with Rydberg excitations and contributions from doubly excited configurations, the correspondence is more difficult to determine. In the following, the singlet and triplet excitation energies and the computed SOMEs of all polyatomic are presented and analysed. Tables with detailed information on the experimental reference data, oscillator strengths for electric dipole-allowed transitions and the wave function composition, as well as the optimised geometries, are given in the SM. Additionally, statistical data for each molecule is shown there.

Diatomic molecules
Spatially degenerate states of diatomic molecules exhibit a first-order spin-orbit splitting. From spectral data involving different substates, the fine-structure parameter A SO can be derived which quantitatively describes the firstorder zero-field splitting at the equilibrium bond distance of the respective state [58]. A SO can be correlated with calculated SOMEs and hence gives valuable reference data for testing theoretical predictions. From the axial symmetry of diatomic molecules, with the z-axis being the internuclear axis, it follows that only the z-component of the spin-orbit Hamiltonian couples the degenerate states. This simplifies the phenomenological SO Hamiltonian to the expression Here, and are projections of the total spatial and spin angular momenta on the internuclear axis, respectively. Our state of interest is the first triplet of symmetry, in which case = = 1. The spin-orbit coupling is introduced perturbatively. If we set up a matrix for the first-order degenerate perturbation theory in the basis of LS coupled states, we see that the fine-structure parameter A SO is exactly equal to the offdiagonal matrix element between the M s = 1 components of the degenerate x and y sublevels, i.e. A SO = 3 x , M s = 1|Ĥ SO z | 3 y , M s = 1 . The first triplet state of symmetry in the considered molecules is a regular state with configuration σ 1 π 1 . In P 2 , it is symmetric with respect to inversion, g , and its configuration is σ 1 g π 4 u π 1 g . For a straightforward comparison with experimental data, adiabatic energies are calculated using the experimental values for the bond lengths. Tables 1 and 2 summarise the calculated and experimental energies and fine-structure parameters of the 1 3 r states, while Tables 3 and 4 present the corresponding statistical data. As a general trend, we notice from the graphs in Figure 2 that all methods underestimate the excitation energies of the 1 3 r state on the average. With correlation coefficients of 0.9856 (DFT/MRCI-R) and 0.9834 (DFT/MRCI), the MRCI methods perform somewhat better than the MR-MP2 methods (0.9629 for MR-MP2(HF) and 0.9658 for MR-MP2(BH-LYP)) when it comes to 3 r excitation energies. The maximum negative deviations are observed for the SiO molecule. A reason may be the unusually strong Si=O bond with at least double bond character (bond length 1.509739Å [46]) in the electronic ground state which poses a challenge for the quantum chemical description.
The correlation plots for the SOMEs of the diatomic molecules are presented in Figure 3. (Note that the experimentally derived zero-field splittings may contain also higher order spin-orbit interactions and contributions from electronic spin-spin coupling which are assumed to be small.) The agreement with the experimental splittings is excellent for all four methods, with equal RMSD of 2 cm −1 , and maximal positive/negative deviations of 3.1/2.6, 3.0/2.7, 2.1/3.5 and 3.8/1.6 cm −1 for Table . Fine-structure splitting of the   (  g for P  ) state with the configuration σ  π  (σ 1 g π 4 u π 1 g ) in selected diatomic molecules, characterised by the presented spin-orbit coupling coefficient A SO (cm −

... Vertical energies
Individual data regarding computed vertical excitation energies of all polyatomic molecules are compiled in             Table 7 contains the RMSD values, maximum positive as well as negative deviations for the complete set of experimentally known excitation energies of electronic states studied in this work. (Please note that the total number of electronic states, for which SOMEs have been determined, is substantially larger.) Very good agreement between the computed and measured excitation energies is observed for all methods ( Figure 4). All correlation coefficients are close to 1. As for the diatomic molecules, the semi-empirical DFT/MRCI methods perform slightly better than the MR-MP2 calculations that also show a broader scattering (Table 7). In contrast to our expectations, MR-MP2 is in somewhat better agreement with experiment on the average when HF orbitals are employed instead of BH-LYP KS orbitals.
If the excitation energies of all calculated states are compared among the different computation methods, the correlation deteriorates, of course. Excitation energies for low-lying Rydberg states have been computed for formaldehyde, thioformaldehyde, furan, thiophene, and nitromethane. For some of the singlet states, experimental reference data are available; in the other cases, CASPT2 energies may be used for comparison. Both DFT/MRCI parameterisations and the two MR-MP2 variants show very good agreement with the reference data. In general, the agreement for the valence-excited states is satisfactory, but in a few cases, larger deviations are recognised ( Figure 5). Inspection of Tables 5 and 6 reveals that the deviations between the two DFT/MRCI parameterisations are all related to doubly excited states of one type or the other. These cases will be analysed in the following on an individual basis. Furthermore, we will address some cases where the different energy splittings between states of equal symmetry in MR-MP2, as compared to DFT/MRCI, result in different wave function compositions, which in turn lead to deviations of the corresponding SOMEs.

Special cases.
Formaldehyde. Being the smallest of the polyatomic molecules, formaldehyde represents the least demanding system for the MR-MP2 methods with regard to the size of the reference space and the first-order space selection threshold. The energy contribution arising from discarded configurations of the first-order space is less than one percent of the total correlation energy, making the results for formaldehyde nearly converged with respect to the selection procedures in the MR-MP2 approach. CASSCF/CASPT2 calculations were carried out for two different active spaces. The first active space where eight electrons were distributed in nine active orbitals comprised two σ MOs, the valence π and π * MOs, the n orbital and the Rydberg 3s and 3p orbitals; in the larger (8,11) calculations, two σ * MOs were added to the active space. The level shift of 0.2 that was used for CASPT2 (8,9) had to be increased to 0.3 in the CASPT2 (8,11) calculations due to problems with intruder states. In addition to our CASPT2 results, high-level ab initio values are available from the literature for comparison [34][35][36][37]59].
For all methods, the calculated excitation energies of the 1 1 A 2 state (Table 5) are in good agreement with the experimentally observed energies of 3.79 [60] and 3.94 eV [61], respectively, for the energy-loss maximum (vapour) of the 1 1 A 2 state. Equally good correspondence (Table 6) with the experimentally measured value of 3.50 eV [60,61] is found for the first triplet state, 1 3 A 2 . The excitation energy of the 1 3 A 1 (ππ * ) state is somewhat underestimated by the DFT/MRCI methods, whereas the results for the MR-MP2(HF) and MR-MP2(BH-LYP) are a bit higher than the energy-loss maximum at 5.82 eV [60] or 5.86 eV [61]. All methods find the Rydberg n → 3s C and n → 3p C excitations in good agreement with experimental data [62]. The DFT/MRCI and MR-MP2 approaches predict similar energies for the valence 4 1 A 1 (ππ * ) state (5 1 A 1 in DFT/MRCI-R) which is experimentally not observed, but suspected to lie in the Rydberg region between 7 and 12 eV [60,61]. Nevertheless, there is one peculiarity with regard to the wave function of that state that gives rise to two outliers when correlating the SOMEs of the DFT/MRCI and DFT/MRCI-R methods. We recognise a pronounced weight of a doubly excited closed-shell configuration (n 2 → π * 2 ) in the DFT/MRCI-R wave function which is not seen among the leading configurations of the DFT/MRCI and MR-MP2 wave functions (Tables S11 and S12 of the SM). The state with leading (n 2 → π * 2 ) term shows up at 9.2 eV in DFT/MRCI-R, whereas DFT/MRCI and both MR-MP2 variants places it at markedly higher excitation energies. The CASSCF/CASPT2 results for the smaller (8,9) active space support the order found by the DFT/MRCI and MR-MP2 methods. The (n 2 → π * 2 ) double excitation exhibits only a small coefficient in the CASSCF (8,9) wavefunction, and an optically very bright ππ * state is obtained with an excitation energy of 9.33 eV after PT2 correction. The state with leading (n 2 → π * 2 ) configuration is found at 10.29 eV in these calculations. Literature values of the vertical 4 1 A 1 (ππ * ) excitation energies range between 9.47 eV for EOM-CCSD [35] over 9.77 eV for CASPT2 [37] and 9.80 for MR-CISD+Q [34] to 10.15 eV for MR-CISD [34]. For the doubly excited n 2 → π * 2 state, MR-CISD(+Q) results are available for comparison [34] that place this state at 10.63 eV (10.54 eV). All these data point toward a problem of the redesigned DFT/MRCI-R in properly describing this closed-shell double excitation of type (n 2 → π * 2 ). The mixture of these high-lying singly and doubly excited configurations is also an issue in the CASSCF/CASPT2 calculations. In the larger active space (8,11) that was primarily chosen to improve the description of the (σ → π * ) excitations, strong configuration interaction of the A 1 symmetric (π → π * ), (n 2 → π * 2 ), and (π → 3p x ) excitations takes place leading to a substantially higher excitation energies for these states than the CASPT2 (8,9) treatment (see Tables S10 and S12 in the SM).
Thioformaldehyde. Substitution of oxygen in formaldehyde by sulphur to obtain thioformaldehyde additionally lowers the energy of the zero-open shell doubly excited 1 A 1 state (n 2 H → π * 2 L ). The effect is very pronounced in the case of the DFT/MRCI-R Hamiltonian which brings this state below the 1 A 1 (π H−1 → π * L ) state and quite close to it ( Table 5). The energetic proximity results in strong wave function mixing so that the 1 A 1 state with leading (54%) (π H−1 → π * L ) term exhibits major contributions (22%) from the (n 2 H → π * 2 L ) double excitation. In addition, contributions from the nearby (n H → 4p x ) excitation are found in the wave function (Table S19 in the SM). In contrast, the admixture is quite small for the standard DFT/MRCI with minor configurations contributing by only 2% (Table S18 in the SM). All methods slightly overestimate the vertical energy of the bright (π H−1 → π * L ) state of thioformaldehyde. With a value of 6.30 eV, MR-MP2(HF) is closest to the CASPT2 (8,11) value of 6.29 eV and the experimentally measured peak maximum which is located at 6.2 eV [63]. No experimental reference data are available for the doubly excited 1 A 1 state with leading (n 2 H → π * 2 L ) term. The original DFT/MRCI parameterisation places this state about 0.7 eV above the CASPT2 value, whereas DFT/MRCI-R places it about 1 eV below the CASPT2 reference. The MR-MP2 methods yield energies in good agreement with CASPT2. Notice, however, that the CASSCF wave function shows configuration mixings very similar to the DFT/MRCI-R Hamiltonian (Table S19 in the SM).
We find also a pronounced difference for the excitation energy of another pair of dark states which are doubly excited with respect to the electronic ground state. The (π 1 H−1 n 1 H π * 2 L ) occupation gives rise to a singlet and a triplet state of A 2 spatial symmetry, 3 1 A 2 and 4 3 A 2 . The energy of the 3 1 A 2 state (7.19 eV) is markedly lower for the standard DFT/MRCI method than for all other methods employed (see Table 5). Since it is also substantially lower than the energy of the corresponding 4 3 A 2 state, we consider this to be a problem of the original DFT/MRCI parameterisation. The energy shift of that configuration has almost no consequence with regard to the SOMEs, however, as there is no close-lying state of equal symmetry with which the wave function could mix.
o-Benzyne. The 2 1 A 1 state of o-benzyne is dominated by a configuration with four open shells which arises from an excitation of the in-plane π orbital and the highest occupied out-of-plane π H orbital to the corresponding antibonding π * and π * L+1 orbitals. In the standard DFT/MRCI calculations, it is located at 4.24 eV. A hint that there might be something wrong comes from the fact that the lowest triplet state with the same leading configuration appears at 5.69 eV. The other methods place this doubly excited 1 A 1 state more than 2 eV higher in energy (Table 5) while the corresponding triplet is shifted upwards by about 0.7 eV (Table 6). In all cases except for DFT/MRCI, it is located well below the respective singlet state. The reason for the problems of the standard DFT/MRCI method in describing doubly excited states with four open shells properly has been traced back to the parameterisation of the effective Hamiltonian [16]. In particular, the parameter that scales the exchange contributions to the diagonal elements depends strongly on the excitation class, the number of open shells, and the multiplicity [1]. Also for the second low-lying singlet state with four open shells, 2 1 B 1 (π 1 π 1 H−1 π * 1 π * 1 L+1 ), good agreement among the DFT/MRCI-R (7.25 eV) and the MR-MP2 methods (both 7.22 eV) is observed, whereas the standard DFT/MRCI value is substantially smaller (5.74 eV).
Nitrobenzene. The questionable description of doubly excited singlet and triplet states with four open shells by the DFT/MRCI method as manifested, for example, in large negative singlet-triplet energy splittings, actually was the reason why we started looking for an alternative method to compute SOMEs of nitrobenzene. Double excitations with four open shells play an important role even in low-lying excited states of this molecule. The standard DFT/MRCI parameterisation associates with them too low energies, and as a consequence, too high contributions in the state vector [17]. Our calculated singlet excitation energies of nitrobenzene can be compared with experimental data for an (n → π * ) transition and an optically bright (π → π * ) transition. In addition, high-level CASPT2 (14,11) excitation energies are available in the literature for comparison [38]. As a general trend, all methods employed in this work place the singlet states at too low excitation energies. Particularly problematic is the description of the B 1 states in the standard DFT/MRCI, where a configuration with spatial occupation (π 1 H−2 π 1 H π * 1 L π * 1 L+1 ) dominates the wave function of the lowest 1 B 1 state. The electron density in the π H−2 and π * L MOs of nitrobenzene is mainly localised on the NO 2 group, whereas π H and π * L+1 are pure benzene MOs (see Figure S6 in the SM). In this way, nitrobenzene exhibits a tendency towards a bi-chromophoric system for which the redesigned DFT/MRCI-R method [16] was developed. In the DFT/MRCI-R wave function, the two CSFs associated with this configuration contribute only with about 0.1% each. The MR-MP2 wave functions have this configuration contributing by less than 2%. The next DFT/MRCI vector (2 1 B 1 ) has, as a dominant, the four open-shell configurations (33%). Its excitation energy of 4.63 eV is much lower than that of the corresponding triplet state (4 3 B 1 ) which is found at 5.88 eV. For that reason, the DFT/MRCI description of the 2 1 B 1 state is regarded as being unphysical. The MR-MP2 methods place the singlet state with leading (π 1 H−2 π 1 H π * 1 L π * 1 L+1 ) configuration at 7.19 eV (MR-MP2(HF)) and 7.00 eV (MR-MP2(BH-LYP)). In the DFT/MRCI-R calculation, a state of that type does not appear among the five lowest roots of 1 B 1 symmetry. Its energy must, therefore, be higher than that of 5 1 B 1 , i.e. >7.3 eV.
Dithiosuccinimide. DFT/MRCI-R finds a low-lying pair of A 1 states (1 3 A 1 at 3.36 eV and 2 1 A 1 at 3.40 eV) with a leading configuration that corresponds to a double excitation from the in-plane lone-pair orbitals on sulphur, i.e. n H−1 in the a 1 irreducible representation (irrep) and n H in the b 2 irrep, to π * L and π * L+1 which correspond to the a 2 symmetric (π * L ) and b 1 symmetric (π * L+1 ) linear combination of the antibonding thiocarbonyl π * orbitals. In addition, π * L+1 exhibits contributions from nitrogen p . The standard DFT/MRCI places these states with leading n 1 H−1 n 1 H π * 1 L π * 1 L+1 term at 4.68 eV (2 3 A 1 ) and 4.71 eV (2 1 A 1 ), about 1.3 eV higher in energy. Since the parameters of the redesigned DFT/MRCI-R Hamiltonian are spin independent, we carried out a calculation for the lowest quintet states, too. Not surprisingly, we find a state with leading n 1 H−1 n 1 H π * 1 L π * 1 L+1 configuration as the 1 5 A 1 state with an excitation energy of 3.24 eV. This means, that the 1 5 A 1 , 1 3 A 1 , and 2 1 A 1 states of dithiosuccinimide, as obtained from DFT/MRCI-R calculations, originate from a rather loosely coupled pair of 3 (n → π * ) excitations. While the multiplet splitting appears to be of similar size in the DFT/MRCI Hamiltonians, the Coulomb repulsion between the pairs is substantially larger in the original DFT/MRCI case. MR-MP2(HF) puts the corresponding triplet at 6.67 eV (4 3 A 1 ) and the quintet at 6.47 (1 5 A 1 ). In the singlet manifold, a state with such electronic structure is not found among the lowest six states of A 1 symmetry. We, therefore, performed an additional MR-MP2(HF) calculation where we solved for 10 roots in A 1 symmetry. In that calculation, the desired singletcoupled wave function is found for the 5 1 A 1 state with an energy of 6.37 eV. At present, it is not clear whether the DFT/MRCI methods, and in particular the redesigned DFT/MRCI-R method, have problems with a balanced description of this particular (n, n → π * π * ) double excitation in dithiosuccinimide or whether the true problem lies rather in the truncation of the MR-MP2 first-order space.
Also, the 3 1 A 1 state of the DFT/MRCI-R at 5.08 eV stems from a double excitation, in this case with two leading closed-shell configurations, i.e. n 2 H π * 2 L and n 2 H−1 π * 2 L . MR-MP2(BH-LYP) finds this state as 4 1 A 1 at 5.26 eV, MR-MP2(HF) at 5.51 eV (MR-MP2(HF) with ten roots in 1 A 1 symmetry yields 5.26 eV), whereas it is the fifth root of 1 A 1 symmetry in the DFT/MRCI calculation with an energy of 6.07 eV. This finding is consistent with the observations made for n 2 → π * 2 double excitations in formaldehyde and thioformaldehyde. The second pair of A 2 states, 2 3 A 2 and 2 1 A 2 , should be uncritical as their wave functions are primarily made up from single excitations with n 1 H π * 1 L+1 as the leading term. Their energies agree well among the DFT/MRCI and DFT/MRCI-R methods (Table 5 and 6). In this case, the MR-MP2 methods add stronger contributions from doubly excited configurations involving the π H−2 orbital (the b 1 symmetric combination of the bonding thiocarbonyl π orbitals) and one of the n orbitals. This brings their energy further down. While the MR-MP2 (HF) energies are in reasonable agreement with the DFT/MRCI results, MR-MP2 (BH-LYP) places these states at much lower energies. A similar observation is made for the second pair of B 1 states, 2 3 B 1 and 2 1 B 1 (see Tables S52 and S53 of the SM) which has a similar composition as the just described with two open shells, but again two n electrons are transferred to π shells. Nitromethane. We performed CASSCF/CASPT2 calculations with various active spaces for this compound in the augmented basis. Unfortunately, we were not able to find an active space yielding a balanced description of the π → π * and Rydberg transitions. Attempts to use a large active spaces, for example (12,11), were unsuccessful because the state-averaged CASPT2 calculations did not converge. The results presented here are for the (8,7) active space which comprise 2 π and 2 n orbitals, 1 π * and 1 σ * MO as well as a 3s orbital. The fact that the n → 3s excitations are stabilised so strongly by the CASPT2 corrections is caused by a mixture of 3s and σ * character in the active orbitals. All methods, in particularly the MR-MP2 methods, predict too low energies for the first two singlet states of nitromethane ( Table 5). The DFT/MRCI-R and CASPT2 excitation energies of these states are practically identical. Both states are of nπ * character and A symmetry, with the electron energy-loss maximum of 1 1 A being measured at 4.25 eV in the gas phase and the absorption maximum of 2 1 A at 4.50 eV in the gas phase [64]. When Rydberg functions are not included in the basis set (Table S42 in the SM), the excitation energy of the first ππ * singlet state (2 1 A ) is overestimated, again in larger measure by the MR-MP2 method. Augmentation of the basis by diffuse functions brings the excitation energy of this state down for all methods and in good agreement with the expeirmental value of 6.25 eV [64] obtained as the maximum of absorption in the gas phase. In the triplet manifold, we observe one state that attracts our special attention. The leading term of the 2 3 A state is a double excitation with two open shells: from each of the in-plane n and n orbitals of the nitro group, one electron is excited to the lowest lying π * MO. The results of the MR-MP2 calculations agree well with the CASPT2 excitation energy, whereas DFT/MRCI and in particular DFT/MRCI-R places this double excitation at significantly lower energies.
Furan and thiophene. CASSCF/CASPT2 calculations of furan distributed six active electrons in nine active orbitals comprising five valence MOs (3 π, 2 π * ) and the 3s and 3p orbitals. For thiophene, one σ * orbital was added to the active space. Very good consensus is achieved among the theoretical methods for furan and the agreement with the available experimental data is good. We, therefore, refrain from discussing this compound further. In thiophene, similar conclusions can be drawn with regard to the 1, 3 (π → π * ) excited states which are found in the A 1 and B 1 irreducible representations for our choice of coordinate system. There is one exception though. MR-MP2(HF) yields an excitation energy of 6.14 eV in the augmented basis while a value of 5.75 eV was obtained in the valence basis. In contrast, DFT/MRCI, DFT/MRCI-R, and MR-MP2(BH-LYP) excitation energies are lowered by about 0.15 eV upon adding diffuse functions to the basis set which brings them into excellent agreement with the experimental value of 5.61 eV determined by electron energyloss spectroscopy [65]. The wave functions show some mixing with (π → 4p y ) excitations, with the diffuseness increasing from DFT/MRCI over MR-MP2(BH-LYP) to MR-MP2(HF). Although barely visible from the wave function composition (Table S30 in the SM), the MR-MP2(HF) state must be considered significantly more Rydberg-like than the wave functions obtained by the other methods. To understand the trends observed for the SOMEs (see below), it is important to notice also that the 1 1 B 2 , 1 3 B 2 , 1 1 A 2 , and 2 3 A 2 states are mixtures of (π → σ * ) and (π → 4p x ) excitations, with the DFT/MRCI and DFT/MRCI-R wave functions being more valencelike than their MR-MP2(HF), MR-MP2(BHLYP), and CASSCF counterparts. Furthermore, we note that the wave function composition of the 3 1 A 1 differs dramatically between the DFT/MRCI and MR-MP2 methods in the augmented basis sets (Tables S30 and S31 in the SM), although the energy is nearly the same for all methods and in good agreement with the experimental value.
Quinoxaline and quinazoline. For these two naphthalene analogues, we observe a general trend that the MR-MP2(HF) singlet excitation energies are markedly lower than their DFT/MRCI counterparts ( Table 5). The difference is even more pronounced for MR-MP2(BH-LYP). There is one exception from this trend, namely for the 1 1 B 2 state of quinoxaline and the 3 1 A state of quinazoline. Both transitions are dominated by (π H → π * L ) excitations, thus representing L a states. In quinoxaline, the L b state which is characterised by a nearly equal mixture of (π H−1 → π * L ) and (π H → π * L+1 ) excitations can be identified unambiguously with the 2 1 A 1 state, whereas the assignment to the 3 1 A state of quinazoline is not so clear-cut because of the lower symmetry of that compound. The reasons underlying these differences are not obvious. One might guess that the first-order interacting spaces in the MR-MP2 calculations were too small and that the PT2 corrections overshooted. In a previous work, DFT/MRCI has been shown to perform very well for the L a and L b states of polyacenes [66]. We, therefore, consider the DFT/MRCI energies trustworthy. As the energy shifts do not lead to substantial variations of the wave function composition, they do not have consequences with regard to the SOMEs.  (Table 5). As a general trend, all DFT/MRCI excitation energies of this compound are higher than the MR-MP2 ones and the maximal negative deviations are 0.53 and 0.76 eV for the MR-MP2(HF) and MR-MP2(BH-LYP), respectively. With respect to the forthcoming analysis of the SOMEs, it is important to note that the energetic separation of the second and third 3 B 1 states that are dominated by linear combinations of (π → π * ) and (π → σ * ) excitations is much smaller in the MR-MP2 calculations compared to the DFT/MRCI methods. The same is true for the second and third 3 A 1 states. Remarkably, DFT/MRCI and DFT/MRCI-R yield very similar excitation energies for the doubly excited 3 1 A 1 state. The leading term in the wave function of this state originates from the closedshell (π 2 H → π * 2 L ) excitation. Bithiophene. The electronic structure of s-trans bithiophene has already been treated in our laboratory using the DFT/MRCI method [68]. Since the same basis and geometry have been used, the vertical excitation energies of S 1 , S 2 , T 1 , T 2 , T 3 and T 4 differ only marginally to the ones presented here due to a different choice of reference space. In the present study, higher lying states are added among which the πσ * ones (3 3 B and 3 3 A) are especially interesting, because of their strong SO interaction with the lower ππ * states. The S 1 and T 1 energies of 3.86 and 2.32 eV, respectively, determined by Siegert et al. [68] in anion photodetachment studies in a jet stream, are 0-0 energies which are close to the adiabatic DFT/MRCI energies as shown by these authors. The computed vertical DFT/MRCI energy of 1 1 B is in better agreement with the experimentally measured absorption maxima in apolar solvents (4.09 eV in dioxane, 4.11 eV in methylcyclohexane) [69]. There is a second, less intensive peak at 5.02 eV in the absorption spectrum in dioxane which can be assigned to the 2 1 B state of bithiophene with corresponding DFT/MRCI energy of 4.95 eV. In order to include both linear combinations of the lone-pair n(S) orbitals, the 12-12-2 space was used in the first iteration. This is of importance only for the MR-MP2 method, since the DFT/MRCI results are not significantly different compared to the calculation with the initial 10-10-2 RAS. The MR-MP2 energies of S 1 are slightly lower than the DFT/MRCI values but all agree well with the experimental 1 1 B absorption energy (Table 5). In contrast, the 2 1 B MR-MP2 energies are significantly lower than the experimentally determined and the DFT/MRCI energies, especially when BH-LYP orbitals are employed. It appears that the first-order interacting space is too small in this case and that the PT2 corrections of the discarded configurations overshoot. Lowering the selection threshold for the inclusion of configurations in the first-order perturbed wave function to 5 × 10 −9 E h increases the MR-MP2(BH-LYP) energy of 2 1 B to 4.29 eV. A further decrease of the threshold to test whether the results are converged was technically not feasible.
Methionine. To the best of our knowledge, there are no experimental data of excitation energies of methionine which could be used for comparison with our vertical energies. DFT/MRCI-R energies are consistently higher than the DFT/MRCI ones by about 0.2-0.3 eV. The striking result for methionine is the huge variation of the MR-MP2 energies regarding the employed type of MOs. MR-MP2(HF) energies are much higher than MR-MP2(BH-LYP), with the difference reaching 1 eV for the 3 3 A state. On the other hand, MR-MP2(BH-LYP) excitation energies are very close to the DFT/MRCI values. Examination of the HF and BH-LYP MOs confirms that virtual orbitals obtained by HF and DFT procedure differ by such an extent that it is difficult to make a correspondence between them. This, of course, also leads to a quite different MR-MP2(HF) and MR-MP2(BH-LYP) state vectors (see Table S55, SM).
Isoalloxazine. Isoalloxazine is the photosensitive core of the flavin family of chromophores. They play an important role in blue-light sensing proteins. It is a heteroaromatic compound with four of its six heteroatoms carrying lone-pair orbitals. The photophysics of this chromophore is determined by the coupling between optically dark nπ * and bright ππ * states and their subtle interactions with the environment. The only available experimental data for comparison are absorption spectra of flavin derivatives in solution. 8-Methyl isoalloxazine in ethanol absorbs with peak maxima at 2.85 eV (2 1 A , ππ * ) and 3.76 eV ( 3 1 A , ππ * ) [70]. Theoretical studies on flavins in various solvent environments point out the small solvatochromicity of the first band (stabilisation with regard to vacuum by 0.06 eV in methanol) and a more pronounced red shift of the second maximum with increasing solvent polarity (stabilisation by 0.21 eV) [12]. Thus, the DFT/MRCI and DFT/MRCI-R vertical excitation energies of the first two ππ * states (Table 5) are considered to be in good agreement with experimental evidences. No experimental data are available for the nπ * states. For these states, we notice that the DFT/MRCI-R energies are systematically higher by approximately 0.1 eV. This finding is in line with the general trend for the DFT/MRCI-R energies of nπ * states observed by Lyskov et al. [16] Isoalloxazine is a challenging molecule for MR-MP2. The method, as implemented in the mrci program, is at its limits here, since the code is not parallelised and allows only for a maximum number of 1000 reference configurations. The latter restriction is severe, because of the large number of active orbitals in isoalloxazine. The computed MR-MP2 excitation energies of isoalloxazine are strongly dependent on the selection threshold, the choice of reference space, and the employed orbital basis, with MR-MP2(HF) being more consistent with the experimental findings. Setting the selection threshold (Esel) for the inclusion of configurations in the first-order space to the standard value of 10 −7 E h leads to energy contributions from the discarded configurations of up to 45% of the total correlation energy in the case of isoalloxazine, while this contribution is lower than 20% for the other molecules (with exception of methionine). Unfortunately, a lowering of the Esel parameter is technically not feasible for the triplet multiplicity because of the large number of CSFs already needed for Esel = 10 −7 E h . In the singlet manifold, the expansion length of the first-order perturbed wave function increases from 122,474,081 CSFs for 10 −7 E h to 348,025,394 CSFs for 10 −8 E h . Test calculations for the singlet states with Esel = 10 −8 E h revealed that the excitation energies of the nπ * states are lowered by up to 0.5 eV in comparison with the standard MR-MP2 calculations.

... Spin-orbit coupling matrix elements
Since there are no experimental data which could be used to judge the quality of our SOMEs for polyatomic molecules, we adopted here the following approach. The matrix elements of the redesigned DFT/MRCI-R and this particular implementation of the MR-MP2 method are presented here for the first time. The standard DFT/MRCI has been successfully used in a number of studies on photophysical and photochemical processes such as intersystem crossing and phosphorescence that depend indirectly on the spin-orbit coupling between the involved states. Absolute values of the calculated SOMEs of all polyatomic molecules are compiled in Table 8. In addition, SOMEs calculated for CASSCF wave functions of formaldehyde, thioformaldehyde, furan, thiophene, and nitromethane in the augmented basis, as well as literature data on formaldehyde [59], have been used for comparison. The phase factors of the SOMEs have been omitted because they depend on the arbitrary phases of the underlying wave functions and MOs. They are relevant only when second-order spin-dependent properties such as phosphorescence probabilities are to be determined [13,57].

General trends.
In order to examine whether the DFT/MRCI-R and the first-order MR-MP2 wave functions provide a reasonable description of the spin-orbit interaction, we have correlated the SOMEs obtained with these methods to the ones calculated with DFT/MRCI. The plots showing this correlation for all SOMEs of polyatomic molecules are presented in Figure 6. Correlation plots for the individual molecules as well as the wave function compositions of the involved states are provided in the SM together with statistical data for each molecule individually and for the complete set of selected SOMEs, 278 in total. Absolute values have been employed when computing the maximum positive and negative deviations of the SOMEs. Normalised RMSD (NRMSD), expressed as percentage, are given as the RMSD normalised by the range, i.e. the maximum value minus the minimum absolute value of the evaluated data.
Excellent agreement is observed between the SOMEs computed with the DFT/MRCI and DFT/MRCI-R wave functions, respectively, with two pairs of outliers that spoil the otherwise good correlation ( Figure S1 of the SM). The outliers are associated with high-lying states of formaldehyde and thioformaldehyde, respectively. Their origin will be discussed in more detail below. The SOMEs obtained for the first-order MR-MP2 wave functions show larger scattering that is even increased when Rydberg functions are added to the AO basis. At first sight, there is not much difference to be found whether HF (Figure S2 of the SM) or KS orbitals ( Figure S3 of the SM) have been employed. In these cases, the largest outliers are found for dithiin, bithiophene, and isoalloxazine.
Special cases. Formaldehyde. SOMEs have been calculated for the CASSCF (8,9) and CASSCF (8,11) wave functions of formaldehyde. In addition, literature values for the spin-orbit coupling in formaldehyde are available for comparison. Formaldehyde has been studied intensively by Langhoff and Davidson by ab initio MRCI methods using the full Breit-Pauli spinorbit Hamiltonian [59]. The SO coupling is relatively strong with SOMEs up to 70 cm −1 and very similar values for all four methods (see Table 8). The MR-MP2 The only exceptions from the good agreement between the DFT/MRCI, MR-MP2, and CASSCF results for the SOMEs of formaldehyde, listed in Table 8, are the matrix elements of the 1 3 A 2 state with the high-lying 4 1 A 1 and 5 1 A 1 states. The leading configuration of the 1 3 A 2 state originates from an (n H → π * L ) (1 3 A 2 ) with respect to the electronic ground state. As discussed at length in Section 3.3.1.2, the 4 1 A 1 and 5 1 A 1 arise from mixtures of the (π H−1 → π * L ) and (n 2 H → π * 2 L ) with weights varying from method to method. A further complication arises due to the energetic proximity of a third 1 A 1 state with Rydberg character which mixes with the valence states. The (n H → π * L ) excitation and the just mentioned valence-excited singlet A 1 configurations are coupled through the large n H |Ĥ SO z |π H−1 and n H |Ĥ SO z |π * L integrals, respectively. The latter is the same integral that dominates the matrix element between the 1 3 A 2 state and the electronic ground state. The Rydberg configurations contribute little to the SOMEs. Depending on the size and phase of the wave function coefficients of the 4 1 A 1 state, SOMEs varying between 22 cm −1 (CASSCF (8,11)) and 57 cm −1 (DFT/MRCI-R) are obtained (Table 8). To clarify this issue further, we can compare our SOMEs with the MRCI values presented by Langhoff and Davidson. These authors found two 1 A 1 valence states, the energetic order of which depends on the chosen MO basis for the MRCI expansion. For the matrix element of the 1 3 A 2 state and the 1 A 1 with leading (π → π * ) configuration, they obtain a value of 41.32 cm −1 which is a bit larger than our SOMEs, save for the DFT/MRCI-R value. For the other SOME involving the 1 3 A 2 state and the 1 A 1 state with leading (n 2 → π * 2 ) configuration, Langhoff and Davidson list a value of 53.08 cm −1 which compares favourably with the corresponding SOME of the DFT/MRCI wave functions (55.0 cm −1 ). The DFT/MRCI-R wave function yields a SOME of 31.5 cm −1 in that case, in accidental coincidence with the CASSCF (8,11) value. The varying weight of the (π → π * ) excitation in the wave functions of the three densely spaced electronically excited singlet states is also reflected in the electric dipole moments for a radiative transition to the electronic ground state (see Table S10 of the SM). Concluding this discussion, we believe that it makes little sense to compare the matrix elements of energetically close-lying interacting electronic states on an individual basis. Rather, their contributions to spin-forbidden transitions such as the phosphorescence of the lowest lying triplet state should be considered as a whole.
Thioformaldehyde. Like in formaldehyde, very good agreement among the methods is observed with respect to the majority of the SOMEs (see Table 8). The above-mentioned large outliers in the DFT/MRCI vs. DFT/MRCI-R correlation of the SOMEs all involve the 2 1 A 1 and 3 1 A 1 states in the valence AO basis which relate to 2 1 A 1 and 5 1 A 1 in the augmented basis due to interjacent Rydberg states. They are coupled with the 1 3 A 2 and 1 3 B 2 states via very large MO integrals, n H |Ĥ SO z |π H−1 , n H |Ĥ SO z |π * L , and σ H−2 |Ĥ SO x |π H−1 , which amount to 273.9, -236.1, and -215.2 cm −1 , respectively, for the BH-LYP MOs in the augmented basis. Such a strong SO interaction makes the above discussed differences in the DFT/MRCI and DFT/MRCI-R wave functions become apparent. The phase factors of the nearly equally large wave function coefficients of the (π H−1 → π * L ) and (n 2 H → π * 2 L ) excitations in the CASSCF leads to an almost perfect cancellation of their contributions to 1 3 A 2 |Ĥ SO z |2 1 A 1 , whereas they add up in 1 3 A 2 |Ĥ SO z |5 1 A 1 (Table 8). Such a cancellation does not happen in 1 3 B 2 |Ĥ SO x |2 1 A 1 or 1 3 B 2 |Ĥ SO x |5 1 A 1 . The leading term of the 1 3 B 2 wave function (σ H−2 → π * L ) connects solely to one of the major contributors to the 2 1 A 1 and 5 1 A 1 wave functions by a single excitation, namely to (π H−1 → π * L ). Indirectly, the 1 3 B 2 |Ĥ SO x |2 1 A 1 and 1 3 B 2 |Ĥ SO x |5 1 A 1 SOMEs, therefore, reflect the weight of that configuration in the 2 1 A 1 and 5 1 A 1 wave functions of the various methods. The SOMEs of the 1 3 B 1 state with the ground state and the first excited singlet state 1 1 A 2 reflect the amount of valence character in the 1 3 B 1 wave function. This state originates mainly from an n → 4s Rydberg excitation. As may be seen in Tables S18 and S19 of the SM, the coefficient of the n H → σ * excitation which is decisive for the magnitude of the spin-orbit interaction is larger for the DFT/MRCI methods compared to MR-MP2 and CASSCF. The size of the 1 3 A 1 |Ĥ SO z |3 1 A 2 SOME is controlled mainly by the weight of the (π 1 H−1 n 1 H π * 2 L ) configuration in the 3 1 A 2 wave function which constitutes a (n H → π * L ) single excitation with respect to the leading (π 1 H−1 π * L ) term of the 1 3 A 1 wave function.
Thiophene. The correspondence between the SOMEs obtained by the two flavours of DFT/MRCI is so close that the differences need not be discussed. Focusing instead on the MR-MP2 approaches, we see that the matrix elements involving π → σ * excitations are systematically smaller for the MR-MP2(HF) method than for MR-MP2(BH-LYP) ( Table 8). The largest SOME is obtained for the second triplet of A 2 symmetry with the ground state. The dominant configuration in the 2 3 A 2 wave function arises from the π H−1 → σ * L+10 excitation, but also substantial contributions from π H−1 to Rydberg excitations are recognised in the wave function expansions (Tables S30 and S31 of the SM). The π H−1 MO exhibits substantial electron density at the sulphur centre. This is true also for the σ * L+10 orbital. As a consequence, the SOC integral π H−1 |Ĥ SO z |σ * L+10 is large, 142.8 cm −1 (BH-LYP orbitals) and 105.2 cm −1 (HF orbitals), indicating that the electron density distributions of the σ * L+10 orbitals differ significantly between the HF and BH-LYP one-particle basis. Most importantly, the electron density at the sulphur centre is lower in the σ * L+10 HF orbital compared to the corresponding BH-LYP orbital. This is the reason for a lower value of π H−1 |Ĥ SO z |σ * L+10 in the case of HF orbitals. The MR-MP2 SOMEs agree better with the CASSCF values in general. The DFT/MRCI approaches tend to exhibit larger valence character in the energetically lower lying states (1 1 B 2 , 2 1 A 1 1 3 B 2 , 1 3 A 1 ) and more Rydberg character in the higher lying ones (2 1 B 2 , 3 1 A 1 , 2 3 A 1 ) (Tables S30 and S31 of the SM), thus explaining the shifts in the SOMEs. For example, the decrease of 2 3 A 2 |Ĥ SO z |2 1 A 1 and concomitant increase of 2 3 A 2 |Ĥ SO z |3 1 A 1 when going from DFT/MRCI to MR-MP2 are due to different valence-Rydberg mixings in the A 1 symmetric singlet states.
Pyranthione. Some of the states in pyranthione are coupled by the largest SOMEs presented in this paper (see Table 8). We should point out that a similarly strong coupling of the 1 3 A 2 state with the GS (127.4 cm −1 [9]) and 2 1 A 1 (126.1 cm −1 [9]) is predicted when HF/MRD-CI wave functions are employed, where MRD-CI stands for ab initio multi-reference singles and double excitation CI with individual configuration selection [71,72]. Very strong coupling could have been foreseen due to the observable absorption of the first triplet state. The correlation between the MR-MP2(BH-LYP) and DFT/MRCI SOMEs is very good, having one of the smallest values of NRMSD, only 3%. MR-MP2(HF) has SOMEs somewhat shifted toward lower values, with a relatively large maximum negative deviation of 17.3 cm −1 , but with still good NRMSD of 6%.
Dithiin. With regard to the treatment of SO interaction, dithiin is the most challenging molecule examined in the present study. The DFT/MRCI and DFT/MRCI-R SOMEs are so close that we refrain from discussing the DFT/MRCI-R values further. In contrast, the MR-MP2 methods show very large NRMSDs values for both variants. Inspection of Table S2 in the SM reveals that the maximum positive and negative deviations for the set of all presented SOMEs both come from dithiin. Large negative deviations arise from the coupling between the 1 3 A and the ground state (Table 8). On the opposite side, there is a stronger coupling of 2 3 A and the ground state in the case of MR-MP2. In Section 3.3.1.2, it was already noted that the energy separation of the 1 3 A and 2 3 A states is smaller in the MR-MP2 cases compared to DFT/MRCI (see also Table 6). Accordingly, the configuration mixing is more pronounced for the MR-MP2 methods (Table S39 in the SM). Due to the distortion from the planar symmetry, a strict characterisation of the MOs as σ , n or π is not possible anymore. This mixing of orbital characters, together with the two sulphur centres located next to each other, leads to very large spin-orbit coupling of all orbitals with electron densities at the sulphur centres. This applies, for example, to the π H−1 and π * L orbitals (see Figure S5, SM) so that the integral π H−1 |Ĥ SO z |π * L has a value of 143.8 cm −1 (BH-LYP orbitals) and 114.7 cm −1 (HF orbitals). In contrast, the π * L+4 orbital has very little electron density at the sulphur atoms and the integrals involving it have very modest values. π H |Ĥ SO z |π * L+4 amounts to 25.1 cm −1 for the BH-LYP orbitals and 18.5 cm −1 for the HF orbitals. Likewise, π H−2 |Ĥ SO z |π * L+4 is 34.5 cm −1 for the BH-LYP orbitals and 31.5 cm −1 for the HF orbitals. Many other configurations, even with small coefficients, contribute to the final values of the 1 3 A|Ĥ SO z |GS and 2 3 A|Ĥ SO z |GS matrix elements. Since the correspondence between the 1 3 A and 2 3 A states in the DFT/MRCI and MR-MP2 treatments is shadowed by their strong mixing, the correlation of their SOMEs leads to outliers on opposite sides (see Figures S2 and S3 of the SM), with the DFT/MRCI matrix elements larger for the coupling of 1 3 A and the GS and smaller in the case of 2 3 A.
The maximal positive deviations in the correlation of DFT/MRCI and MR-MP2 SOMEs of dithiin are introduced by the coupling between triplet states, 2 3 B and 3 3 B ( Table 8). The state vectors are dominated by the negative and positive linear combinations arising from the (π H−2 → π * L ) and (π H → σ * L+1 ) excitations. As the two leading configurations in 2 3 B and 3 3 B are either the same or differ by more than one excitation, they do not interact by the here employed effective one-electron, mean-field SO Hamiltonian. Thus, the SO coupling is achieved by the involvement of configurations with minor coefficients in the state vectors. The 2 3 B and 3 3 B states are no exception from the usual trend that the dominant configurations have smaller coefficients in the MR-MP2 wave functions compared to DFT/MRCI. In the present case, the substantial admixture of minor configurations results in a stronger SO coupling of the MR-MP2 wave functions. Especially important are excitations from the n H−4 orbital which are more pronounced in the 3 3 B state calculated by MR-MP2 (see Table S39 of the SM). The n H−4 orbital has most of its electron density located at the two sulphur centres (see Figure S5, SM) and SO coupling involving n H−4 is very strong. For example, the configurations obtained by the (π H−2 → π * L ) and (n H−4 → π * L ) excitations are strongly coupled via the n H−4 |Ĥ SO z |π H−2 integral which has values of 193.7 and 203.1 cm −1 for the BH-LYP and HF orbitals, respectively. Such large values of SOMEs between molecular orbitals are the reason why dithiin is a challenging case for the SO treatment, especially when comparing different theoretical methods, since all differences in calculated wave functions are amplified.
Bithiophene. Bithiophene has the smallest correlation coefficient of only 0.3280 for the MR-MP2(HF) and 0.5416 for the MR-MP2(BH-LYP), suggesting very poor correlation between MR-MP2 and DFT/MRCI SOMEs. Indeed, the graphs presented in Figures S2 and S3 of the SM show the occurrence of both positive and negative outliers. The outliers are in the region of stronger spinorbit interaction, suggesting the involvement of the πσ * states. Note that bithiophene, like dithiin, is not planar but exhibits a trans gauche conformation in the electronic ground state. The πσ * states, 3 3 A and 3 3 B, can therefore mix with the close-lying 4 3 A and 4 3 B states and thus also have some ππ * character. The mixing is stronger for the MR-MP2 triplet wave functions (compare Table S41 in the SM) so that the pair of SOMEs (x and y components) involving 3 3 B has larger values for DFT/MRCI (more πσ * character) than for MR-MP2. The 4 3 B has more ππ * character in the case of DFT/MRCI and thus lower matrix elements (Table 8). A similar observation is made for 3 3 A and 4 3 A. This rotation of characters between the closelying triplet states in MR-MP2 leads to correlated pairs of positive and negative outliers compared to DFT/MRCI. Again, DFT/MRCI and DFT/MRCI-R values agree nearly perfectly ( Figure S1 of the SM).
Nitromethane. The agreement between SOMEs is quite good among the DFT/MRCI and MR-MP2 treatments (Table 8). We ascribe the deviations for the corresponding CASSCF SOMEs to the insufficient active space employed in the augmented basis (see Section 3.3.1.2).
Nitrobenzene. Despite the substantially different wave function compositions of the 1 1 B 1 state, good agreement between the DFT/MRCI and MR-MP2 SOMEs are found. Maximum deviations are around 5 cm −1 and NRMSDs are approximately 4%. The correlation between DFT/MRCI and DFT/MRCI-R SOMEs is even better. Why is that? The reason is tracked back to the small size of the spin-orbit integrals involved in the coupling of the 1 1 B 1 state with the low-lying 3 A 1 , 3 A 2 , and 3 B 2 states. At the planar ground-state geometry, the coupling between 1 1 B 1 and 1 3 A 1 is El-Sayed forbidden because both states represent (π → π * ) excitations. The origin of the small values for 1 3 A 2 |Ĥ SO |1 1 B 1 and 1 3 B 2 |Ĥ SO |1 1 B 1 is more subtle. Their leading configurations differ by an El-Sayed allowed single excitation from π H to n H−3 ( 3 A 2 ) and n H−4 ( 3 B 2 ). However, the π H orbital has very little electron density on the nitro group (see Figure S6 in the SM), whereas n H−3 and n H−4 represent predominantly the negative and positive linear combinations of the oxygen lonepair orbitals. In addition, the n H−4 MO exhibits density along the C-N σ bond. Because of its r −3 dependence, the spin-orbit operator is near-sighted and has large matrix elements only for electronic wave functions located at the same centre. As a consequence, the spin-orbit integrals π H |Ĥ SO |n H−3 and π H |Ĥ SO |n H−4 have very small values. The photophysically relevant SOMEs of the 1 3 B 1 and 1 3 A 2 states with the electronic ground state as well as the coupling between the first excited singlet state (1 1 A 2 ) and 1 3 B 1 are free from complications due to four open-shell configurations. Hence, in nitrobenzene, the SOMEs are seen to be more robust with regard to the wave function quality than the energies. This observation is confirmed by good agreement between the DFT/MRCI and MR-MP2 SOMEs. Maximum deviations are around 5 cm −1 and NRMSDs are 4%.
Methionine. The large variation of the MR-MP2 energies with regard to the employed type of MOs seems to have no effect on the SO interaction since the net coupling between the states is almost equal for the MR-MP2(HF) and MR-MP2(BH-LYP) wave functions (Table 8). While the SOMEs of states connected by an (n → π * ) excitation, 2 3 A|Ĥ SO x |1 1 A for example, are nearly equal for all methods, we observe systematically lower MR-MP2 SOMEs for states that are coupled by a π → σ * excitation such as 1 3 A|Ĥ SO y |1 1 A . Like in thiophene, DFT/MRCI and DFT/MRCI-R wave functions of the lower lying excited states tend to have more valence and less Rydberg character than the first-order perturbed wave functions of the MR-MP2 methods. Overall, the agreement between the DFT/MRCI and MR-MP2 SOMEs is satisfactory, especially considering the difference in energies between DFT/MRCI and MR-MP2(HF), with NRMSD of 7% and 5% for MR-MP2(HF) and MR-MP2(BH-LYP), respectively.
Isoalloxazine. The agreement between the DFT/MRCI and DFT/MRCI-R SOMEs (RMSD of 1 cm −1 , NRMSD of 5%) is satisfactory, considering the high density of states and the medium size of the selected SOMEs (ranging between 2 and 23 cm −1 ). In contrast, the correlation of DFT/MRCI and MR-MP2 SOMEs is rather poor (see graphs in Figures S2 and S3 of the SM). We could not find any pattern in the scattering depicted in the graphs, as is it not caused by a simple 2 × 2 rotation of wave functions. Examination of the DFT/MRCI and MR-MP2 state vectors (Table S57 in the SM) reveals that the nπ * states are heavy mixtures of excitations originating from lone-pair orbitals centred on N and O. Since the molecular orbitals located at O contribute with larger spin-orbit integrals, the size of the computed SOME is very sensitive with respect to wave function composition. In an effort to improve the results, we expanded the initial RAS space to 16-11-2 to include all the leading orbitals in the desired states. To this end, the number of excited states was reduced to five per irreducible representation, but Esel still could not be lowered. The computed excitation energies differ by up to 0.2 eV from the results of the corresponding calculations with an initial 10-10-2 RAS, but the changes of the SOMEs are relatively small (compare Table 8 and Table S58 in the SM). Thus, the bad correlation between DFT/MRCI and MR-MP2 SOMEs is retained. After all, it will be difficult to carry out MR-MP2 calculations on isoalloxazine that are converged with respect to the choice of first-order interacting space.

Conclusion
In this work, we have presented a thorough comparison of electronic excitation energies and spin-orbit matrix elements (SOMEs) of a representative set of molecules, obtained with different electronic structure methods. Among them are two variants of the DFT/MRCI method, the original parameterisation by Grimme and Waletzke [1] and the newer parameterisation of a redesigned DFT/MRCI-R Hamiltonian by Lyskov et al. [16]. While SOMEs of DFT/MRCI wave functions have been successfully employed in many photophysical studies, SOMEs of DFT/MRCI-R wave functions are presented here for the first time. Furthermore, we have tested the performance of the MR-MP2 method [33] with respect to spin-orbit coupling. Herein, a truncated first-order perturbed wave function is employed where the inclusion of configurations in the first-order interacting space is controlled by a selection threshold. Two variants have been employed, MR-MP2(HF) based on Hartree-Fock orbitals and MR-MP2(BH-LYP) based on Kohn-Sham orbitals using the BH-LYP density functional.
For the diatomic molecules, our quantum chemical results could be validated with respect to experimental data. For the polyatomic molecules, this was possible only to a limited extent. In particular, experimental data on optically dark states are missing for comparison. For some of the experimentally unknown states and their coupling matrix elements, we therefore carried out reference calculations at the CASSCF/CASPT2 levels.
Overall, we find very satisfactory agreement between the excitation energies and the SOMEs obtained with the four approaches. The correlation between the experimentally known adiabatic excitation energies of the first 3 state and the two variants of DFT/MRCI is very good (correlation coefficient above 0.98). Both MR-MP2 approaches seem to underestimate these energies somewhat but the correlation is still considered good (correlation coefficient above 0.96). The agreement between the experimentally derived fine-structure splitting constant A SO and the computed SOME is excellent (correlation coefficient >0.997) for all four approaches. For the experimentally known excitation energies of the polyatomic test molecules, correlation coefficients of 0.9806 (MR-MP2(BH-LYP)), 0.9781 (MR-MP2(HF)), 0.9883 (DFT/MRCI), and 0.9910 (DFT/MRCI-R) are found. This correlation includes both valence and Rydberg excitations. With regard to the uncertainties arising from a spread of experimental values, the influence of solvent effects, and the more systematic question how well vertical excitation energies can be compared with peak maxima of absorption spectra, we are very satisfied with this correspondence.
Since our comparison of SOMEs also involves higher lying electronically excited states, the agreement among the four approaches deteriorates somewhat, with the correlation coefficients still extending 0.95. The correlation plots show some pronounced outliers that can be grouped into several types. The most innocuous one arises from rotations among the wave functions of close-lying electronic states in molecules with low spatial symmetry. These types of outliers are characterised by a good correlation of the corresponding excitation energies. An example is the 3 3 A and 4 3 A pair of states of bithiophene. The different mixing of πσ * and ππ * in MR-MP2 as compared to DFT/MRCI leads to pairs of positive and negative outliers. Dithiin is a similar case with σ − π mixing. Isoalloxazine also falls into this group although it is planar. Here, several nπ * states are involved where the contribution of the lone-pair excitations on the oxygen and nitrogen centres to the wave function varies substantially among the states. In several thio compounds (thioformaldehyde, thiophene, pyranthione, methionine), we find systematically lower SOMEs between (π → σ * ) and (π → π * ) excited states in the MR-MP2 approaches compared to DFT/MRCR and DFT/MRCI-R. The analysis is complicated by the fact that the first-order perturbed wave functions of the MR-MP2 expansions often do not exhibit a few dominant terms plus many terms with small coefficients. In contrast to DFT/MRCI, the major wave function contributions are spread over many configurations with medium-sized coefficients. Since the deviation is more pronounced when HF orbitals are employed, we ascribe the main difference in the SOMEs to an orbital effect. The density distribution in the lowest lying σ * HF MO is more diffuse than in the corresponding KS MO.
Other reasons for large outliers are more serious as they involve the appearance of doubly excited configurations in the wave functions of (low-lying) electronic states. The large outliers observed in the correlation of the SOMEs of DFT/MRCI-R, on the one hand, and DFT/MRCI, on the other hand, arise from the interaction of the 2 1 A 1 and 5 1 A 1 states in thioformaldehyde with the 3 A 2 state that stems from a (n → π * ) single excitation. In particular, a double excitation of the n 2 → π * 2 type gains a much larger weight in the DFT/MRCI-R wave function of 2 1 A 1 so that this becomes the leading configuration, whereas a (π → π * ) excitation dominates the 2 1 A 1 wave functions of the other methods. Comparison with CASPT2 energies point toward a problem of the redesigned DFT/MRCI-R in properly describing this closed-shell double excitation of (n 2 → π * 2 ) type. For such a small molecule, also the MR-MP2 method represents a good alternative because the reference and first-order interacting spaces can be converged.
Some of the different wave function compositions are not reflected in the SOMEs. In o-benzyne and nitrobenzene, for example, large negative singlet-triplet splittings (meaning that among states with equal spatial composition, the singlet has a substantially lower energy) point towards an unbalanced description of doubly excited configurations with four open shells by the original DFT/MRCI method. This problem is remedied by the redesigned DFT/MRCI-R Hamiltonian. However, since the spin-orbit integrals are close to zero, the changes in the wave function have almost no effect on the SOMEs. Closing the discussion on doubly excited states, a comparison of the higher lying electronic states of dithiosuccinimide suggests that it might be worthwhile to reinvestigate the performance of the DFT/MRCI-R method on n, n → π * π * and n 2 → π * 2 double excitations. The difference between the two parameterisations arises predominantly from two integrals, i.e. o 1 o 2 ||o 1 o 2 and v 1 v 2 ||v 1 v 2 where o denotes an occupied and v a vacant orbital. When validating the DFT/MRCI-R approach, we computed energy profiles for a simultaneous rotation of the CH 2 end groups about the carbon-carbon double bonds of s-trans-butadiene [16]. Using this example, it was shown by us that such types of integrals play an important role for (π 1 π 2 → π * 1 π * 2 ) excitations, causing a large change of the adiabatic potential energy surfaces of the first excited states. Moreover, the good agreement with the present MR-MP2 results for the 2 1 B 1 state of obenzyne corroborates our conclusions that DFT/MRCI-R describes the interaction of π electrons properly. In cases in which electrons are annihilated from n orbitals (which are in general much more compact than π orbitals) and promoted to π orbitals, the change of correlation energy associated with the electron motion becomes more significant compared to a (π 1 π 2 → π * 1 π * 2 ) transition. Thus, the doubtful states in thioformaldehyde and dithiosuccinimide can be treated more accurately by addressing the correlation-exchange potential of the underlying density functional.
In cases of doubt, the use of MR-MP2 is a good alternative. Herein, the choice of orbitals seems to play a minor role. On the average, MR-MP2 energies are in slightly better agreement with experiment when Hartree-Fock orbitals are employed instead of Kohn-Sham (BH-LYP) orbitals, whereas the opposite is found for SOMEs. In connection with MR-MP2 calculations, it is recommended to solve for more roots than the minimum number, thereby improving the first-order interacting space. Otherwise, states with large perturbation corrections are easily missed, in analogy to CASSCF/CASPT2 treatments. For the photophysically important electronic states, SOMEs agree well with the corresponding coupling matrix elements from DFT/MRCI and DFT/MRCI-R calculations. The mapping of states among the methods becomes more and more difficult as we move up in energy. Concomitantly, the matching of the SOMEs deteriorates in these cases. Since MR-MP2 is computationally substantially more demanding than DFT/MRCI-R, the latter method is generally preferable. However, if doubly excited configurations appear in the wave function with large weights, we recommend the use of the MR-MP2 approach as a control that the semi-empirical DFT/MRCI-R approach is functioning properly.

Disclosure statement
No potential conflict of interest was reported by the authors.