Accurate prediction of core-level spectra of radicals at density functional theory cost via square gradient minimization and recoupling of mixed configurations

Linear-response time dependent density functional theory (TDDFT)^1–3 is very widely used to model electronic excited states of chemical species. TDDFT is an appealing approach as it is computationally inexpensive [O(N³⁻⁴) scaling vs number of basis functions N], nearly black-box, and able to simultaneously compute a large number of excited states. However, the lack of explicit orbital relaxation renders it unsuitable for describing excitations that involve substantial reorganization of electron density, such as charge transfer^3,4 or Rydberg excited states.^5,6 The excitation of core electrons, in particular, involves a substantial relaxation of the core-hole (and an accompanying reorganization of valence electron density), which leads to substantial errors in excitation energies predicted by TDDFT with standard functionals. It is consequently not unusual to blueshift TDDFT core-level spectra by ∼10 eV for alignment with experiment^7–12 (though the qualitative nature of transitions is typically reasonably predicted). Some specialized short-range corrected functionals specifically trained to predict core-level spectra¹³ tend to fare better,^14–20 but the very strong sensitivity of TDDFT excitation energies on delocalization error^3,21 is troubling (as even small perturbations could have a disproportionate impact on relative peak positions).

In contrast, linear-response based wave function theories, such as equation of motion coupled cluster singles and doubles (EOM-CCSD),^22–29 tend to systematically overestimate core-excitation energies^25,30–32 due to the lack of explicit orbital relaxation, often necessitating empirical redshifting by 1–2 eV for alignment with experiment.^25,26,31,32 Encouragingly, however, the use of different core-valence separation³³ (CVS) schemes has been observed to reduce the magnitude of the shift required.^29,34,35 A flavor of second-order extended algebraic diagrammatic construction [ADC(2)-x;³⁶ specifically, CVS-ADC(2)-x⁸] is also often employed to calculate core-level spectra.^8,37–39 The accuracy of CVS-ADC(2)-x, however, owes a great deal to fortuitous cancellation between various sources of error,^8,40 and the performance actually worsens when third-order ADC is employed.⁴⁰ At any rate, the higher computational cost of these wave function theories [O(N⁶) for both EOM-CCSD and ADC(2)-x] and slower basis set convergence renders them impractical for large molecular systems or extended materials, relative to computationally inexpensive DFT approaches. Nonetheless, the development of lower-scaling approximations to these wave function based methods is expected to broaden their applicability considerably.^25,41

In contrast to these linear-response based protocols, state-specific orbital optimized (OO) methods have been much more successful for accurate prediction of core-level spectra even within the DFT paradigm.^42–46 The main difficulty with these methods is the potential for “variational collapse” of the target excited state down to the ground state or another excited state, as it is challenging to optimize excited state orbitals (by virtue of excited states typically being saddle points of energy). The maximum overlap method (MOM)^47,48 was developed to address this problem for repeated Fock matrix diagonalization based methods like DIIS,⁴⁹ though convergence failures and variational collapse (via slow drifting of orbitals) are not always prevented.^50,51 More recently, some of us have proposed a square gradient minimization (SGM)⁵¹ based direct minimization approach that appears to be robust against both modes of MOM failure. SGM has been employed in conjunction with the spin-pure restricted open-shell Kohn–Sham (ROKS) method^52,53 to predict highly accurate (<0.5 eV error) core-level spectra of closed-shell molecules⁴⁵ at local DFT cost (using the modern SCAN⁵⁴ functional). It is also worth noting that there exist linear-response methods that incorporate partial OO character through relaxed core-ionized states, like Static Exchange (STEX)⁵⁵ or Non-orthogonal Configuration Interaction Singles (NOCIS),^56–58 though such treatments are wave function based and ∼1 eV error remains common due to the lack of dynamic correlation.

Stable open-shell molecules are fairly uncommon in nature, and there is consequently a scarcity of static experimental spectra for such species. However, open-shell systems are omnipresent in chemical dynamics experiments (either as fragments or excited states of closed-shell molecules) where transient x-ray absorption spectroscopy (XAS) is often employed.^59–62 It is consequently useful to have cheap and reliable theoretical techniques capable of modeling core-level spectra of such species. The highly accurate ROKS method is however not applicable to most open-shell systems, as it is explicitly designed for singlet states with one broken electron pair. In fact, open-shell systems pose additional challenges for many of the methods described above, as a spin-pure treatment of excited states necessitates inclusion of some double excitations^57,58,63 even for states that conventionally appear to be single excitations breaking one electron pair. This is not too difficult for wave function approaches, as shown by the extended CIS (XCIS⁶³) and open-shell NOCIS^57,58 methods. However, it is not at all straightforward to achieve this within TDDFT, which has no route for describing double excitations within the widely used adiabatic approximation.^3,64,65 It is tempting to believe that missing such configurations would not be particularly significant if the unpaired electrons interact only weakly, but the failure of TDDFT in describing excited state single bond dissociations despite the unrestricted reference state being reasonable⁶⁶ indicates some cause for caution.

In this work, we apply OO excited state DFT in conjunction with SGM to study single core-excitations of open-shell systems. This entails investigation of excitations to both singly occupied levels (which can be well described by single determinants, in principle) and completely unoccupied levels (which result in intrinsically multiconfigurational states). We present a scheme for recoupling multiple configurations to obtain an approximate doublet state for the latter class of excitations and demonstrate the utility of this protocol by considering the C K-edge spectra of the allyl radical, the O K-edge of CO⁺, and the N K-edge of NO₂. We also discuss general principles for reliably using these techniques to predict core-excitation spectra. Overall, we demonstrate that highly accurate DFT results can be obtained via orbital optimization with the modern local SCAN functional at low computational cost, similar to behavior observed for closed-shell systems. Low error can also be achieved via cam-B3LYP, TPSS, and ωB97X-D3 functionals (albeit at a somewhat higher asymptotic cost for the hybrid functionals).

