Theoretical approaches for treating non-valence correlation-bound anions Free

Bound anionic states of molecules can be classified as valence-bound or non-valence-bound depending on the nature of the interaction between the excess electron and the molecule. In valence-bound anion states, the excess electron is bound to the molecule by short-range interactions and occupies an orbital localized in the valence region. By contrast, in non-valence-bound anion states, the binding of the excess electron is dominated by long-range electrostatics or dispersion-type correlation effects or a combination of these two effects, and the excess electron occupies an orbital that has relatively little weight in the valence region.^1–13 Dipole-bound anions are a well-known class of non-valence anions for which electrostatic interactions dominate the electron binding. Within the Born-Oppenheimer approximation, a molecule has to have a dipole moment in excess of 1.625 D to possess a dipole-bound anion, whereas the inclusion of non-Born-Oppenheimer corrections causes the critical moment to depend on the moment of inertia and to assume a value closer to 2.5 D.^1,2,11 Examples of systems for which long-range correlation effects dominate the binding of the excess electron are sufficiently large Xe_n clusters,³ the s-type anion states of C₆₀,⁷ and large acenes.⁹ In addition to these limiting cases, there are systems, e.g., succinonitrile⁸ and certain structures of (NaCl)_n clusters,¹⁰ for which both electrostatics and long-range correlation effects play an important role in the binding of the excess electron. We refer to non-valence anions for which inclusion of long-range correlation effects is essential for electron binding as non-valence correlation-bound (NVCB) anions.

For non-valence anions for which the binding of the excess electron is dominated by long-range electrostatic interactions, both the Hartree-Fock (HF) and Koopmans’ theorem (KT)¹⁴ approximations give a bound anion, provided that sufficiently flexible basis sets are employed. For these species, quantum chemistry methods such as coupled-cluster singles and doubles plus perturbative triples [CCSD(T)]¹⁵ and, in some cases, even second-order Möller-Plesset perturbation theory (MP2)¹⁶ can accurately predict the electron binding energies (EBEs) and other properties of the anions. By contrast, NVCB anions, by definition, are unbound at the HF level of theory. In fact, when flexible basis sets are employed for such systems, the HF wave function for the excess electron system actually describes the neutral system plus an electron in a discretized continuum orbital. As a result, electronic structure methods such as MP2 and CCSD(T) that start from the HF approximation are expected to fail to bind the excess electron when flexible basis sets are employed. Thus, it is of considerable interest to establish which electronic structure methods are suitable for describing NVCB anions. In this work, we consider a (H₂O)₄ cluster with zero net dipole moment and tetracyanoethylene (TCNE) as model systems for accessing the performance of various theoretical methods for describing NVCB anions. We also consider the case of CO₂⁻ for which the inclusion of correlation effects is essential for the binding of the excess electron as the bending potential of the anion approaches that of the neutral. Although, the wave function of the ground electronic state of CO₂⁻ acquires considerable non-valence character at these geometries, the nature of the correlation contribution appears to be different from that of typical NVCB anions.

The theoretical methods considered include HF, MP2, second-order algebraic diagrammatic construction (ADC(2)),¹⁷ CCSD, CCSD(T), various equation-of-motion (EOM) coupled cluster methods,^18–22 orbital-optimized MP2 (OO-MP2),^23,24 Brueckner coupled-cluster doubles with perturbative triples (BCCD(T)),²⁵ and the G₀W₀ Green’s function method.²⁶ In the EOM, ADC(2), and G₀W₀ methods, the electron binding energy (EBE) is calculated directly. For the remaining methods, the energies of the neutral and anionic systems are calculated separately, with spin-unrestricted methods being employed for the anions, and the EBEs are obtained from the difference of the energies of the neutral and anionic species. These are referred to as delta (Δ) methods.

The ADC(2) and G₀W₀ methods are Green’s function methods, with the former employing a second-order approximation to the self-energy. Unlike the CCSD(T) and MP2 methods, the BCCD(T) and OO-MP2 methods employ orbitals that are optimized in the presence of correlation effects. The BCCD(T) method optimizes the orbitals so as to eliminate single excitations to all orders in the inter-electron interaction, while the OO-MP2 method optimizes the orbitals in the presence of the second-order correlation effects. The EOM approximations considered are EOM-CCSD(2) (here referred to as EOM-MP2),^18,19 EOM-CCSD,²⁰ EOM-CCSDT,²¹ and EOM-CCSD(T)(a)*.²² In the EOM-MP2, EOM-CCSD, and EOM-CCSDT methods, the energy of the neutral system is calculated using the MP2, CCSD, and CCSDT²⁷ methods, respectively, and the resulting excitation amplitudes (doubles for EOM-MP2, singles and doubles for EOM-CCSD, and doubles, singles, and triples for EOM-CCSDT) are used to perform a similarity transform of the Hamiltonian which is then used to carry out a CI calculation on the anion. In the EOM-MP2 and EOM-CCSD methods, the CI calculations include all symmetry-allowed one-particle (1p) and two-particle-one-hole (2p1h) configurations, and in the EOM-CCSDT method, the 3p2h configurations are included as well. The EOM-CCSD(T)(a)* method includes the effects of triple excitations in a non-iterative manner and, thus, is computationally considerably less demanding than full EOM-CCSDT calculations.

