Comparative study of methodologies for calculating metastable states of small to medium-sized molecules

Electron-induced reactions are ubiquitous in a myriad of chemical, biological, and astronomical processes.^1–8 High energy environments are associated with free electrons, and these can interact with molecules leading to the formation of states that facilitate chemical reactions through new pathways. Attachment of an electron to a neutral molecule yields an anion, which may be more (positive electron affinity) or less (negative electron affinity) stable than the parent neutral molecule. The latter case leads to a metastable state, which is unstable and can undergo autodetachment. In addition to that, stable anions usually do not support electronically bound excited states and as a result, excited states of anionic species are often metastable. The finite lifetimes associated with such species make their spectroscopy difficult, and therefore, it is essential to investigate them theoretically.^9,10 Since the anion resonance states are metastable, they are embedded within the infinite set of continuum states that converge to the neutral molecule and an electron at an infinite separation limit. The rate of electron-detachment is, therefore, dependent on the coupling between the resonance and continuum states. The associated total molecular wavefunction thus becomes non-L²-integrable and thus violates the boundary conditions required for Hermitian quantum mechanics.¹¹ The traditional Hermitian methods of electronic structure theory fail to describe such states, and some modifications to their foundation are needed.¹⁰ 

To date, several molecular systems that decay via dissociative electron attachment (DEA) have been investigated using a variety of theoretical methods.^12–15 Scattering calculations are a well-known example that contains a traditional foundation in the earlier studies of metastable resonances.¹⁶ Scattering theory has a firm theoretical foundation, and it involves the calculation of scattering cross sections. Commonly used methods are R-matrix theory, complex Kohn method, Schwinger multichannel method, etc.^17–19 Often, calculations based on scattering theory are devoid of appropriate correlation effects owing to the difficulty in its implementation, are computationally expensive, and are met with numerous convergence issues.²⁰ In addition, scattering calculations do not provide chemical insights since they cannot be readily extended to calculate potential energy surfaces and do not provide information about Frank-Condon factors, Dyson orbitals, etc.

Methods based on electronic structure theory are also a popular choice and exploit the existing formalism of quantum chemistry, without the need for invasive redevelopment.²¹ Methods based on non-Hermitian quantum mechanics have proven to be quite useful in uniquely identifying resonances, which are supported by a rigorous mathematical framework. The basic idea behind the non-Hermitian approach is revealed by solving the Schrödinger equation for the standard molecular Hamiltonian with outgoing boundary conditions, which gives rise to a non-Hermitian eigenproblem from which the resonance is obtained as a solution with Siegert energies,²² 

Here, E_R is the position of the resonance, while Γ is the width associated with the lifetime τ = ℏ/Γ. The non-Hermitian approach provides a way for extending quantum-chemistry methods developed for bound states to resonances. The complex energies are obtained by performing the analytic continuation of the existing Hermitian Hamiltonian into the complex plane. Complex scaling is one such method, where the electronic coordinates are scaled with a complex number, consequently rotating the string of pseudocontinuum states into the complex plane about the respective thresholds, thereby revealing the resonant states with square-integrable wavefunctions.^22–25 This procedure works very well in identifying atomic resonances, but it is problematic to apply to molecular systems.²⁵ Complex absorbing potentials (CAP) and complex basis functions (CBFs) represent an improvement on the complex scaling methods for computing molecular resonances. The application of CAPs to electronic resonances was put on a solid mathematical foundation by Riss and Meyer.²⁶ To date, many mainstream electronic structure methods [e.g., SAC-CI, DFT, ADC(2), XMCQDPT2, MRCI, and EOM-CCSD] have been augmented with a complex absorbing potential.^27–33 CAPs will be discussed in greater detail in Sec. II. The use of complex basis functions to determine molecular resonances is another rigorous approach, which was demonstrated for the first time by McCurdy and Rescigno.³⁴ The rigorous justification of CBFs as a finite basis approximation to exterior complex scaling was explored in the work of Morgan and Simon³⁵ and validated numerically by Ben-Asher and Moiseyev.³⁶ The use of complex scaled basis functions has recently been revitalized, with its efficient implementation at higher-level correlated theories (MP2 and EOM-EA-CCSD) for calculating shape and Feshbach resonances in small polyatomic systems.^37–39 

