A study of core-excited states of organic molecules computed with the generalized active space driven similarity renormalization group

X-ray absorption spectroscopy (XAS) is a valuable technique for probing the electronic and nuclear structure of molecules as it provides complementary information to other spectroscopies.^1,2 In XAS experiments, absorption of a photon excites a core electron to a valence or continuum orbital, providing information about unoccupied states and the local geometric structure. The development of x-ray transient absorption spectroscopy (XTAS)^3–5 in very recent years has allowed the high-resolution, time-resolved detection of molecules during chemical reactions. This new technique has been applied to investigate various processes, including bond dissociation,^6,7 ring-opening reactions,⁸ intersystem crossing,^9,10 and non-adiabatic dynamics.^11–14 

The rapid development of x-ray spectroscopies has increased the need for accurate and efficient electronic structure theories.^15,16 Many computational approaches have been developed for predicting core-excited states. Methods based on density functional theory^17–36 are widely used in the simulation of core-excited states due to their low computational cost. On the other hand, systematically improvable (and generally more expensive) wave function methods have been proposed.^10,37–65 The majority of these approaches are based on a single-reference formalism and often express core-excited states starting from the ground state wave function. As a consequence, when the initial state is not well-described by a single electron configuration (e.g., a molecule far from its equilibrium geometry), errors in the wave function propagate to the excited states, resulting in a poor description of potential energy surfaces. To overcome the limitations of single-reference approaches, many multireference (MR) methods have been developed for treating core-excited states, including multi-configurational self-consistent-field (MCSCF) theory,^66–69 MR algebraic diagrammatic construction,⁷⁰ MR configuration interaction,^71–73 and MR coupled cluster theories.^74–77 

We recently developed a state-specific MR approach to compute core-excited states using general active space self-consistent field (GASSCF) reference wave functions and treating dynamical correlation effects with the driven similarity renormalization group (DSRG).⁷⁸ By simulating the vibrationally resolved XAS of diatomic molecules using a full-valence active space, we have shown that this method can accurately describe core-excited states in the bond-dissociation region, even for systems like CO⁺ and $N2+$ that display significant spin coupling effects. Our work also showed the importance of treating dynamical electron correlation beyond second-order perturbation theory, considering third-order and nonperturbative truncation schemes that include single and double substitutions. In this work, we examine various ways to extend the applicability of the GASSCF-based DSRG approach to medium-sized molecules and evaluate transition intensities. First, due to limitations in the size of the generalized active space (GAS) that can be treated efficiently, we investigate how to extend computations to large molecules by approximating the full-valence active space. Second, we employ a state-averaged (SA) extension of the MR-DSRG⁷⁹ to evaluate intensities and deal with problematic conjugated organic systems, where both core and delocalized π orbitals may be near degenerate and can cause root flipping problems. To test the performance of this MR-DSRG approach, we use the x-ray absorption benchmark of organic molecule (XABOOM) set, an XAS benchmark set containing 40 molecules, most of which are organic compounds.⁸⁰ 

This article is structured as follows: in Sec. II, we briefly summarize the GAS and state-averaged MR-DSRG approaches. In Sec. III, we provide details of the computations performed in this study, including the approximations we applied for larger systems and the choice of the DSRG flow parameter used. In Sec. IV, we compare vertical transition energies and relative splittings computed with GAS-DSRG against those from the frozen-core core-valence separated equation-of-motion coupled cluster theory with singles and doubles (fc-CVS-EOM-CCSD)^51,81 and experiment. We also discuss several systems for which MR-DSRG predictions deviate significantly from the fc-CVS-EOM-CCSD ones. In Sec. V, we summarize our results and provide a perspective for future applications of the MR-DSRG to XTAS simulations.

This section summarizes the generalized active space multireference driven similarity renormalization group (GAS-DSRG) theory we introduced in our previous study on core-excited states of diatomic molecules.⁷⁸ The combination of the state-averaged GASSCF and DSRG theory allows us to account for static and dynamical electron correlation in the ground and core-excited states. We will focus on the main features of the state-averaged scheme and the evaluation of intensities, two new aspects introduced in this study. The details of the state-averaged implementation of DSRG can be found in Ref. 79. An outline of the state-averaged GAS-DSRG computational scheme is shown in Fig. 1. Compared to the state-specific scheme, which targets one state at the time, the state-averaged formalism can obtain a manifold of solutions in one computation, avoiding root flipping problems.

We first model the core-excited states using a set of state-averaged GASSCF zeroth-order reference states Ψ_α, α = 1, 2, …, m, which account for static correlation in the core and valence orbitals (with all states built from a unique set of orbitals). Each reference state is a generalized active space self-consistent field wave function of the form