Excitations from the core to singly occupied molecular orbitals (SOMOs) of open-shell systems result in states representable via a single Slater determinant, as there is no change in the number of unpaired electrons. The simplest approach for modeling such states is Δ Self-Consistent Field (ΔSCF),^42,47,67,68 where the non-Aufbau solution to the Hartree–Fock⁶⁹ or Kohn–Sham⁷⁰ DFT equations is converged via an excited state solver like SGM or MOM. The resulting excited KS determinant would not necessarily be exactly orthogonal to the ground state determinant, but this is generally of little concern since KS determinants are fictitious entities useful for finding densities, and thus, there exists no requirement that ground and excited state determinants be orthogonal. Nonetheless, a significant (>0.1, for example) squared overlap between the ground and excited state configurations would be concerning, but we have not observed such occurrences in our investigations and do not believe them to be likely without at least partial variational collapse of the core-hole.

The principal dilemma for such states is choosing between spin-restricted and unrestricted orbitals for ΔSCF. Unrestricted orbitals are typically more suitable for DFT studies on open-shell systems, though some functionals are known to yield atypically unphysical behavior in certain limits away from equilibrium.⁷¹ On the other hand, restricted open-shell (RO) orbitals artificially enforce a spin-symmetry that does not exist in radicals. As will be shown later (in Table II), the use of unrestricted orbitals appears to systematically lower the core-excitation energies (via extra stabilization of the core-excited state relative to the ground state). The best functionals for predicting spectra of closed-shell species yield lower errors for radicals when unrestricted orbitals are employed, and we thus recommend the use of unrestricted orbitals over RO orbitals for radicals. RO orbitals however should be employed for closed-shell systems (via ROKS or related methods),⁴⁵ on account of the existence of spin-symmetry in such species.

Multiconfigurational DFT is a difficult challenge even outside the unique challenges of TDDFT for double excitations, as the Kohn–Sham (KS) exchange-correlation energy is defined for a single determinant reference. KS-DFT target states therefore should be single determinants, and directly recoupling them via configuration interaction (CI) would result in double counting of some electron–electron interactions through both the functional and the CI off-diagonal terms. This is quite undesirable, making modeling such states fairly challenging.

One very reasonable solution is to note that single determinants with both α and β unpaired electrons are mixtures of different spin-states, and the highest spin-state within that ensemble can be well approximated by a single determinant by merely making all unpaired spins point in the same direction. Approximate spin-projection (AP)⁷² can consequently be applied to remove this high spin contribution from a spin impure mixed determinant. This approach should be sufficient when there are only two eigenstates that significantly contribute to the mixed configuration, as is the case for single excitations out of closed-shell molecules (where only the singlet and triplet states contribute). ROKS in fact utilizes this very feature to ensure spin-purity. ROKS employs a mixed configuration that has one unpaired α spin and one unpaired β spin (which has energy E_M) and a triplet configuration that has both unpaired spins as α (which has energy E_T). The use of RO orbitals forces the mixed configuration to be exactly halfway between the singlet and triplet, indicating $EM=ES+ET2$⁠, where E_S is the true singlet energy. ROKS consequently optimizes the purified singlet energy E_S = 2E_M − E_T.

Things are however substantially more challenging for doublet states. A mixed configuration with two unpaired α electrons and one unpaired β electron is a mixture of three states—two doublets and a quartet. The quartet contribution can be easily removed using an AP protocol similar to ROKS, but disentangling the two doublet energies is nontrivial.

Looking at the pure wave function based CI approach however offers some hints as to how to proceed. If we consider restricted open-shell configurations with three unpaired electrons occupying three spin-restricted orbitals (labeled 1, 2, and 3, respectively), eight possible configurations exist. Spin-inversion symmetry in the absence of magnetic fields however indicates that only four provide unique information:

1. $Q=↑↑↑$⁠: All three spins are α. This is the pure quartet with energy E_Q.

2. $M1=↓↑↑$⁠: Only the spin at orbital 1 is β. This is a mixed configuration with energy $EM1=EQ+K12+K13$⁠, where K_pq is the exchange interaction $pq|qp$ between an electron in orbital p and another in orbital q. The inversion of the spin in orbital 1 relative to the quartet leads to a loss of exchange stabilization between this orbital and the other two, leading to the energy going up by K₁₂ + K₁₃.

3. $M2=↑↓↑$⁠: Only the spin at orbital 2 is β. Consequently, $EM2=EQ+K12+K23$⁠.

4. $M3=↑↑↓$⁠: Only the spin at orbital 3 is β. Consequently, $EM3=EQ+K13+K23$⁠.

Having the single determinant energies $EQ,EM1,EM2,EM3$ is sufficient to uniquely solve for the exchange interactions K_pq with $K12=EM1+EM2−EQ−EM32$ and so on. This is quite useful as the off-diagonal CI coupling elements are $MiHMj=−Kij$ from Slater–Condon rules for double excitations.⁶⁹ This indicates that the knowledge of the single determinant energies is sufficient for solving the CI problem. With this, we find the eigenvalues of H within the subspace spanned by $M1,2,3$ to be

The first eigenvalue corresponds to the quartet within the $MS=12$ subspace (which is a linear combination of all three configurations with equal weights). The other two correspond to the energies of the two possible doublet states.