In an alternative approach, resonances can be obtained by analytic continuation of eigenvalues from the stabilization graph into the complex plane.^40–45 In this approach, a specified parameter is varied until an eigenfunction is identical to the resonance wavefunction localized in the interaction region of the potential. This leads to, what is known as, the “stabilization” of a given resonance. The energy of the resonance remains invariant to any further changes in the specified parameter, once stabilization has been achieved. The continuum and bound-state-like solutions interact as the scaling parameter is varied, which eventually manifests itself in the form of avoided crossings. Simons demonstrated that one can obtain an estimate of the desired lifetime directly from the stabilization graph highlighting a connection with the complex coordinate rotation method for which a satisfactory mathematical basis exists.⁴² Later, Thomson and Truhlar⁴⁶ used a single eigenvalue level from the stabilization graph, including both the stable part and the avoided crossing region, which was analytically dilated into the complex plane. McCurdy and McNutt suggested an analytic continuation method,⁴⁷ which was first developed by Jordan⁴⁸ using the branch structure near avoided crossings. When using electronic structure calculations, a common variable parameter is the exponents of the diffuse functions of the basis set. Other potential candidates for appropriate scaling parameters could involve variations in the induced external dipole or quadruple field, since such fields form a harmonic-type potential around the continuum wavefunction. Charge stabilization is another formalism in which the charge on the nuclei spanning the orbital with the extra electron is increased in infinitesimal amounts, and the generated curve is extrapolated to reveal the electron affinity of the metastable state.⁹ In the charge stabilization and other similar approaches, the introduction of a perturbation potential transforms resonance states into bound states. Analytic continuation can provide information about the Siegert energies.^49,50

In the present manuscript, we benchmark the efficiency of the CAP method against the orbital stabilization method (OSM). Our goal is to see how efficiently these methods can be used for medium-sized molecules (about 15 atoms) and to compare the results at the same underlying electronic structure theory. Initially, each method will be used to compute the metastable resonance anion states for the highly symmetric small systems: N₂, CO, and HCN, which have been studied extensively and there is a wealth of information to compare. Then, we will extend to the high symmetry and larger benzoquinone and finally to the lower symmetry uracil molecule. N₂, CO, and para-benzoquinone have been the subject of previous studies both with stabilization and CAP-based methods.^{29,43,51–54} When it comes to small molecules such as CO, N₂, and HCN, it is easy to perform the calculations with larger basis sets to assess the improvements in results as we extrapolate it to the complete basis limit. Uracil, on the other hand, is a more complicated system and one whose resonances are of particular interest. Hence, it serves as an excellent system to see how these methods perform for the larger molecules. This work has only considered the lowest shape resonances of all systems for consistency and for more robust comparisons, since the lowest resonances for these systems have been studied before experimentally and theoretically. The calculation of higher-lying resonances requires the calculation of more states, which makes the task more expensive, so we have left that out of the scope of the current work.

This section will describe the theoretical methods that are used in the present benchmark.

As outlined briefly in the Introduction, the stabilization method makes use of unmodified mainstream electronic structure methods by introducing a parameter, which ultimately identifies metastable resonances upon the infinite set of continuum states. In the OSM, additional diffuse functions are added and they are scaled by a scaling parameter (henceforth α). The α parameter is varied, progressively increasing the relative size of the augmented component of the basis functions. The electronic energy of the ground and lowest few excited states of the anion is then computed as a function of α using electronic structure methods. In our present case, the electronic structure theory of choice is the EOM-EA-CCSD method,⁵⁵ as implemented in the Q-Chem computational package,⁵⁶ as it is a very accurate method for computing electron attachment energies. Once calculated, the many electronic states show distinct variations in energy as a function of α, depending on whether the states are true resonance or continuum states. The resonances are almost independent of α, while the continuum states show large fluctuation upon variation in α.²¹ Avoided crossings between the resonance and continuum states reveal the resonance-to-continuum coupling strength. Such critical points therefore represent regions in which analytical continuation can be used in order to compute the resonance positions and lifetimes (width). There have also been methods that use regions away from the avoided crossing to apply analytic continuation,⁵⁷ but in this work, we will only use the avoided crossings. Using analytic continuation, the complex stationary points and energies can be located using various forms of the Generalized Padé Approximant (GPA). The GPA was introduced by Jordan for studying the temporary anion states of molecular systems.^48,58 In the present study, the GPA is used in the form shown in Eq. (2), where the real energy is expanded as a quadratic polynomial,^48,58,59

The quantities P, Q, and R in Eqs. (1) and (2) are themselves polynomial expansions of the scaling parameter (α) with coefficients p_i, q_j, and r_k, respectively,

Equation (2) is henceforth denoted as the (i, j, k) General Padé Approximants. The indices i, j, and k are held equivalent, giving rise to 11, 14, and 17 unknowns for the (3, 3, 3), (4, 4, 4), and (5, 5, 5) approximations, respectively. These unknowns are found by fitting to the calculated EOM-EA-CCSD energies. Higher order expansions can be used as well, when more than two states are involved in the avoided crossings.

Once P, Q, and R are computed, the location of the complex stationary points (α^*) may be determined by calculating the roots of the following equation:^43,59,60