where |Φ_μ⟩ is a determinant that satisfies the GAS restrictions and $Cαμ$ is the corresponding coefficient. The GAS approach partitions the orbitals into the core, active, and virtual subspaces as in CAS, but further splits the active space into multiple subspaces (GASn, with n = 1, 2, …). By imposing restrictions on the maximum/minimum number of electrons occupying each GASn space, specific electronic states can be targeted. For core-excited states, the GAS1 space comprises the core orbitals from which electrons are removed, and GAS2 include selected occupied and unoccupied valence orbitals. Note that we describe different electronic states with different groups of GASSCF solutions (herein referred to as ensembles). For example, the ground electronic state is described by a GAS reference, restricting GAS1 to be fully occupied. Other GAS subspaces can be introduced to include the effect of semi-core levels, an extension not explored in this work. In these GASSCF computations, orbital rotations between two different GAS subspaces are frozen to prevent the collapse of a core-excited state to the ground state, a widely used approximation.^68,82–84 Furthermore, we freeze the number of electrons within each GAS in both the ground and core-excited state calculations. An approach alternative to the one used in this work is to employ one set of reference GASSCF orbitals for both the ground and core-excited states. In our experience, we found that the resulting orbitals are inadequate to describe states with such different electronic structures. The resulting transition energies typically show absolute errors of the order of 2 eV or more and, therefore, are too inaccurate for simulating XAS.

After generating the GASSCF reference states, we include the missing dynamical electron correlation via the state-averaged multireference DSRG.⁷⁹ To this end, for each ensemble of states, we build an effective Hamiltonian $(H̄)$⁠, which is later diagonalized in the space of ensemble GAS determinants. This state-averaged formalism starts by considering an ensemble of GASCI states ({Ψ_α}) described by the density operator $(ρ̂)$⁠,

where m is the number of states, ω_α is the weight of each state, and $∑α=1mωα=1$⁠.

In the state-averaged version of the DSRG, we build an effective Hamiltonian via the following unitary transformation:

where the operator $Â=T̂−T̂†$ is anti-Hermitian and expressed in terms of a generalized form of the coupled cluster excitation operator $T̂$⁠. In the state-averaged DSRG, the goal of this unitary transformation is to fold dynamic correlation effects into $H̄$ and achieve a partial decoupling of the ensemble density matrix from the components that involve excited configurations. These components of $H̄$ are referred to as non-diagonal and are indicated with $H̄N$⁠. When this decoupling is achieved exactly (i.e., $H̄N=0$⁠), one can obtain the eigenvalues for a manifold of states by diagonalizing $H̄$ in the space of GAS determinants that enter each ensemble. To obtain $H̄$⁠, we employ a form of operator normal ordering in which expectation values computed with respect to the density matrix $ρ̂$ are zero.^85,86

To avoid numerical instabilities in solving for the condition, $H̄N=0$⁠, caused by small denominators (related to the intruder state problem), the DSRG amplitudes are obtained by solving the following set of many-body equations:

where the off-diagonal term $H̄N(s)$ is gradually driven to zero by the source operator, $R̂(s)$⁠, in such a way that when s → ∞, we recover the condition $H̄N=0$⁠. The source operator that enters in Eq. (4) depends on the number s, referred to as the flow parameter. This quantity controls the magnitude of the amplitudes that enter the operator $Â$⁠, and its presence imparts an s-dependence onto $H̄$⁠. The choice of the flow parameter used in the GAS-DSRG computations is discussed in Sec. III B. For the full definition of the source operator and the state-averaged DSRG equations, we refer the reader to Ref. 79. The only modification applied compared to the original state-averaged formalism is setting to zero amplitudes that correspond to excitations between different GAS subspaces.

Once the DSRG amplitudes are obtained, the state-averaged DSRG energy is computed by diagonalizing $H̄(s)$ in the space of GAS determinants

where $Eα′$ and $|Ψα′〉$⁠, with α = 1, …, m, are relaxed energies and references, respectively.

In this work, we focus mainly on perturbative approximations of the MR-DSRG^4,87 for core-excited states. The zeroth-order Hamiltonian is the diagonal component of the generalized Fock matrix, and we work in a semi-canonical basis so this operator is diagonal in the core, valence, and individual GAS subspaces. Truncation of the MR-DSRG using this partitioning leads to second and third-order DSRG multireference perturbation theories (DSRG-MRPTn, n = 2, 3).^4,87 The availability of efficient implementations of these DSRG perturbative methods⁸⁸ (and that require at most the three-body reduced density matrices of the GAS references) allows us to easily perform DSRG-MRPT2/3 calculations of multiple core-excited states of relatively large molecules such as DNA nucleobases. For the calculations on ozone, we also consider a nonperturbative approximation to the DSRG termed LDSRG(2), which truncates $Â$ to one- and two-body substitution operators.⁸⁹ Like the parent MR-DSRG theory based on CASSCF references, the GAS-DSRG approach is size-consistent and size-extensive, with the former condition being satisfied if the GAS reference is multiplicatively separable. The GAS-DSRG excitation energies are also size-intensive, i.e., the transition energies of a supermolecular system composed of noninteracting fragments are identical to those of the isolated fragments. We verified the size intensivity of GAS-DSRG at the level of second-order perturbation theory by considering the lowest O K-edge transition in CO. The energy and intensity of a supermolecular system containing CO and a Ne atom at a finite distance were found to be different from those of isolated CO, but converge to those for an isolated CO as the distance between these fragments goes to infinity.

