Non-orthogonal configuration interaction with single substitutions for the calculation of core-excited states

Due to the increasing accessibility of synchrotron radiation, x-ray absorption spectroscopy (XAS) methods have become far more powerful and are able to study an increasing number of systems. However, the theoretical calculation of core-excited states, which are often needed to identify XAS peaks, still has space to improve. An ideal core-excited state calculation is one which, while still efficient, considers the interactions between the different possible core excitations and correctly treats electron localization on atoms of the same type. If possible, it should also be systematically improvable.

There are many different approaches to calculating core-excited states. Early methods used multiple scattering X_α,¹ static exchange (STEX),² and transition potential calculations,³ but several methods have evolved since then. One such method is configuration interaction singles (CIS),⁴ which is based on a Hartree-Fock (HF) reference for an excited state that considers all single-electron excitations. However, CIS calculations do not allow for orbital relaxation,⁵ or include dynamic correlation,⁶ meaning that their results are not particularly accurate. Dynamic correlation via CIS(D)^7,8 offers only limited improvement.⁹ 

Another method is time-dependent density functional theory (TDDFT), which calculates poles of the DFT response function from a CIS-like equation. While it allows for only very limited orbital relaxation, it does include differential electron correlation effects. Another problem is that a TDDFT potential does not have the correct long-range Coulomb tail due to the self-interaction problem, meaning that excitation energies corresponding to states that use this tail are not accurate,¹⁰ so this method has been most widely used for valence-excited states. These methods can scale from N² to N⁴ depending on how they are implemented.¹¹ Significant effort, nonetheless, has been devoted to modifying standard density functionals to improve the performance of TDDFT for predicting core excitation energies.^12–16 In particular, in contrast to the success of long-range corrected functionals¹⁷ for valence Rydberg and charge-transfer excited states, short-range corrected functionals¹⁸ improve core-excited states.

Another set of methods is delta-self-consistent field (Δ-SCF) and Δ-DFT, which scale N² − N³ depending on the implementation. In these methods, the excited state energy is obtained by subtracting a ground-state from an excited-state HF or Kohn-Sham energy calculation.¹⁸ Since the excited state is computed separately and explicitly, this method allows for orbital relaxation. Δ-SCF does not include electron correlation effects, while Δ-DFT does through its functional.¹⁹ Another crucial component of this method is the functional used, as different functionals produce different results for different types of calculations. Some shortfalls of these methods are that Δ-DFT suffers from the same self-interaction problem as TDDFT and that both Δ methods are not spin pure, which is important when simulating XAS spectra. Additionally, they do not consider static correlation in the case of molecules with several of the same atom.

A different method is algebraic diagrammatic construction (ADC),²⁰ which scales as N⁵ for ADC(2), the most common variant used.^21,22 This approximation to the polarization propagator comes from Green’s function theory and it has been derived with the diagrammatic perturbation expansion of the polarization propagator using the Møller-Plesset (MP) partition of the Hamiltonian. For core-excited states, there exists the core-valence separation (CVS) addition to the ADC scheme,^23–27 which uses the decoupling of the space of core-excitations from the valence in the ADC matrix. Orbital relaxation effects are included indirectly in ADC(2), the second-order Møller-Plesset perturbation (MP2)-corrected version. Electron correlation is also included, resulting in a relatively high degree of accuracy.²⁸ 

A final set of methods to consider are those that come from equation-of-motion coupled-cluster singles and doubles (EOM-CCSD) methods,²⁹ specifically those for electronically excited states (EOM-EE) and electron affinities (EOM-EA), which scale as N⁶. In these types of calculations, target states are excitations from a reference state, and they are constructed using an excitation operator. For EOM-EE, the reference state is the closed-shell CC ground state, and the excitation operator conserves the number of α and β electrons. While EOM-EA and EOM-EE have the same reference state (typically the CCSD ground state), they differ in that EOM-EA uses electron attaching excitation operators, while EOM-EE uses electron exciting operators. Both methods include correlation effects,³⁰ but while EOM-EE does not fully include orbital relaxation,³¹ EOM-EA does, as it is used with a core-ionized reference.³² 

There have been several methods developed recently that have also made important strides toward the computation of high-energy states. One such example is the DFT/ROCIS work of Roemelt et al.³³ that was specifically designed for L-edge XAS and combines DFT with CI through the use of Kohn-Sham orbitals. Other work has been done by Coriani et al.³⁴ on damped coupled-cluster linear response functions to extend this method from UV to XAS spectra.