We propose that the same approach be employed for recoupling DFT configurations, with the KS energies of configurations $M1,2,3$ being employed instead of the HF ones used in the wave function theory approach. The risk of double counting should be greatly reduced as the effective off-diagonal elements are found directly from the KS energies vs Slater–Condon rules. Indeed, the off-diagonal elements should no longer be viewed as exchange interactions but rather effective spin–spin coupling elements. The entire approach is basically equivalent to solving for the eigenstates of the effective Ising like Hamiltonian $H′=−2J12S→1⋅S→2−2J13S→1⋅S→3−2J23S→2⋅S→3$ for three interacting spins, where the couplings J_ij are obtained from DFT (and are equivalent to the exchange interactions K_ij if HF is used as the functional). Such approaches have been used within broken-symmetry DFT to calculate spin coupling constants of transition metal species with reasonable accuracy,^73–80 and it is hoped that similar behavior will transfer over. Furthermore, equivalent logic for the case of two unpaired spins yields ROKS, which is known to be quite accurate for singlet states with one broken electron pair.^51,53,81 These known instances of successful behavior encourage us to believe that this protocol is worth exploring. We also note that Eqs. (2) and (3) were reported in Ref. 68 without an explicit description of the derivation, but these have not been actually applied to core-level spectroscopy (or any excited state problem) to the best of our knowledge.

Having obtained E_2,3 as spin-purified energies, we next seek to determine how to obtain the optimal orbitals. It is tempting to directly optimize E_2,3 in a manner analogous to ROKS, but we have elected not to do so at present. This optimization is nontrivial due to the nonlinear nature of the energy expression (vs the simpler form for ROKS). In addition, the derived equation is only precisely true for restricted open-shell orbitals, while Sec. II A seems to suggest unrestricted orbitals are optimal. We therefore look to AP-ΔSCF^67,68 for singlet excited states for inspiration, where the mixed determinant and triplet determinants are individually optimized (resulting in two sets of orbitals) and the singlet energy is simply computed as 2E_M − E_T from the individually optimized energies, instead of optimizing a single set of orbitals as in ROKS. The resulting energies however are often not dramatically different from ROKS,⁵³ and so we choose to follow a similar protocol here to determine if there is sufficient utility in this route for recoupling mixed determinants to justify optimizing a single set of unrestricted orbitals for computing the doublet energies. We consequently optimize $Q$ and $M1,2,3$ individually and compute E_2,3 from those optimized energies.

One rather inconvenient detail is that individually optimized $M1,2,3$ configurations would thus not be strictly orthogonal to each other due to slight differences in the orbitals. However, we do not consider any non-orthogonality derived terms arising from mixed configurations, as the KS determinants are fictitious constructs. On a more practical note, we ensure low overlap via providing restricted open-shell quartet orbitals as the initial guess for SGM optimization of the mixed determinants. The initial guesses are thus orthogonal, and orbital relaxation to the closest stationary point (which SGM is supposed to achieve) in unrestricted space should not lead to significant non-orthogonality for cases where this model of three unpaired electrons is a good approximation. Further details about initial guesses are enumerated in Sec. IV.

The magnitude of the transition dipole moment between the ground and excited states is essential for computing oscillator strengths (and thus relative intensities in computed spectra). The fictitious nature of the KS determinant (which represents a wave function of noninteracting electrons subjected to a fictitious potential) is a significant obstacle here as it implies there is no rigourous route for computing transition dipole moments. However, treating the KS determinants as real wave functions might be a reasonable approximation for computing this quantity, in the hope that the KS determinants (or superpositions thereof) would have a reasonably large overlap with the true wave functions to make this exercise worthwhile. Indeed, spectra computed via this route show fairly good agreement with experiment (as can be seen from previous work⁴⁵ by some of us, for instance). Such a protocol can (and should) account for non-orthogonality between ground and mixed determinants as it is fairly simple to compute NOCI dipole matrix elements.⁸² 

There are some additional factors to consider for the recoupled multiconfigurational states. The wave function inspired approach indicates that transition dipole moments should be computed via a linear combination of the transition dipole moments of individual determinants, as weighted by their coefficients in the eigenvectors corresponding to E_2,3. The effect of non-orthogonality between mixed determinants $M1,2,3$ on eigenvector coefficients is neglected here both because such terms are relatively small (because the mixed determinants have fairly low overlap with each other) and because it is not straightforward to calculate these effects. The decision to not consider this form of non-orthogonality does not appear to have any significant deleterious impact, as shown by the spectra presented later.

The other important factor to consider is that the analysis in Sec. II B found off-diagonal coupling elements directly from the energies $EM1,2,3$ and thus did not account for phases of $M1,2,3$⁠. These phases however are critical for estimating transition dipole moments and thus must be obtained somehow. A protocol for estimating these phases via the formally “quartet” state is supplied in the  Appendix.

The relative scarcity of experimental XAS data for radicals leaves us with a fairly small dataset of 17 excitations for assessing the performance of single determinant ΔSCF. The precise statistical values here are thus less reliable than those obtained in Ref. 45 from 40 excitations out of closed-shell molecules, but general qualitative trends can be drawn even from this restricted amount of data. The experimental excitation energies for all the C K-edge excitations (save allyl and CO⁺) were measured by some of us, via radicals obtained from the photodissociation of the corresponding iodide.^62,83 These values should have an uncertainty of ±0.1 eV, although vibrational excitations induced by photodissociation could shift the values somewhat. However, the resulting excitation energy for CH₃ agrees well with vibrationally resolved spectra obtained from radicals generated from flash pyrolysis.⁸⁴ Furthermore (as can be seen from Table I), the experimental shifts between the C K-edge of the allyl radical (obtained by authors of Ref. 85 on cold radicals generated via flash pyrolysis) and other C K-edges are very well reproduced by theoretical methods, suggesting that any vibrational excitation induced effect was small overall. A full Frank–Condon analysis could prove useful in quantifying any such effect but was not pursued at present.