In this study, we also extend our approach to evaluate the transition oscillator strength between ground and core-excited states. We approximate the transition oscillator strength using the relaxed GAS states, where the ground state wave function employs orbitals optimized for core-excited states. This approximation captures the bulk of the transition oscillator strength but neglects orbital relaxation effects in the ground state and contributions from dynamical correlation. In practice, for each ensemble of states, we perform a state-averaged GASSCF computation in which, in addition to the m target core-excited states, we include a state with fully occupied GAS1(core) and 0 weight (ω₀ = 0). We then approximate the transition dipole moment between the ground state (α = 0) and core-excited state α (d_0α) as

where $μ⃗pq=⟨ϕp|μ⃗|ϕq⟩$ is the molecular orbital dipole moment matrix and $⟨Ψα|âqp|Ψβ⟩$ is the transition one-body reduce density matrix. Note that transition energies are, instead, computed using a fully optimized ground state GAS reference.

To benchmark the accuracy of GAS-DSRG theories on core-excited states, we employ the XABOOM benchmark set,⁸⁰ which contains 1s → π* transitions for a total of 40 molecules of different sizes. We use molecular geometries optimized at the frozen-core MP2/cc-pVTZ level of theory (taken from Ref. 80). Vertical transition energies computed with GAS-DSRG are computed according to a procedure similar to our previous study.⁷⁸ First, ground-state molecular orbitals are obtained using restricted Hartree–Fock (RHF) and the cc-pVQZ basis.⁹⁰ These orbitals are then separately optimized using a state-averaged GASSCF for the ground and core-excited electronic states.

For larger systems, it is necessary to introduce approximations to reduce the time and memory cost of the GAS-DSRG calculations. To reduce storage costs, density fitting of the two-electron integrals was applied on molecules with more than five atoms, using the cc-pVQZ-JKFIT⁹¹ and cc-pVQZ-RIFIT⁹² auxiliary basis sets for state-averaged GASSCF and MR-DSRG, respectively. Furthermore, we did not account for relativistic effects. Following the XABOOM study,⁸⁰ we also used the cc-pVQZ basis, which lacks the flexibility to treat core-valence correlation.

As in our previous treatment,⁷⁸ the core orbitals pertinent to a target transition are included in GAS1, while the active valence orbitals are all in GAS2. To denote the choice of active space and the constraints on the number of electrons, a GAS occupation is denoted by $(no1o,ne1e;no2o,ne2e;⋯)$⁠, where $nom$ and $nem$ are the number of orbitals and electrons in the m-th GAS, respectively. For molecules with two heavy (non-hydrogen) atoms, GAS2 spans the full valence orbitals (2s + 2p). For molecules with three heavy atoms, only the 2p orbitals of the heavy atoms are included in GAS2, except for the highly symmetric O₃ and N₂O molecules for which full valence orbitals are used. For systems with four or more heavy atoms, GAS1 consists of the 1s core orbitals of the heavy atoms and GAS2 consists of (2p)π orbitals. The selection of the π orbitals is straightforward for molecules with C_s or higher symmetry. However, for systems with C₁ symmetry, especially the nucleobases with pseudo-planar equilibrium geometries, the π orbitals must be manually selected. For C₄F₆ in the benchmark set (C₁ symmetry), GAS2 contains only six orbitals around the LUMO/HOMO levels.

Most of the GAS-DSRG calculations are performed averaging over several electronic states. This is mainly to avoid variational collapse in the GASSCF optimization due to the small energy gap between core-excited states. Moreover, it is crucial to ensure a consistent determinantal composition of the electronic state(s) before and after diagonalization of the DSRG effective Hamiltonian. By default, all the N_c × N_π singly excited states, where N_c is the number of cores and N_π is the number of low-energy unoccupied π orbitals, are averaged. However, for some of the molecules, we observe root flipping after the DSRG perturbative treatment, especially at the DSRG-MRPT2 level, a problem we address by increasing the number of roots averaged. For larger systems with a high density of states (especially for the carbon and nitrogen K-edge calculations of nucleobases), we perform separate computations on atoms with different chemical environments (split-core method), whereby GAS1 only consists of one core orbital for each calculation. The active spaces and states averaged for all GAS-DSRG calculations are given in Table S1 of the supplementary material.