It is also important to discuss non-orthogonal configuration interaction (NOCI) methods for excited states.^35–38 NOCI uses different Hartree-Fock solutions as the basis for the CI calculation, which renders the basis non-orthogonal. This allows the basis to more effectively span the relevant part of the configuration space than an orthogonal calculation.³⁹ Since the creation of this method, further work has been done to increase the accuracy of NOCI calculations, for example, the spin-flip NOCI (SF-NOCI)⁴⁰ and NOCI-MP2⁴¹ methods. In SF-NOCI, the wavefunction is defined as a linear combination of independently relaxed Slater determinants obtained from all possible spin-flipping excitations, and it improves on several less-desirable elements in the original NOCI calculation process.⁴⁰ NOCI-MP2 was designed to add dynamic correlation to NOCI. The formalism has recently been improved to permit size-consistency, and an efficient implementation has been developed.⁴¹ 

In this paper, we present a novel method for the calculation of single and multi-reference core-excited states: non-orthogonal configuration interaction singles or NOCIS. This method considers the interactions not only between singly-excited states of different core orbitals but also between localized excited states for a molecule with two or more of the same atom. While NOCI is a general framework, where the appropriate determinants are chosen on a case-by-case basis, NOCIS is designed as a black-box method that includes a precisely defined set of determinants. NOCIS allows for orbital relaxation by using as its basis the electron-attached states of independently optimized, core-ionized references, providing, as we shall demonstrate, an impressive improvement over the standard orthogonal CIS method. Thus, it can efficiently produce all core-excitation energy values for a given molecule with a high degree of accuracy. It is important to note that this method is identical to the static exchange (STEX) approximation method⁴² for a single reference, but its invocation of non-orthogonality for multi-reference states renders it unique from other average-of-configuration methods such as 4C STEX,⁴³ which still use orthogonal orbitals.

The NOCIS wavefunction is constructed from R different core-ionized determinants corresponding to core electron removal from each identical atom of interest (e.g., R = 3 for C in C₃H₄O) and restricted open-shell Hartree-Fock (ROHF) optimization of these cation determinants. Electrons are then reattached to the V virtual orbitals (different for each cation determinant), and the resulting core-excited configuration is projected against the closed-shell ground state determinant to ensure that the core-excited states are strongly orthogonal to the ground state. The amplitudes $tia$ are unknowns to be variationally determined,

Here ∣Φ_i⟩ is the ROHF determinant with its core hole on atom i filled. $P^$ acts by projecting against the neutral ground state reference. Since Eq. (1) is a linear trial function, application of the variational principle gives a generalized eigenvalue problem for the core-excited state energies E_k and amplitudes t_k,

We can at this stage note that in addition to being variational, the NOCIS core-excited state energies are also size-consistent, just as for the standard orthogonal CIS energy. For a supersystem composed of two non-interacting molecules, M₁ and M₂, the NOCIS core excitation energies on M₁ are unaffected by the presence of M₂. This property is essential to correctly predict substituent effects on core excitation energies.

The algorithm for a NOCIS calculation is as follows:

1. A self-consistent field (SCF) calculation is run on the ground state of the given molecule.

2. If more than two core orbitals are to be considered, a Boys Localization⁴⁴ is performed on those orbitals (we restrict ourselves to K shells in this discussion though generalization is possible).

3. The molecular orbital (MO) coefficient matrices are rearranged to reflect the removal of an electron from the core orbital of interest.

4. A restricted open-shell Hartree-Fock (ROHF) optimization is performed on this cation state without freezing the core orbital.

5. Steps 3 and 4 are repeated as necessary for each relevant core orbital.

6. The diagonal blocks of the H, S, and S² matrices are constructed for each i (determinants in a block are orthonormal).

7. The off-diagonal elements of the matrices are constructed using NOCI.

8. The generalized eigenvalue problem is solved, the SCF ground energy is subtracted from the resulting energies to obtain the excitation energy of each state, and the ⟨S²⟩ value for each state can be obtained.

The number of core orbitals of a particular type defines the number of references in the calculation (e.g., a 1s calculation for C₂H₂ would have 2 references). Figure 1 shows a visual example of a completed NOCIS matrix. As NOCIS calculations have only been implemented so far for closed-shell molecules, those are what are discussed here. For NOCIS, it is necessary to consider both alpha and beta ionization, but, for closed shell reference, it is only necessary to perform the optimization for either alpha or beta.