We only consider a relatively small number of density functionals as a combination of large experimental uncertainty (typically 0.1 eV) and limited number of data points would make precise rankings of many functionals meaningless. We think it is more useful to investigate the performance of some representative functionals and see if they are sufficiently accurate to justify wider use. We therefore consider the following functionals from various rungs of Jacob’s ladder:²¹ 

1. Rung 1 (local spin-density approximation/LSDA^86–88): not considered due to very large errors found in Ref. 45.

2. Rung 2 (generalized gradient approximation/GGA): BLYP^89,90 and PBE.⁹¹ 

3. Rung 3 (meta-GGA): TPSS⁹² and SCAN.⁵⁴ 

4. Rung 4 (hybrids): B3LYP⁹³ and PBE0⁹⁴ (global hybrids); cam-B3LYP,⁹⁵ ωB97X-D3,⁹⁶ and ωB97X-V⁹⁷ (range separated hybrids).

5. Rung 5 (double hybrids): not considered due to significant computational expense.

The Hartree–Fock (HF) wave function method is also considered in order to determine the impact of neglecting correlation entirely. The choice of functionals here was guided by both a desire to compare against closed-shell results reported earlier⁴⁵ and a wish to examine the behavior of classic, minimally parameterized functionals like B3LYP.

Table I presents the excitation energies calculated using the chosen approaches (using spin-unrestricted orbitals), along with statistical measures of error like the root mean squared error (RMSE). None of the density functional methods deviate from experiment by more than 1 eV, which is in sharp contrast to the typical behavior of TDDFT with the same functionals.⁴⁵ Even HF has only <2 eV error despite complete absence of correlation. We specifically observe that the BLYP, TPSS, SCAN, B3LYP, cam-B3LYP, and ωB97X-D3 functionals yield 0.3 eV or lower RMSE and do not deviate by more than 0.5 eV from the experimental reference values. The good performance by local functionals like BLYP, TPSS, and SCAN is quite impressive as these functionals are much more computationally efficient than hybrids. Of the trio, the performance of only SCAN has been characterized for closed-shell systems,⁴⁵ where it was also found to be similarly accurate. We consequently focus on the performance of SCAN in Secs. III B–III D of this work as good performance in both the closed and open-shell limits is critical for prediction of transient x-ray absorption spectroscopy. However, we believe that good performance can be obtained from many functionals considered in this work (as partially demonstrated in the supplementary material). Interestingly, PBE and PBE0 perform surprisingly poorly, especially relative to BLYP and B3LYP, respectively.

Table I furthermore shows that the small errors for many functionals are mostly systematic, which appears to suggest that the change in excitation energy between two species (say between methyl and tert-butyl, for instance) would be reproduced fairly accurately by most functionals. This is also in principle true for TDDFT, although the massive (∼10 eV) errors in the individual excitation energies mean that even a relatively small variation in absolute error could have significant impact on relative peak positions (made more likely by very high sensitivity of TDDFT results to delocalization error^3,21). Most functionals (including SCAN) appear to systematically overestimate energies, while PBE and PBE0 systematically underestimate (which might be the reason for their poor overall performance). The inclusion of relativistic effects¹⁰⁹ (which systematically increase excitation energies by binding core electrons more tightly) would therefore degrade the performance of many functionals, while improving the performance of PBE and PBE0. The atom specific relativistic corrections for C, N, and O are however quite small¹⁰⁹ (0.1 eV–0.3 eV) and therefore are often neglected in studies (such as by the SRC functionals trained for TDDFT spectrum prediction,¹³ which has these effects implicitly baked into what is fundamentally a nonrelativistic theory). The impact of incorporating these corrections on the errors of various models is provided in the supplementary material, which shows that the RMSE of functionals (other than PBE and PBE0) goes up by 0.1 eV at most, suggesting that this is not a major issue in practice. We also note that HF systematically overestimates excitation energies by ∼1 eV due to missing correlation, which indicates that simple models for dynamical correlation (such as perturbative approaches^69,110) might be adequate for substantially lowering error, albeit at higher computational cost than DFT. HF however has a strong propensity to spuriously spin-contaminate Slater determinants, and the performance of perturbative corrections to HF references could consequently be greatly degraded.^110,111

We also consider whether there is any benefit of using restricted open-shell orbitals over unrestricted orbitals. Table II indicates that the use of unrestricted orbitals systematically lowers excitation energies by ∼0.1 eV relative to restricted open-shell results. This consequently indicates that the use of RO orbitals instead of U would degrade the performance of most of the studied functionals (as they systematically overestimate with U orbitals) and improve the behavior of PBE and PBE0. Indeed, Table II shows that both RO-PBE0 and RO-SCAN have the same RMSE of 0.4 eV. This potentially argues that RO-PBE0 is perhaps preferable to USCAN as the small relativistic corrections further improve the RO-PBE0 RMSE to 0.2 eV (while degrading USCAN’s RMSE to 0.4 eV, as shown in the supplementary material). However, we believe that SCAN with unrestricted orbitals is still the preferred route, even aside from greater asymptotic computational efficiency. Open-shell systems tend to often arise in transient absorption experiments starting from closed-shell species, and so it is important to use an approach that is effective at predicting the spectra for both types of systems. PBE0 is perceptibly inferior to SCAN when it comes to closed-shell systems⁴⁵ (irrespective of inclusion of relativistic effects), and the two are fairly close in predictive ability for open-shell systems, making SCAN with unrestricted orbitals the preferred choice. We also note that a comparison between aug-cc-pCVTZ and aug-cc-pCVQZ results shows that a small part (∼0.1 eV) of the systematic overestimation predicted by SCAN for Table I values stems from basis set incompleteness (as shown by a comparison in the supplementary material), similar to behavior of closed-shell species.⁴⁵ 