Both the EOM and ADC(2) methods are free of the “collapse” to the continuum problem for which the HF method and the MP2 and coupled cluster methods employing HF orbitals are susceptible. The OO-MP2 and BCCD(T) methods can also avoid the collapse onto the continuum in calculations on the anion, provided suitable initial guess orbitals are employed. For the BCCD(T) calculations, orbitals from a HF calculation on the cation were used for the starting guess for the anion. Thus, when employing a large basis set with very diffuse functions, the singly occupied orbital of the initial guess wave function for the anion actually corresponds to a Rydberg orbital, which is more localized in the molecular region than is a discretized continuum orbital. For the OO-MP2 calculations, orbitals from a B3LYP^28–30 density functional theory calculation were used as the initial guess orbitals.

The EOM, MP2, BCCD(T), and CCSD(T) calculations were carried out using the CFOUR code,³¹ the ADC(2) and OO-MP2 calculations were carried out using the PSI4 code,³² and the G₀W₀ calculations were performed using Turbomole.³³ The ADC(2) algorithm, as implemented in PSI4, actually calculates excitation energies rather than EBEs. However, this was “tricked” into calculating EBEs by choosing as the ground state a “dianion” with two electrons occupying a very diffuse “continuum” orbital. The electron binding energy of the NVCB anion state is then recovered as an excitation energy (which is of negative sign) in the ADC(2) calculations. The EOM-CCSDT and EOM-CCSD(T)(a)* calculations were also carried out using the electronic excitation implementations of the algorithms together with the “continuum” orbital trick. The basis sets employed in the various calculations are described in Secs. III–V.

The (H₂O)₄ model cluster employed in this work was taken from Ref. 34 and is depicted in Fig. 1. By design, it has a zero net dipole moment, with the water molecules arranged so that each monomer points an OH group toward the cluster interior, which results in a favorable electrostatic interaction for the excess electron near the cluster center. The cluster is treated as two rigid dimers with OO separations of 3.461 Å, the distance (R) between which is varied. This arrangement is not intended to correspond to an experimentally observable structure for (H₂O)₄⁻, but rather is designed as a model for which one can vary the relative importance of correlation and non-correlation contributions to the EBE. In particular, starting at large R and progressing to shorter R, we expect correlation effects to grow in importance relative to the sum of the other contributions to the EBE. Calculations using the various theoretical methods considered are carried out for R values ranging from 2.5 to 9.5 Å.

The majority of the calculations on the (H₂O)₄ cluster were carried out using an aug-cc-pVDZ+7s7p Gaussian-type orbital basis formed by supplementing the standard aug-cc-pVDZ basis set^35,36 with a set of seven diffuse s and seven diffuse p primitive Gaussian functions located at the center of mass of the cluster. The exponents of the supplemental Gaussians were chosen to be in geometric ratios with values ranging from 0.025 to 0.000 025 and 0.022 to 0.000 022 for s and p, respectively. To help elucidate the role of these diffuse functions in describing the NVCB anion, calculations were also carried out with the aug-cc-pVTZ basis set,^35,36 which lacks supplemental diffuse functions. Also, in order to establish the sensitivity of the EBE to the basis set employed, EOM-CCSD(T)(a)* calculations were also carried out for the (H₂O)₄ model with R = 4.2 Å and using an aug-cc-pVTZ+6s6p6d basis set, with the supplemental functions again being centered at the center-of-mass of the cluster.

We first consider the results obtained with the aug-cc-pVDZ+7s7p basis set. Figure 2(a) reports the EBEs from the KT, ΔHF, ΔMP2, ΔCCSD, and ΔCCSD(T) calculations, and Fig. 2(b) reports the EBEs from the ADC(2), ΔOO-MP2, ΔBCCD(T), and various EOM methods using this basis set. At the largest separation considered, R = 9.5 Å, the anion is calculated to be bound by ∼35 and 40 meV in the KT and ΔHF approximations, respectively. Here, and elsewhere in this paper, we use the convention that bound anions correspond to positive EBE values. As R decreases, the KT and ΔHF EBEs decrease, going to zero at R ∼ 4.2 Å. For R $≲$ 4.2 Å, the LUMO from the HF calculation on the neutral molecule has collapsed onto an approximate continuum orbital and the HF wave function of the excess electron system has collapsed onto an approximate continuum state. In contrast to the KT and HF results, the EBEs obtained from the various methods that recover correlation effects (with the exception of MP2 which will be discussed below) increase as R decreases from 9.5 to 4.2 Å. However, while the ADC(2), OO-MP2, EOM, and BCCD(T) methods continue to bind the excess electron for R values $≲$4.2 Å, the EBEs calculated using the CCSD and CCSD(T) methods collapse to zero once the HF method ceases to bind the electron.

