Benchmarking the completely renormalised equation-of-motion coupled-cluster approaches for vertical excitation energies

The vertical excitation energies for a comprehensive test set of about 150 singlet excited states of 28 medium-sized organic molecules computed using two variants of the completely renormalised (CR) equation-of-motion (EOM) coupled-cluster (CC) method with singles, doubles, and non-iterative triples, abbreviated as δ-CR-EOMCCSD(T), and the analogous two variants of the newer, left-eigenstate δ-CR-EOMCC(2,3) approach are benchmarked against the previously published CASPT2, CC3, and EOMCCSDT-3 results, as well as the suggested theoretical best estimate (TBE) values. The δ-CR-EOMCC approaches are also used to identify and characterise about 50 additional excited states, including several states having substantial two-electron excitation components, which have not been found in the previous work and which can be used in future benchmark studies. It is demonstrated that the non-iterative triples corrections to the EOMCCSD excitation energies defining the relatively inexpensive, single-reference, black-box δ-CR-EOMCC approaches provide significant improvements in the EOMCCSD data, while closely matching the results of the iterative and considerably more expensive CC3 and EOMCCSDT-3 calculations and their CASPT2 and TBE counterparts. It is also shown that the δ-CR-EOMCC methods, especially δ-CR-EOMCC(2,3), are capable of bringing the results of the CC3 and EOMCCSDT-3 calculations to a closer agreement with the CASPT2 and TBE data, demonstrating the utility of the cost-effective δ-CR-EOMCC methods in applications involving molecular electronic spectra. We show that there may exist a relationship between the magnitude of the triples corrections defining δ-CR-EOMCC approaches and the reduced excitation level diagnostic resulting from EOMCCSD.