Upon optimization of the stationary points, the computed values are then substituted back into Eq. (2) in order to obtain the complex energies.

Often many avoided crossings are obtained for the same resonance. In the past, we used all of these avoided crossings to obtain the complex energies.⁶¹ The avoided crossings for lower values of α, however, should be more accurate since the basis set is more diffuse and can represent the resonance better. For this reason, in this work, we will use the avoided crossing located at the smaller values of α as our more accurate estimate. Results from the other avoided crossings will also be presented in order to validate this assumption.

The underlying foundation of the CAP method is the augmentation of a standard Hamiltonian (H) with a complex potential. The resulting CAP-augmented total Hamiltonian is then given by the following equation:

In the current work, W of Eq. (5) represents a quadratic potential^{28,29,51,62–64} given by Eq. (6) and contains the conditions given below,

In Eq. (7), $rj0$ indicates the j = x, y and z Cartesian coordinates associated with the size of the box. The CAP is therefore dependent on four parameters (i.e., the three Cartesian coordinates associated with the box size and the CAP strength η). In the complete one-electron basis, the exact energy of a given metastable resonance, in the complex plane, is obtained as lim_η→0E(η). Conceptually, this indicates that the resonance position is stabilized with an infinitesimally small CAP, at which the exact energy of the resonance is represented in the infinite basis set limit. The inevitable use of finite basis sets in electronic structure theory, however, requires the computation of a series of energies by varying values of η in order to identify an optimal value for the η at which the resonance position is stabilized (henceforth η_opt).^26,65 A common way in which η_opt is obtained is by finding the minimum of the logarithmic velocity, as given by the following equation:

The introduction of CAP makes the results sensitive to the parameters used, specifically the size of the box used, since the resonance wavefunction in the interior region will be somewhat unphysically disturbed. Within this representation, the energetic positions and widths of a given resonance have been shown to be very sensitive to the CAP onset.^28,51,65 Krylov and co-workers have recently provided an alternative approach to finding η_opt.⁵¹ This involves augmenting the total complex eigenvalues (will be referred to as zeroth-order from now on) (E^R and E^I) with the trace of the one-particle density matrix (γ) times W, as given in the following equations:

The resulting first-order corrected complex eigenvalues (U^R and U^I) display stationary behavior asymptotically at large η values, beyond the stabilization point. In this case, we find that trajectories can be harder to analyze. Alternatively, one can obtain the minima in the derivatives of U^R(η) and U^I(η).

In our present benchmark manuscript, we are using the recent implementation of CAP into the stable EOM-EA-CCSD method (see Sec. II C for more details) in the Q-Chem computational package,⁵⁶ which ultimately computes U^R and U^I. The uncorrected (or zeroth-order) affinities and widths are obtained from the turning point of the η-trajectory (E^R vs E^I), which, as mentioned earlier, corresponds to the minimum of the logarithmic velocity. The corrected (or first-order) affinities and widths are obtained from the minimum of $dURdη$ and $dUIdη$⁠, respectively. In the literature, both the turning points of trajectories and the minima employed here have been used.^51,52 In our work, we found that it is not always easy to interpret the turning point of corrected trajectories(U^I vs U^R) by mere visual inspection, while it is easier to use the minima of $dUdη$⁠, so we are employing this approach here.

Using the Q-Chem computational package,⁵⁶ the equation-of-motion coupled-cluster singles and doubles method, specifically the electron attachment version (i.e., EOM-EA-CCSD),⁵⁵ was used in all calculations in the present manuscript. The EOM-EA-CCSD method is capable of producing reliable theoretical estimates of attachment energies, and it can be systematically improved to obtain the exact solution. It is size-extensive, and an apt description of dynamical and static correlation is included in one computational step. These methods are also devoid of any system-dependent parameterization.^29,66 In all systems, the CAP method and OSM were used and benchmarked against one another. It should be noted that for CAP-EOM-EA-CCSD calculations, the CAP is added at the HF stage. Details of the methodology for each molecule are given below, while the choice of basis set and box size is summarized in Table I.

For N₂ and CO, Dunning’s augmented correlation-consistent basis set of triple-ζ quality, aug-cc-pVTZ, was used. Extra sets of 3s and 3p diffuse functions were added to the atoms of each molecule—giving the total basis: aug-cc-pVTZ+[3s3p]. This basis set was the same as that used by Jagau et al.⁵¹ We also used the same ground state equilibrium geometries used in that work.⁵¹ The size of the boxes in this case was also taken directly from the literature⁵¹ since we are trying to reproduce the CAP results and compare them to the OSM. We report the results at three box sizes in the supplementary material and retain the results for only one box size (centered at ⟨R²⟩) in the main paper.