The MP2/aug-cc-pVDZ+7s7p calculations considerably underbind the excess electron for the (H₂O)₄ model at all R values considered. Also, the MP2 method, unlike the CCSD and CCSD(T) methods, weakly binds the excess electron for R values between about 3.6 and 4.2 Å. The binding by the MP2 calculations of the excess electron over this distance range is fortuitous and is a consequence of the fact that although the singly occupied orbital of the HF wave function has collapsed onto a discretized continuum solution, it does have some weight in the molecular region. As a result, the MP2 method gives a slightly larger correlation correction for the excess electron system than for the neutral cluster, leading to very weak binding of the excess electron. Coupled cluster doubles calculations (without singles corrections) also weakly bind the excess electron for R values slightly less than 4.2 Å, leading us to conclude that for these geometries, the single excitations in the CCSD and CCSD(T) calculations stabilize the neutral molecule more than the excess electron system.

For large R values, the OO-MP2, BCCD(T), CCSD(T), and various EOM methods give similar EBE values, while the ADC(2) and CCSD methods give weaker binding due to their less complete treatment of relaxation in response to correlation effects. In progressing from large R to R = 4.2 Å, there is a growing discrepancy between both the EOM-MP2 and EOM-CCSD values of the EBEs and the EOM-CCSDT and EOM-CCSD(T)(a)* values of the EBEs, with the latter two methods giving larger values of the EBE (by about 10% at R = 4.5 Å). As R approaches 4.2 Å, the EBEs calculated using the OO-MP2, BCCD(T), and CCSD(T) methods all remain close to the EOM-CCSDT results, which we take as the benchmark values. By contrast, as noted above, the CCSD and CCSD(T) calculations collapse onto the neutral plus an approximate continuum orbital for R values less than 4.2 Å.

For R $≲$ 4.2 Å, the ADC(2) method continues to give weaker binding than do the OO-MP2, EOM, and BCCD(T) methods. In addition, at these short R values, the EOM-MP2 and EOM-CCSD methods bind the excess electron more strongly than does the ADC(2) method but less strongly than do the EOM-CCSDT, EOM-CCSD(T)(a)*, BCCD(T), and OO-MP2 methods, all of which give EBEs in close agreement with one another. At R = 3.5 Å, the EBEs from the EOM-MP2 and EOM-CCSD methods are about 15% smaller than the EOM-CCSD(T)(a)* and EOM-CCSDT values. These results lead us to conclude that for the (H₂O)₄ model at geometries at which the HF method fails to bind the excess electron, the OO-MP2 and BCCD(T) methods are more successful at describing relaxation in response to correlation effects than are the EOM-MP2 and EOM-CCSD methods.

It is instructive to compare, at different R values, the LUMO from the HF calculations on the neutral cluster and the singly occupied natural orbital (SONO) from the EOM-CCSD calculations on the anion, both calculated using the aug-cc-pVDZ+7s7p basis set. Figure 3 plots these two orbitals along the axis perpendicular to the cluster center for R = 2.5, 3.5, 4.5, and 8.5 Å. For each of these R values, the charge distribution associated with the SONO from the EOM-CCSD calculations is largely localized in a region within 10 Å of the center of the cluster, with the extent of localization increasing slightly as R decreases. For R = 8.5 Å, the LUMO from the HF calculations is qualitatively similar to the SONO from the EOM-CCSD calculations, but as R decreases, the LUMO from the HF calculations becomes progressively more diffuse than the EOM-CCSD SONO. For R = 4.5 Å, the LUMO, while still bound, is far more radially extended than the EOM-CCSD SONO, and for R = 2.5 and 3.5 Å, the LUMO has little weight in the vicinity of the cluster, consistent with its collapse onto an approximate continuum function. In fact, at R = 3.5 Å, it is the sixth lowest unoccupied orbital from the HF calculations of the neutral cluster that most closely resembles the SONO from the EOM-CCSD calculations.

The results presented above were obtained using the aug-cc-pVDZ+7s7p basis set. We now examine the sensitivity of the EBE of the (H₂O)₄ model to further expansion of the basis set. This analysis is done at R = 4.2 Å using the EOM-CCSD(T)(a)* method. The resulting EBE obtained using the aug-cc-pVTZ+6s6p6d basis set is 208 meV, which is about 20% greater than that obtained using the aug-cc-pVDZ+7s7p basis set. We also carried out EOM-CCSD(T)(a)* calculations using the aug-cc-pVDZ+6s6p6d and aug-cc-pVTZ+7s7p basis sets, and, based on the results of these calculations, we conclude that roughly half of the 20% increase in EBE noted above comes from the improvement of the description of the charge distribution of the excess electron upon inclusion of the diffuse d functions and about half is due to the improved description of the response of the valence electrons of the water molecules in going from the aug-cc-pVDZ+7s7p to the aug-cc-pVTZ+6s6p6d basis set. This is consistent with the dominant correlation effects being long-range dispersion-like interactions between the excess electron and the valence electrons of the water molecules. The value of the EBE obtained from the EOM-CCSD(T)(a)*/aug-cc-pVTZ+6s6p6d calculations is expected to be close to the exact value, i.e., that which would be obtained from a full CI calculation in the complete-basis-set limit. In particular, we note that in an earlier publication, where we considered the (H₂O)₄ model at a single geometry (R = 5.3 Å), we found that the EBE increased by only 5 meV in going from the aug-cc-pVTZ+6s6p6d to the aug-cc-pVQZ+6s6p6d.³⁷ 