Introduction
Modern developments in quantum chemistry (QC) have provided versatile computational tools, whose application in areas of physical and biological sciences offers quantitative results as well as useful insights, helping to predict, verify, and understand complicated experimental observations and measurements. Given the plethora of QC methods, accompanied by the emergence and development of available program suites, it is essential to critically evaluate the performance and accuracy of existing QC approximations as well as the formulation of new and improved ideas. A few examples of the databases that exist for benchmarking a broad range of ground-state properties of ab initio and density functional theory (DFT) methodologies are the Gaussian Gn test sets [1][2][3][4][5][6] for atomisation energies, ionisation potentials, electron affinities, proton affinities, and enthalpies of formation, the S22 [7], S26 [8], and S66 [9] training sets for non-covalent interactions, the ATcT [10][11][12] benchmark set for thermochemical data, the DBH24 [13] test set for barrier heights, the W4 [14] training set for total atomisation energies, the 3dBE70 [15] benchmark set EOMCCSD approach [59][60][61][62] for the geometries and selected spectroscopic properties of a variety of doublet radicals, comparing its performance with experiment and other CC approaches.
It is often stated that the basic EOMCCSD [48][49][50] and linear-response CCSD [41,42] approximations, which construct the excited-state information on top of the CCSD [77][78][79][80] ground state and which are characterised by the relatively inexpensive iterative CPU steps that scale as n 2 o n 4 u (n o and n u are the numbers of occupied and unoccupied orbitals, respectively, used in the correlated calculations) or N 6 with the system size N , provide an accurate description of excited states dominated by one-electron transitions. However, EOMCCSD is not robust enough to obtain a quantitative description of many of such states, especially when larger polyatomic species are examined (cf., e.g., [64,[81][82][83]; for a thorough evaluation of a number of EOMCC methods, including EOMCCSD, illustrating the same, see [17,19,20,22,24,25]). It also fails to describe states with significant double excitation contributions. One can address these shortcomings by including triple excitations, as is done in the EOMCC approach with singles, doubles, and triples (EOMCCSDT) [84][85][86]. While the full treatment of triple excitations substantially improves the description of excited electronic states (see, e.g., [70,71,[84][85][86][87][88][89][90][91]), it is also accompanied by a steep increase in the iterative CPU time scaling characterising the EOMCCSDT approximation (n 3 o n 5 u or N 8 ), limiting its applicability to systems with up to a dozen or so correlated electrons and smaller basis sets. Thus, if one is to make use of the EOMCC methodologies in accurate calculations of molecular electronic spectra, EOMCC schemes which can account for the effects of triples in an approximate, cost-effective, and yet reliable manner need to be employed.
There are several ways of incorporating triple excitations in the EOMCC and linear-response CC formalisms without running into the prohibitive computational costs of full EOMCCSDT. For example, one can select the dom-inant triply excited components of the cluster operator T that defines the underlying ground-state CC wave function and three-body components of the linear excitation operator R μ in the EOMCC wave function ansatz through the use of active orbitals, as in the EOMCCSDt method [84,85,92] (see [72,91] for reviews; cf., also [93][94][95][96][97][98][99] for the related ground-state approaches). While this allows for full-EOMCCSDT-quality results at the cost of EOMCCSD times a small prefactor, the approach is no longer a black box as one has to select the active orbitals.

Theory
In the δ-CR-EOMCCSD(T) and δ-CR-EOMCC (2,3) methods developed in [63] and [64], we correct the vertical excitation energies obtained with EOMCCSD for the leading triples effects extracted from the MMCC considerations. Thus, if ω μ = E μ − E 0 represents the vertical excitation energy from the ground state | 0 to the excited state | μ , the δ-CR-EOMCCSD(T) and δ-CR-EOMCC (2,3) values of ω μ can be given the following general form: where ω (CCSD) μ is the EOMCCSD excitation energy obtained by diagonalising the similarity-transformed Hamiltonian H (CCSD) = e −T 1 −T 2 H e T 1 +T 2 in the space of singly and doubly excited determinants, | a i and | ab ij , respectively, with T 1 and T 2 representing the one-and two-body clusters obtained with CCSD, and δ μ is the triples correction to ω and the only difference between δ-CR-EOMCCSD(T) and δ-CR-EOMCC (2,3) is in the formulas for abc μ,ij k and M ij k μ,abc . In the δ-CR-EOMCCSD(T) approximation, which is based on the more general CR-EOMCCSD(T) considerations described in [63], M ij k μ,abc are the generalised moments of the EOMCCSD equations corresponding to projections of these equations on the triply excited determinants | abc ij k : M ij k μ,abc = abc ij k |H (CCSD) (R μ,0 + R μ,1 + R μ,2 )| , (3) where R μ,0 , R μ,1 , and R μ,2 are the zero-, one-, and two-body components of the linear excitation operator R μ defining the EOMCC wave function ansatz | μ = R μ e T | resulting from the EOMCCSD calculations, in which T = T 1 + T 2 and R μ = R μ,0 + R μ,1 + R μ,2 (R μ,0 is defined as r μ,0 1, where r μ,0 provides the weight of the reference determinant in the EOMCCSD wave function and 1 is the unit operator). The corresponding amplitudes abc μ,ij k that multiply moments M ij k μ,abc to produce the δ μ correction are calculated using the following expression: where | (CCSD) μ = (R μ,0 + R μ,1 + R μ,2 )e T 1 +T 2 | is the EOMCCSD wave function of excited state μ and |˜ μ = {R μ,0 + (R μ,1 + R μ,0 T 1 ) Here,R μ, 3 represents the approximate form of the threebody component of R μ given bỹ R μ,3 = i < j < k a < b < cr ij k μ,abc a a a b a c a k a j a i , are the corresponding triple excitation amplitudes and a p (a p ) is the creation (annihilation) operator associated with spin orbital p. In the most complete δ-CR-EOMCCSD(T) treatment, defining variant ID of it and abbreviated, following [63], as δ-CR-EOMCCSD(T),ID, the D abc μ,ij k denominator entering Equation (7) is calculated as whereH (CCSD) n is the n-body component ofH (CCSD) . In the simplified IA version, abbreviated as δ-CR-EOMCCSD(T),IA, we replace the Epstein-Nesbet form of D abc μ,ij k (Equation (8)) by the Møller-Plesset-style expression: where ε's are the Hartree-Fock (in our case, RHF) spinorbital energies.
As mentioned in Introduction, one of the biggest advantages of the δ-CR-EOMCCSD(T) and δ-CR-EOMCC (2,3) approaches, in addition to black-box character, is their relatively low computational cost compared to EOMCCSDT. In analogy to the EOMCCSD(T) [51], EOMCCSD(T) [52], and similar methods and their linear-response CCSDR(3) counterpart [46,47], the most expensive steps of δ-CR-EOMCCSD(T) and δ-CR-EOMCC (2,3) calculations scale as n 2 o n 4 u in the iterative EOMCCSD part and n 3 o n 4 u in the determination of the non-iterative triples corrections δ μ . Clearly, this is a lot less than the iterative n 3 o n 5 u steps of EOMCCSDT. Furthermore, unlike in EOMCCSDT, one does not have to store triply excited amplitudes in δ-CR-EOMCCSD(T) and δ-CR-EOMCC (2,3) calculations, since the non-iterative corrections δ μ rely on one-and two-body components of the T, R μ , and L μ operators obtained with EOMCCSD. There are other ways of reducing prohibitive costs of EOMCCSDT calculations, including the EOMCCSDT-3 [52] and CC3 [44][45][46][47] approaches that provide reference information in this work, but we must keep in mind that EOMCCSDT-3 and CC3 are considerably more expensive than δ-CR-EOMCCSD(T) and δ-CR-EOMCC (2,3). In the case of EOMCCSDT-3 and its EOMCCSDT-1 [51] predecessor, one uses MBPT-type arguments, similar to those exploited in the ground-state CCSDT-n considerations, to eliminate terms in the EOM-CCSDT eigenvalue problem that cause the n 3 o n 5 u scaling, reducing the resulting CPU operation count to a n 3 o n 4 u level. However, the n 3 o n 4 u steps of EOMCCSDT-3 are iterative and one still has to diagonalise the similaritytransformed Hamiltonian in the space containing triply excited determinants to obtain the corresponding excitation energies, i.e., one has to store triply excited amplitudes in the EOMCCSDT-3 iterations, so that computer costs of EOMCCSDT-3 calculations are much larger than in the case of δ-CR-EOMCCSD(T) and δ-CR-EOMCC (2,3). Similar remarks apply to the CC3 model, which is a linearresponse analogue of EOMCCSDT-3 using iterative n 3 o n 4 u steps in determining excitation energies. It is, thus, interesting to examine how accurate our δ-CR-EOMCCSD(T) and δ-CR-EOMCC (2,3) approaches are compared to their more expensive EOMCCSDT-3 and CC3 counterparts, which are known to provide results of the full EOMCCSDT quality as long as the excited states of interest do not have substantial double excitation character.

Benchmark molecules and computational details
The database set of 28 organic molecules proposed in [17] is composed of 7 unsaturated aliphatic hydrocarbons, 11 aromatic hydrocarbons and heterocycles, 6 carbonyl compounds, and 4 nucleobases, all shown in Figure 1. The  main Tables I and II of Schreiber et al. [17] contain a total of 221 electronically excited states, namely 149 singlet and 72 triplet excitations. The Supporting Information to [17] provides data on 22 additional singlet states, but this study, particularly in its statistical error analysis in Section 4.3, focuses on the 149 singlet excitations listed in Table I of [17]. Our EOMCC calculations have produced 54 additional singlet excited states in the energy range covered by Table I of [17], including five states that can be found in the Supporting Information to [17] and 49 states that have not been considered in the earlier benchmark work [17][18][19][20][21][22][23][24][25]. We provide information about these 54 additional singlet excitations in our tables as well, but we do not include them in our statistical error analyses, since the EOMCCSDT-3 and CC3 reference data are not available for them.
In determining the corresponding vertical excitation energies, we considered the ground-state equilibrium geometries taken from Schreiber et al. [17], which were optimised at the Møller-Plesset second-order perturbation theory (MP2) level using the 6-31G * [148] basis. In addition, to examine the effect of geometry on the calculated vertical excitation energies, we also re-optimised the ground-state geometries of the same set of molecules at the higher CR-CC(2,3),D level [149,150], consistent with the δ-CR-EOMCC (2,3),D approximation, using the TZVP [151] basis, which was used in the vertical excitation energy calculations reported in [17][18][19][20][22][23][24][25] and which is used in the EOMCC computations discussed in the present study. These additional CR-CC(2,3),D geometry optimisations were carried out using the parallel coarse-grain finite-difference model available in the CIOpt program suite [152,153], which we interfaced with the CR-CC(2,3) routines [149,150,154] available in the GAMESS package [155,156]. All of the δ-CR-EOMCCSD(T),IA, δ-CR-EOMCCSD(T),ID, δ-CR-EOMCC (2,3),A, and δ-CR-EOMCC (2,3),D computations reported in this work and the underlying EOMCCSD calculations were performed using the CC/EOMCC routines developed in [63,76,149,157], available in GAMESS as well. In all of the correlated calculations reported in this work, core electrons were kept frozen and spherical components of d basis functions were employed throughout.

Geometries and their effect on the calculated vertical excitation energies
Comparing the MP2/6-31G * and CR-CC(2,3),D/TZVP geometries, which can be found in the supplemental data, we see that the corresponding bond lengths and bond angles, which differ by 0.005Å and 0.2 degrees, respectively, on average, are in excellent agreement with each other. The largest differences in the bond lengths and angles occur in the nucleobases, namely cytosine (0.019Å) and adenine (0.6 degrees). For the dihedral angles, the largest difference between the MP2/6-31G * and CR-CC(2,3),D/TZVP results, of 0.3 degrees, occurs for cyclopropene. Using the above two sets of geometries, we computed the 203 vertical excitation energies for the 28 molecules comprising the database of Schreiber et al. [17], which are presented in Table 1 for the MP2/6-31G * geometries and Table 2 for the CR-CC(2,3),D/TZVP geometries. Comparing the two sets of vertical excitation energies, we can see that they are in very good agreement for each of the EOMCC approaches employed in this work. For example, if we examine the correlation plot for the vertical excitation energies corresponding to the 149 singlet excited states listed in Table I of Schreiber et al. [17] and computed at the δ-CR-EOMCC (2,3),D/TZVP level using the MP2/6-31G * and CR-CC(2,3),D/TZVP geometries, shown in Figure 2, we observe that their correlation coefficient is 0.9994 and maximum energy difference (MaxE) is 0.18 eV. The corresponding mean unsigned error (MUE) Figure 2. Correlation plot of the δ-CR-EOMCC (2,3),D/TZVP excitation energies (in eV) computed using the MP2/6-31G * and CR-CC(2,3),D/TZVP geometries. and mean signed error (MSE) values are 0.05 and 0.04 eV, respectively. This similarity of the δ-CR-EOMCC (2,3),D excitation energies obtained at the MP2/6-31G * and CR-CC (2,3),D/TZVP geometries remains virtually the same if  we include the entire set of 203 excitations listed in Tables 1  and 2 in our calculations. Thus, we can safely rely on the MP2/6-31G * geometries in discussing the relative performance of various EOMCC methods in this study. For this reason, much of our assessment of the δ-CR-EOMCC approaches examined in this work is based on the results obtained with the MP2/6-31G * geometries. The same geometries were used in the previous method assessments using the database of Schreiber et al. [17], reported in [17][18][19][20][21][22][23][24][25], which helps us in making judgements on the performance of our δ-CR-EOMCC methods, particularly when compared with the previously published EOMCCSDT-3 [22,24], CC3 [17,19,21,25], and CASPT2 [17,21] data.

Remarks on the calculated vertical excitation spectra for molecules in the database of Schreiber et al. [17]
As already mentioned, the 203 singlet excitation energies of the 28 molecules constituting the benchmark set of Schreiber et al. [17] determined in the present study are collected in Tables 1 and 2. We compare the EOMCCSD, δ-CR-EOMCCSD(T),IX, and δ-CR-EOMCC (2,3),X (X = A, D) results, obtained in this work using the MP2/6-31G * (Table 1) and CR-CC(2,3),D/TZVP (Table 2) geometries, with the previously reported CASPT2 [17,21], CC3 [17,19,21,25], EOMCCSDT-3 [22,24], and theoretical best estimate, TBE-1 [17] and TBE-2 [21], values. In addition to dominant orbital excitations, the nature of each computed excited state (singly, doubly, or somewhere in between) is characterised using the reduced excitation level (REL) diagnostic [74] resulting from the EOMCCSD calculations, which is defined such that REL ≈ 1 implies a one-electron transition and REL close to 2 indicates a doubly excited state (see Equation (26) in [74]). Typically, states dominated by single excitations, for which REL is close to 1, are characterised by smaller triples corrections. These corrections become larger for states having a more substantial double excitation character, where REL is much bigger than 1. In other words, there is a certain degree of proportionality between REL and the magnitude of the triples corrections to EOMCCSD excitation energies. Although not a main topic of this study, we examine this proportionality in the appendix, showing its potential usefulness in predicting electronic excitation spectra of δ-CR-EOMCC quality only on the basis of EOMCCSD calculations.
The main part of the original benchmark set of Schreiber et al. [17], presented in Table I of this reference, contains 149 π → π * , n → π * , and σ → π * vertical singlet excitations, the majority of which have REL values close to 1, with only relatively few (45 or ∼30%) having non-negligible, Table 1. Vertical excitation energies (in eV) for singlet states of all molecules in the test set of Schreiber et al. [17] using the geometries optimised at the MP2/6-31G * CC3 results taken from [17], except for those for the nucleobases which were reported in [25]. h EOMCCSDT-3 values taken from [24].
j Additional excited states not included among the 149 singlet excitations listed in Table I of [17]. States in this category that have been characterised in the Supporting Information to [17] are designated as n 1 X, where X is the irreducible representation and n indicates the state number within symmetry X (n is underlined to distinguish from the 149 states included in Table I of [17], which are labelled as n 1 X). States that have not been found in [17] and other prior benchmark work [18][19][20][21][22][23][24][25] are designated as 1 X, i.e., without the state number in front of the term symbol. In ordering these additional states, we In [17], these states were reported as having doubly excited character. However, as determined in the present study, they are predominantly one-electron transitions (REL for 2 1 state to the EOMCCSD and CC3 roots at 10.54 and 10.45 eV, respectively, producing large discrepancies, of more than 1 eV, with the CASPT2 and TBE-1 values reported in [17], of 9.31 and 9.3 eV, and the CC3/aug-cc-pVQZ value of 9.29 eV obtained in [21]. We have found another EOMCCSD root of the 1 A 1 symmetry at 9.77 eV, which results in δ-CR-EOMCC values in the 9.37-9.57 eV range, in very good agreement with the above CASPT2, TBE-1, and CC3/aug-cc-pVQZ data. Thus, the new EOMCC root that we have found must be the 2 1 q CC3/aug-cc-pVQZ results taken from [21]. r The CC3/aug-cc-pVTZ result reported in [21]. Our δ-CR-EOMCC excitation energies for the 2 1 A 1 (π → π * ) state of acetone and the corresponding CC3/TZVP, EOMCCSDT-3/TZVP, and TBE values reported in [17,22,24] suggest that the CC3/aug-cc-pVTZ root at 8.90 eV found in [21] does not represent the 2 1 A 1 (π → π * ) excitation. Based on our δ-CR-EOMCC analysis, the CC3/aug-cc-pVTZ root at 8.90 eV corresponds, most likely, to the 1 A 1 (π → π * /n → σ * ) state that has not been considered in the previous benchmark studies [17][18][19][20][21][22][23][24][25]. For this reason, in comparisons of our δ-CR-EOMCC and EOMCCSD data for the 2 1 A 1 (π → π * ) state of acetone with a TBE-2 result we use the TBE-1 value taken from [17] instead. s Schreiber et al. [17] missed the EOMCCSD root of the 1 A symmetry at 7.52 eV, which corresponds to the CC3 root at 7.24 eV found in [25] and which matches the TBE values reported in [17,21] of about 7.4 eV. As a result, they incorrectly interpreted the higher-energy EOMCCSD and CC3 roots at 8.52 (in our calculations, 8.50) and 8.27 eV, respectively, as the 2 1 A (π → π * ) state. To avoid relabeling of all 1 A (π → π * ) states, we designate the state at 8.50 (EOMCCSD) and 8.27 (CC3) eV as the 1 A (π → π * ) excitation and assign the 2 1 A (π → π * ) state to the EOMCCSD root at 7.52 eV and CC3 root found in [25] at 7.24 eV, in perfect agreement with the δ-CR-EOMCC values that range from 7.09 to 7.31 eV when the MP2/6-31G * geometry is employed. t Taken from [25]. The EOMCCSD excitation energies for the 2 1 A (n → π * ) and 3 1 A (n → π * ) states of uracil are 6.50 and 7.69 eV, respectively, not the other way around, as reported in [17]. Table 2. Vertical excitation energies (in eV) for singlet states of all molecules in the test set of Schreiber et al. [17] using, in the case of the EOMCCSD and δ-CR-EOMCC results, the geometries optimised at the CR-CC(2,3),D/TZVP level. a Molecule    f CC3 results taken from [17], except for those for the nucleobases which were reported in [25]. g EOMCCSDT-3 values taken from [24]. h CASPT2 and TBE-1 excitation energies taken from [17], with updates from [21], and TBE-2 values taken from [21].
i Additional excited states not included among the 149 singlet excitations listed in Table I of [17]. States in this category that have been characterised in the Supporting Information to [17] are designated as n 1 X, where X is the irreducible representation and n indicates the state number within symmetry X (n is underlined to distinguish from the 149 states included in Table I of [17], which are labelled as n 1 X). States that have not been found in [17] and other prior benchmark work [18][19][20][21][22][23][24][25]  These two TBE-1 values had to be taken from [161], since Schreiber et al. [17] provided no information what the best estimates for them might be.
n This state is the only almost purely doubly excited state among the 149 singlet excited states included in Table I of [17]. The large discrepancies between our  Table I of Schreiber et al. [17], i.e., the 1 1 B 3g (n 2 → π * 2 ) state of s-tetrazine, is a truly multi-reference-type two-electron transition with REL ≈ 2.
Thus, most of the singlet excited states found in [17] have a predominantly one-electron excitation nature. The doubly excited 1 1 B 3g (n 2 → π * 2 ) state of s-tetrazine is excluded from our statistical error evaluations discussed in Section 4.3, since we have no access to the corresponding CC3 and EOMCCSDT-3 reference data, while the existing CASPT2 and TBE values seem unreliable (see footnote o to Table 1 and the discussion later). In other words, in our overall statistical error analyses in Section 4.3, we start with the 148 singlet excited states from Table I of [17], as described above, and then, in each comparison of our δ-CR-EOMCC data with another method, we use as many states from this set as computed with this other method.

Linear polyenes: ethene, E-butadiene, all-E-hexatriene, and all-E-octatetraene
Comparing the δ-CR-EOMCC vertical excitation energies of the seven states of ethene, E-butadiene, all-Ehexatriene, and all-E-octatetraene, shown in Table 1, with their CASPT2 and TBE-2 counterparts from [17] and [21], Schreiber et al. [17] state that the 2 1 A g (π → π * ) states of E-butadiene, all-E-hexatriene, and all-E-octatetraene have large contributions from two-electron transitions. As a result, they and other authors [17,19,20,[23][24][25] remove these states from some of their statistical error evaluations, particularly when CC vs. CASPT2 comparisons are made, but we will not do this in the present work. The REL values of about 1.2 characterising the 2 1 A g (π → π * ) states of E-butadiene, all-E-hexatriene, and all-E-octatetraene indicate that these three states are predominantly one-electron transitions, with relatively small contributions from double excitations, which higher and more robust EOMCC levels, such as our δ-CR-EOMCC(2,3) approaches, should be able to handle quite well.

Unsaturated cyclic hydrocarbons: cyclopropene,
cyclopentadiene, and norbornadiene For this group of molecules, the results of our δ-CR-EOMCC calculations for the nine states included in Table I of Schreiber et al. [17] are of a similar quality as the CASPT2 and TBE-2 data [17,21], with MaxE and MUE values ranging from 0.16 to 0.42 eV and 0.08 and 0.21 eV, respectively. Compared to the CC3 [17] and EOMCCSDT-3 [22,24] calculations, the δ-CR-EOMCC results for the nine states of unsaturated cyclic hydrocarbons included in Table I of [17] are characterised by MaxE values in the range of 0.16-0.32 eV and MUEs ranging from 0.09 to 0.28 eV.
Among the various excited electronic states in this group of molecules, the 1 1 B 2 (π → π * ) state of cyclopentadiene is of particular note, as it has been the focus of several theoretical studies where the location of this state has been disputed by various high-level single-and multireference QC approaches (see [158,159] and references therein). Our δ-CR-EOMCC vertical excitation energies for this state agree very well with the CASPT2 and TBE-1 or TBE-2 (both equal to EOMCCSDT in this case [17,159]) values of 5.51 and 5.55 eV, respectively, with the δ-CR-EOMCC(2,3),D approach yielding the experimental band maximum, which is, after correcting for vibronic interactions, 5.43 eV [159]. Both the CC3 and EOMCCSDT-3 results place this state ∼0.2-0.3 eV higher in energy compared to our δ-CR-EOMCC values, suggesting that the δ-CR-EOMCC excitation energies, which almost perfectly match the TBE/EOMCCSDT and experimental data for this state, may be more accurate.
Schreiber et al. [17] mention that the valence excited states of cyclopentadiene, 2 1 A 1 (π → π * ) and 3 1 A 1 (π → π * ), contain significant contributions from double excitations. Their argument is based on the ∼0.3-0.4 eV lowering when going from EOMCCSD to CC3. However, as can be seen in our tables and as shown in previous studies [64,[81][82][83], it is not unusual to observe errors of this magnitude in EOMCCSD calculations for singly excited states. Indeed, the REL values characterising the 2 1 A 1 (π → π * ) and 3 1 A 1 (π → π * ) states of cyclopentadiene, of about 1.153 and 1.055, respectively, when the MP2/6-31G * geometry is employed, show that these states are predominantly oneelectron transitions and that double excitations do not significantly contribute to their wave functions. Thus, in analogy to the 2 1 A g (π → π * ) excitations in linear polyenes, there is no reason to exclude them from the statistical error analyses of various methods, as has been done in some earlier studies [17,19,20,[23][24][25]. The 2 1 A 1 (π → π * ) and 3 1 A 1 (π → π * ) states of cyclopentadiene are included in our statistical error analyses in Section 4.3.

Aromatic hydrocarbons: benzene and
naphthalene In analogy to the unsaturated cyclic hydrocarbons, the δ-CR-EOMCC results for the 4 excited states of benzene and 10 excited states of naphthalene included in Table I of Schreiber et al. [17] are in generally good agreement with the CASPT2 and TBE-2 results reported in [17,21]. This is particularly true for the δ-CR-EOMCC (2,3),D calculations, which are characterised by the MUE and MaxE values of 0.17 and 0.39 eV, respectively, when compared to CASPT2, and 0.20 and 0.40 eV, respectively, when compared to TBE-2. The remaining three δ-CR-EOMCC approaches have very similar MUE values relative to the CASPT2 and TBE-2 data, of 0.17-0.29 and 0.20-0.27 eV, respectively, but the corresponding MaxE values are somewhat higher, especially in the δ-CR-EOMCCSD(T),IA case, where MaxE relative to CASPT2 is 0.87 eV and relative to TBE-2 0.86 eV. Comparing to the iterative triples CC3 and EOMCCSDT-3 approaches exploited in [17,22,24,25], we see that the non-iterative δ-CR-EOMCCSD(T) corrections produce MUEs of 0.18 and 0.16-0.18 eV, respectively. The δ-CR-EOMCC(2,3) approaches give similar deviations from CC3 and EOMCCSDT-3 (all in the range of ∼0.2-0.3 eV).
The 14 excited states of benzene and naphthalene considered in this work are predominantly one-electron transitions. In particular, the four states of naphthalene, namely 2 1 A g (π → π * ), 1 1 B 1g (π → π * ), 3 1 A g (π → π * ), and 3 1 B 3u (π → π * ), in contrast to their characterisation in [17] as having strong contributions from doubly excited configurations, are shown to be predominantly one-electron transitions, with REL values ranging from 1.109 to 1.177. Therefore, unlike in the earlier ab initio benchmark studies [17,19,20,[23][24][25], we do not exclude them from the statistical error analyses in Section 4.3, as there are no reasons to do so.

Heterocycles: furan, pyrrole, imidazole, pyridine,
pyrazine, pyrimidine, pyridazine, s-triazine, and s-tetrazine This group is the largest subgroup among the 28 molecules in the database of Schreiber et al. [17]. As shown in Tables 1 and 2, our δ-CR-EOMCC calculations have identified a total of 87 excited states for the molecules in it. Among them are 66 states listed in Table I of [17], the 2 1 A g (n 2 → π * 2 ) state of s-tetrazine that can be found in the Supporting Information to [17], and 20 other additional states found in this work, but not considered in the prior benchmark studies [17][18][19][20][21][22][23][24][25]. As already alluded to above and as further elaborated on below, one of the states of s-tetrazine, namely 1 1 B 3g (n 2 → π * 2 ), included in Table I of Schreiber et al. [17], is an almost pure two-electron transition, which is excluded from our statistical analyses due to the absence of sufficiently many reliable data for it. If we compare the results of our δ-CR-EOMCC calculations for the remaining 65 states of the heterocycles considered here and listed in Table I of [17], for which the CASPT2, CC3, and EOMCCSDT-3 reference data are available [17,22,24,25] and which are of predominantly single excitation nature, we can see a great deal of consistency among the various theoretical values. Indeed, the MUEs characterising the δ-CR-EOMCC calculations for these states relative to the corresponding CASPT2, CC3, and EOMCCSDT-3 excitation energies are 0.17-0.24, 0.14-0.27, and 0.15-0.36 eV, respectively. Unsurprisingly, the corresponding MaxE values are larger, particularly when one compares the δ-CR-EOMCCSD(T),IA and CASPT2 results, where MaxE is 1.16 eV, but the δ-CR-EOMCC(2,3) approaches, especially variant D, are very effective in reducing them to an acceptable level. The largest deviation between the δ-CR-EOMCC(2,3),D and CASPT2/CC3/EOMCCSDT-3 excitation energies of the 65 states of heterocycles included in Table I of Schreiber et al. [17], other than the 1 1 B 3g (n 2 → π * 2 ) doubly excited state of s-tetrazine excluded from our statistics, is 0.58 eV. Similar remarks apply to a comparison of our δ-CR-EOMCC excitation energies with the TBE-2 data reported in [21]. The MUE and MaxE values characterising the differences between the δ-CR-EOMCC and TBE-2 excitation energies of the 65 states of heterocycles included in Table I of [17], other than the twoelectron 1 1 B 3g (n 2 → π * 2 ) transition in s-tetrazine, are 0.13-0.20 and 0.40-0.50 eV, respectively. In other words, if we limit ourselves to the excited states of heterocycles dominated by one-electron transitions, the δ-CR-EOMCC approaches, especially the triples corrections defining the δ-CR-EOMCC (2,3) approximations, are capable of providing the CC3/EOMCCSDT-3-quality data, while matching the results of the multi-reference CASPT2 calculations and TBE-2 values. An interesting question arises if the above observations apply to the doubly excited 1 1 B 3g (n 2 → π * 2 ) state of s-tetrazine excluded from the above error analysis and disregarded in all of the previously published CC/EOMCC work [17,[19][20][21][22][23][24][25].
Among the 21 additional states of the heterocycles considered in this work, six π → σ * excitations in pyrrole, five states of pyridine representing n → π * , π → σ * , and n → π * /π → σ * excitations, five states of s-triazine of the n → σ * and n → π * types, and one n → π * excitation in s-tetrazine are predominantly one-electron transitions, as demonstrated by their REL ≈ 1 values. As a result, all triples corrections to EOMCCSD are quite small (∼0.1-0.4 eV) and all δ-CR-EOMCC approaches provide similar values. This is especially true when we compare the δ-CR-EOMCC (2,3),A and δ-CR-EOMCC (2,3),D data, which agree to within 0.1 eV. Such an agreement typically implies that the δ-CR-EOMCC(2,3) excitation en-ergies are reasonably well converged. Thus, in the absence of other high-level CC/EOMCC data for the 21 additional states of the heterocycles considered in this work, we can treat our δ-CR-EOMCC (2,3) values for the 17 singly excited states among them as reference data for future studies.
Among all the heterocycles listed in Tables 1 and 2, s-tetrazine is the only molecule that has doubly excited states in the prescribed energy range, which are characterised by REL close to 2. One of them is the 1 1 B 3g (n 2 → π * 2 ) excitation, which we have already discussed above. The remaining doubly excited states of s-tetrazine, which we have found, are the 2 1 A g (n 2 → π * 2 ) state, which was also found by Schreiber et al. [17] (cf. the Supporting Information to [17]), and the n 2 → π * 2 excitations of the 1 B 3g , 1 B 1g , and 1 B 1u symmetries that have not been considered before. As in the previously discussed 1 1 B 3g (n 2 → π * 2 ) state, the excitation energies for all of these doubly excited states lie quite high in the EOMCCSD spectrum. However, when triples are properly accounted for, using, for example, our δ-CR-EOMCC (2,3),D approach, we see a significant lowering of the EOMCCSD excitation energies. For example, for the 2 1 A g (n 2 → π * 2 ) state found by Schreiber et al. [17], the more robust triples correction of δ-CR-EOMCC (2,3),D lowers the EOMCCSD energy of 11.47 to 5.04 eV. This result is in reasonable agreement with the CASPT2 excitation energy of 4.55 eV reported in the Supporting Information to [17], to which our comments above on the underestimation of excitation energies by CASPT2 still apply. For the 1 B 3g (n 2 → π * 2 ) two-electron transition that has not been considered before, the EOMCCSD excitation energy of 15.10 eV is lowered by almost 8 eV when the δ-CR-EOMCC (2,3),D method is employed, resulting in 7.31 eV, and for the remaining two doubly excited states of s-tetrazine, namely 1 B 1g (n 2 → π * 2 ) and 1 B 1u (n 2 → π * 2 ), the δ-CR-EOMCC (2,3),D corrections to the EOMCCSD values of 14.31 and 15.04 eV, respectively, result in the final excitation energies of 7.96 and 8.07 eV, respectively. As discussed above, the δ-CR-EOMCC(2,3),D method has been demonstrated to provide an accurate description of the challenging electronically excited states dominated by two-electron transitions. Thus, we expect our δ-CR-EOMCC (2,3),D estimates for all five doubly excited states of s-tetrazine considered in this work to be reasonable, possibly serving as reference data in future benchmark studies.

Aldehydes and ketones: formaldehyde, acetone, and p-benzoquinone
While this subgroup is much smaller than the previous one, with only 3 molecules and 13 one-electron excitations reported in Table I of Schreiber et al. [17], it is of interest, as we have found several additional excited states in the prescribed energy range, including 3 states of p-benzoquinone characterised in the Supporting Information to [17] (the 1 1 B 2g (n → π * ), 2 1 B 2g (n → π * ), and 1 1 B 2u (π → π * ) states) and 11 other states that have not been considered before. Among the other states, two are dominated by n, π → π * 2 two-electron transitions, which we elaborate on later.
Comparing to the other high-level CC/EOMCC methods, we can see that our δ-CR-EOMCC results for this group of molecules are in good agreement with the CC3 [17] and EOMCCSDT-3 [22,24] data. The MUE values characterising the deviations of the δ-CR-EOMCC excitation energies from the CC3 and EOMCCSDT-3 results are in the range of 0.11-0.29 eV. The corresponding MaxE values range from 0.32 to 0.52 eV, when the δ-CR-EOMCC energies are compared to CC3, and from 0.26 to 0.47 eV, when we replace CC3 by EOMCCSDT-3. Similar remarks apply to comparisons of the δ-CR-EOMCC excitation energies characterising the 13 states of formaldehyde, acetone, and p-benzoquinone listed in Table I of Schreiber et al. [17] with their CASPT2 and TBE-2 counterparts reported in [17,21]. In this case, the MUE values characterising the differences between the δ-CR-EOMCC and CASPT2 data are 0.16-0.24 eV. For the differences between the δ-CR-EOMCC and TBE-2 excitation energies, we obtain 0.14-0.18 eV. Once again, the corresponding MaxE values are somewhat higher, but the δ-CR-EOMCC (2,3),D approach reduces them to an acceptable level, especially when compared with CASPT2, where we obtain 0.34 eV. In comparing the δ-CR-EOMCC and TBE-2 excitation energies for formaldehyde, acetone, and p-benzoquinone, we treat one of the 12 states for which such comparison can be made, namely the 2 1 A 1 (π → π * ) state of acetone, differently than the remaining states. As explained in footnote r to Table 1 (cf. also footnote q to Table 2), the TBE-2 value for the 2 1 A 1 (π → π * ) state of acetone reported in [21] is a result of an incorrect assignment of the CC3/aug-cc-pVTZ root at 8.90 eV. Based on our δ-CR-EOMCC results for acetone, we think that the CC3/aug-cc-pVTZ root at 8.90 eV found in [21] corresponds to the 1 A 1 (π → π * /n → σ * ) state of acetone, which has not been considered in the previous benchmark studies [17][18][19][20][21][22][23][24][25], not to the 2 1 A 1 (π → π * ) excitation. Indeed, we observe an excellent agreement between our δ-CR-EOMCC excitation energies for the 1 A 1 (π → π * /n → σ * ) state of acetone, which range from 8.77 to 8.99 eV when the MP2/6-31G * geometry is employed, and the CC3/aug-cc-pVTZ value of 8.90 eV reported in [21]. We now move to the 14 additional states of formaldehyde, acetone, and p-benzoquinone, which are not included in Table I of Schreiber et al. [17] and which we have found in our calculations. As already alluded to above, three of these additional states are the 1 1 B 2g (n → π * ), 2 1 B 2g (n → π * ), and 1 1 B 2u (π → π * ) excitations in p-benzoquinone described in the Supporting Information to [17], where the authors characterised them using CASPT2. As shown in Table 1, our best δ-CR-EOMCC (2,3),D results for these three predominantly singly excited states agree with the corresponding CASPT2 data reported in the Supporting Information to [17] to within ∼0.3-0.4 eV. We can view this as a good agreement given the fact that CASPT2 tends to underestimate excitation energies.
For the remaining 11 states that have not been considered before, we have to make the following two comments. The first one deals with the confusion regarding the 2 1 A 1 (π → π * ) and 1 A 1 (n → π * ) states of formaldehyde (see footnotes p and q in Table 1). Schreiber et al. [17] have interpreted the EOMCCSD and CC3 roots at 10.54 and 10.45 eV, respectively, as the 2 1 A 1 (π → π * ) state, producing large discrepancies, of more than 1 eV, with the CASPT2 and TBE values reported in [17], of 9.31 and 9.3 eV, respectively, and the CC3/aug-cc-pVQZ value of 9.29 eV obtained in [21]. It is hard to understand such a discrepancy given the fact that all excited states of formaldehyde considered in [17] are almost pure single excitations. After thorough examination of the problem, we have found another EOMCCSD root of the 1 A 1 symmetry at 9.77 eV, which results in the δ-CR-EOMCC values that range from 9.37 eV in the δ-CR-EOMCC (2,3),D case to 9.57 eV when the δ-CR-EOMCCSD(T),IA approach is employed, in very good agreement with the above CASPT2, TBE, and CC3 values, particularly when our best δ-CR-EOMCC (2,3),D result of 9.37 eV is considered. Thus, the new EOMCC root that we have found in this work must be the 2 1 A 1 (π → π * ) state, whereas the 1 A 1 state at 10.54 eV in the EOM-CCSD calculations and at 10.45 eV in the CC3 calculations, interpreted in [17] as the 2 1 A 1 (π → π * ) excitation, is the next state of the 1 A 1 symmetry, designated in our Tables 1 and 2 as 1 A 1 (n → π * ).
Our second comment regarding the additional states of aldehydes and ketones found in this work deals with the 1 B 3u (n, π → π * 2 ) and 1 B 2g (n, π → π * 2 ) doubly excited states of p-benzoquinone. In analogy to the previously discussed two-electron transitions in s-tetrazine, these two challenging multi-reference states cannot be properly charactersed by EOMCCSD. When our best δ-CR-EOMCC(2,3),D approach is used, we see the excitation energies of the 1 B 3u (n, π → π * 2 ) and 1 B 2g (n, π → π * 2 ) states of p-benzoquinone lower by about 5 eV compared to EOM-CCSD, giving 6.70 and 6.92 eV, respectively. The energy lowering provided by the remaining three δ-CR-EOMCC approaches, particularly in the δ-CR-EOMCCSD(T),IA case, is not as large, but this is consistent with our earlier experiences with the δ-CR-EOMCC methods. It would be desirable to perform the full EOMCCSDT or active space EOMCCSDt calculations for p-benzoquinone to verify if our best δ-CR-EOMCC (2,3),D estimates of the vertical excitation energies of the 1 B 3u (n, π → π * 2 ) and 1 B 2g (n, π → π * 2 ) doubly excited states of p-benzoquinone are accurate.

Amides: formamide, acetamide, and
propanamide Unlike in the previous group of aldehydes and ketones, all of the excited states of the three amides considered in this work, including 9 states listed in Table I of Schreiber et al. [17] and 14 additional states found in our calculations, but not considered in the prior benchmark studies [17][18][19][20][21][22][23][24][25], are dominated by one-electron transitions. As a result, there is a very good agreement between our δ-CR-EOMCC excitation energies and the previously published [17,21,22,24] CASPT2, CC3, EOMCCSDT-3, and TBE-2 values for the subsets of nine (CASPT2, CC3), eight (EOMCCSDT-3), and six (TBE-2) states of formamide, acetamide, and propanamide, for which the latter data are available. Indeed, the MUE values characterising the deviations of the δ-CR-EOMCC excitation energies for the nine states of the three amides included in Table I of [17] from the corresponding CASPT2 data, reported in [17] as well, range from 0.25 to 0.30 eV, with the largest deviation (MaxE) of 0.56 eV observed when the δ-CR-EOMCCSD(T),IA and CASPT2 results are compared with each other. As in the case of the other molecular groups examined in this work, the δ-CR-EOMCC (2,3),D approach improves the agreement with CASPT2, reducing the above MaxE value of 0.56 eV to 0.35 eV. Similar remarks apply to a comparison of the δ-CR-EOMCC results with the available [21] TBE-2 values. In this case, the MUE values characterising the differences between the δ-CR-EOMCC and TBE-2 data are 0.23-0.27 eV and, once again, the δ-CR-EOMCC (2,3),D approach produces the smallest MaxE value of 0.42 eV. The agreement between our δ-CR-EOMCC excitation energies for formamide, acetamide, and propanamide, and the previously published reference values improves even further when the δ-CR-EOMCC results are compared with the CC3 and EOMCCSDT-3 data reported in [17,22,24]. The MUE and MaxE values characterising the differences between our δ-CR-EOMCC results and their CC3 counterparts reported in [17] are 0.08-0.16 and 0.14-0.27 eV, respectively. They are 0.05-0.24 and 0.12-0.28 eV, respectively, when the δ-CR-EOMCC and EOMCCSDT-3 excitation energies are compared with each other.
As has been the case with the other subgroups of molecules examined in this work, we have found several additional excited states of formamide, acetamide, and propanamide that have not been considered in the previous benchmark studies [17][18][19][20][21][22][23][24][25], with acetamide having the richest variation of the excitation types. Of the 14 additional excited states we have found in the prescribed energy range, there are one π → π * excitation in formamide (which we will further comment on later), two n → π * , three n → σ * , three π → σ * , and one π → π * /n → σ * excitations in acetamide, and one n → π * and three π → σ * excitations in propanamide. As already mentioned, all of these additional excited states are predominantly one-electron transitions. Thus, as in the case of the previously discussed 9 states of formamide, acetamide, and propanamide considered in the prior benchmark studies [17][18][19][20][21][22][23][24][25], the triples corrections to the EOMCCSD excitation energies resulting from our δ-CR-EOMCC calculations for the 14 additional states of the three amides found in this work are rather small (∼0.1-0.3 eV). Our best δ-CR-EOMCC (2,3),A and δ-CR-EOMCC (2,3),D values, in analogy to the aldehydes and ketones, agree almost perfectly, differing by ∼0.05 eV. Thus, in the absence of other high-level CC/EOMCC results for the 14 additional excited states of formamide, acetamide, and propanamide found in the present study, our δ-CR-EOMCC(2,3) excitation energies can be used as accurate reference data in future benchmark calculations or applications involving these three molecules.
The 2 1 A (π → π * ) state of formamide considered in the previous benchmark work [17][18][19][20][21][22][23]25] and the extra 1 A (π → π * ) state of the same molecule found in the present study require an additional comment. As explained in footnote s to Table 1, Schreiber et al. [17] have incorrectly assigned the EOMCCSD and CC3 roots at 8.52 (8.50 in our calculations) and 8.27 eV, respectively, to the 2 1 A (π → π * ) transition, where the corresponding CASPT2, TBE-1, and TBE-2 values are 7.39, 7.39, and 7.35 eV, respectively. The 2 1 A (π → π * ) state of formamide is an almost pure one-electron transition (REL = 1.096 when the MP2/6-31G * geometry is employed) and so differences between the CASPT2 and CC3 results on the order of 1 eV are unlikely. Our EOMCCSD calculations show that there is another root of the 1 A symmetry and dominated by π → π * excitations, designated in Tables 1 and 2 as 1 A (π → π * ) and located at 7.52 eV, which has the δ-CR-EOMCC excitation energies ranging from 7.09 to 7.31 eV, in excellent agreement with the above CASPT2 and TBE values. We have, thus, reassigned the previously reported [17] EOM-CCSD and CC3 roots at 8.52 (in our calculations, 8.50) and 8.27 eV, respectively, to the higher energy 1 A (π → π * ) state, while interpreting the EOMCCSD root at 7.52 eV as the 2 1 A (π → π * ) state. Kánnár and Szalay [25] have found a CC3 root of formamide of 1 A symmetry and π → π * character at 7.24 eV, lending support to the above state reassignment.

Nucleobases: cytosine, uracil, thymine, and
adenine For this final group of molecules, there are a total of 31 excited states listed in Table I of [17]. Our δ-CR-EOMCC calculations have also found two additional states in the prescribed energy range that have not been considered before (two π → π * excitations in adenine). As shown in Tables 1 and 2, all of these states are dominated by one-electron transitions. The 31 states included in Table I of Schreiber et al. [17] have been characterised by CASPT2 and 19 of them have been assigned the TBE values of the corresponding excitation energies [17,21]. Due to prohibitive compu-tational costs of the high-level iterative triples CC3 and EOMCCSDT-3 calculations, none of the initial benchmark studies [17][18][19][20][21][22][23][24] using the database of [17] have provided the CC3 and EOMCCSDT-3 data for nucleobases. The CC3 reference data for 22 of the 31 excitations of nucleobases listed in Table I of Schreiber et al. [17] have finally become available in [25]. All of this means that in assessing performance of the δ-CR-EOMCC approaches in calculations for nucleobases, we are somewhat limited, although comparisons with the available CASPT2 [17], CC3 [25], and TBE [17,21] data allow us to make some useful observations. Compared with the CASPT2 excitation energies for all 31 states included in Table I [17], it is encouraging to observe the very good agreement between our best δ-CR-EOMCC (2,3) calculations and the considerably more demanding CC3 computations.

Statistical error analysis for the entire benchmark set
We now turn to the examination of the overall quality of the δ-CR-EOMCCSD(T),IA, δ-CR-EOMCCSD(T),ID, δ-CR-EOMCC (2,3),A, and δ-CR-EOMCC (2,3),D excitation energies obtained in this work, as compared to the previously published CASPT2 [17,21], TBE-2 [21], CC3 [17,25], and EOMCCSDT-3 [22,24] data and the underlying EOMCCSD results, using a larger set of 148 singlet excited states of 28 molecules considered in this work, taken from Table I of Schreiber et al. [17], as a starting point for the corresponding statistical error evaluations. As explained in the beginning of Section 4.2, the list of 148 singlet excited states used to initiate the statistical error analyses discussed in this section has been obtained by excluding the doubly excited 1 1 B 3g (n 2 → π * 2 ) state of s-tetrazine, for which the CC3 and EOMCCSDT-3 reference data are unavailable and the existing CASPT2 and TBE excitation energies reported in [17,21] unreliable, from the 149 singlet excitations collected in Table I of [17]. As a result, the vast majority of states entering the statistical error analysis presented in this section are dominated by one-electron transitions, with only 45 out of the above 148 states having 1.1 < REL < 1.2, three states having REL ≥1.2, and no state having REL greater than 1.225. As pointed out in Section 4.2, we have identified several additional states with significant contributions from double excitations, including, in addition to the above 1 1 B 3g (n 2 → π * 2 ) state of s-tetrazine, 6 other states with REL close to 2 that can be found among the 54 extra roots outside the set of 149 states listed in Table I of Schreiber et al. [17] obtained in our EOMCC calculations, but none of these additional states can be considered in the overall error evaluation, since we do not have access to the appropriate high-level reference data which would allow us to assess performance of the δ-CR-EOMCC approaches in this case. On the basis of our past experiences with the δ-CR-EOMCC calculations [63,70,71,[74][75][76]107,111,127,140,143] and the previously discussed comparison of the δ-CR-EOMCC (2,3),D and NEVPT2 [23] results for the doubly excited 1 1 B 3g (n 2 → π * 2 ) state of s-tetrazine (cf. footnotes o and n to Tables 1  and 2, respectively), we may expect the results of our overall best δ-CR-EOMCC (2,3),D calculations for the excited states with REL close to 2 identified in this work to be accurate to within ∼0.2 − 0.3 eV, but we would need to perform additional high-level EOMCCSDT/EOMCCSDt or MRCI calculations, which are beyond the scope of the present study, to verify such a statement.
The results of our statistical error analyses involving various ab initio methods considered in this work and the 148 singlet excited states of 28 molecules taken from Table  I of Schreiber et al. [17], as described above, or their appropriate subsets for which the relevant data are available, as further elaborated on below, are summarised in Tables 3-6.  Table 3 compares the MUE, MSE, and MaxE values characterising the differences between the excitation energies obtained in our EOMCCSD and δ-CR-EOMCC calculations, the CC3 calculations reported in [17,25], and the EOMCCSDT-3 calculations reported in [22,24] from the corresponding CASPT2 data given in [17], with updates provided in [21]. Table 4 does the same for comparisons of the EOMCCSD, δ-CR-EOMCC, CC3, and EOMCCSDT-3 data with their TBE-2 counterparts taken from [21]. Tables 5 and 6 compare the MUE, MSE, and MaxE values characterising the differences between the excitation energies resulting from our EOMCCSD and δ-CR-EOMCC calculations with the CC3 (Table 5) and EOMCCSDT-3 (Table 6) data reported in [17,22,24,25]. To make sure that our statistical analyses of errors are as systematic and as fair as possible, we focus on the δ-CR-EOMCCSD(T),IX and δ-CR-EOMCC (2,3),X (X = A, D) excitation energies, and their EOMCCSD counterparts summarised in Table 1, i.e., those obtained using the TZVP basis set and the MP2/6-31G * geometries, since the reference CASPT2, CC3, and EOMCCSDT-3 results reported in [17,21,24,25] and Table 3. Statistical error analyses of the various CC/EOMCC data, including EOMCCSD, δ-CR-EOMCCSD(T),IA (δ-CR(T),IA), δ-CR-EOMCCSD(T),ID (δ-CR(T),ID), δ-CR-EOMCC (2,3),A (δ-CR(2,3),A), δ-CR-EOMCC (2,3),D (δ-CR(2,3),D), CC3, and EOMCCSDT-3, from their CASPT2 counterparts.       collected in Table 1 as well have been obtained using the TZVP basis and the geometries resulting from the MP2/6-31G * optimisations. Clearly, our δ-CR-EOMCC excitation energies would change somewhat if we used basis sets larger than TZVP, but comparing the results obtained with various methods where different calculations use different basis sets would make our comparisons less systematic. The same applies to comparisons of vertical excitation energies using different ground-state geometries in different excited-state calculations, although comparing vertical excitation spectra at the ground-state geometries obtained with the same method as used in excited-state calculations for each of the ab initio approaches analysed in this study would be interesting. As explained in Section 4.1 (cf. Figure 2), the effect of ground-state geometries on the vertical excitation spectra resulting from the δ-CR-EOMCC calculations is small and unlikely to have an effect on our main conclusions regarding performance of the δ-CR-EOMCC methods relative to other approaches, but the topic of nuclear geometries used in the various computations for the 28 molecules comprising the benchmark set of Schreiber et al. [17] is worth further exploration. The effect of the basis set on the δ-CR-EOMCC calculations summarised in Tables 1 and 2 is interesting too, and we plan to examine it in the future. We should also point out that while in the original benchmark set presented in Table I of [17] there are a total of 149 excited states, including the aforementioned 148 excitations that provide the basis for the statistical error analyses discussed in this section, we do not have access to all 148 vertical excitation energies for every approach providing the reference data for judging the δ-CR-EOMCC methods. Thus, in making comparisons of the δ-CR-EOMCC and CC3 results, we have to limit ourselves to the subset of 139 excited states, for which the CC3 results are available [17,25]. When comparing the δ-CR-EOMCC and EOMCCSDT-3 excitation energies, we limit ourselves to the subset of 115 states, for which the EOMCCSDT-3 results are available [22,24]. Comparisons of our δ-CR-EOMCC excitation energies with the corresponding TBE-2 values are limited to the subset of 102 states, for which the latter values are available [21]. The only comparison that can utilise all 148 excited states, ob-tained by excluding the 1 1 B 3g (n 2 → π * 2 ) state of s-tetrazine from the 149 states listed in Table I of Schreiber et al. [17], in determining the corresponding MUE, MSE, and MaxE values is that involving the δ-CR-EOMCC and CASPT2 approaches, where the relevant CASPT2 data can be found in [17], with updates in [21]. Information about the numbers of excited states included in the statistical error analyses involving various methods considered in this work can be found in Tables 3-6.
Examining the error values in Tables 3-6 and the correlation plots shown in Figures 3-22, we can immediately see that all four δ-CR-EOMCC approaches considered in this work provide significant improvements in the EOM-CCSD excitation energies relative to the CASPT2, TBE-2, CC3, and EOMCCSDT-3 data. In the case of com- parisons with CASPT2 and TBE-2 (Tables 3 and 4 and Figures 3-12), we see steady error decreases when going from EOMCCSD, through δ-CR-EOMCCSD(T),IA, δ-CR-EOMCCSD(T),ID, and δ-CR-EOMCC (2,3),A, to   δ-CR-EOMCC (2,3),D, with the δ-CR-EOMCC (2,3),D approach providing the most accurate description which seems better than that provided by CC3 and EOMCCSDT-3, whereas comparisons with CC3 and EOMCCSDT-3 (Tables 5 and 6 and Figures 13-22) indicate a more uniform performance of all four δ-CR-EOMCC methods, although, as further elaborated on later, one may also argue that the overall agreement with the CC3 and EOMCCSDT-3 data continues to be best in the case of the δ-CR-EOMCC (2,3) triples corrections. Let us discuss these observations some more, starting with comparisons of the various CC/EOMCC approaches with the CASPT2 and TBE-2 data.
As shown in Table 3, the MUE, MSE, and MaxE values relative to CASPT2 decrease from 0.50, 0.49, and 1.63 eV, respectively, when the EOMCCSD results are considered, to 0.18, −0.06, and 0.53 eV, when the EOMCCSD excitation energies are corrected for triples using the δ-CR-EOMCC (2,3),D approach. Variant A of the δ-CR-EOMCC(2,3) method, which is characterised by the MUE, MSE, and MaxE values relative to CASPT2 of 0.17, 0.00, and 0.68 eV, respectively, is essentially equally accurate from the point of view of these three error measures, but the triples corrections of δ-CR-EOMCCSD(T),IA and δ-CR-EOMCCSD(T),ID, although offering significant improvements in the EOMCCSD results, are not as effective as their biorthogonal δ-CR-EOMCC(2,3) counterparts. This becomes particularly evident when the MSE and MaxE values relative to CASPT2 characterising the various δ-CR-EOMCC approaches in Table 3 are examined. Similar improvements in the accuracy of the calculated excitation energies relative to CASPT2, when going from EOMCCSD, through δ-CR-EOMCCSD(T),IA, δ-CR-EOMCCSD(T),ID, and δ-CR-EOMCC (2,3),A, to δ-CR-EOMCC(2,3),D, can be seen when analysing the correlation plots in Figures 3-7. Indeed, we observe a systematic decrease in the number of the outlier states whose energies differ from the CASPT2 data by more than 0. 75     are so substantial that one can regard them as competitive with the considerably more expensive iterative CC3 and EOMCCSDT-3 calculations. Indeed, as shown in Table 3, the MUE, MSE, and MaxE values relative to CASPT2 characterising the δ-CR-EOMCC (2,3),D excitation energies are smaller than those characterising the CC3 and EOMCCSDT-3 computations. This is especially true in the latter case, suggesting that the linear-response CC3 approach is somewhat more accurate than its EOMCCbased EOMCCSDT-3 analogue, when both are compared to CASPT2, and that our δ-CR-EOMCC (2,3),D calculations are even more accurate. We would not be surprised by such a statement if our statistical error evaluation involved doubly excited states, since the MMCC-based CR-EOMCC approaches are generally more robust than the perturbative CC3 and EOMCCSDT-3 approximations in applications involving quasi-degenerate excited states (cf., e.g., [63,70,71,74,75,160]), but it is interesting to see that similar might be true when the excited states of interest are of a predominantly single excitation character.
It is encouraging to observe the improvements in the overall accuracy of the EOMCCSD, CC3, and EOMCCSDT-3 calculations offered by the δ-CR-EOMCC(2,3),D approach, when compared to the CASPT2 and TBE-2 excitation energies, and the improvements in the accuracy of the EOMCCSD results by all four δ-CR-EOMCC methods examined in this study, but one would also like to know how accurate the various non-iterative triples δ-CR-EOMCC corrections are when compared with their iterative CC3 and EOMCCSDT-3 counterparts. This is examined in Tables 5 and 6 and Figures 13-23. The MUE and MSE values shown in Table 5 and the MUE, MSE, and MaxE values shown in Table 6 may create an impression that the δ-CR-EOMCCSD(T),IA and δ-CR-EOMCCSD(T),ID approaches are somewhat more accurate, when compared to the CC3 and EOMCCSDT-3 results, than their biorthogonal δ-CR-EOMCC (2,3),A and δ-CR-EOMCC(2,3),D counterparts, but given the fact that the MUE and MSE values characterising the average differences between the δ-CR-EOMCC and CC3/EOMCCSDT-3 excitation energies are all very small this might be a somewhat misleading conclusion. Indeed, when we look at the correlation plots in Figures 13-17, we can see a fairly systematic decrease in the number of the outlier states whose energies differ from the CC3 data by more than 0.50 eV when going from EOMCCSD, through the δ-CR-EOMCCSD(T) approximations, to the δ-CR-EOMCC (2,3) approaches, from 18 in the EOMCCSD case, to 3 and   (2,3),D approach, but knowing that the error distribution relative to the CC3 data is the narrowest in the δ-CR-EOMCC (2,3),D case and keeping in mind that the CC3 approach is not necessarily more accurate than the δ-CR-EOMCC (2,3),D method, especially when the excited states in question have some contributions from double excitations, makes us believe that the overall best approach among the four types of triples corrections to EOMCCSD excitation energies considered in this study is δ-CR-EOMCC (2,3),D, with variant A of δ-CR-EOMCC(2,3) offering a similar performance as long as the excited states of interest are not dominated by two-electron transitions, which is the case here. As already alluded to above, using the 1 1 B 3g (n 2 → π * 2 ) state of s-tetrazine as an example, and as shown in our past work [75,76], variants D of the CR-EOMCC(2,3) and δ-CR-EOMCC (2,3) methodologies are generally more robust than other CR-EOMCC/δ-CR-EOMCC corrections when the doubly excited states are considered. Most of the above remarks related to statistical comparisons of the δ-CR-EOMCC and CC3 excitation energies apply to the analogous comparisons of the δ-CR-EOMCC and EOMCCSDT-3 data. With the cut-off threshold of 0.50 eV used in the examination of the EOM-CCSD and δ-CR-EOMCC vs. CC3 and EOMCCSDT-3 results, we cannot say as much about the relative performance of the various non-iterative triples δ-CR-EOMCC corrections relative to EOMCCSDT-3 as in the case of the analogous comparisons with CC3, since all four corrections work equally well, producing no outliers, but we continue to observe a steady increase in the R 2 correlation factor when going from EOMCCSD, through δ-CR-EOMCCSD(T),IA, δ-CR-EOMCCSD(T),ID, and δ-CR-EOMCC(2,3),A, to δ-CR-EOMCC (2,3),D, from 0.9860 when the EOMCCSD and EOMCCSDT-3 excitation energies are compared to 0.9885, 0.9893, 0.9904, and 0.9904, when we compare the EOMCCSDT-3 data with their δ-CR-EOMCCSD(T),IA, δ-CR-EOMCCSD(T),ID, δ-CR-EOMCC(2,3),A, and δ-CR-EOMCC (2,3),D counterparts, respectively. This is reflected in the error distribution curves shown in Figure 23, where we can see that the error distribution relative to the EOMCCSDT-3 data becomes increasingly narrower as we go from EOMCCSD, through δ-CR-EOMCCSD(T),IA, δ-CR-EOMCCSD(T),ID, and δ-CR-EOMCC(2,3),A, to δ-CR-EOMCC (2,3),D. Once again, based on the MUE and MSE values collected in Table 6, one might argue that δ-CR-EOMCCSD(T),IA and δ-CR-EOMCCSD(T),ID approaches are more effective in reproducing the EOMCCSDT-3 results than their biorthogonal δ-CR-EOMCC (2,3),A and δ-CR-EOMCC (2,3),D counterparts, but knowing that the error distribution relative to the excitation energies obtained in the EOMCCSDT-3 calculations is the narrowest in the δ-CR-EOMCC (2,3),D case and keeping in mind that the EOMCCSDT-3 method is not necessarily more accurate than the δ-CR-EOMCC (2,3),D approach, particularly when the excited states in question have some contributions from double excitations, reinforces our belief that the overall best approach among the four types of triples corrections to EOMCCSD excitation energies investigated in this work is δ-CR-EOMCC (2,3),D, with δ-CR-EOMCC (2,3),A offering similar accuracies as long as the excited states of interest are not dominated by two-electron transitions.
By comparing our δ-CR-EOMCCSD(T),IA, δ-CR-EOMCCSD(T),ID, δ-CR-EOMCC (2,3),A, and δ-CR-EOMCC(2,3),D excitation energies for the 148 excited states of 28 molecules taken from Table I of Schreiber et al. [17], as described above, or their appropriate subsets for which the relevant reference CASPT2, TBE-2, CC3, and EOMCCSDT-3 data are available, we have shown that the non-iterative triples corrections to the EOMCCSD excitation energies defining the relatively inexpensive, singlereference, black-box δ-CR-EOMCC approaches provide significant improvements in the EOMCCSD data, while closely matching the results of the iterative and considerably more expensive CC3 and EOMCCSDT-3 calculations and their CASPT2 and TBE counterparts, typically to within ∼0.1 − 0.2 eV, i.e., to within intrinsic errors of the CC3, EOMCCSDT-3, CASPT2, and TBE estimates. We have also demonstrated that the δ-CR-EOMCC methods, especially the most robust δ-CR-EOMCC (2,3),D approach that works well for singly as well as doubly excited states, are capable of bringing the results of the CC3 and EOMCCSDT-3 calculations to a closer agreement with the CASPT2 and TBE data, demonstrating the utility of the cost-effective δ-CR-EOMCC methods in application involving molecular electronic spectra. This has allowed us to conclude that the overall best balanced approach among the four types of triples corrections to EOMCCSD excitation energies investigated in this work is δ-CR-EOMCC (2,3),D, with δ-CR-EOMCC(2,3),A offering similar accuracies as long as the excited states of interest are not dominated by two-electron transitions. We have reached these conclusions by performing a variety of full and partial statistical error analyses and examining the suitably designed correlation and error distribution plots.
We have also used the four δ-CR-EOMCC approaches considered in this study to identify and accurately characterise 54 additional singlet excited states in the energy range covered by Table I of Schreiber et al. [17], including five states that can be found in the Supporting Information to [17] and 49 states that have not been considered in the earlier benchmark work [17][18][19][20][21][22][23][24][25], which can be used in future benchmark studies. The aforementioned 1 1 B 3g (n 2 → π * 2 ) state of s-tetrazine, listed in Table I of [17], which could not be used in our overall statistical error analyses due to the absence of the reliable benchmark data to judge our δ-CR-EOMCC results, and 6 other states among the 54 states outside the set of 149 states listed in Table I of [17] are almost pure two-electron transitions, which many QC methods have problems with, but we have provided arguments, based on the successful track record involving various CR-EOMCC or δ-CR-EOMCC calculations, including quasi-degenerate excited states dominated by double excitations [63,70,71,[74][75][76]107,110,111,114,121,127,140,143] and the comparison of our best δ-CR-EOMCC (2,3),D excitation energies for the 1 1 B 3g (n 2 → π * 2 ) state of s-tetrazine with the recently published NEVPT2 data [23], that our δ-CR-EOMCC calculations for the doubly excited states found in this work and other additional states that have not been considered in the prior work [17][18][19][20][21][22][23][24][25] are accurate to within ∼0.2 − 0.3 eV. We have suggested full EOMCCSDT, active-space EOMCCSDt, or accurate MRCI calculations for all of the additional excited states found in our calculations to verify if our assessment of the accuracy of the δ-CR-EOMCC calculations for these extra states is correct. We have also suggested that there may exist a relationship between the REL diagnostic introduced in [74], resulting from EOMCCSD calculations, and the magnitude of δ-CR-EOMCC triples corrections to EOMCCSD excitation energies, which might serve as a potentially useful tool in predicting electronic excitation spectra of δ-CR-EOMCC quality only on the basis of EOMCCSD data.

Supplemental data
Supplemental data for this article containing symmetry unique Cartesian coordinates for the ground-state geometries of the 28 molecules comprising the benchmark set of Schreiber et al. [17] resulting from the MP2/6-31G * and Table A1. Mean signed errors (MSE), mean unsigned errors (MUE), and maximum energy differences (MaxE) characterising the extrapolated δ-CR-EOMCC (2,3),D excitation energies determined using Equation (A1) for the triples corrections δ μ relative to their true values.  Table I of Schreiber et al. [17]. b 48 singlet excited states obtained by removing six doubly excited states with REL ≈ 2 from the 54 additional excitations calculated in this work. c All of the excitations described in footnotes a and b. d All electronic excitations calculated in this work including seven doubly excited states with REL ≈ 2 found in our calculations.
the 54 additional excitations discussed in Sections 3 and 4.2 that are, in most cases, predominantly one-electron transitions, we obtained the MSE, MUE, and MaxE results of −0.12, 0.12, and 0.30 eV, respectively, which, given the extreme simplicity of Equation (A1), can be viewed as an encouraging result. When we apply the same empirical formula to the entire set of 203 excited states calculated in this work, which includes several doubly excited states with REL close to 2, the quality of the fit represented by Equation (A1) worsens, but it is still not bad given the nature of the procedure discussed here. Most likely, we would have to perform additional δ-CR-EOMCC calculations to produce a considerably larger number of data points involving doubly excited states to come up with a more accurate representation of the dependence of δ μ on REL for states dominated by two-electron transitions. It is also possible that states with REL close to 2 have to be fitted to a different functional representation than that represented by Equation (A1). Nevertheless, in order to see if Equation (A1) (or similar equations relating δ μ to REL that we may propose in the future) can be used in practice to generate the δ-CR-EOMCC-level information on the basis of EOMCCSD calculations for organic species outside the set of 28 molecules examined in this work, which may have a mixture of singly and doubly excited states in their electronic spectra, we have used it to determine the vertical excitation spectrum of azulene, examined by us in [140], where we have discovered the previously unknown doubly excited state in the higher energy part of the spectrum near the ionisation threshold. The results are shown in Table A2. As one can see, our simple empirical formula for δ μ , added to the EOMCCSD excitation energies, reproduces the true δ-CR-EOMCC (2,3),D data for azulene to within 0.06 eV on average. This includes the doubly excited S 9 state with REL of about 1.9, which has been the main objective of the study reported in [140].  Figure A1. Correlation plot for the triples corrections δ μ obtained in the δ-CR-EOMCC (2,3),D/TZVP calculations for the 148 singlet excited states of 28 molecules taken from Table I of [17] and discussed in the main text vs. the corresponding REL values determined at the EOMCCSD/TZVP level. The red line in the plot is a result of the least-squares fit.
All of this is telling us that there may exist approximate relationships between the REL diagnostic obtained in the EOMCCSD computations and the magnitude of δ-CR-EOMCC triples corrections to EOMCCSD excitation energies (most likely, other similar corrections as well), which might allow us to predict δ-CR-EOMCC-level spectra for organic molecules of the types included in the database of Schreiber et al. [17] and related species to within 0.1-0.2 eV on average and which could, given the simplicity of such a procedure, become useful tools in practical applications. We have to keep in mind that δ-CR-EOMCC and similar corrections to EOMCCSD excitation energies have computer costs that scale as n 2 o n 4 u in the iterative EOMCCSD part and n 3 o n 4 u in the noniterative steps related to the determination of triples corrections. If one can develop empirical formulas for the triples corrections of the type discussed in this appendix, we may avoid the most expensive n 3 o n 4 u steps and generate excitation spectra of δ-CR-EOMCC-quality with n 2 o n 4 u steps of EOMCCSD. Given the fact that δ-CR-EOMCC methods, especially δ-CR-EOMCC (2,3), are accurate enough to reproduce other high-level EOMCC and non-CC calculations, the procedure discussed in this appendix may help us develop empirical, but useful tools for determining reasonably accurate molecular electronic spectra in large molecules that are outside the reach of δ-CR-EOMCC and similar approaches.