For HCN, the ground state geometry of the neutral molecule was optimized at MP2 theory with Dunning’s augmented correlation-consistent triple valence basis set, aug-cc-pVTZ. Calculations for OSM and CAP were once again done using EOM-EA-CCSD with aug-cc-pVTZ+[3s3p] and cc-pVTZ+[2s5p2d] basis sets, with the additional set of diffuse functions added on heavy atoms in both basis sets. The “s” and “d” functions in the latter basis set were scaled with a factor of 2, and the “p” functions were scaled with a factor of 1.5. The aug-cc-pVTZ+[3s3p] basis set has been used by Krylov to describe resonances in small molecules, and the cc-pVTZ+[2s5p2d] was used by Ehara et al.⁶⁷ in their paper to describe anionic resonances in the HCN anion with SAC-CI. The parameters of the box size, in this case, were taken from ⟨R²⟩, which has been proposed as a default value for box size in CAP calculations.²⁹ 

In the case of benzoquinone, the anionic ground state geometry was optimized at the EOM-EA-CCSD/aug-cc-pVDZ level of theory. The OSM and CAP EOM-EA-CCSD calculations were carried out using Dunning’s correlation consistent cc-pVTZ basis set with extra sets of three “s” functions added to ghost atoms 1.45 Å above and below the plane of carbon atoms, similarly to the functions used by Kunitsa and Bravaya.⁵² The resulting basis set is designated as cc-pVTZ+[24s]. Two box sizes have been tested, one is centered around the ⟨R²⟩ expectation value at the ground state of the radical anion, while the other one is larger and is taken from the work of Kunitsa and Bravaya.⁵² The results of the larger box size are reported in the main paper, while the results from ⟨R²⟩ are reported in the supplementary material.

The ground state geometry of neutral uracil was optimized at the B3LYP/cc-pVTZ level of theory. For calculations on the metastable anion of uracil using EOM-EA-CCSD, Dunning’s augmented correlation-consistent double zeta valence basis set was used and additional basis functions of s, p, and d character were added on the heavy atoms atom in an even tempered manner resulting in aug-cc-pVDZ+[spd]. In addition to that, we have also computed the resonances using the cc-pVDZ+[2s5p2d] basis set with the 2s5p2d diffuse functions added to heavy atoms.²⁷ The set of “s” and “d” functions was scaled by a factor of 2, and the “p” was scaled by a factor of 1.5. A box size of ⟨R²⟩ is used in this case as well.

Dyson orbitals were calculated for the resonant states of the listed molecules using CAP-EOM-EA-CCSD as implemented in QChem,⁵⁴ and the orbitals are summarized as well.

We start by exploring the resonances of anions of the diatomic molecules N₂ and CO as test cases since they have been studied extensively in the past. Figure 1 shows the stabilization curves associated with the ground ²Π_g state of $N2−$ and the ground ²Π state of CO⁻, computed using EOM-EA-CCSD. The avoided crossings are indicated by xi. The top part of Tables II and III lists the theoretical resonance positions and widths obtained from these curves, which are compared to experimental values.^68,69 Each AC has been analyzed using three GPA polynomial expansions, as indicated by (n, n, n). Within the same avoided crossing, the different GPA expansions show minor variations in the positions and widths, indicating a robustness of the resonances with respect to GPA expansions. Results for the different avoided crossings are also shown in Tables II and III in order to highlight the differences. There is a ≈0.2 eV and ≈0.1 eV (for CO and N₂, respectively) difference in the positions when considering the different avoided crossing. The associated resonance widths show larger variations, especially for CO (0.5 eV). The width associated with the avoided crossing at the smallest value of α (more diffuse functions) however shows the best agreement with the experiment in both cases, as one would expect. This agrees with our choice of using the first AC for more accurate results.

Comparative computations, using CAP-EOM-EA-CCSD, are listed at the bottom of Tables II and III and in Fig. 2. Figure 2 shows the first-order corrected energies (real and imaginary) as a function of changing CAP strength for CO and N₂ anions, respectively. In addition to that, the zeroth order real energies are also plotted against the imaginary energies (η-trajectories) to locate the turning point, which corresponds to the stabilized values for the resonance. The optimum value of the CAP strength, at which the resonant states of CO⁻ and $N2−$ are stabilized, is summarized in Table III of the supplementary material (denoted as SI-III, and similarly denoting all figures and tables in the supplementary material). Figures SI-1 and SI-3 in the supplementary material show plots of the logarithmic velocity $ηdEdη$ to indicate where it is minimum. We use these plots to confirm the turning point of the trajectory. In addition, plots of ⟨R²⟩ are also shown in Figs. SI-1 and SI-3. A stabilized value of this property is another indication that a true resonance is located.