To assess the impact of these approximations and the choice of the active space, we ran benchmark computations on two systems, N₂O and O₃. The active spaces tested for O₃ are shown in Fig. 2. The “standard” choice of active spaces includes all the 1s core orbital in GAS1 and all the valence orbitals in GAS2, leading to a (3o, 5e; 12o, 19e) GAS. The “π only” (PO) treatment includes only π orbitals in GAS2 , reducing the GAS to the much smaller (3o, 5e; 3o, 5e). In both cases, a total of eight states are averaged. The “split core”(SC) strategy performs two different state-averaged calculations; the first one focuses on the central oxygen core excitations and uses the GAS (1o, 1e; 12o, 19e), averaging two states. The second one considers core excitations from the terminal oxygens and uses the GAS (2o, 3e; 12o, 19e), averaging over six states. We can also apply the PO and SC approximation together (PO + SC). Similar approximations are considered in the N K-edge calculation of N₂O in which the “standard” choice of active space is (2o, 3e; 12o, 17e). The transition energies and oscillator strengths of excitations from the 1s core orbital of the central atom (denoted by the subscript C) and the terminal atom (denoted by the subscript T) computed with different active spaces are given in Table I. We also reported the energy difference, ΔE = E_C − E_T, and the ratio of oscillator strengths, f_T/f_C, between the two transitions, which are perhaps better criteria for measuring the accuracy of core-excited state calculations.

Several observations can be made from this initial test. First, the effect of partitioning and excluding orbitals from the active space is much more substantial in O₃ than in N₂O. Under all approximations, ΔE values for N₂O deviate less than 0.3 eV from the “standard” results. For O₃, at the DSRG-MRPT2 level, the deviation from the “standard” results is already about 1 eV from both the PO and SC approximations, increasing significantly (to 4.14 eV) when PO and SC are combined. A similar observation can also be made for the oscillator strength ratio. This significant difference is most likely due to the orbital near-degeneracies in the O₃ molecule, which will be discussed in Sec. IV C 1. Our second observation is that the deviations of ΔE from the standard calculation at the DSRG-MRPT3 level are less than 0.40 eV for all approximations, an error significantly smaller than the aforementioned deviation in DSRG-MRPT2. In addition, compared to the PO approximation, the SC approximation has a stronger impact on the energy splitting between two core-excited states, though the number of determinants in the SC approximation is generally much closer to the standard one.

We tested the sensitivity of transition energies and dipole moments with respect to the DSRG flow parameter, s, for molecules with different sizes. For valence excitations, Wang et al.⁹³ conducted a study of 280 transitions and showed that the optimal choice of flow parameter is 0.35/2 $Eh−2$ for SA-DSRG-MRPT2/3. In this work, three molecules (ethene, benzene, and naphthalene) are selected as representatives for benchmarking the choice of s, as they share the same elements but vary in size. The active space for ethene [(2o, 3e; 12o, 11e)] consists of all the valence orbitals, while the active spaces for benzene [(6o, 11e; 6o, 7e)] and naphthalene [(10o, 19e; 10o, 11e)] only include core and π orbitals. The computed transition energies and oscillator strengths of the C 1s → π* transitions are shown in Table II. These were computed with different MR-DSRG treatments of dynamical correlation and flow parameter values in the range 0.25–4 $Eh−2$⁠.

The dependence of oscillator strengths on s is not very significant, but the dependence of the transition energies varies with system size. The transition energies for ethene vary by up to 0.6/0.08 eV for DSRG-MRPT2/3, indicating a more pronounced dependence on s than previously observed for $N2+$ (0.1 eV).⁷⁸ The variation in transition energy as a function of s is even larger for benzene (∼3.5 eV) and naphthalene (∼4 eV for three transitions). On the other hand, the relative splittings of three states of naphthalene are affected less by the s parameter, with the splitting between the highest and lowest transitions ranging from 0.90 to 1.35 eV at the DSRG-MRPT2 level and from 0.86 to 0.94 eV at the DSRG-MRPT3 level.

An analysis of these data shows that this strong s-dependence of the transition energies is mainly observed for small values of s (0.125–0.25 $Eh−2$⁠), when important DSRG amplitudes are suppressed and increases with system size due to the extensivity of dynamical correlation. We also note that the absolute energies of the ground and core-excited state show different dependences on s, with the latter requiring larger s values to converge. For example, at the DSRG-MRPT2 level, the difference between the 1s → π* core-excited state of benzene computed with s = 1 and $0.5Eh−2$ is 29.2 mE_h, but the same difference is only 4.3 mE_h for the ground state. In practice, we find that using s = 1 $Eh−2$ is a good compromise for the absolute transition energy for all our reported calculations. In benzene and naphthalene, this value leads to a shift in absolute transition energies no more than 0.3/0.1 eV for DSRG-MRPT2/3 calculations compared to s = 2 $Eh−2$⁠.