The spectrum of $NH2+$ has been experimentally characterized¹⁰¹ but was not considered in Table I as the two possible excitations to singly occupied levels are unresolved experimentally (assuming the radical cation is in the ³B₁ ground state). We used ΔSCF to compute the two transitions separately and report them in Table III. These transitions have nearly the same oscillator strength, and thus, their average should roughly correspond to the experimental peak. The ROKS results for the lowest lying singlet ¹A₁ excited state are also reported, in case it contributes to the experimental spectrum as well. Figure 1 presents the representative case of the SCAN functional with other methods yielding similar spectra.

The computed average triplet excitation energies in Table I agree fairly well with experiment, especially for good performers like SCAN, B3LYP, or ωB97X-D3. However, the values are somewhat redshifted, in stark contrast to the general behavior seen in Table I. One possible explanation for this would be blueshifting of the experimental spectrum due to the presence of singlet $NH2+$⁠, since this state absorbs fairly strongly (roughly twice the oscillator strength than the individual triplet transitions) at a slightly higher energy than the triplet, pushing the overall center of the band to higher energies (as hinted by Fig. 1). However, the computed triplet excitation average and the experimental maximum are not too far from each other (0.5 eV for SCAN), so it is not entirely impossible for DFT errors to be the sole reason behind the discrepancies. ωB97X-V, for instance, gives quite good agreement with experiment, without needing to invoke the singlet state.

Having explored the utility of ΔSCF in predicting excitation energies to the SOMO, we next seek to investigate the utility of the theory described in Secs. II B and II C for predicting the full core-excitation spectrum. The recoupling approach described therein is expected to be most effective for excitations to unoccupied valence orbitals, as then all three unpaired spins (in the core, SOMO and valence excited levels) will be interacting strongly. The scarcity of experimental spectra to compare against is again a problem and restricts us to only a few data points. Fortunately, the allyl radical has an experimentally characterized spectrum⁸⁵ that is dominated by excitations to the unoccupied π^* LUMO orbital, making it an excellent example for determining the utility of our recoupling approach, relative to simply using mixed configurations alone.

Figure 2 compares the performance of the orbital optimized methods in reproducing the C K-edge spectrum of the allyl radical. The performance of fc-CVS-EOM-CCSD and TDDFT with the specialized short-range corrected SRC2-R1¹³ functional is also considered. All three DFT methods are reasonable at predicting the lowest energy allowed excitation (from the terminal C atoms to the SOMO, the corresponding transition from the central C atom being symmetry forbidden), though all systematically overestimate by approximately 0.5 eV, resulting in the computed peak aligning with the vibrational fine structure of the experimental band. This is potentially indicative of some multireference character of this excited state, though it is difficult to draw firm conclusions from density functional data alone (especially, since it is possible to get better agreement via a functional that systematically underestimates 1s → SOMO excitation energies, like PBE0). It is however worth noting that fc-CVS-EOM-CCSD is spot on for this excitation, without any need for empirical translation of spectrum [as can be seen from Fig. 2(d)].

Figure 2(c) lays bare the failure of TDDFT at predicting excitations to the LUMO as the peak positions are completely off. This is not a pecularity of the SRC2-R1 functional but rather a failure of the TDDFT family of methods, as translated TDDFT spectra from other functionals yield a similarly poor picture (as shown in the supplementary material). Figure 2(d) also shows that fc-CVS-EOM-CCSD is unable to yield a qualitatively better spectrum than TDDFT, further highlighting the inadequacies of linear-response methods for this system. It is somewhat interesting that the inclusion of double excitations in fc-CVS-EOM-CCSD did not lead to any significant improvement over TDDFT (which is restricted to single excitations alone). The qualitative failure of both linear-response methods is likely a consequence of both spin-contamination and lack of orbital relaxation. The explicit inclusion of triple excitations should ameliorate both issues, but the significant computational expense of full EOM-CCSDT would dramatically constrain practical use.

The SCAN based orbital optimized approaches fare better, with both spin-contaminated mixed determinants and the recoupling approach yielding roughly qualitatively correct behavior. However, Fig. 2(b) shows that the mixed determinant approach fails to accurately predict the energy of the higher energy central C to LUMO transition, underestimating it by an electronvolt. This substantially damages the quality of the predicted spectrum, by making this peak appear in an area where none are present experimentally.