We now turn our attention to the calculations using the aug-cc-pVTZ basis set, with the main results being reported in Fig. 4. For this basis set, the anion is calculated to be unbound by 200 to 550 meV in the KT and Hartree-Fock approximations over the entire range of R values considered. The failure of the KT and ΔHF calculations using the aug-cc-pVTZ basis set to bind the excess electron to the (H₂O)₄ cluster at large R values is a consequence of the artificial confinement (due to the lack of very diffuse functions in the basis set) of the orbital occupied by the excess electron, which, in turn, causes the kinetic energy contribution to be too high. For R values ranging from 2.5 to about 5.5 Å, all of the methods that include correlation corrections bind the excess electron when employing the aug-cc-pVTZ basis set. However, this binding is fortuitous as it is a consequence of the unrealistically large correlation contributions to the EBE that result from the artificial confinement of the singly occupied orbital in the vicinity of the water monomers. Thus, it is not surprising that the R dependence of the binding energy is not properly reproduced in the calculations with the aug-cc-pVTZ basis set. In fact, none of the methods considered give a bound anion for R > 5.5 Å, in contrast to the calculations carried out with the aug-cc-pVDZ+7s7p basis set.

The most important take-away message from the calculations on the NVCB anion of the (H₂O)₄ cluster model is that in order to properly describe this system, it is essential to employ a theoretical method that adequately accounts for the coupling of orbital relaxation and dispersion-like correlation effects and to use a basis set that includes s and p basis functions considerably more diffuse than those in the standard aug-cc-pVXZ basis sets. The fact that the EBEs from the OO-MP2 calculations are nearly identical to those from the EOM-CCSDT calculations indicates that for this cluster, even a second-order treatment of the dispersion-like correlation effects suffices for describing the NVCB anion, provided that orbitals optimized in the presence of these correlation effects are employed. In the case of the ADC(2) method, this relaxation is accomplished through the off-diagonal terms in the self-energy. When the off-diagonal coupling is neglected, the EBE from the ADC(2) calculations using the aug-cc-pVDZ+7s7p basis set goes to zero with decreasing R even more rapidly than the MP2 value of the EBE, and the diagonal ADC(2) method fails to bind the excess electron for R values less than 4.2 Å. The non-self-consistent G₀W₀ method, which also neglects off-diagonal self-energy terms, is found to give results similar to those obtained using the diagonal ADC(2) method. The EBEs obtained with the ADC(2), diagonal-ADC(2), and G₀W₀ methods are compared with the corresponding EOM-CCSDT results in Fig. 5.

The success of the BCCD(T) method, in contrast to the CCSD(T) method for R values for which the HF method does not bind the excess electron, is particularly noteworthy. Because a Slater determinant of Brueckner orbitals is that with the maximum overlap with the true wave function, the Brueckner procedure results in a singly occupied orbital for the anion with a charge distribution close to that of the SONO from the EOM calculations. Figure 6 compares the EBEs obtained using the difference of the energy of the anion and the neutral cluster both calculated using Slater determinants of Brueckner orbitals (and here referred to as ΔSD-BO) with the corresponding ΔHF results. For large R, the ΔSD-BO calculations bind the excess electron albeit less strongly than the ΔHF method, but for R values $≲$6 Å, they give a negative EBE. At R = 2.5 Å, the shortest distance considered, the anion is predicted to be unbound by 220 meV in the ΔSD-BO method. These results are consistent with the singly occupied orbital from the BCCD calculations on the anion being slightly more localized than the HF LUMO at large R and with the former becoming more localized and the latter more delocalized as R decreases. We also see from Fig. 6 that the difference between the ΔBCCD(T) and ΔSD-BO values of the binding energy grows in magnitude from 69 meV at R = 9.5 Å to 350 meV at R = 2.5 Å, consistent with the growing importance of the correlation contribution to the EBE as R decreases.

We also applied the self-consistent polarization model of Ref. 34 to calculate the electron binding energy of the (H₂O)₄ cluster at R = 4.2 Å. With the model potential approach, the net EBE can be decomposed into kinetic energy, electrostatics, repulsion, and polarization contributions, where the kinetic energy contribution results from localization of the excess electron, the electrostatics contribution results from the interaction of the electron with the three point charges employed on each water monomer, the repulsion contribution results from the pseudopotential used to prevent overattraction to the positively charged H atoms and to effectively account for charge penetration and exchange effects, and the polarization contribution describes the correlation effects between the excess electron and the valence electrons. The results of the energy decomposition are summarized in Table I. At R = 4.2 Å, the model potential gives an electrostatic contribution to the binding energy of 709 meV. However, the sum of the kinetic energy and repulsive contributions is even greater in magnitude, with the sum of these three contributions being −104 meV, roughly comparable to the ΔSD-BO value of the EBE at R = 4.2 Å. The polarization contribution to the EBE is 312 meV, giving a net binding of 208 meV, which is in close agreement with the EOM-CCSDT and ΔBCCD(T) results.