To investigate the accuracy of x-ray absorption transition energies predicted with the GAS-DSRG approach, we first compare the DSRG-MRPT2 and DSRG-MRPT3 results with these from fc-CVS-EOM-CCSD,⁵¹ the most accurate theoretical results in the XABOOM set.⁸⁰ We note, however, that works that have compared the performance of EOM-CCSD to higher-order coupled cluster methods⁵⁶ show that triple excitations can shift core-excitation energies by more than 1 eV and that quadruple excitations are required to achieve sub-eV accuracy. The fc-CVS-EOM-CCSD method applies the frozen-core and core-valence separation, which provides an efficient and robust method for evaluating core-excited state energies. This method uses coupled cluster amplitudes that exclude core orbitals (fc) and restricts the EOM operator to excitations that involve at least one core orbital (CVS). In comparison, in the current DSRG approach, we also enforce the core-valence separation (since the GASSCF solutions that represent core-excited states enforce occupation restrictions on the core orbitals), while we do not apply the frozen-core approximation since the core orbitals are treated either as active or restricted, and therefore, their occupation is modified by the operator $Â(s)$⁠. For simplicity, we denote the fc-CVS-EOM-CCSD method as EOM-CCSD in the following sections.

Table III reports error statistics for the relative energy difference between these theories for each of the C, N, and O K-edge excitation energies. As it can be seen, the mean signed deviations (MSD) (PT2 − CC/PT3 − CC) are all negative and relatively close to the mean absolute deviation (MAD) of each set. Figure 3(a) shows the distribution of computed transition energies differences with respect to EOM-CCSD excitation energies. For the C K-edge transitions, absolute energies computed with the DSRG-MRPT3 and EOM-CCSD show good agreement. Most of the transitions are less than 0.5 eV apart from each other, and the MAD between the two sets is 0.43 eV. The DSRG-MRPT2 results show a much larger deviation (MAD = 1.32 eV) from EOM-CCSD. On the other hand, the N and O K-edge show a somehow different behavior, where the deviation between the two DSRG theories is generally smaller than the deviation from EOM-CCSD. We also computed oscillator strengths, and plotted the distribution of the f_0α ratios between the three theories in Fig. 3(b). The transition oscillator strength ratio between EOM-CCSD and DSRG-MRPT2/3 can vary dramatically, likely due to our approximate way of computing f_0α within the state-averaged GAS framework. Computing transition dipole moments using a state-averaged GAS with a nonzero weight for the ground state may be a way to improve this approximation. However, preliminary computations show that the agreement with EOM-CCSD results does not improve systematically. Two other potential sources of error in our approximate oscillator strengths come from ignoring dynamical correlation effects and constraining the ground and excited states to use the same orbitals. Improving upon these approximations will be pursued in future research.

The theoretical transition energies are also compared with available experimental values,^13,94–117 which are absent from the XABOOM set. Transition energy distributions are shown in Fig. 4 and error statistics are given in Table III. The DSRG-MRPT3 results best match the experimental transitions, with the MAD from the experiment being only 0.32 eV, a much smaller value than for DSRG-MRPT2 (0.76 eV) and EOM-CCSD (0.70 eV). Moreover, DSRG-MRPT2 (MSE = −0.74 eV) and EOM-CCSD (MSE = 0.70 eV) are more likely to underestimate and overestimate transition energies, respectively. The overestimation of EOM-CCSD for valence transition energies has been observed previously,^118,119 and Fransson et al.⁸⁰ also reported that EOM-CCSD with relativistic corrections overestimates the CO C K-edge, HCN N K-edge, and CO O K-edge by 0.3, 0.7, and 1.2 eV when compared to the experiment.

The excellent agreement between DSRG-MRPT3 and experiment may disguise some error cancellation due to the complexity of comparing these x-ray absorption data, as it has been illustrated in the XABOOM study.⁸⁰ Some common assumptions made in the theoretical calculation of x-ray absorption spectra have been assessed in our previous study on diatomic molecules.⁷⁸ For example, the difference between the vertical and adiabatic transition energies of the CO 1s → π* excitation is 0.6 eV at the MR-LDSRG(2) truncation level. Scalar relativistic effects contribute to shifts in the order of 0.2 eV for C, N, and O K-edge transitions, and core-valence correlation basis sets can lead to absolute energy differences as large as 0.5 eV. The existence of all of these uncontrolled factors suggests comparing, instead, the energy difference between transitions, a topic discussed in Sec. IV B.