The recoupling approach shifts this peak to the appropriate location and predicts a spectrum in excellent agreement with experiment [as can be seen from Fig. 2(a)]. Indeed, Table IV shows that the peaks predicted by recoupled SCAN agree better with experiment than MCSCF calculations reported in Ref. 85 (though not too much should be inferred from this single data point). This good performance is not unique to SCAN alone, as several other functionals yield similar spectra in both the recoupled and mixed regimes (as shown in the supplementary material). Specifically, we find that recoupled cam-B3LYP, PBE0, and TPSS give good predictions for the 1s → LUMO portion of the spectrum, while BLYP and PBE yield rather poor performance even after recoupling. SCAN and cam-B3LYP appear to give the best performance, while some of the higher energy peaks with PBE0 and TPSS are somewhat redshifted with respect to the experimental spectrum. This supports our decision of selecting SCAN as the principal functional for the manuscript, despite BLYP and TPSS having the same computational scaling and slightly lower RMSE for excitations to SOMO (as shown in Table I). The poor qualitative performance by BLYP and PBE also serves as a potential warning against attempting to use GGAs for prediction of core spectra, despite BLYP’s excellent behavior for excitations to SOMO. Ultimately, the recoupling scheme cannot correct for deficient physics in the mixed configuration energies and a poor choice of functional could lead to poor results. Nonetheless, it is encouraging to see that all “advanced” functionals (Rungs 3 and 4) tested yield a reasonable spectrum after recoupling.

Overall, this example seems to suggest that orbital optimized approaches have an edge over TDDFT/EOM-CCSD when it comes to predicting core-excitation spectra of radicals. Furthermore, recoupling spin-contaminated mixed configurations to yield approximate doublets appears to not degrade performance and leads to some improvements. The overall accuracy of recoupled SCAN in predicting the spectrum of allyl certainly appears to hint at the efficacy of using this approach for XAS studies of large carbon based polyradical systems, such as the ones that might arise in soot formation during combustion.¹¹² 

We next consider the rather challenging case of the CO⁺ radical cation, whose experimental spectrum has been characterized very recently.⁹⁸ We focus on the O K-edge as the two doublet states corresponding to the 1s → LUMO excitation are experimentally well resolved, unlike the C K-edge (where the vibrational fine structure of the lower energy excitation overlaps with the higher energy one).

Figure 3 presents the orbital optimized SCAN spectrum (both recoupled and mixed), along with those from translated TDDFT and fc-CVS-EOM-CCSD. There are three peaks in all cases: the 1s → SOMO excitation (lowest in energy) and the two doublets arising from 1s → LUMO excitations. We observe that the linear-response approaches yield a fairly poor picture. Both TD-SRC2-R1 [Fig. 3(c)] and EOM-CCSD [Fig. 3(d)] need to be redshifted by ∼2 eV to align the 1s → SOMO peak with experiment (vs the orbital optimized DFT spectra, which needs no such translation). The translated spectra are nonetheless greatly compressed relative to experiment, and the relative intensities of the two 1s → LUMO peaks are incorrect. This is not merely a consequence of spin-contamination as Fig. 3(b) shows that SCAN using mixed configurations does better at reproducing the overall shape of the spectrum, despite having quartet contamination as well. The lack of orbital relaxation thus appears to be the critical factor that compromises the performance of TDDFT and EOM-CCSD for this system.

Figure 3(b) however also shows that SCAN with mixed configurations has too small a spacing between the two 1s → LUMO doublets (the two highest energy peaks). Figure 3(a) demonstrates that recoupling fixes this problem (and correctly reduces the intensity of the highest energy peak), yielding a spectrum that is in decent agreement with experiment. The spacing between the two highest energy peaks remains somewhat small (2.8 eV) vs experiment (∼3.4 eV but the unresolved broadness of the experimental second peak makes this hard to pinpoint). Other DFT functionals similarly underestimate this splitting (to varying extents), while reproducing the general shape of the spectrum (as can be seen from the supplementary material). Nonetheless, it is undeniable that the spectrum quality is greatly improved by recoupling. We also note that the NOCIS method^56–58 (which performs linear-response atop orbitals relaxed for the core-ionized state and is spin-pure in a manner analogous to XCIS⁶³) yields spectra that are in excellent agreement with experiment (as shown in the supplementary material), further demonstrating the utility of orbital relaxation and configuration recoupling, in an unambiguous, wave function based manner. At any rate, the qualitative failure of TDDFT and EOM-CCSD seems to argue for the use of methods with explicit orbital relaxation and configuration recoupling (like the scheme presented here or NOCIS) for the computation of core-level spectra of open-shell systems, irrespective of whether the computed spectra is translated or not.

NO₂ is another rare open-shell system with a known experimental high resolution core-level spectrum,¹⁰² by virtue of being quite stable for a radical. It is isoelectronic with allyl, although the SOMO is not a π^* orbital (but is rather a σ orbital mostly localized on N). The spectrum is nonetheless dominated by the transitions to the SOMO and the π^* LUMO levels. However, some Rydberg states have also been characterized, indicating that it could serve as an example to demonstrate whether our approach is balanced at predicting both valence and Rydberg excitations simultaneously.

Figure 4 compares the experimental spectrum at the N K-edge with those predicted via DFT (employing the doubly augmented d-aug-cc-pCVTZ basis to properly converge Rydberg states). The valence regime spectrum in Fig. 4(a) shows that all methods get the qualitative form right, though the 1s to π^* LUMO transition is somewhat redshifted by all methods. The success of TDDFT here stands in contrast to the failure observed for the valence regime of the allyl radical, although the different symmetry of the SOMO (σ vs π^*) may contribute to this. Recoupled SCAN performs better than mixed configuration SCAN for the second excitation, by removing the quartet contribution to the energy. This blueshifts the 402.3 eV excitation energy predicted by the mixed configuration approach to 402.9 eV, which is much closer to the experimentally observed peak at 403.3 eV. This disagreement is not particularly small (and is in the opposite direction to the systematic overestimation exhibited by SCAN for excitations to the SOMO), but the recoupled DFT method gives best agreement with experiment.