It is important to note that both the BCCD(T) and polarization potential approaches use orbitals that have been optimized in the presence of correlation effects. It is more traditional to define correlation corrections with reference to the HF energy, in which orbitals are calculated using orbitals determined in the absence of correlation effects. However, for NVCB anions, the HF wave function collapses onto a discretized continuum solution, and we believe it is more instructive to perform the energy decomposition using orbitals optimized in the presence of correlation effects.

As our second test case, we consider the bending potential of the ground state anion of CO₂. It is known from electron scattering measurements that vertical electron capture into the valence π* orbital of CO₂ results in a temporary anion located energetically about 3.7 eV above the neutral molecule.³⁸ It is also known that CO₂⁻ has a valence-like ²A₁ anion state with a minimum energy structure with an OCO angle of about 138° at which it lies energetically about 1.5 eV below the neutral molecule at the same geometry and about 0.5 eV above the neutral molecule at its minimum energy linear structure.^39–42 The temporary anion of CO₂⁻ formed by vertical electron capture splits into ²B₁ and ²A₁ components as the molecule is bent, and, in the standard textbook explanation, the ²A₁ component of the π* temporary anion state correlates with the bound valence-type ²A₁ anion state with OCO angle ∼138°.^39–42 In this picture, the bending potential of the lowest energy anion state of CO₂⁻ would be expected to be similar to that of the isoelectronic ²A₁ state of NO₂. However, theoretical studies of Vanrose and co-workers,⁴³ Sommerfeld,³⁹ and Sommerfeld et al.⁴² have shown that this textbook picture of CO₂⁻ is not correct. In particular, as shown by Sommerfeld, EOM-CCSD calculations employing a large basis set with several sets of diffuse functions predict that as the OCO angle increases from ∼138°, the potential energy curve of the anion first rises slightly in energy, having a maximum near 150°, and then descends in energy, approaching the potential energy curve of the neutral molecule, ceasing to be bound near an OCO angle of 154°. Moreover, as the potential energy curve of the anion approaches that of the neutral, the wave function of the anion becomes much more diffuse than that associated with a typical valence state.

In Sommerfeld’s study of CO₂⁻, the CCSD(T), EOM-CCSD, ADC(2), and ADC(3)⁴⁴ methods were considered. In the present investigation, we extend this work to examine also the performance of the OO-MP2, EOM-CCSD(T)(a)*, EOM-CCSDT, and BCCD(T) methods for this challenging system. We employ an ANOTZ+3s3p basis set formed by augmenting the ANOTZ basis set⁴⁵ on each atom with three diffuse s and three diffuse p functions, the exponents of which were taken from Ref. 39. The CO distance is allowed to vary with the OCO angle, as was done in Ref. 39. In particular, CO₂⁻ in its bent structure and CO₂ in its linear structure were optimized using the CCSD method, and a linear interpolation of these two geometries was used to determine the CO distances at other values of the OCO angle.

In our calculations on CO₂⁻, we discovered that for OCO angles between 147° and 154°, it was possible to obtain two different HF solutions: one corresponds to the valence state and the other to the neutral molecule plus an electron in a discretized continuum orbital. One HF solution was obtained by using orbitals from calculations at smaller angles as the initial guess for calculations at larger angles, and the other by using orbitals from calculations at larger OCO angles as the initial guess for the calculations at smaller angles. Figure 7(a) depicts for the 147°–154° range of angles the anion potentials obtained using the valence-like HF orbitals, while Fig. 7(b) shows the corresponding potentials obtained using orbitals from the HF solution that collapsed onto the discretized continuum. The potential for the neutral molecule calculated using the various theoretical methods is also reported. From Fig. 7(a) it is seen that the HF solution that corresponds to the valence state rises in energy with increasing bending angle up to 154° (beyond that it collapses onto the neutral plus discretized continuum orbital). By contrast, the HF potential corresponding to the discretized continuum solution is essentially parallel to the potential for the neutral, lying about 0.04 eV above it. (The reason it lies slightly above the neutral in energy rather than on top of it is a consequence of the basis set employing only three sets of supplemental diffuse functions on each atom.)