Resonance positions and widths obtained from these plots are shown in Tables II and III. In CAP, the box-size is a user-defined parameter that requires optimization, and the choice of the size of the box causes variation in the results. We performed the CAP calculations at three different box sizes taken from previous studies, as indicated in the methodology. Previous work has recommended a consistent use of a box size corresponding to ⟨R²⟩ of the molecule, so we mainly report these results here. The results for the remaining box sizes are summarized in Tables SI-I and SI-II of the supplementary material. We emphasize, however, that the standard deviations with respect to the box size are moderate (0.01–0.03 eV for the position of resonances taken from corrected energies, and 0.02–0.06 eV for the widths of resonances, again taken from corrected energies). Focusing on the results reported in Tables II and III for the ⟨R²⟩ box, we can compare the corrected and uncorrected results. The difference in positions is 0.13 eV and 0.03 eV, while the corresponding difference in widths is 0.16 and 0.14 eV for CO and N₂, respectively. Based on previous work, we will use the corrected results as our optimum ones to compare to the OSM. It should be noted that for the same box size and basis set, our results are somewhat different from Krylov and co-workers because we use the U vs η plots to obtain resonances, while they used the trajectories.

Resonances for both $N2−$ and CO⁻ have been studied extensively in the past with a variety of methods including OSM and CAP.^{29,33,51,70–82} Our purpose in including these results is to benchmark our approaches comparing to the previous results and to compare OSM and CAP using exactly the same electronic structure theory. Falcetta et al.⁴⁵ have examined the performance of orbital stabilization methods using a variety of electronic structure methods and applied to various molecules, including N₂. In that work, the EOM-CCSD method and the aug-cc-pVTZ basis set were used as well, but they only used p extra diffuse functions. Our results are very similar to that work. They also used larger basis sets, aug-cc-pVQZ, or aug-cc-pV5Z, which made some small improvements to the results lowering the position by 0.1 eV and the width by 0.07 eV.

Krylov and co-workers have used CO and N₂ in their CAP studies. They showed that the resonance positions are heavily dependent upon the basis set.⁵¹ Specifically, they showed that increasing the number of additional “s” and “p” diffuse functions on the aug-cc-pVTZ basis set leads to a substantial lowering of the first-order resonance position of the ²Π state of CO⁻, from 2.517 eV to 2.081 eV,^29,51 approaching the experimental value. In contrast, the studied box-sizes made only minor variations to the width of the resonance at a given basis set. They found that the optimal box-size is centered around the ⟨R²⟩ expectation value of the total ground-state electronic wavefunction of the neutral molecule. Landau and Moiseyev have applied a different approach to remove the perturbations due to CAP using analytical continuation with Padé approximants. They demonstrated this approach to the resonances of CO and N₂ anions.⁸³ Their approach gave similar results to Krylov and co-workers when the same electronic structure and basis set were used, and improved results with larger basis sets.

To compare the two methods, OSM and CAP, with each other, we pick the OSM values with the (5, 5, 5) expansion and the x1 avoided crossing, and the first-order corrected CAP values at the box size corresponding to ⟨R²⟩. For CO, the difference in position is 0.20 eV, while the difference in width is 0.29 eV. For N₂, the position differs by 0.05 eV, while the width differs by 0.26 eV. The two methods give different results, even though the underlying theory is the same.

Recently, White et al.³⁹ have explored the resonances in small molecules, including N₂ and CO anions, using CBFs and a variety of electronic structure methods, including EOM-EA-CCSD. This provides an opportunity to compare all three methods (although they used a somewhat different basis set than ours). All three methods, OSM, CAP, and CBFs, give very similar results for the positions of the resonances in N₂ and CO anions, with a standard deviation of 0.02–0.1 eV. The widths however have much larger deviations, 0.17–0.29 eV, highlighting once again the difficulty in obtaining accurate widths. The widths are the largest when using CBFs and the smallest when using CAP.

With both the CAP and OSM methods, the returned theoretical resonance positions over estimate the experimental value. This is likely a manifestation of the inadequate size of the basis set. As already discussed above, use of larger basis sets improved the agreement with the experiment in CO. The inclusion of correlation beyond singles and doubles is also a non-negligible factor.^84,85 Jagau showed that the inclusion of triple excitations, in some cases, tends to correct the resonance parameters with the same magnitude as that of first-order correction.⁸⁵ It must be emphasized that a direct comparison of Siegert widths to experimental ones is slightly precarious owing to its sensitivity to nuclear positions.⁸⁶ The current results only examine the electronic contributions to the positions and widths, while the experimental values include vibrational contributions.