The Rydberg regime depicted in Fig. 4(b) however shows somewhat surprising behavior. It was tempting to believe that the weak coupling between the excited electron and the other unpaired electrons would lead to good performance by all methods. However, TDDFT absolutely fails to reproduce the spectrum in this regime, significantly blueshifting the experimental peak at 408.9 eV–410.0 eV. On the other hand, the mixed configurations are quartet contaminated and are thus slightly redshifted from their optimal location. Our recoupling protocol eliminates this problem, giving excellent agreement with experiment. It is also worth noting that the recoupled approach appears to predict the shape of the curve better than individual mixed configurations, indicating that the protocol described in Sec. II C was reasonably effective. This is however ultimately only one data point, and the comparison against more high resolution experimental spectra would be useful in validating our observation. We therefore hope that spectra of more open-shell species in the Rydberg regime will be available in the near future. We note that high energy spectra for $N2+$ and CO⁺ have been very recently reported,^98,100 but the Rydberg region appears to also contain a large number of doubly excited states with significant multiconfigurational character (involving more than three orbitals), that DFT based methods are unlikely to successfully model. This is less likely to be the case for neutral species.

The proposed protocol for recoupling mixed configurations appears to yield improved agreement with experiment, relative to simply using the two individual mixed configurations that correspond to single excitations. Nonetheless, it entails individual optimization of four configurations per excitation (⁠$Q,M1,2,3$⁠) to get two doublet state energies. We subsequently recommend the following protocol for ensuring maximum agreement between these configurations and minimizing computational cost:

1. Optimize unrestricted KS ground state orbitals.

2. Use these orbitals as initial guesses to optimize RO orbitals for the ground state.

3. Using the RO ground state orbitals as the initial guess, optimize the RO orbitals for the core-ionized state via SGM. This decouples the relaxation of the core-hole from the rest of the computations.

4. Using the RO core-ionized orbitals as the initial guess, optimize RO orbitals corresponding to the desired quartet state with SGM. The core-ionized orbitals can thus be computed only once and repeatedly utilized for multiple excitations. Furthermore, the unoccupied orbitals for the core ionized state are much more representative of the optimized orbitals for the excited electron than canonical ground state orbitals.

5. Using the RO core-excited quartet orbitals as initial guesses, find the unrestricted orbitals for the quartet $Q$ and mixed configurations $M1,2,3$ with SGM.

Steps 1–3 also apply for excitations to the SOMO level, followed by the use of the RO core ionized orbitals to initialize the excited state optimization for the core to SOMO excited configuration. They also apply for the computation of core-excitations in closed-shell species via ROKS. We believe that the RO energies themselves are not particularly useful for radicals, but the RO orbitals act as useful intermediates to prevent the alpha and beta spatial orbitals from differing prior to the last optimization step (step 5). The RO orbital space in fact is much more tightly constrained, and SGM is faster at those optimizations in practice. Difficult convergence cases in general could also be addressed via converging to the same state with a different (ideally, cheaper) functional and using the resulting orbitals as initial guesses.

The following three additional points regarding orbital optimized core-excitation calculations in general (for both closed and open-shell systems) are worth noting as well:

1. The use of a localized core-hole is absolutely critical for systems where there are symmetry equivalent atoms (like the terminal carbons of allyl). Delocalized core-holes lead to substantial delocalization error^21,114 driven underestimation of energy, as shown in Ref. 45. The localization of core orbitals can be achieved via explicit localization, or via weak electric fields that break symmetry. The mixed basis strategy described in the next point also leads to core orbital localizing symmetry breaking.

2. It is absolutely essential to use at least a triple zeta level basis with split core functions (like cc-pCVTZ) at the local site of the core-excitation. The core-hole would otherwise not be able to adequately relax, and energies be systematically overestimated.⁴⁵ However, a smaller basis can be used for all other atoms, with cc-pVDZ being adequate in our experience⁴⁵ (though even smaller bases could potentially be fine). This mixed basis strategy helps bring down the computational cost considerably as well, as the overall computation cost is comparable to a double zeta basis DFT ground state calculation per iteration, though excited state orbital optimization does often require many more iterations than ground state computations.

3. Many core-excited states possess significant Rydberg character. A good description of these states necessitates the presence of diffuse functions in the basis, and even double augmentation is sometimes necessary [such as the NO₂ spectrum presented in Fig. 4, where singly augmented aug-cc-pCVTZ blueshifts the Rydberg peaks in Fig. 4(b) by 0.2 eV]. This is easily the most onerous basis set requirement for such calculations but is functionally unavoidable for any electronic structure method seeking a correct description of Rydberg states.

We have investigated orbital optimized density functional approaches to study core-excitation spectra of open-shell systems by employing the SGM approach for averting variational collapse. The lack of gas-phase experimental data proves to be a hindrance for assessing the performance of these methods, but existing data show encouraging behavior. We first find that several density functionals like SCAN, TPSS, BLYP, B3LYP, cam-B3LYP, and ωB97X-D3 can be employed to predict excitation energies corresponding to 1s to SOMO transitions in radicals, with RMSE at or below 0.3 eV. The 1s → SOMO transitions are, however, not very challenging excitations as they do not result in a change in the total number of unpaired electrons and thus can be well approximated by single Slater determinants.