The singly occupied orbital for the two HF solutions of the anion at OCO angles of 147° and 154° is shown in Fig. 8 as is the SONO from the EOM-MP2 calculations. The valence-type HOMO of the anion becomes more extended as the OCO angle increases from 147° to 154°, but the increase is relatively modest. Interestingly, the SONO from the EOM calculations is somewhat more extended than the HOMO from the HF calculations on the valence-type state, and it grows in extent with the increase in the OCO angle from 147° to 154° more than the HF HOMO. Over the 147°–154° range, the bending potentials of CO₂⁻ obtained from the two sets of CCSD calculations are fundamentally different. In particular, over this range of angles, the CCSD calculations starting from the valence-like HF orbitals do not collapse onto the “continuum” solution even when the continuum solution lies lower in energy and the CCSD calculations starting from the “collapsed” HF wave function do not give the valence-like solution even when the valence-like solution lies lower in energy. In other words, there is a range of angles for which the single excitations in the CCSD wave function are not sufficient for overcoming the limitations of the HF starting point.

Figure 9 reports the calculated bending potentials of CO₂ and CO₂⁻ for the entire range of angles considered. The bending potential for the neutral CO₂ molecule was generated at the CCSD(T) level of theory, and for the anion, results are reported for the MP2, CCSD(T), EOM-CCSD, ADC(2), OO-MP2, and BCCD(T) methods, with the anion potentials being generated by subtracting the EBE calculated with each method from the CCSD(T) value of the energy of the neutral at the same geometry. (Results obtained using higher-order EOM methods will be presented below.) In the case of the MP2 and CCSD(T) calculations, the results shown are based on orbitals from the valence-like HF solution. The EOM-CCSD and BCCD(T) methods give similar bending potentials for the anion. In both cases, the anion bending potential is very flat over the 135°–148° range, with a small barrier near 150°, and then descends slightly in energy, intercepting the neutral potential near an OCO angle of 156°. The bending potentials of CO₂⁻ from the CCSD and CCSD(T) calculations based on the HF valence-type wave function for the anion agree closely with those from the EOM-CCSD and BCCD(T) calculations up to an angle of about 151° but then rise in energy rather than descending as found in the EOM-CCSD and BCCD(T) calculations as the OCO angle is further increased, collapsing onto the neutral plus continuum for angles greater than ∼154°. The MP2 potential, on the other hand, lies about 0.3 eV above the EOM-CCSD potential for the whole range of angles for which the former is bound.

We note that for OCO angles greater than about 154°, the CCSD(T) potential reported by Sommerfeld³⁹ descends in energy rather than rising as found in the present study. This difference in behavior is due to the different basis sets used in the two studies: whereas Sommerfeld started with a valence-type basis set of Dunning⁴⁶ contracted to 5s3p1d, the more flexible augmented ANOTZ basis set contracted to 5s4p3d2f was used in this work.

The OO-MP2 potential starts out about 0.3 eV too high in energy at an OCO angle of 135° and “rapidly” descends in energy crossing the neutral potential at essentially the same angle as that found with the EOM-CCSD and BCCD(T) methods. The ADC(2) bending potential lies too low in energy at all angles and crosses the neutral potential near an OCO angle of 159°. Although the MP2, OO-MP2, and ADC(2) methods all fail to adequately describe correlation effects important for the binding of the excess electron to CO₂ at angles for which the anion is bound, the latter two methods do predict a rapid energy descent as the anion potential approaches that of the neutral. This is a consequence of their including orbital relaxation in response to correlation effects. The theoretical characterization of CO₂⁻ is complicated by the fact that the anion has significant multiconfigurational character at small OCO angles where the anion is valence-like and by the rapid growth of non-valence character as the anion potential approaches that of the neutral. The former is primarily responsible for the inadequacy of the MP2, OO-MP2, and ADC(2) methods for describing the bending potential of the anion at small OCO angles.

Figure 10 reports as a function of the bending angle the EBE of CO₂ calculated using the EOM-MP2, EOM-CCSD(T)(a)*, and EOM-CCSDT methods as well the EOM-CCSD method, already considered above. It is seen from this figure that for OCO angles between 135° and 154°, the EOM-MP2, EOM-CCSD and EOM-CCSDT methods give essentially the same EBE values, while the EOM-CCSD(T)(a)* method slightly overbinds the excess electron for small values of the OCO angle.

The results presented above indicate that there is a fundamental difference in the nature of the (H₂O)₄ and CO₂ anions near the crossing points of the anion and neutral potential energy surfaces. In the case of the (H₂O)₄ cluster, the singly occupied orbital of the HF wave function of the anion rapidly expands in radial extent as the crossing point is approached from the large R side. By contrast, for CO₂, two different HF solutions can be found for OCO angles ranging from 147° to 154°. One of these is valence-like and the other corresponds to the neutral molecule plus an excess electron in an approximate continuum orbital. The HF bending potential of the valence-like anion state crosses the potential of the neutral near an OCO angle of 149° and continues to rise up to an angle of about 154°, after which it collapses to the lower-energy “continuum solution.” Most importantly, the wave function of the valence-like HF solution for CO₂⁻ undergoes only a moderate growth in its extent as the crossing point is approached (from the bound state side). Moreover, while the SONO from the EOM calculations on (H₂O)₄⁻ does not change significantly in its extent as R passes through the HF crossing point, the SONO from the EOM calculations on CO₂⁻ (shown in Fig. 8) evolves from being relatively compact (valence-like) at an OCO angle of 135°, to being much more extended for values of the OCO angle greater than 149°. For the (H₂O)₄ model system, the anion evolves from being largely electrostatically bound at large separation between the two (H₂O)₂ dimers to being correlation bound at short separation between the dimers, with correlation effects causing the anion wave function to contract at R values where the anion is only weakly bound in the HF approximation. By contrast, the CO₂⁻ ion evolves from being valence-like for values of the OCO near 135°, but with correlation effects being responsible for the more diffuse character of the anion function near the crossing point of the anion and neutral bending potentials.