From our results, we evaluated the energy splittings between transitions for all the molecules for which we compute more than one core-excited state from any type of core atom. For molecules that have more than two transitions, we calculate all the splittings between adjacent transitions. A total of 47 energy splittings are obtained and the error statistics with respect to experimental and theoretical (EOM-CCSD) results are shown in Table IV. Due to the relatively small number of splittings for N/O K-edge transitions, we only report aggregated statistics. For each of these splittings, we also computed the ratio between their transition oscillator strengths. All splittings and oscillator strength ratios are given in Table S2 of the supplementary material.

Compared to the absolute transition energy, the energy splittings given in Table IV are less sensitive to different levels of theory. For the 34 splittings where an experimental value is available, the statistics show no significant difference among theories. This is partly due to the limited resolution of experiments, where some splittings are simply unresolved. As most systems in the XABOOM set have close-shell ground states, this relative good agreement is expected. However, there are several systems where large deviations in the energy splittings (⁠$>$0.40 eV) can still be observed among the theories. We will discuss these systems in Sec. IV C.

Our computations show significant discrepancies between the GAS-DSRG and EOM-CCSD predictions for ozone. The ground state of ozone is composed of two closed-shell electron configurations, the first one corresponding to the occupation $(1a1)2(2a1)2(1b2)2…(6a1)2(1a2)2$ and the second one corresponding to the doubly excited determinant $(1a2)−2(2b1)2$⁠, where we use a compact notation to express the occupation pattern with respect to the dominant electron configuration. Our GAS-DSRG computations show good agreement for the splittings between the two 1s → π* transitions (2a₁ → 2b₁ and 1a₁ → 2b₁), predicted to be 5.16 and 5.14 eV by DSRG-MRPT2/3. However, the splitting predicted by EOM-CCSD (∼4.5 eV) significantly underestimates the experimental value¹²³ (estimated to be 5.55 eV). This large deviation is likely due to the open-shell character of the ground state of ozone, which has been studied theoretically in the low-lying states,^124–133 core-excited states,^51,121 and core-ionized states.^70,121 To further analyze this discrepancy, we extended our GAS-DSRG computations to all transitions in the x-ray absorption spectrum of O₃ and compared them with those measured experimentally.^101,123,134 We focus, in particular, on the more recent studies of Rosenqvist et al.¹⁰¹ and Stranges et al.,¹²³ which have addressed deficiencies in earlier measurements.¹³⁴ We also applied the MR-LDSRG(2) method, the most accurate nonperturbative truncation scheme available. The active space used in GAS-DSRG calculations is the “standard” active space in Fig. 2.

According to the experimental assignment,^101,123 four core-to-valence transitions in the x-ray absorption region have significant oscillator strength: 1a₁ → 2b₁(π*), 2a₁ → 2b₁(π*), 1b₂ → 7a₁(σ*), and a weak 2a₁ → 7a₁(σ*) transition that overlaps with the 1b₂ → 7a₁(σ*) transition. Including these 1s → σ* transitions requires averaging 16 states. The experimental spectrum and theoretical simulations shifted to match the experiment are shown in Fig. 5. For comparison, we have also plotted existing theoretical XAS spectra of ozone from fc-CVS-EOM-CCSD,⁵¹ MS-RASPT2,¹²¹ and EOM-CC3.¹²⁰ We note that the work of Stranges et al. assigns excitation energies relative to the first line of 5.55 and 6.50 eV.¹²³ All the GAS-DSRG methods predict a similar triplet feature, with some differences in the relative shifts.

For absolute transition frequencies, it is not surprising to see that DSRG-MRPT2 requires a large blue shift (+1.3 eV) to match the experiment (a trend discussed in Sec. IV). However, a large red shift (−1.0 eV) is also necessary for DSRG-MRPT3, while MR-LDSRG(2) requires only a small shift (−0.3 eV). For all the DSRG simulated spectra, the strong transitions [1a₁ → 2b₁(π*) and 1b₂(2a₁) → 7a₁(σ*)] agree well with the experiment. The third 2a₁ → 2b₁(π*) transition is predicted to be 6.4–6.7 eV above the first line, in agreement with the experimental assignment (6.5 eV);¹²³ however, the intensity appears to be overestimated, as in the experimental spectrum this peak is not as prominent as in our simulations. There might be several reasons for the observed discrepancy between GAS-DSRG and the experimental spectrum; however, a more detailed investigation of this molecule is beyond the scope of this work. Going beyond EOM-CCSD, the EOM-CC3¹²⁰ predictions reproduce well the transition frequencies for all the peaks in the experiment and requires a very small shift (−0.1 eV). Recent MS-RASPT2 results¹²¹ yield a spectrum similar to the EOM-CC3 one; however, the splitting between the two lowest transitions appears to be overestimated (6.8 eV vs the experimental value 5.55 eV). A GAS-DSRG test calculation using the same active space and basis set used in the MS-RASPT2 study gives transition energies that are within 0.21 eV from our values. Despite the good agreement that can be reached with the excitation energies, none of the theoretical computations satisfactorily reproduce the intensity profile of the ozone XAS spectrum.