Higher excitations entail breaking of electron pairs and thus are natively multiconfigurational. These states therefore cannot be described by single determinants, although somewhat reliable results can at times be obtained from symmetry broken mixed determinants in the limit of weak coupling between unpaired spins (analogous to how unrestricted HF/DFT being effective for single bond dissociations in closed-shell species). For a more general accuracy, we present a CI inspired approach for self-consistently recoupling these single determinant mixed configurations with unpaired spins to yield approximately spin-pure results corresponding to multiconfigurational doublet states. The performance of this approach is compared against that of using unrecoupled mixed determinants alone and TDDFT/fc-CVS-EOM-CCSD for the core-level spectra of the allyl radical and CO⁺ at the O K-edge. The N K-edge spectrum for NO₂ is also studied with both orbital optimized DFT and TDDFT. We find that the recoupling scheme leads to no degradation of performance and, in fact, consistently improves upon results obtained by merely using single mixed determinants (significantly so for the O K-edge of CO⁺). It is nonetheless worth appreciating that unrecoupled determinants often yield fairly reasonable answers by themselves, especially relative to TDDFT/EOM-CCSD for the allyl radical and the O K-edge of CO⁺. Our work therefore shows promise in using orbital optimized DFT approaches for predicting core-level spectra of radicals, where high accuracy can be obtained even from local functionals like SCAN, at low computational cost. Available evidence also appears to argue for recoupling mixed configurations, although this is roughly computationally twice as expensive (as four configurations need to be optimized as opposed to only two). The O K-edge of CO⁺ also seems to suggest that our recoupling scheme somewhat underestimates doublet–doublet splitting in the strong coupling limit. More experimental spectra for open-shell systems (involving transitions to unoccupied valence orbitals) would, however, be immensely useful in fully characterizing the limitations of the recoupling approach. We will consequently continue to attempt to validate this approach via comparison to experiment as new data arise.

In the future, we will also seek to develop approaches that optimize a single set of unrestricted orbitals for recoupling mixed configurations vs separately optimizing all four relevant states. This should reduce the computational cost of such calculations substantially and enhance their utility. It would also be useful to generalize the recoupling approach to higher spin states like triplets, where there are more spins to recouple and a correspondingly larger number of coupling constants. Work along these directions is presently in progress.

All calculations were performed with the Q-Chem 5.3¹¹⁵ package. Local exchange-correlation integrals were calculated over a radial grid with 99 points and an angular Lebedev grid with 590 points. Experimental geometries (from the NIST database¹¹⁶) were used whenever possible, with MP2¹¹⁰/cc-pVTZ¹⁰⁶ optimized geometries being employed in their absence. The plots labeled “mixed” only used the two mixed configurations corresponding to single excitations from the ground state, as the third configuration is technically a double excitation that would not usually be considered due to formally zero (and in practice, typically small) oscillator strength. All TDDFT calculations employed the Tamm–Dancoff approximation.^3,117–119

See the supplementary material for additional spectra for the allyl radical and CO⁺, raw data, and molecular geometries.

The data that support the findings of this study are available within the article and its supplementary material.

D.H., K.J.O., and M.H.-G. were supported by the Director, Office of Science, Office of Basic Energy Sciences, of the U.S. Department of Energy under Contract No. DE-AC02-05CH11231, through the Atomic, Molecular, and Optical Sciences Program of the Chemical Sciences Division of Lawrence Berkeley National Laboratory. E.A.H., Z.Y., and S.R.L. were funded by the Gas Phase Chemical Physics program, which operates under the same DOE-OS-BES contract number. We also thank Marta L. Vidal for helpful discussions regarding EOM-CCSD calculations.

The phase convention chosen for $M1,2,3$ in Sec. II B ensures that the off-diagonal elements of the coupling matrix are −J_ij, where the couplings J_ij are given by

However, the determinants $Mi′$ obtained from orbital optimization can differ from this ideal phase. Specifically, DFT can yield $Mi′=piMi$⁠, where p_i = ±1. This has no implication for the energies but will affect properties like the transition dipole moment for which the relative phases of the configurations matter (as these properties depend on off-diagonal elements and are computed from $Mi′$ vs the idealized $Mi$⁠).

The easiest route for phase finding seems to be via the quartet state, which has the eigenvector $131313T$ in the $Mi$ basis. This state should formally have zero transition dipole moment and thus could be employed to compute relative phases.

Specifically, let $Mi′$ have transition dipole moments $μ→i=0μ→^Mi′$ against the ground state determinant $0$⁠. Without loss of generality, we can set the phase p₁ of $M1′$ to 1 (as only relative phases matter). Then, the transition dipole moment of the ostensibly quartet state is $μ→Q=μ→1+p2μ→2+p3μ→33$⁠. Consequently, the signs of p_2,3 should be chosen to minimize this quantity. In practice, this protocol is often simplified on account of one of the three mixed determinants being a formal double excitation (⁠$M2$ if the orbitals are ordered by energy), which would typically have very low transition dipole moments (though generally not exactly zero on account of non-orthogonality between the ground and mixed determinant orbitals). The phase estimation problem here is thus often just finding whether p₃ (say) should be 1 or −1 for $μ→Q$ to be smallest.

In fact, this is essentially an internal consistency check for determining the impact of neglecting non-orthogonality between mixed determinants and the overall quality of the optimized orbitals as this “quartet” transition dipole moment $μ→Q$ should be at least an order of magnitude smaller (and hopefully even less) than the largest transition dipole moment corresponding to the two doublet states, after finding optimal phases. The oscillator strength scales as square of the transition dipole and thus any spurrious “quartet” peak stemming from neglect of non-orthogonality would be at least a hundred times weaker than the strongest doublet peak, and thus, the quality of the spectrum will be preserved.

As an example, let us consider the N 1s → π^* transition in NO₂. The orbital optimized determinants $M1,2,3′$ we obtained had transition dipole moments (after ignoring terms smaller than 10⁻⁴)

$μ→Q$ is minimized if p₃ = 1 as then the dipoles will mostly cancel each other.