The change in the charge distribution of the CO₂⁻ with bending angle appears to result from an avoided crossing between the valence-like state and a non-valence anion state. Three candidates for the non-valence state (or states) can be envisioned: (1) a NVCB anion, which, over a range of OCO angles near 149° and in the absence of mixing with the valence state, would lie slightly below the neutral in energy; (2) a virtual state,^43,47 and (3) the manifold of continuum levels. We note that an avoided crossing between valence states and a NVCB anion state has been demonstrated in a theoretical study of C₆F₆⁻.⁴⁸ NVCB anion states and virtual states are closely related in that a virtual state turns into a NVCB anion state if the potential experienced by the electron is made slightly more attractive. With respect to the possibility of a NVCB anion of CO₂, we note that for an OCO angle of 149°, the calculated (HF level) dipole moment and average dipole polarizability of the neutral molecule are 1.05 D and 14.3 a.u.,³ respectively. Based on our experience with a wide range of NVCB anions, we would not expect a molecule with a polarizability this small, even with a dipole moment of 1.05 D, to support a NVCB anion in which long-range dispersion interactions are primarily responsible for the binding of the excess electron. Rather, we favor the third possibility, namely, that near OCO angles of 149°, the bound anion state from strong mixing of the valence state and the continuum. We suggest, therefore, that it may be useful to broaden the definition of NVCB anions to include such species.

Although our focus in this section has been on the performance of various theoretical methods for describing the anion of CO₂ near the crossing point of the anion and neutral bending potentials, it is interesting to reflect on the implications of the current results for the interpretation of experiment. In this context, we note that CO₂⁻ has been observed mass spectroscopically upon collision of high energy alkali atoms with CO₂.⁴⁹ This study gave, for CO₂, an adiabatic electron affinity of −0.60 ± 0.2 eV, which is consistent with the results of the present EOM and coupled cluster calculations. However, there is the issue that the calculated bending potential has only a very small (∼0.036 eV or 290 cm⁻¹) barrier before it bends downward near an OCO angle of ∼150°, and with such a small barrier, the bent valence-type anion would be expected to autoionize in a short time compared to the experiment if the zero-point level were above or near the top of the barrier. We have calculated, using the Gaussian 09 program,⁵⁰ anharmonic vibrational frequencies of bent CO₂⁻ using the second-order vibrational perturbation theory (VPT2)⁵¹ approach with the MP2 electronic structure method in conjunction with the ANOTZ+3s3p basis set. These calculations give an anharmonic frequency of 600 cm⁻¹ for the bending vibration, which places the zero-point vibrational level slightly above the barrier. However, we note that the bending potential of CO₂⁻is very flat, and it is likely that the VPT2 calculations considerably overestimate the value of the frequency of the bending mode, and it seems highly likely that the zero-point level is below the barrier.

In the two examples considered above, there are no valence anion states energetically below the NVCB anion state. However, as we demonstrated for C₆₀, molecules with bound, valence anion states can also possess NVCB anion states.⁷ Accurately describing the NVCB anion states of such systems can be especially challenging since high-order correlation effects can be important. To illustrate this, we consider TCNE which has a valence anion state bound by 3.2 eV⁵² as well as several low-lying valence-like temporary anion states.⁵³ Calculations to search for a possible ²A_g NVCB anion of TCNE were performed using the aug-cc-pVDZ and aug-cc-pVTZ basis sets augmented with 3s, 3s3p, and 3s3p3d sets of supplemental diffuse functions centered at the midpoint of the central CC bond. The calculations indicate that TCNE does have a NVCB anion, with the EBEs obtained using the various theoretical methods summarized in Table II. The SONO from the EOM-MP2/aug-cc-pvDZ+7s7p calculations on the ²A_g NVCB anion state of TCNE is depicted in Fig. 11.