An analysis of the wave function of the core-excited states of ozone validates the need for a multireference treatment. Table V lists the leading contributions for the three core-excited states, 2a₁ → 2b₁(π*), 1a₁ → 2b₁(π*), and 1b₂ → 7a₁(σ*), computed at the MR-LDSRG(2) level of theory. The most important contribution to all three core-excited states is due to singly excited configurations with respect to the leading ground state determinant. However, these appear with a weight (54%–65%) smaller than the dominant ground state determinant (∼90%), indicating large configuration mixing in the core-excited states. The second largest contribution to the core-excited states is due to a single excitation out of the $(1a2)−2(2b1)2$ configuration, which appears with a large weight (∼16%–23%). The rest of the significant configurations are double excitations (one core-to-valence and one valence-to-valence), and add up to around 20% of the wave function probability. This significant multi-determinantal character of the core-excited states of ozone is consistent with our previous discussion in Sec. III A, where the 1s → π* transition frequencies showed a strong dependence on the choice of the GAS. All these observations suggest that modeling the core-excited state of ozone requires an accurate treatment of static and dynamical correlation.

Glyoxylic acid (HCOCOOH) is another system where we observe significant discrepancies between different theories. Specifically of interest are the two 4a′ → 4a″ and 5a′ → 4a″ carbon K-edge transitions, denoted by the orbitals involved. These orbitals are shown in Fig. 6, and the computed frequencies and oscillator strengths of the two transitions are given in Table VI. Unlike the case of ozone, we observe a significant underestimation of the splitting between transitions at the DSRG-MRPT2 level (1.63 eV for DSRG-MRPT2 vs ∼2.10 eV for DSRG-MRPT3 and EOM-CCSD). We also notice underestimation of the 4a′ → 4a″ transition oscillator strength by DSRG-MRPT2, as the ratio of the oscillator strengths is 0.51 for DSRG-MRPT2, while 0.95 for DSRG-MRPT3 and 0.92 for EOM-CCSD.

A detailed analysis of the DSRG-MRPT2 wave function shows the 4a′ → 4a″ excited state being primarily composed of the $(4a′)−1(4a′′)1$ configuration (84.1%) plus a small contribution from the $(4a′)−1(5a′′)1$ configuration (4.3%). However, DSRG-MRPT3 shows a more pronounced multi-determinantal wave function (58.2% and 24.7%), much closer to the composition of the GASSCF reference (38.2% and 36.1%). This is one of several cases where dynamical correlation from second- and third-order perturbation theories give opposite contributions to the determinantal composition of a state, with pyridine and pyridazine being two other examples.

For glyoxylic acid, we have traced this problem to the largest two-body amplitude (0.126), observing that when s is reduced to 0.5 $Eh−2$⁠, its value is reduced to 0.087 and the predicted energy splitting (1.91 eV) and an intensity ratio (0.74) between the two 1s → π* transitions agree well with DSRG-MRPT3 and EOM-CCSD. The DSRG-MRPT2 wave function is also more multi-determinantal with the $(4a′)−1(4a′′)1$ and $(4a′)−1(5a′′)1$ configurations contributing 73.4% and 13.6%, respectively. On the other hand, DSRG-MRPT3 is more robust to changes in s, predicting an energy splitting (2.14 eV) and an intensity ratio (0.98) similar to the corresponding values computed with $s=1Eh−2$⁠. The strong dependence on flow parameters indicates that, despite regularization of the small denominators, the DSRG-MRPT2 can show a more pronounced s dependence due to intruders. Interestingly, the split core approximation can also help mitigate the effect of intruder states. In the split-core calculations (denoted as SC in Table VI), the largest amplitude (⁠$s=1Eh−2$⁠) for the 4a′ → 4a″ transition reduces to 0.097 at the DSRG-MRPT2 level of theory. This treatment leads to a 1.72 eV energy splitting and an oscillator strength ratio of 0.60 between the two transitions, closer to the DSRG-MRPT3 (⁠$s=0.5Eh−2$⁠) values. A similar trend has also been found for 1,2-benzoquinone, for which the energy splittings and oscillator strength ratio evaluated with s = 0.5 and 1 $Eh−2$ are given in Table S3 of the supplementary material.