In the aforementioned peach blocks of Fig. 1, NOCI and CIS calculations produce the same results. These diagonal blocks represent particle states, which are formed by adding an electron to a determinant, and thus they are simply the $vv$ block of the Fock matrix for ionized states. In the basis of these particle states, the matrix elements of H are then the matrix elements of the Fock operator, with the addition of the excited ion energy on the diagonal elements.

The S matrix is simply the identity because these matrix elements are between configurations constructed from the same set of orthonormal MOs. This orthonormality also makes the S² matrix simple to compute.

The main substance of NOCI can be found in Refs. 39 and 45, but the theory is summarized here. The intention is to build and diagonalize the Hamiltonian in a basis of non-orthogonal Slater determinants or basis states. This renders the equation a generalized eigenvalue problem. The mathematical method for computation of the overlap and Hamiltonian between two generally non-orthogonal basis states is described in detail by Thom and Head-Gordon.⁴⁵ 

To obtain the orbitals, we use the orbitals of Amos and Hall,^46,47 alternately called Löwdin-Paired orbitals.⁴⁸ After these biorthogonal states have been obtained, generalized Slater-Condon rules^45,46 are used to compute the H matrix elements and the generalized eigenvalue problem is solved.

One further element of note is that the paper of Sundstrom and Head-Gordon³⁹ contains an error in the derivation of the S² matrix elements. This was corrected by Levine and Head-Gordon.⁴⁹ 

The code for NOCIS was implemented in a development version of Q-Chem 5.0⁵⁰ utilizing the Armadillo linear algebra library.⁵¹ All molecular geometries were obtained from Curtiss et al.⁵² The ROHF optimizations were performed using the Newton Minres algorithm,⁵³ and the NOCI code was a modified version of the previous code developed by Sundstrom.³⁹ The maximum overlap method (MOM) is used to keep the core-hole from collapsing.⁵⁴ 

The approximate scaling of this code is N⁴. Each NOCI element requires a 2-electron integral build, which is a source of computational bottleneck. One solution, which has been implemented, is to calculate these integrals in batches instead of calculating them individually, but this also causes a significant increase in memory use. The code can be further accelerated by noting that many matrix elements share a common orbital basis, but this remains a further project. However, even in its current unoptimized state, NOCIS does not use any large intermediates that require disk storage, an advantage over ADC and EOM-CCSD methods, where the unknowns are as large as fourth-rank tensors.

The NOCIS calculations were tested using both ROHF and UHF methods for the cation optimization. Although the UHF method provided slightly more accurate energies, they were not spin-pure. Since S² values are vital to identifying XAS states, we decided to move forward with the spin-pure ROHF calculations.

We first assess the basis set convergence of NOCIS calculations on the 3 lowest allowed singlet core-excited states of CO (1s C and 1s O) and NH₃. The results are shown in Table I for the standard sequence of Dunning correlation consistent basis sets⁵⁵ and the sequences augmented with diffuse and/or tight functions. The use of diffuse functions, as expected, becomes more important for states higher than the lowest, which contain some Rydberg character in the particle levels that become partly occupied in the transition. Use of the CV augmentation, on the other hand, helps accelerate basis set convergence of the core hole. The largest improvement, however, is when both augmented and core-polarized basis sets are used. This is reasonable, as including core correlated functions as well as diffuse functions should result in a more accurate result for a core excitation. Looking more in-depth at the aug-cc-pcVXZ basis sets, the balance of error reduction and efficiency suggests that a good basis set for both efficiency and accuracy is aug-cc-pCVTZ. These conclusions are generally consistent with the recent basis set assessment of Fouda and Besley.⁵⁶ 

Having selected the aug-cc-pCVTZ basis, we next consider the comparison of the NOCIS method to three other methods, as shown in Table II, which shows the results of CIS, Δ-SCF, Δ-DFT, TDDFT, and NOCIS calculations expressed as deviations from experimental values. As expected, the CIS method performs very poorly due to its lack of orbital relaxation and correlation, with a RMS error of 19.64 eV. The Δ-SCF method, which includes orbital relaxation but not dynamic correlation, is a large improvement on CIS methods, with a RMS error of 6.55 eV. Notably, Δ-SCF calculations perform worse in cases such as N₂ where there is a problem of static correlation as a result of core orbitals with the same energy. Static correlation will be discussed in more depth later.

