Communication: Avoiding unbound anions in density functional calculations Free

Anions and radicals are important for many applications including environmental chemistry,^1–3 semiconductors,^4,5 fullerene chemistry,^6–10 charge transfer,¹¹ and solar cells.^12,13 Recently, electron affinities of biological species become of great interest, especially in studies of low-energy electron DNA damage.^14–20 Low-energy electrons cause single-strand breaking, double-strand breaking, and supercoil loss in DNA even below the DNA ionization potential. The electron affinity of DNA bases and base-pairs is important in determining damage mechanism.

Density functional theory (DFT) has become a standard method for electronic structure calculations in chemistry, and the standard functionals can be applied with standard basis sets to calculate electron affinities. The results are excellent, with mean absolute errors (MAE) below about 0.2 eV.²¹ However, there is a theoretical fly in the computational ointment: Inspection of the orbital energies shows that the HOMO of the anion is usually positive. This implies that, in principle, the calculation is unconverged.²² If a sufficiently large basis set had been used, a fraction of the additional electron would ionize²³ and the HOMO drops to zero. This is due to the self-interaction error that all the standard density functional approximations suffer from. This error is especially large for anions, because of their additional electron. Self-interaction error is caused by using approximate exchange functionals, i.e., inexact exchange, in DFT calculations. As a consequence, it produces an exchange-correlation potential that incorrectly decays exponentially in the asymptotic region instead of decaying as −1/r. For atomic anions, a large positive barrier appears in the Kohn-Sham potential (see Fig. 1 of Ref. 24) resulting in positive HOMO resonances. These metastable states are occupied and artificially bound by moderate basis sets (MBS), and so produce a positive HOMO. Because the positive barrier is often very wide, the total energy appears converged unless extreme basis functions are used.

There have been strong discussions about this issue.^22,24–27 Users find reasonable results with MBS for most cases and ignore the positive HOMO. Purists consider all such calculations as unconverged and so their results are suspect.²² Pragmatists will report results with the standard methods, but attach a caveat emptor footnote.²¹ The paradox has recently been addressed in several papers,^24,27 which explain how accurate results can come from such unconverged calculations but also suggest an alternative procedure that avoids the dilemma: Evaluate the density functional total energies on Hartree-Fock (HF) densities for both the neutral and anion, and calculate the electron affinity from the total energy difference. We refer to this method as HF-DFT.

Because HF treats exchange exactly and ignores correlation, all self-interaction is cancelled,^28,29 and its HOMOs are bound, even for anions. Electron affinity calculations for atoms and their anions show excellent results with either method,^24,27 with MAEs about 0.1 eV, about half of that for ionization potentials.

In the present work, we test the new procedure for the adiabatic electronic affinity of the molecules in the G2-1 data set³⁰ excluding CN.^31,32 All calculations are performed with TURBOMOLE 6.2.³³ The functionals used in our DFT calculations are PBE,³⁴ B3LYP,^35–37 and PBE0.³⁸ Here, we consider calculations with hybrid functionals such as B3LYP as being approximate DFT. This is because one can calculate these within a Kohn-Sham scheme, treating them as orbital-dependent, using for example the local Hartree-Fock technique.^39,40 In such a case, their Kohn-Sham potentials capture a fraction of the exact long-range decay (the fraction of exact-exchange mixing) and so are also incorrect, and suffer from the same conceptual problems. But their total energies are indistinguishable from those of a HF-style calculation, which is how they are commonly calculated and the prescription we follow here. Since it is well-known that HF may give poor densities when using small basis sets,⁴¹ to overcome this problem, we use Dunning's augmented correlation-consistent pVTZ (AVTZ) basis set.^42,43 Structures of neutral molecules and molecular ions are optimized with self-consistent Kohn-Sham-DFT and non-scaled zero-point vibrational correction with the same functional is added.³² For HF-DFT calculations, unrestricted HF calculations are performed on both neutral and anion DFT optimized structures. Based on these HF orbitals, the same functional used in the structure optimization was selected for energy evaluation. We distinguish such calculations by HF-XC, where XC indicates the exchange-correlation approximation used.

In Fig. 1, we plot calculated versus experimental electron affinities showing just how good the overall agreement is.³² In Table I, we report results for both methods for all molecules in the G2-1 set, but with averages excluding CN. Averages are reported in MAE and mean of errors (ME). We find excellent results once again for almost all molecules, with either conventional DFT with MBS or the HF-DFT method with MAEs again about 0.1 eV. We find HF-PBE yields the best results overall but differences are slight and not significant.

In Fig. 2, we plot errors versus electron affinities including also the results from atoms of Refs 24 and 27. Notice that switching the density from self-consistent to HF always either reduces the electron affinity (sometimes increasing the error) or increases it by no more than 0.02 eV. In the case of PBE, the self-consistent electron affinities are mostly too large, a systematic error inherited from the local density approximation (LDA, sometimes called VWN(Ref. 44)). This is reflected in the fact that the ME ∼ MAE on the scale of the MAE in Table I. The reduction in electron affinities on using the HF density, which leads to subsequent reduction in MAE and large reduction in ME, suggests that this is largely a self-interaction error in the density not the energy functional.