Figure 3 shows the Dyson orbitals describing the shape resonance in CO⁻ and $N2−$⁠. Both the real and the imaginary part are shown, and they agree with similar results shown by Jagau and Krylov.⁵⁴ The real Dyson orbitals have a similar shape to the valence π orbitals that we expect the electron to be attached, validating the characterization of these solutions as resonances, rather than pseudocontinuum states. As discussed by Jagau and Krylov, the symmetry of the Dyson orbital can also provide information about the wavefunction of the outgoing electron.^54,87 Experiments have shown that the outgoing electron in $N2−$ is of pure d-character (as shown in the symmetry of the corresponding Dyson orbital), while the outgoing electron in CO⁻ has also some p-wave contributions.^68,88 In summary, the Dyson orbitals provide useful characteristics of the resonances.

HCN, a prototypical triatomic molecule, has been the subject of many studies throughout the years. It is known to exist in high abundance in space, especially the interstellar medium. Recent studies have highlighted the importance of HCN species as the principal source of carbon and nitrogen in the prebiotic formation of nucleobases because of its rich presence in the early atmosphere on the Earth.⁹⁰ 

The stabilization curves obtained in this work for HCN are shown in Fig. 4. The resonances obtained are reported in Table IV. For the OSM calculations performed using the aug-cc-pVTZ+[3s3p] basis set, the position of the resonance and the width show sensitivity to the GPA expansion unlike the cc-pVTZ+[2s5p2d] basis, where both the position and width of the ²Π resonance remain well precise for different orders of expansions of the GPA polynomial as seen from the numbers summarized in Table IV. The results obtained from the latter basis set are suggestive of its better performance compared to the former.

In the CAP studies, only one box size was used corresponding to the expectation value of ⟨R²⟩. Figure 5 shows the variation in zeroth-order and first-order complex energies with CAP strength. In addition, the zeroth-order real energies are plotted against corresponding imaginary energies for the estimation of trajectory-stabilized resonance position and width. The optimum value of CAP strength corresponding to the stabilized resonance is reported in Table SI-IV. The variation in absolute values of the first-order correction and ⟨R²⟩ with CAP strength are summarized in Figs. SI-5 and SI-6. Table IV reports the resulting resonances. The difference between the first and zeroth-order positions is between 0.06 and 0.14 eV, while the difference seen in the widths is 0.15–0.3 eV, depending on the basis set used.

Using what we believe are the best estimates from OSM and CAP, as we did for the diatomics, we see that the differences in the position between the two methods are 0.29 eV and 0.18 eV at aug-cc-pVTZ+[3s3p] and cc-pVTZ+[2s5p2d] basis sets, respectively, while the difference in widths is 0.56 eV at aug-cc-pVTZ+[3s3p] and 0.67 eV at cc-pVTZ+[2s5p2d]. Although the positions are not very different, the two methods predict very different widths. As there is no experimental value for the widths, it is hard to estimate the correct value. All calculated values are relatively large, indicating a broad resonance with short lifetime. In such cases, it is expected that the uncertainty in the calculations will be large. As the width increases, one would expect a need for very diffuse basis sets to describe it accurately, increasing the sensitivity to the basis set. Practically, in the OSM, when the width is large the avoided crossings become less well defined, making the analysis less accurate.

Electron transmission studies of selected cyano compounds by Burrow and co-workers found the π^* resonance position of the HCN anion at 2.26 eV.⁸⁹ Ehara et al., using the CAP/SAC-CI methodology and the cc-pVTZ+[2s5p2d] basis set, have predicted the resonance position for the HCN anion in the range of 2.51–2.52 eV and the width being 1.09–1.10 eV.⁶⁷ Their calculations employ a smooth Voronoi potential,⁹¹ unlike the conventional box potential used in this work. Both our OSM and CAP results are qualitatively similar to these CAP/SAC-CI results. More specifically, the CAP results differ by 0.02 eV in position and by 0.4 eV in width. The calculated position obtained with the cc-pVTZ+[2s5p2d] basis set, and CAP is closer to the experimental value, being overestimated by about 0.4 eV compared to the experimental one. However, when using OSM, the aug-cc-pVTZ+[3s3p] gives better agreement with the experiment. From these observations, it is difficult to assess which basis set is better since the performance seems to be linked to the method for extracting resonances as well. There is no experimental value for the width reported in the literature, so we cannot compare our values.

Benzoquinone is a high symmetry aromatic molecule, which has recently attracted attention from the negative-ion community.^92–94 Quinones, in general, act as excellent electron carriers, which is evident from their roles in important metabolic processes such as respiration and photosynthesis.^95–97 Naturally occurring quinones and their synthetic derivatives are also being considered as potential drug candidates because of their effectiveness in suppressing tumor growth and hence are of pharmacological interest.⁹⁸ The simplest of the quinones, para-benzoquinone (p-BQ), has been the subject of plenty of theoretical as well as experimental studies.^99–106 The shape ²A_u resonance of the p-BQ anion obtained at various levels of theory is tabulated in Table V. Cheng et al., using the orbital stabilization method, have studied the shape and core-excited resonances of the p-BQ anion within the framework of density functional theory and Koopmans’ theorem (denoted as SKT).⁵³ The R-matrix method was also used in both the Static-Exchange plus Polarization (SEP) and close-coupling (CC) models. The more sophisticated CC approach gave a position of 1.1 eV and a width of 0.03 eV. The theoretical predictions in previous works underestimate the experimental position by 0.36–0.52 eV.