For each basis set considered, the EOM-CCSD(T)(a)* and EOM-CCSD values of the EBE are in close agreement, and in the assessment of the performance of the various approaches, the results of the EOM-CCSD(T)(a)* calculations will be used as the benchmark reference. The EOM-CCSD(T)(a)* EBEs obtained with the aug-cc-pvDZ+3s3p and aug-cc-pvTZ+3s3p basis sets are 63 and 101 meV, respectively. The BCCD(T) calculations give EBE values for the NVCB anion of TCNE very close to the corresponding EOM-CCSD(T)(a)* results when the same basis set is employed for the two methods. By contrast, the EOM-MP2 OO-MP2 and ADC(2) methods all give EBE values much larger than the EOM-CCSD(T)(a)* values, with the EOM-MP2 method overestimating the EBE by 20%-30% and the OO-MP2 and ADC(2) methods by factors of 2.7 and 4.7, respectively. These results indicate that the inclusion of high-order correlation effects is essential for obtaining an accurate estimate of the EBE of the NVCB anion of TCNE. We have calculated the dipole polarizability of TCNE using the Perdew-Burke-Ernzerhof (PBE) density functional method⁵⁴ in conjunction with the aug-cc-pVTZ basis set. The resulting average polarizabilities in the uncoupled and coupled Kohn-Sham approximations are 229 and 99 a.u.,³ respectively. The large reduction of the polarizability due to the coupling of excitations is consistent with screening effects playing an important role in the interaction of the excess electron with TCNE, which, in turn, could explain the large overestimation of the EBE as calculated with the OO-MP2 and ADC(2) methods.

One surprising result of our calculations on the NVCB anion of TCNE is that the mid-bond supplemental diffuse p functions play a negligible role in the electron binding. This indicates that the excitations of the excess electron contributing to the dispersion interaction between it and the valence electrons of the molecule are dominated by excitations into virtual orbitals with a large component of the diffuse aug-functions in the basis sets.

As TCNE supports both a valence-bound (²B_3g) anion and a NVCB (²A_g) anion, it provides an opportunity for a direct comparison of the properties of these two types of anion states. We note that EOM-CCSD/aug-cc-pVDZ+3s3p3d calculations predict the valence-type anion to be bound by 2.91 eV, as compared to the 79 meV EBE for the NVCB anion. For both states, the natural orbital occupied by the excess electron has an occupation number close to 0.985. (All other natural orbitals have occupation numbers greater than 1.9 or less than 0.1). Figure 12 reports the densities of these natural orbitals as a function of r_neartest which is defined as

where $r→$ is the position of the excess electron and $R→α$ is the position of nucleus α. The corresponding expectation value 〈r_nearest〉 provides a measure of the extent of the orbital.¹² 

The densities of the EOM-CCSD natural orbitals associated with the excess electron of the ground state valence-bound anion and of the NVCB anion of TCNE are shown in Fig. 10, and, as expected, the difference is striking. In the valence-bound state, the excess electron is localized close to the molecular framework, with the maximum of the density distribution function occurring at 0.70 Å and 〈r_nearest〉 being only slightly larger at 0.83 Å. By contrast, the excess electron of the NVCB state extends much further out and shows a broad distribution with a long-range tail. The maximum of its distribution occurs at 2.4 Å, while 〈r_nearest〉 is 4.4 Å, emphasizing the smeared-out character of a NVCB electron.

In contrast to valence and dipole-bound anions, non-valence correlation bound anions have been much less studied. Theoretical studies from our group and others have predicted the existence of NVCB anions for several clusters and molecules.^{7–10,34,48,55} The picture that has emerged from the theoretical studies is that any system for which the electrostatics are not sufficiently attractive to give electron binding but which has a sufficiently high polarizability will possess one or more NVCB anion states. Although there is a paucity of experimental data on NVCB anions of gas-phase molecules and clusters, we note that NVCB anions are the finite system analogs of image potential states of metallic and graphitic surfaces, which have been characterized experimentally and theoretically.^56–58 In both the finite and extended systems, the binding of the excess electron results from long-range dispersion interactions between the excess electron and the other electrons of the system.

Theoretical characterization of NVCB anions requires the use of methods that allow for orbital relaxation effects in response to the dispersion-like correlation effects. We have demonstrated that the equation-of-motion coupled cluster and algebraic diagrammatic Green’s function direct methods and also the orbital-optimized MP2 and Brueckner coupled cluster delta methods are all able to treat such anion states. However, depending on the nature of the neutral molecule or cluster, lower-order methods such as OO-MP2, ADC(2), and even EOM-MP2 may have significant errors in the EBEs. The need to include higher-order correlation effects is especially important when considering NVCB anions of systems with small HOMO/LUMO gaps and/or for which screening effects are important. This is evident from the results presented here for the NVCB anion of TCNE and from the previously reported results for the NVCB anion of C₆₀. While methods such as BCCD(T) and EOM-CCSD, particularly when used with the large, flexible basis sets required for adequate characterization of NVCB anions, are computationally prohibitive for large molecules and clusters, it is possible to design one-electron models that can accurately describe these species at low computational cost.^34,55,59

K.D.J. acknowledges support from NSF Grant No. CHE1362334, and T.S. acknowledges support from NSF Grant No. CHE1565495. The calculations were carried out on computers in the University of Pittsburgh’s Center for Research Computing. We thank Professor D. Sherrill for helpful discussions concerning orbital optimized methods and Professor John Stanton for helpful discussions about EOM-CCSD(T)(a)* calculations.