The Δ-DFT methods, which include both orbital relaxation and dynamic correlation, perform even better, as expected, with RMS errors of 2.62 eV and 1.97 eV, respectively, for B3LYP and SRC1-R1. The TDDFT methods perform extremely well, with an RMS error of 1.77. This is likely due to the fact that the SRC1-R1 functional⁷² was specifically designed for core excitations in TDDFT. This result shows the importance of correlation in obtaining an accurate core-excitation energy. However, TDDFT is extremely functional dependent, and performing it with other functionals produces results that vary considerably in accuracy.

The NOCIS calculations, which include orbital relaxation and static correlation but not dynamic correlation, perform remarkably well, even when compared to the DFT methods mentioned above, outperforming both Δ and TDDFT in the root-mean-squared error (RMSE). One particular area of strength for this calculation is its ability to account for static correlation, specifically excitations from core orbitals of the same type, as in C₂H₂. This is accomplished through Boys localization⁴⁴ on the relevant core orbitals and the non-orthogonal interactions between them. The addition of this element significantly improved the calculation results, often by several eV. In the cc-pCVDZ basis set, the errors on C₂H₂, N₂, CO₂ O, and F₂ were reduced by 6.53, 8.03, 10.24, and 7.61 eV, respectively, relative to non-localized NOCIS calculations. The NOCIS results together with Δ-SCF and standard CIS give interesting insight into the relative roles of orbital relaxation and configuration interaction in core excitation energies. It is clear that both are essential in cases where there are multiple equivalent nuclei, although the orbital relaxation effect is dominant in cases where only a single nucleus is core-excited (e.g., H₂CO). Cases where Δ-SCF performs better than NOCIS may well be due to fortuitous error cancellation. A further note is that none of these methods include relativistic effects, and thus the larger the atom, the more poorly the method performs, as shown in the F₂ and HF calculations.

Table III features NOCIS calculations compared with ADC calculations from Ref. 28. While ADC performed better with the CO C 1s, NOCIS performed better for the CO O 1s and the NH3 N 1s, which is consistent with the fact that ADC has been shown to overestimate energies for heavier atoms.²⁸ This shows that NOCIS is at least competitive with ADC and, with its smaller scaling, more efficient. Additionally, in some smaller basis sets, ADC does not provide the second and third excited states for these molecules, while NOCIS provides reasonable estimates for these states even in smaller basis sets. The fact that NOCIS, a method without dynamic correlation, is competitive with ADC-CVS(2), which includes it, shows the value of treating orbital relaxation via NOCI. Formulating a dynamic correlation to NOCIS looks to be potentially very worthwhile.

We also verified these k-edges by computing oscillator strengths in a method identical to those of the STEX method.⁴² It is worth noting that this has been done in an early version of STEX in the aug-cc-pVTZ basis set.⁷³ For NOCIS, in all cases, the lowest singlet excitation energy produced the largest oscillator strength, confirming our results. The results are shown in Table IV. The NOCIS results are competitive with or better than the ADC results and vastly out-perform the CIS results. We believe that this once again illustrates the importance of orbital relaxation.

This paper describes and benchmarks a new NOCI method specifically designed for core-excited states. Unlike the standard NOCI, this non-orthogonal CIS method is well-defined, using a precise number of configurations that depend only on the molecular composition and the basis set. This method compares favorably to all current methods of calculating core excitations in both efficiency and accuracy. However, the implementation is still preliminary, and some work is needed before we can claim a truly efficient method for the robust identification of XAS peaks. In this vein, a further effort for this project is to construct the output in such a way that it can be used to identify spectroscopic peaks. This will be done by implementing oscillator strength calculations and explicitly connecting orbital transitions to their energies.

There is also further theoretical and computational work to be done on NOCIS. As previously mentioned, it does not contain dynamic correlation, so a possible next step would be implementing something like NOCI-MP2⁴¹ in the construction of the matrix elements. A different direction is to include relativistic effects in NOCI calculations so that it can better calculate core-excitations for larger atoms such as silicon, which is beyond the capability of non-relativistic core-excited calculations. Computationally, the code can be sped up by implementing batching for the 2-electron integrals, as described in Sec. II, which reduces the prefactor and thus the total computation time.

Overall, NOCIS is a very promising method for calculating core excitations. Even without dynamic correlation, it still has a high level of accuracy, and it has been shown to deal well with static correlation. With multi-state return (in contrast to Δ methods) and competitive scaling, NOCIS provides a convenient path to obtain accurate calculations that can be used to simulate an XAS spectrum.

K.J.O. thanks the members of the Head-Gordon group for fruitful discussions, particularly Daniel Levine. This work was supported by the Director, Office of Science, Office of Basic Energy Sciences of the US Department of Energy under Contract No. DE-AC02-05CH11231.