In our work, we employed both OSM and CAP to obtain the ²A_u resonance. Figure 6 shows the stabilization curves obtained at the EOM-EA-CCSD/cc-pVTZ+[24s] level of theory. There are two avoided crossings corresponding to the ²A_u resonance, but we were only able to apply the AC for one of them. The results are shown in Table VI. The position of the resonance is found to be at 0.86 eV and the width 0.008–0.022 eV (depending on GPA). The position is in good agreement with many of the previous results shown in Table V and is smaller than the experimental value (Cooper et al.¹⁰²) by 0.5 eV. There is a great variation in the values obtained for the width in previous studies, and the OSM results we obtain are within this wide range. The position is very stable with the various GPA expansions, but the width changes from 0.008 eV in the (3, 3, 3) GPA to 0.022 eV for (5, 5, 5). We expect the value of 0.022 eV obtained from the larger GPA to be more accurate.

CAP based EOM-EA-CCSD calculations are shown in Table VI based on the trajectories and complex energies shown in Fig. 7. The box size presented here is taken from the work of Kunitsa and Bravaya.⁵² An additional box size centered around the ⟨R²⟩ expectation value at the ground state of the radical anion has been used, and the results are shown in the supplementary material. In the case of benzoquinone, the additional basis functions are located away from the plane of carbon atoms. The optimum η-value for the stabilized resonance is reported in Table SI-VI. These results reveal that the zeroth and first order positions in this case are almost identical and range between 0.90 and 0.99 eV depending on the box size. These values are in close agreement with the previously obtained positions. The width is more sensitive to both the box size, and whether the corrected energy is used to obtain it. The zeroth order width is 0.012–0.033 eV, whereas the corrected width is half of that at 0.005–0.016 eV. Our results are almost same as those reported by Kunitsa and Bravaya since we used the same parameters, except for an anion geometry optimized at somewhat different levels.

Comparing OSM and CAP, we see that the position differs by 0.1 eV, while the width by 0.005 eV. Interestingly, these differences are smaller than what we saw in the smaller molecules. The width is in reasonably good agreement with the experimental value, while the position from both CAP and OSM is smaller than the experiment, with the CAP result being somewhat more accurate.

Uracil is the largest system we examine in this work. Our ultimate goal is to determine which methods can be used for systems such as uracil. The resonances of uracil may be very important in determining the mechanisms for radiation damage in nucleosides, and our group has shown much recent interest.^61,108 Experiments have shown that low-energy electron irradiation (5–20 eV) leads to dissociative electron attachment (DEA)—yielding CO and H products.^109–112 Such fragmentation represents the instability of a biologically relevant system upon electron irradiation, and it is, therefore, paramount to be able to model the intrinsic mechanisms responsible for such decay. Recent work in our group showed that the DEA of uracil involves a multitude of anion electronic states with reach nonadiabatic coupling en route to CO and H elimination—forming an imidazole derivative.¹⁰⁸ In this study, the importance of effectively computing the excited states of the anion, as well as the ground state, was particularly portrayed. In other recent work, we have computed the lowest three anion resonance positions and widths of uracil using the OSM method.⁶¹ Here, we will focus on comparisons between OSM and CAP for obtaining the lowest π^* resonance in uracil.

Figure 8 shows the EOM-EA-CCSD energy as a function of the basis set scaling parameter α for two basis sets: aug-cc-pVDZ+[spd] and cc-pVDZ+[2s5p2d]. The stabilization curves yield two well-defined resonances represented by avoided crossings, which may be analyzed further using AC, coupled with GPA methods. Table VII presents the stabilized resonance positions and widths for the first lowest shape resonance using the avoided crossings and the various GPA polynomial expansions. The computed resonance energies overestimate the experimentally measured electron affinity—highlighting, again, that the basis set is too small to fully reproduce the true resonance energy. We also examine the basis set dependence by comparing the aug-cc-pVDZ+[spd] basis set with the cc-pVDZ+[2s5p2d], which was also used by Ehara et al.²⁷ 