This study benchmarks a state-averaged version of the generalized active space driven similarity renormalization group (GAS-DSRG) approach for computing core-excited electronic states using the XABOOM set.⁸⁰ The core-excited states are first modeled by a state-averaged generalized active space self-consistent field reference, followed by a treatment of dynamical electron correlation via the state-averaged MR-DSRG approach. The state-averaged formalism is essential for modeling the high density of core-excited states in large organic systems and helps avoid root-flipping issues that arise in state-specific methods. Several approximations to the full valence active space were considered to truncate the determinant space of the largest molecules in the XABOOM set. The effect of the two approximations that split the orbital space and restrict the valence space to only the π orbitals are tested on the core-excited states of N₂O and O₃. Both approximations affect the energy splittings by less than 0.3 eV in N₂O. Still, these approximations can lead to large energy errors in the case of O₃ since the partitioning of the determinant space is inconsistent with the strong coupling between core-excited states observed in the full valence space. We also investigated the choice of the flow parameter s, finding that $s=1Eh−2$ is an appropriate value for GAS-DSRG calculations on organic systems.

We report the 1s → π* vertical transition energies and oscillator strengths of 40 molecules from second- and third-order DSRG treatments (DSRG-MRPT2/3) and find that these values are in good agreement with EOM-CCSD results.⁸⁰ Compared to existing experimental data, DSRG-MRPT3 transition energies agree better with absolute transition energies. At the same time, DSRG-MRPT2 and EOM-CCSD under and overestimate the transitions energies, showing a mean signed deviation of −0.74 and +0.69 eV, respectively. We also carefully analyzed ozone and glyoxylic acid, two molecules for which theoretical predictions show significant deviations. In addition to its higher accuracy, the third-order GAS-DSRG approach shows several other advantages compared to DSRG-MRPT2 in computations of core-excited states. DSRG-MRPT3 is more likely to avoid any inconsistency in the order of eigenstates before and after the relaxation step, and it is less sensitive to intruders in systems like glyoxylic acid.

The present work also highlights some remaining challenges in describing core-excited methods with multireference methods. First, limits on the size of the active space that can be treated with the GAS approach may affect the accuracy of computations on molecules larger than those in the XABOOM set. However, this restriction may also be relevant to special cases, like ozone, where a GAS larger than the full-valence space applied in this work may be necessary to reconcile the remaining differences between the experimental and theoretical spectra. Other challenges that should be targeted in future research include simplifying the selection of GAS orbitals, improving algorithms for the optimization of GASSCF reference states, and extending the formalism to improve the estimation of oscillator strengths. Finally, we note that the approach pursued in this work obtains core-excitation energies from separate computations on the ground and core-excited states. As such, one of its strengths is that it can accurately describe core-excited states, while one of its weaknesses is that it works best to describe a small number of states.⁷⁸ In contrast, methods that start from a correlated ground state and directly target the energy differences—including time-dependent DFT, EOM methods, and Green’s function approaches—can directly obtain a large number of states in one computation and are formally more appealing. Therefore, a desirable extension of the present work would be to formulate a direct approach to target core-excited states. Examples of multireference formalisms that explore this direction can be found in Refs. 70, 75–77, and 135.

Despite these challenges, the GAS-DSRG method is well-suited for the simulation of x-ray spectra of transient species, including applications to XTAS of molecules during chemical reactions. Further theoretical extensions of the approach to account for valence-excited states and their coupling to core-excited states would also allow simulating pump-probe UV/x-ray experiments to track electron and nuclear dynamics in photochemical processes. Other interesting future directions include the application of the GAS-DSRG method to the simulation of x-ray photoelectron spectroscopy and resonant-inelastic x-ray scattering.

See the supplementary material for (1) active spaces and number of states used for the GAS-DSRG approach of each molecule; (2) energy splittings and oscillator strength ratios for molecules with two or more transitions from the same element; (3) energies and oscillator strengths for 1,2-benzoquinone computed with different flow parameter values; and (4) GAS-DSRG transition energies, available experimental transition energies, and oscillator strengths of all molecules (xaboom_dsrg_data.xlsx).

The authors would like to thank Alexander Paul and Henrik Koch for providing the EOM-CC3 ozone results reported in this paper. This research was supported by the U.S. National Science Foundation under Award No. CHEM-1900532.

The authors have no conflicts to disclose.

Meng Huang: Conceptualization (equal); Data curation (equal); Investigation (equal); Methodology (equal); Validation (equal); Visualization (equal); Writing – original draft (equal); Writing – review & editing (equal). Francesco A. Evangelista: Conceptualization (equal); Funding acquisition (equal); Methodology (equal); Supervision (equal); Visualization (equal); Writing – review & editing (equal).

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