Charge-transfer in time-dependent density-functional theory via spin-symmetry-breaking
Long-range charge-transfer excitations pose a major challenge for time-dependent density functional approximations. We show that spin-symmetry-breaking offers a simple solution for molecules composed of open-shell fragments, yielding accurate excitations at large separations when the acceptor effectively contains one active electron. Unrestricted exact-exchange and self-interaction-corrected functionals are performed on one-dimensional models and the real LiH molecule within the pseudopotential approximation to demonstrate our results.
Although time-dependent density-functional theory (TDDFT) has had resounding success in predicting accurate excitation spectra in a wide variety of systems RG84 (); TDDFTbook (), difficulties still plague its application to certain areas. The problem of charge-transfer (CT) excitations has drawn especially significant attention in recent years DWH03 (); T03 (); GB04bNGB06 (); M05cMT06 (); TTYY04p (); SKB09 (); HIG09 (), due to its relevance for biomolecules, molecular conductance, solar cell design; these are systems for which TDDFT would be particularly attractive due to its favorable system-size scaling. However, it has been challenging to find a satisfactory universal solution to the CT problem: ab initio approaches based on modeling the exact kernel appear impractical, while practical approaches tend to involve empirical parameters. Here we present a new approach to calculate CT excitations in TDDFT for certain cases, based on spin-symmetry-breaking. We show that accurate excitations are obtained when the acceptor is an effectively one-electron system, e.g. an element in Group 1 of the periodic table treated in the pseudopotential approximation, and justify why this is so. For large separations, the leading order behavior is captured solely from Kohn-Sham (KS) orbital energy differences. Results are given for model systems and for the LiH molecule, and suggest a type of Koopmans’ concept for one-electron systems.
The usual approximations in TDDFT notoriously underestimate CT excitations between fragments at large separation . To leading order in , the exact answer for the lowest CT frequency is:
where is the ionization energy of the donor, is the electron affinity of the acceptor and is the first electrostatic correction between the now charged species. (Atomic units are used throughout). It is well-understood why TDDFT severely underestimates CT DWH03 (); T03 (); GB04bNGB06 (); M05cMT06 (): In TDDFT, the first step is to compute the Kohn-Sham (KS) orbital energy differences between occupied () and unoccupied () orbitals, . In a second step, these frequencies are corrected to the true excitations via the Hartree-exchange-correlation kernel, , which shifts and mixes the KS excitations within a matrix formulation. The kernel is a functional of the ground-state density , with matrix elements . For CT excitations, the vanishing spatial overlap at large separations between occupied donor and unoccupied acceptor orbitals sitting on different nuclei means that the TDDFT predictions for CT excitations reduce to the KS orbital energy difference, when using usual semi-local functional approximations for . With approximate ground-state functionals the highest occupied molecular orbital (HOMO) energy, , underestimates the true ionization energy, while the lowest unoccupied molecular orbital (LUMO), , lacks relaxation contributions to the electron affinity. The last few years have seen many methods to correct the underestimation of CT excitations, e.g. Refs. TTYY04p (); SKB09 (); HIG09 (); most modify the ground-state functional to correct the approximate KS HOMO’s underestimation of , and mix in some degree of Hartree-Fock, and most, but not all SKB09 (); HIG09 () determine this mixing via at least one empirical parameter. Fundamentally, staying within pure DFT, both the relaxation contributions to and the tail in Eq. 1 come from , which must exponentially diverge with fragment separation GB04bNGB06 (); M05cMT06 (). Worse, in the case of open-shell fragments, not covered by most of the recent fixes, additionally the exact is strongly frequency-dependent M05cMT06 ().
The major reason for the awkward kernel structure in the case of open-shell fragments lies in the KS ground-state description: the HOMO (and LUMO) are delocalized over the whole molecule, quite distinct from the Heitler-London-like nature of the true wavefunction. In either the case of the exact or semi-local approximations, their orbital energy difference tends to zero as the molecule is pulled apart, and so must be responsible for the entire CT energy M05cMT06 (). It is long-recognized that this static correlation error MCY09 () is the root of the problem of poor ground-state energies, studied extensively in molecules like H GL76p (), and that a simple way out is to allow the system to break spin-symmetry. An unrestricted calculation with an approximate functional run on a diatomic molecule, leads, at a critical internuclear separation, to the spin-polarized solution obtaining the lowest ground-state energy. Albeit having incorrect spin-symmetry, accurate ground-state energies are achieved essentially because the KS description is rendered to have one electron on each atom.
Although less discussed, the same physics applies for heteroatomic molecules composed of open-shell fragments TMM09 (), and suggests that symmetry-breaking could be a means to obtain its CT excitations. If the exact functional were used, the lowest-energy state remains correctly spin-unpolarized at any , but at the cost of stark step and peak features in the bonding region, and strong-frequency-dependent structure in , difficult for approximations to capture. If instead, correct spin-symmetry is imposed on any existing approximate density-functional, the ground-state displays unphysical fractional charges at large , delocalized HOMO and LUMO orbitals, and again poor CT energies M05cMT06 (); TMM09 ().
The following examples show that remarkably accurate CT excitations can indeed be obtained from TDDFT via spin-symmetry breaking when the acceptor contains effectively one active electron; for large enough separations, these are contained in simply the bare KS excitations. We present first one-dimensional models that enable us to compare with highly accurate numerically-exact solutions (computed using a Runge-Kutta differential equation solver as implemented in the octopus code octopus ()), and to then analyze and understand in detail why the CT in symmetry-broken TDDFT is accurate via the underlying potentials. For the DFT calculations we study the exact-exchange (EXX) and local spin density approximation(LSD) CSS06 () with self-interaction correction (SIC). We shall see their correct long-range behavior yields the correct -dependence at large separations. We stress that this long-rangedness must be combined with symmetry-breaking in order to get accurate CT excitations: restricted calculations using these functionals fail CGGG00 (). The functionals are implemented in octopus within the KLI approximation to OEP KLI (). After studying the models, we then turn to the real LiH molecule. In our first model, the nuclear potentials are represented by short-range wells, at separation :
where the strengths au, and au. We place two electrons interacting soft-Coulombically via
into this heteroatomic molecule. The unrestricted KS calculation yields a spin-unpolarized solution until a separation of about a.u. (coinciding with an avoided crossing in the potential energy surfaces TMM09 ()) where it begins to break symmetry: the eigenvalues for the up and down spin in the ground-state begin to separate and slowly approach those of the isolated wells as the molecule is pulled further apart and each electron settles in a different well. In contrast to Ref. KMK08 (), localization of orbitals is in fact achieved within KLI: for two electrons KLI is equivalent to full OEP, but even for systems of more than two electrons (not shown here) KLI yielded spin-polarized localized orbitals.
Figure 1 plots the lowest orbital excitation energies of the model Eq. (2)-(3). The KS energy differences, especially those of EXX, capture the exact CT excitations throughout with remarkable accuracy!
Why this is so can be seen by studying the underlying KS potentials. Consider first the limit , where symmetry-breaking has placed, say, the spin-up(down) electron in the left(right) well in the ground state. The left well is the donor for the following discussion. Fig. 2 plots the KS potential for the -spin:
and its components at a separation of a.u. In Eq. 4, refers to the atomic acceptor(donor) potential (i.e. the right(left) well), and is the Hartree potential generated by density . As , in the vicinity of the donor and , so : the system is essentially a one-electron system in this limit and local excitations of the -electron are, correctly, just excited states of the donor. In the vicinity of the acceptor, where only a -electron lives, in the limit , and
The Hartree potential generated by the -electron and a small correlation contribution (the two dips in the red curve) result in a net upward shift of : the unoccupied -electron states living in the right-hand-well are shifted up in energy compared with those of the -electron (i.e. those of ). An approximate affinity level is thereby induced in the right-hand-well. So the bare KS orbital energy-differences yield accurate CT excitations: first, the HOMO for the -electron is the lowest orbital in the left-hand-well, for which , the exact ionization potential, due to Koopmans’ theorem Koopmans (), since both EXX and SIC are exact for one electron. Second, and more significantly, the LUMO approximates the affinity level of the right-hand-well, , sensing the presence of the -electron, i.e. Hartree-correlation relaxation contributions to the affinity are already incorporated at the bare KS level. We call this eigenstate of Eq. 5, and the corresponding state of , for which entirely analogous analysis holds, the “induced affinity” levels of the right and left atoms respectively. In the limit of infinite separation, they converge onto the lowest state for the unoccupied spin of the isolated one-electron atom. A key point is that and are different: in any restricted calculation where these are the same, unoccupied levels are excitations (of the same -electron system), not affinities.
How well does the induced affinity level approximate the true affinity generally? DFT theory tells us that the exact affinity of an -electron atom is where the middle expression is computed from total energy differences while the third expression is computed from the highest occupied KS eigenvalue of the relaxed -electron system (see eg. GKKG00p ()). Let us then consider the KS potential of a two-electron atom. Because the density is localized, an unrestricted calculation yields a spin-unpolarized result. Denoting the two-electron ground-state density ,
where is the nuclear potential and we have noted that, for exact exchange, . First, neglecting correlation, Eq. 6 is very close to Eq. 5 if , i.e. if, when a second electron is added to a well in which there is already one electron, there is little density relaxation. This is the case in the model example above, since the dominant part of the energy of the electrons is from the external potential. In such a case, Eqs. 5 and 6 then imply that the induced affinity level approximates the true affinity well. It will be a lower bound (i.e. ), because electron repulsion leads to being a little weaker than . This was borne out in all the EXX results of different models we considered. Correlation tends to raise the induced affinity, sometimes bringing it higher than the true affinity: certainly using LSD-SIC, forms a deeper negative well than . Shortly we will discuss examples in which the density relaxation is important, so that is not very close to , and there the affinity level is not such a good approximation to the true affinity; consequently the charge-transfer excitations are not as accurate.
The arguments above hold only for one-electron acceptors: if the acceptor already has an electron of the transferring spin in it, then excitations of that spin are the usual constant-number excitations of TDDFT, not approximate affinity levels. The donor however may contain any number of electrons: similar models that have, for example, three electrons in one well and one in another again showed excellent CT excitations from the former to the latter, under spin-symmetry-breaking.
We now extend the discussion following Eq. 4 to the case of an -electron donor and 1-electron acceptor, at finite but large separation. For ease of notation, assume again the acceptor carries a -electron in the molecular ground state. First consider Eq. 4 near the -electron donor. For external (nuclear) potentials that decay Coulombically, the (-electron), , where is the KS potential of the donor atom. At the acceptor, (noting that cancels , while ). So, to leading order in , for the -electron, and where is the induced affinity level of the acceptor in the infinite separation limit. Therefore,
as in Eq. 1, with and approximated by the induced affinity level. (For short-ranged potentials as in Eq. 2, the arguments lead instead to and ). Therefore, a long-ranged exchange-correlation as in EXX or SIC, once symmetry-broken, yields good CT excitations, from just its bare KS orbital energies, as a function of , for large . Note the importance of correct asymptotics of the functional used for correct -dependence, as well as for accurate ionization potentials and induced affinities.
For an accurate CT asymptote the density relaxation upon the addition of an electron must be small. But even when density relaxation effects are significant, symmetry-breaking can still be useful as we now explain. In a practical sense, the infinite-separation limit itself is not so much a problem for TDDFT, because total ground-state energy differences computed from DFT LFB10 () can often yield reliable values for and in Eq. 1. Rather it is intermediate but large distances that are the challenge, where CT energies deviate from the asymptotic formula Eq. 1. Our symmetry-breaking approach can capture these deviations, going beyond Eq. 1. The procedure is to compute the (symmetry-broken) KS HOMO and LUMO energy difference, but, when density-relaxation is large, shift it by
where and are computed from total ground-state DFT energy differences. In this way, asymptotically, the curves approach Eq. 1 accurately, but for intermediate to large distances, contain correct physical deviations from Eq. 1 due to polarization. To illustrate this, we now consider soft-Coulomb nuclear wells:
and Fig. 3 takes and . The more diffuse densities of these wells makes them more polarizable so deviations from Eq. 1 are more evident, as shown in Fig. 3: at intermediate separations, the (exact) CT energies fall shy of the asymptotic Eq. 1, shown as the black curve, due to the local polarization of the CT state towards the positive charge at the other nucleus. After applying the shift of Eq. 8 the unrestricted SIC results approach the exact results well, and capture this attractive shift; this holds also for CT in the other direction (not shown). (The blue curve is the asymptote for the unrestricted SIC prediction, but also shifted according to Eq. 8).
So far, we have discussed CT excitations obtained from bare KS excitations alone and argued why these work so well, as demonstrated by the model examples. The second step of TDDFT is to apply to correct the KS excitations towards the exact ones; there are both “diagonal” terms which shift each KS excitation, as well as “off-diagonal” ones that mix them. For the systems so far discussed we expect both these effects are small, because (i) the diagonal term involves overlap of the occupied and unoccupied orbitals in the excitation, which vanish exponentially at large distances, and (ii) there is little mixing with other excitations in the system. Mixing and shifting of KS excitations will be important at small and intermediate distances, leading to further deviations from the asymptotic Eq. 1, especially for real molecules, given their higher density of states.
Turning now to a real molecule, LiH: Using a pseudopotential for the Li atom renders it an effectively one-electron atom, and our method captures CT in both directions. In Fig. 4 we plot the lowest potential energy surfaces of the LiH molecule, computed with unrestricted SIC, with the Troullier Martins pseudopotential coded in octopus octopus (), compared to the highly-accurate configuration-interaction (CI) calculations of Ref. ADD09 (). The induced affinity level of H is 0.0726H while that computed from ground-state DFT energy-differences is 0.0264H, closer to the experimental (0.0277H), so we have applied the shift of Eq. 8 for . As accurate energies are unavailable for CT from H to Li, we do not show this curve here. The excellent agreement of the unrestricted approach with CI can be contrasted with restricted SIC calculations, whose collapse at smaller is due to the near-degeneracy of the HOMO and LUMO mentioned earlier; the latter leads to convergence difficulties for larger separations. The symmetry-broken SIC predicts the separation at which there is a crossing between the ionic curve and the Li(3s) curve very accurately, although it appears more as a direct crossing rather than an avoided one.
In summary, we have shown that symmetry-breaking is a simple non-empirical way to obtain CT excitations from KS orbital energies alone, for acceptors that contain effectively one electron, and explained why. Strikingly good results were obtained for model systems as well as for real molecule LiH and further studies are underway. Applying the TDDFT kernel should improve the accuracy at intermediate distances, capturing mixing of CT and local excitations, while these corrections will vanish asymptotically.
Conceptually, the observation that for one-electron systems, the levels of the unoccupied spin approximate affinity levels can be interpreted in an extended Koopmans’ sense. Koopman’s theorem in DFT states that exactly, while generalized Koopman’s theorem applies to Hartree-Fock where and , leading to the use of hybrid functionals for CT, mentioned earlier TTYY04p (); SKB09 (). Although the LUMO in exact DFT represents an excitation of the -electron system, rather than the -electron one, our results show that when in spin-DFT, the levels of the unoccupied spin can be interpreted in a generalized Koopmans’ sense, as they approximate affinity levels.
We acknowledge support from MEC (FIS2007-65702-C02-01), ACI-promociona (ACI2009-1036), Grupos Consolidados UPV/EHU del Gobierno Vasco (IT-319-07), e-I3 ETSF project (Contract No. 211956), the National Science Foundation (CHE-0647913), the Cottrell Scholar Program of Research Corporation and a grant of computer time from the CUNY High Performance Computing Center under NSF Grants CNS-0855217 and CNS-0958379.
- (1) Time-Dependent Density Functional Theory eds. M.A.L. Marques et al., (Springer, Berlin, 2006).
- (2) E. Runge and E. K. U. Gross, Phys. Rev. Lett. 52, 997 (1984).
- (3) A. Dreuw, J. Weisman, and M. Head-Gordon, J. Chem. Phys. 119, 2943 (2003).
- (4) D. Tozer, J. Chem. Phys. 119, 12697 (2003).
- (5) O. Gritsenko and E.J. Baerends, J. Chem. Phys. 121 655, (2004); J. Neugebauer, O. Gritsenko and E.J. Baerends, J. Chem. Phys. 124, 214102 (2006).
- (6) N.T. Maitra, J. Chem. Phys. 122, 234104 (2005); N. T. Maitra and D.G. Tempel, J. Chem. Phys. 125, 184111 (2006).
- (7) Y. Tawada et al., J. Chem. Phys. 120, 8425 (2004); O. A. Vydrov et al. J. Chem. Phys. 125, 074106 (2006); Y. Zhao and D. G. Truhlar, J. Phys. Chem. A. 110, 13126 (2006); M.A. Rohrdanz, K.M. Martins, and J.M. Herbert, J. Chem. Phys. 130, 054112 (2009); Q. Wu and T. van Voorhis J. Chem. Theor. Comp. 2, 765 (2006); J. Autschbach, ChemPhysChem 10, 1757 (2009).
- (8) T. Stein, L. Kronik and R. Baer, J. Am. Chem. Soc. 131, 2818 (2009).
- (9) A. Hesselmann, M. Ipatov, A. Görling, Phys. Rev. A. 80, 012507 (2009);
- (10) P. Mori-Sanchez, A. J. Cohen, W. Yang, Phys. Rev. Lett. 102, 066403 (2009).
- (11) O. Gunnarson and B. I. Lundquist, Phys. Rev. B. 13, 4274 (1976); J. Perdew, A. Savin and K. Burke, Phys. Rev. A. 51, 4531 (1995).
- (12) D. G. Tempel, T. J. Martinez, N.T. Maitra, J. Chem. Th. Comput. 5, 770 (2009).
- (13) M. Casula, S. Sorella, and G. Senatore, Phys. Rev. B. 74, 245427 (2006).
- (14) http://www.tddft.org/programs/octopus; M.A.L. Marques, A. Castro, G. F. Bertsch, A. Rubio, Comput. Phys. Commun. 151, 60 (2003).
- (15) M.E. Casida et al., J. Chem. Phys. 113, 7062 (2000).
- (16) J. B. Krieger, Y. Li, and G. J. Iafrate, Phys. Rev. A. 46, 5453 (1992).
- (17) T. Körzdörfer, M. Mundt, and S. Kümmel, Phys. Rev. Lett. 100, 133004 (2008).
- (18) T. Koopmans, Physica 1, 104 (1993); O. V. Gritsenko and E. J. Baerends, J. Chem. Phys. 117, 9154 (2002).
- (19) T. Grabo, T. Kreibich, S. Kurth, E.K.U. Gross, in Strong Coulomb Correlations in Electronic Structure Calculations: Beyond the Local Density Approximation, ed. V.I. Anisimov, (Gordon and Breach, 2000).
- (20) D. Lee, F. Furche, and K. Burke, J. Phys. Chem. Lett. 1, 2124 (2010).
- (21) M. Aymar, J. Deiglmayer, O. Dulieu, Can. J. Phys. 87, 543 (2009).