The OSM results [from (5,5,5) GPA expansions] obtained using both basis sets once again overestimate the experimental value by 0.43–0.47 eV, similar to what we have seen in our previous examples. The positions and widths computed at the aug-cc-pVDZ+[spd] basis set remain consistent with respect to expansions of the GPA, while the same results obtained with the cc-pVDZ+[2s5p2d] basis set seem to vary for different expansions of the GPA. For the (4, 4, 4) and (5, 5, 5) GPA results obtained with the cc-pVDZ+[2s5p2d] basis, the calculated positions agree closely with the ones obtained with the aug-cc-pVDZ+[spd] basis. There is a difference of 0.04 eV in the computed widths at (5, 5, 5) expansion of the GPA between both basis sets.

We also obtained the electronic resonances for uracil at the CAP-EOM-EA-CCSD level of theory. Results are shown in Fig. 9 and Table VII for both basis sets. We used a box size centered around the expectation value of ⟨R²⟩ at cc-pVDZ+[2s5p2d] and chose a slightly larger box size for the aug-cc-pVDZ+[spd] basis, where the x coordinate of the box is at ⟨R²⟩ of the system and the other coordinates are similar to the larger box used for benzoquinone. We chose a slightly larger box size for the latter basis due to some convergence issues encountered while using the optimum ⟨R²⟩ box size. The first-order corrected electron affinity is 0.729 eV with aug-cc-pVDZ+[spd] basis and 0.718 eV with cc-pVDZ+[2s5p2d] basis, overestimating the experimental value by ≈0.5 eV. Figures SI-13 and SI-14 summarize the variation of logarithmic velocity and ⟨R²⟩ as a function of CAP strength. The Dyson orbitals corresponding to the stabilized anionic state are also presented in Fig. 3. As discussed earlier, the real Dyson orbital corresponds to the valence π^* orbital where the electron is attached.

Comparison between OSM and CAP shows that the positions obtained from the two methods differ by 0.03–0.08 eV, while the widths 0.01–0.04 eV, depending on the basis set used. Overall, again the differences between the two methods are smaller here than what we observed for the smaller molecules. The above analysis demonstrates the comparable performance of both OSM and CAP methods, provided that the basis set is good enough.

In this manuscript, we have attempted to study the performance of two methods for computing metastable anion resonances. Specifically, we have used CAP and OSM to assess the successes and shortcomings of each method in modeling the lowest resonances of diatomic and triatomic molecules (CO, N₂, and HCN) and the feasibility to which each method extends to medium-sized molecules (benzoquinone and uracil). The results show that both CAP and OSM are adequate for computing the ground state anionic resonance positions and widths of diatomic and triatomic molecules, without a clear winner when it comes to agreement with the experiment. The widths obtained from the two methods differ by up to 0.7 eV, especially when attempting to describe wide resonances. In general, it is more challenging to describe resonances with large widths and short lifetimes. It must be emphasized that the experimental widths are marred by numerous assumptions making them less reliable than the experimental positions. Therefore, theory may be more suited to provide more reliable widths, if a higher-level treatment encompassing high-degree of correlation is used.

Extrapolation of both methods to medium-sized molecules can be done, although the size of the basis sets has to be reduced. Interestingly, our results show that the two methods have better agreement with each other for these larger molecules. The positions differ by less than 0.1 eV, while the widths differ by less than 0.05 eV.

As we have already seen, both methods require the calculations of energy profiles with respect to a primary parameter. In the OSM, it is the exponent of the diffuse Gaussian functions and in CAP based methods, the energies are obtained over a range of CAP strengths. For the calculations on the anionic resonance in HCN, the OSM required the calculations of some ≈20 points and the CAP based EOM-EA-CCSD method required the calculation of ≈40 or more points. The calculations for this particular example with OSM/cc-pVTZ+2s5p2d took a cumulative time of 5 h in wall time [3 h in central processing unit (CPU) time], whereas the estimated cumulative time for CAP-EOM-EA-CCSD/cc-pVTZ+2s5p2d was 71 h in wall time (71 h in CPU h as well). It is obvious that while the CAP method works well for small molecular systems, it is computationally demanding to use in larger systems with large basis sets. On the other hand, while the OSM is based on a weaker theoretical basis than CAP, it shows similar results to CAP without obvious failures and could be used as a cheaper alternative in larger systems. The stabilization method has also been recently used to compute the polyatomic complex potential energy surfaces.^84,114 A strong advantage of CAP, though, is the ease at which it can be extended to compute the molecular properties of metastable states, thereby assisting in its characterization, and the inclusion of analytic gradients and other tools has opened up whole new avenues for the exploration of potential energy surfaces, exceptional points, and crossing regions between neutral and metastable anionic states.^115–117 

See the supplementary material for additional calculations as outlined in the text.

The authors acknowledge support from the National Science Foundation under Grant No. CHE-1800171. This work used the Extreme Science and Engineering Discovery Environment (XSEDE),¹¹⁸ which was supported by the National Science Foundation (Grant No. ACI-1548562).