To analyze this effect in more detail, we note that the LDA generally underestimates the magnitude of the exchange-correlation energy. Because the electronic energy of any species (with fixed nuclear positions) is convex in N, the number of electrons, the ionization energy increases as N decreases. The greater the net charge, Z − N, the more compact the density and the less such as a uniform gas. Thus, the error of LDA typically increases with net charge leading to an overestimate of the ionization energy for given N. This applies to both neutrals and anions. Since PBE is designed to reproduce the systematic behavior of LDA while improving its accuracy, it likewise has a smaller, systematic overestimate of both ionization potentials and electron affinities.

Generally, HF densities are not as accurate as those of approximate density functional calculations,⁴¹ because the latter include correlation albeit approximately. However, for anions this is formally untrue, since a truly self-consistent approximate DFT calculation in the infinite basis limit will lose some fraction of an electron to the continuum, so that the bound density will not even integrate to the correct number of electrons. In the MBS approach, this does not occur, one expects a remnant of the self-interaction error to make the anionic density more diffuse nonetheless. This self-interaction error is about as large as one can get in a single molecule. Thus, for this particular case, the error due to lack of correlation in the HF density must be weighed against the self-interaction error in the MBS approximate DFT density. Figure 6 of Ref. 24 shows that, for Li⁻, the errors in the densities are comparable with either method especially in the tail region. Our calculations here clearly show that the calculation with HF densities generally reduces the PBE electron affinity and so improves the calculated electron affinity. We suggest that this is due to the improved accuracy of the HF anionic densities, the self-interaction error being greater than that due to lack of correlation. (More precisely, it is the change in density between the anion and the neutral that is more accurate in HF.) A similar trend was found for transition-state barriers.⁴⁵ 

On the other hand, B3LYP is a hybrid functional with empirical parameters. It already cancels some self-interaction error and has smaller MAE. But the ME is much smaller than the MAE showing that its errors have random signs, i.e., much less systematic than those of PBE. Inserting the HF density does not improve MAE, and even worsens ME. To check our interpretation of the effect of the hybrid we applied another functional, the non-empirical hybrid PBE0(Ref. 38) and found results with the same trends, but higher MAEs (1.4 eV with or without HF densities).³² 

We also show the HOMO energies for the anions, in both HF and self-consistent DFT calculations. All species except Cl₂ have positive HOMO in the DFT calculations indicating their unbound nature. A sufficiently large basis²⁶ would reduce this value, but this effect may not be noticeable with any standard basis set. On the other hand, although the HF HOMOs are negative, they are not an accurate guide to the true electron affinities. Even anion HOMO levels calculated from density functional calculations with correct asymptotic decay suffers from the same problem.⁴⁶ Using Koopman's theorem, one may estimate electron affinity from the HOMO energy of anion but both relaxation and correlation effects are so large that Koopman's theorem is unhelpful here. Self-interaction corrections have been shown to improve orbital energies but not energy differences for the functionals.^47,48 Since the HF density is calculated upon geometries optimized from DFT, the resulting HF-DFT energy will not typically be a minimum in the HF-DFT potential energy surface. This shows further improvement may be made in HF-DFT by development of potential energy surface scan and optimization techniques.

A useful tool for understanding these effects is the electron affinity density

where n₀(r) is the charge density of the neutral and n₋₁(r) is that of the anion. Figure 6 of Ref. 27 shows this for the Cl atom and anion, and how the HF electron affinity density is more compressed than that of using MBS and standard functionals. We plot the cross sections of electron affinity densities of NH along the molecular axis for different methods. In Fig. 3, the electron affinity densities of self-consistent Kohn-Sham DFT and HF calculations are plotted for NH. NH, which has the largest deviation between the electron affinity error of PBE and HF-PBE, the electron affinity density of the two is nearly identical. Nonetheless, the self-consistent density is more diffuse than the HF density due to the electron leakage in the anion as shown in Fig. 4. To confirm that this is not an artificial effect from the geometry difference in neutral and anion, we present the local vertical electron detachment densities, which are electron affinity densities using the same geometry for the anion and neutral.³² 

In summary, we have shown that the methods proposed in Refs. 24 and 27 work just as well for the small molecules of the G2-1 data set, so long as the HF densities are not too far from the exact densities. The results are equally good with HF-PBE as with B3LYP, perhaps slightly better. However, we found no case where the limited basis set approach fails. Presumably, the unphysical barrier holding the additional electron in for atomic anions²² is sufficiently large for our molecules that standard basis sets show no sign (other than a positive HOMO) that the state being calculated is a resonance rather than an eigenstate. We recommend that the HF-DFT method be applied more broadly for electron affinity calculations, especially for cases where DFT with MBS is believed to be inaccurate.

We thank Professor Fillip Furche and Dr. Donghyung Lee for fruitful discussions. This work was supported by the global research network (Grant No. NRF-2010-220-C00017) and E.S. acknowledges the support under the national research foundation (Grant Nos. NRF-2010-0016487 and NRF-2010-0017172). K.B. acknowledges support under National Science Foundation (NSF) CHE-0809859 and M-C thanks the fellowship of the BK 21 program from MOEHRD.