Direct extraction of excitation energies from ensemble density-functional theory

TL;DR

This paper introduces a new ensemble density-functional theory approach that directly computes excitation energies as simple corrections, avoiding complex calculations required by traditional TDDFT, and captures double excitations.

Contribution

It presents an exact formula for excitation energies from a specific ensemble, enabling direct extraction without extensive unoccupied state calculations or multiple ensemble iterations.

Findings

SEHX approximation yields results comparable to TDDFT for atoms.

The method captures double-excitations missed by adiabatic TDDFT.

Avoids the need for many unoccupied orbitals or self-consistent ensemble calculations.

Abstract

A very specific ensemble of ground and excited states is shown to yield an exact formula for any excitation energy as a simple correction to the energy difference between orbitals of the Kohn-Sham ground state. This alternative scheme avoids either the need to calculate many unoccupied levels as in time-dependent density functional theory (TDDFT) or the need for many self-consistent ensemble calculations. The symmetry-eigenstate Hartree-exchange (SEHX) approximation yields results comparable to standard TDDFT for atoms. With this formalism, SEHX yields approximate double-excitations, which are missed by adiabatic TDDFT.

Figures

Tables

Table 1. Table 1: Exact and approximate singlet Δ ​ ω I Δ subscript 𝜔 𝐼 \Delta\omega_{I} (in mH) of the 1D two-electron contact-interaction Hooke’s atom with λ = 0.2 𝜆 0.2 \lambda=0.2 . The dressed TDDFT results are calculated with an exchange-only frequency-dependent kernel Maitra  et al. ( 2004 ) .

		$\Delta {\omega }_I$
			TDDFT	DEC		TDDFT
$I$	${\omega }_{KS}$	exact	AEXX	SEHX	exact	dressed
Singles
1	962	38	39	39	38	39
3	1953	47	30	30	48	49
5	2948	52	25	27	51	54
Doubles
2	1923	41	–	58	41	39
4	2915	49	–	77	49	47

Equations

$$E_I(\{{\mathtt{w}}_i\})=\sum _{i=0}^I{g}_i{\mathtt{w}}_i{E}_i,\sum _{i=0}^I{g}_i{\mathtt{w}}_i=1,$$

$$n^{\text{ens}}(\mathbf{r})=\sum _{i=0}^I{\mathtt{w}}_i{\tilde{n}}_i(\mathbf{r}),$$

$${\mathtt{w}}_{i\ne I}=(1-g_I\mathtt{w})/M_{I-1},{\mathtt{w}}_I=\mathtt{w},(\mathrm{GOK})$$

$${\omega }_I=\frac{1}{g_I}{\frac{d{E}_I^{\mathrm{GOK}}}{d\mathtt{w}}|}_{{\mathtt{w}}_I}+\sum _{i=1}^{I-1}\frac{1}{M_i}{\frac{d{E}_i^{\mathrm{GOK}}}{d\mathtt{w}}|}_{{\mathtt{w}}_i},$$

$${\omega }_I=M_I{E}_I(\mathtt{w}=M_I^{-1})-M_{I-1}{E}_{I-1}(\mathtt{w}=M_{I-1}^{-1}),$$

$${\mathtt{w}}_0=\frac{1-\mathtt{w}(M_I-g_0)}{g_0},{\mathtt{w}}_{i\ne 0}=\mathtt{w}.(\mathrm{GOKII})$$

$${\omega }_I=\frac{1}{g_I}\left[{\frac{d{E}_I^{\mathrm{GOKII}}}{d\mathtt{w}}|}_{{\mathtt{w}}_I}-{\frac{d{E}_{I-1}^{\mathrm{GOKII}}}{d\mathtt{w}}|}_{{\mathtt{w}}_{I-1}}\right].$$

$$\Delta\omega_{I}=\frac{1}{g_{I}}\frac{d}{d\mathtt{w}}\Bigg{|}_{\mathtt{w}=0}\left(E^{\rm GOKII}_{\scriptscriptstyle\rm{XC},I}-E^{\rm GOKII}_{\scriptscriptstyle\rm{XC},I-1}\right),~{}~{}~{}~{}({\rm DEC})$$

$$\Delta {\omega }_I^{\text{SEHX}}=H_I/g_I-H_0/g_0.$$

$$H_i^{\mathrm{orb}}=\frac{1}{2}\int \frac{{\mathrm{d}}^{3}r{\mathrm{d}}^3{r}^{\prime }}{\left|\mathbf{r}-{\mathbf{r}}^{\prime }\right|}\mathrm{tr}\left\{\mathbf{V}(\mathbf{r},{\mathbf{r}}^{\prime })\cdot {\mathbf{Q}}_i\right\},$$

$$H_i^{\mathrm{dens}}=-\int d^{3}r{v}_{\mathrm{HX}}(\mathbf{r}){\tilde{n}}_i(\mathbf{r}).$$

$$\hat{H}=\frac{1}{2}\sum _{i=1}^2\left(-\frac{{\partial }^2}{\partial x_i^2}+x_i^2\right)+\lambda \delta (x-x^{\prime }),$$

Full text

The Hohenberg-Kohn (HK) theorem Hohenberg and Kohn (1964); Levy (1979, 1982); Lieb (1983) of ground-state density-functional theory (DFT) Hohenberg and Kohn (1964); Kohn and Sham (1965) has several parts. The most-used in practice is the establishment of an exact density functional, $F[n]$, whose minimum yields the exact ground-state density and energy of a given system. Almost all practical calculations use the Kohn-Sham (KS) scheme Kohn and Sham (1965) to minimize $F$ with an approximation to the small exchange-correlation contribution, $E_{\scriptscriptstyle\rm XC}[n]$. In fact, many properties of interest in a modern chemical or materials calculation can be extracted from knowledge of the ground-state energy as a function of nuclear coordinates, or in response to a perturbing field.

However, except under very special circumstances, most optical excitation frequencies cannot be deduced. Hence there has always been interest in extending ground-state DFT to include such excitations. Moreover, another part of the HK theorem guarantees that such frequencies (and all properties) are indeed functionals of the ground-state density. In recent years, linear-response time-dependent DFT (TDDFT) Runge and Gross (1984); Casida (1996); Marques et al. (2012); Ullrich (2012); Ullrich and Yang (2014) has become a popular route for extracting low-lying excitation energies of molecules, because of its unprecedented balance of accuracy with computational speed Maitra (2016). For significantly sized molecules, more CPU time will be expended on a geometry optimization than a single TDDFT calculation on the optimized geometry.

However, while formally exact, TDDFT with standard approximations is far from perfect. If the unknown exchange-correlation (XC) kernel of TDDFT is approximated by its zero-frequency (and hence ground-state) limit, no multiple excitations survive Maitra (2016). While a useful work-around exists for cases where a double is close to a single excitation Maitra et al. (2004); Cave et al. (2004), there is as yet no simple and efficient general procedure for extracting double excitations within adiabatic TDDFT Elliott et al. (2011).

Ensemble DFT (EDFT) Theophilou (1979); Gross et al. (1988a) applies the principles of ground-state DFT to a convex ensemble of the lowest $M$ levels of a system, for which a KS system can be defined Gross et al. (1988b). EDFT is formally exact, but practical calculations require approximations, and initial attempts yielded disappointing results Oliveira et al. (1988). Accuracy is greatly improved when so-called “ghost interactions” between distinct states are removed from the approximations Gidopoulos et al. (2002). EDFT remains an active research area because, being variational, it should not suffer from some of the limitations of standard TDDFT. Recent strides by Pernal and coworkers Pastorczak et al. (2013); Pastorczak and Pernal (2014), Fromager and coworkers Senjean et al. (2015); Deur et al. (2017), and others attempt to create a useful practical alternative to TDDFT, but the difficulty remains in finding accurate low-cost approximations. EDFT usually requires running several different self-consistent ensemble calculations to extract several low-lying excitations.

Here we (a) derive a formula from EDFT to correct a KS orbital energy difference into an exact excitation energy, without doing any self-consistent ensemble calculations, (b) argue that its computational cost should typically be less than either standard TDDFT or EDFT, (c) calculate this correction using the symmetry-eigenstate Hartree-exchange (SEHX) approximation Pribram-Jones et al. (2014); Yang et al. (2014); sup for atoms, demonstrating its accuracy relative to standard TDDFT, and (d) show that SEHX estimates double excitations.

EDFT is a formally exact and variational excited-state method Theophilou (1979); Gross et al. (1988a, b). Let $E_{i}$ be the electronic energy levels, $i=0,1,...$, each with degeneracy $g_{i}$. Construct an ensemble from positive convex weights $w_{i}$, letting $I$ be the maximum non-zero weight. The weights are not variational parameters. Then, from the foundational theorems, the ensemble energy

[TABLE]

is a functional of the ensemble density

[TABLE]

where $\tilde{n}_{i}({\bf r})$ is the sum of all densities in the $i$th multiplet (so that $\int d^{3}r\,\tilde{n}_{i}({\bf r})=g_{i}N$, with N being the number of electrons), and can be found via a minimization, so long as the weights are monotonically non-increasing. Applying the same conditions to a fictitious system of non-interacting electrons with the same weights, one can define a KS system whose ensemble density matches the interacting one. Defining energy components in the usual way, only the XC contribution needs to be approximated to perform an ensemble DFT calculation. Since the ensemble energies depend linearly on the weights (at least, in the exact theory), one can easily deduce transition frequencies.

Infinitely many ensembles can be realized, but the GOK ensemble from the original work Gross et al. (1988b) is particularly useful and popular, in which all weights are the same except for the highest multiplet, i.e.,

[TABLE]

where $M_{I}$ is the number of states up to and including the $I$-multiplet, and $\mathtt{w}\leq M_{I}^{-1}$ to preserve convexity. When $\mathtt{w}=M_{I}^{-1}$, the weights are all equal (an equiensemble). In general, the corresponding ensemble density must be found by self-consistent solution of the ensemble KS equations, for the given weights. The excitation energy of the $I$ multiplet can only be isolated by performing self-consistent calculations for all lower multiplets. The excitation energy from the ground state to the $I$th multiplet is Gross et al. (1988b)

[TABLE]

requiring $I+1$ self-consistent calculations, including the ground-state KS calculation, where the density is held fixed when the derivative is taken.

The weights defined by Eq. (3) are also a linear interpolation between two consecutive equiensembles, containing $M_{I-1}$ and $M_{I}$ states. Thus $\omega_{I}$ can also be calculated via

[TABLE]

which requires only two self-consistent calculations. However, if one needs all excitation energies up to $\omega_{I}$, $I+1$ self-consistent calculations are still needed. The computational costs of Eqs. (4) and (5) are much higher than TDDFT with the Casida equation Casida (1996).

Now we reintroduce an alternative one-parameter ensemble, in which all states have weight $\mathtt{w}$ except the ground state:

[TABLE]

We say reintroduce, as this ensemble was mentioned in a footnote in Ref. Gross et al. (1988b), although never applied (as far as we know). However, we can show (see supplemental material sup ) that the excitation energy using Eq. (6) has a much simpler formula than using Eq. (4):

[TABLE]

Despite the simplicity, in general one still needs to do $I+1$ calculations to get all excitation energies. However, unlike Eq. (3), the set of weights defined by Eq. (6) is now a linear interpolation between the ground state and the equiensemble of $M_{I}$ states. Now, $\mathtt{w}=0$ recovers the ground state, not an equiensemble with one less multiplet. A further simplification is made by noting that the EDFT formalism is valid even as $\mathtt{w}\to 0$. Setting $\mathtt{w}_{I}=\mathtt{w}_{I-1}=0$ in Eq. (7) and defining $\Delta\omega_{I}=\omega_{I}-\omega_{I}^{\text{KS}}$, where $\omega_{I}^{\text{KS}}$ is the KS orbital energy difference, yields

[TABLE]

where $E^{\rm GOKII}_{\scriptscriptstyle\rm{XC},I}$ is the GOKII ensemble XC energy functional Gross et al. (1988b) containing up to the $I$th multiplet. This is a direct ensemble correction (DEC) to the KS transition frequencies.

Equation (8) is the central formal result of this work. Because all elements of the right-hand side are evaluated on the ground-state density, this correction is a formally exact ground-state density functional for correcting KS transitions into physical transitions. If approximated by an explicit density functional, it could be evaluated at no noticeable additional cost to a standard ground-state DFT calculation. Compared with the cubic scaling of the TDDFT linear response equations Casida (1996), Eq. (8) is vastly more efficient. On the other hand, TDDFT yields both transition frequencies and oscillator strengths, as well as dipole overlap matrices. In addition, linear response TDDFT can yield spatially resolved response functions, once perturbations different from a long-wavelength electric field are allowed. In future work, we will explore what else, beyond transition frequencies, might be extracted in a manner similar to Eq. (8).

There is an infinite number of excited-state ensembles. Even if we consider only those that interpolate between the ground state and the equiensemble, Eq. (6) is not the only choice. The exact ensemble functional yields the same result in any ensemble, but approximations yield different results for different ensembles. A DEC expression is a particularly simple route to excitation energies.

Eqs. (3) and (6) are identical for a simple bi-ensemble, the ensemble of a non-degenerate ground and first excited states. Studies of $\mathtt{w}=0$ bi-ensembles have been carried out previously Levy (1995), as well as calculations of the first excitation energy Pribram-Jones et al. (2014); Yang et al. (2014). Thus the DEC of Eq. (8) can be viewed as a generalization of such results to an arbitrary excitation.

The exact $E_{\scriptscriptstyle\rm XC}^{\rm GOKII}$ of Eq. (8) can be obtained numerically for simple cases Pribram-Jones et al. (2014); Yang et al. (2014), but in practical calculations $E_{\scriptscriptstyle\rm XC}^{\rm GOKII}$ must be approximated. In general, the $E_{\scriptscriptstyle\rm XC}^{\rm GOKII}$ of Eq. (8) must account for the state ordering and differences in multiplet structure between the real and KS systems, which poses a challenge for the development of approximations.

SEHX Pribram-Jones et al. (2014); Yang et al. (2014) is an explicit orbital-dependent ensemble-density functional generalization of the exact-exchange approximation (EXX) of ground-state DFT, whose full expression is given in the supplemental materialsup . Using the energy decomposition of Nagy Nagy (2001, 2002), SEHX constructs the combined Hartree-exchange energy from an ensemble sum over spin- and spatially-symmetrized multi-determinant KS wavefunctions, removing “ghost interactions” and approximating the ensemble discontinuity Levy (1995), and yielding good results in the GOK ensemble Pribram-Jones et al. (2014); Yang et al. (2014). Inserting SEHX into Eq. (8), all the contributions from excitations below $I$ cancel, yielding an approximation that depends only on the difference between a contribution from the $I$th multiplet and the ground state:

[TABLE]

Here $H_{i}=H_{i}^{\rm orb}+H_{i}^{\rm dens}$, where

[TABLE]

and

[TABLE]

${\bf V}$ is a matrix containing products of KS orbitals, ${\bf Q}_{i}$ is a matrix containing orbital occupation factors and symmetrization coefficients of KS determinants (see supplemental material sup ), and $v_{\scriptscriptstyle\rm HX}({\bf r})$ is just the ground-state Hartree-exchange potential. As our tests are on atoms and ions, we use the KLI approximation Krieger et al. (1992) for $v_{\scriptscriptstyle\rm HX}({\bf r})$ to obtain more accurate orbital energies than those from semilocal approximations Petersilka et al. (2000). We denote calculations with Eq. (9) as DEC/SEHX. Unlike Eq. (8), Eq. (9) depends only on the ground and excited states in question, so the state ordering problem is bypassed and calculation is highly efficient. The ordering-independency of Eq. (9) is due to SEHX, yielding reasonable excitation energies even if the KS state ordering is different from the real one. On the other hand, the approximate state ordering might not be correct.

To illustrate the performance of DEC/SEHX, we calculated excitation energies of small atoms. To see exclusively the effect of the excitation method Kim et al. (2013), we use the exact KS potential and energies for the He and Be atoms Umrigar and Gonze (1994a, b). We compare with TDDFT using the adiabatic local density approximation (ALDA) Ullrich (2012). For simplicity, we use the Tamm-Dancoff approximation (TDA) Hirata and Head-Gordon (1999) in TDDFT calculations, and we checked to make sure that the results only change slightly with and without TDA. The results are shown in Figs. 1 and 2. More results for atoms are available in the supplemental material sup .

Figure 1 shows the He results. These are all single excitations (as all doubles in He are in the continuum). The DEC/SEHX gives results that are qualitatively similar to those of standard TDDFT. In fact, the mean absolute errors are typically about 30% smaller, despite the lack of approximate correlation in the DEC calculation.

Figure 2 shows the results for Be, again with the exact KS potential. For single excitations, the results are quantitatively similar to those of He, again with DEC errors being noticeably smaller than their TDDFT/ALDA counterparts. But in DEC we can also calculate the double excitations, which are completely absent from any adiabatic TDDFT calculation. We note that the double excitations are less accurate than their single counterparts, but since there are only two, this might be incidental. The supplemental material sup gives many more atomic calculations, using approximate ground-state KS potentials, showing the strong sensitivity of both DEC and TDDFT to the KS levels in atoms.

To better understand the performance of DEC/SEHX for the double excitations, we turn to a much simpler model problem that was designed to study precisely this question. Consider two fermions in a 1d harmonic potential with a contact interaction Maitra et al. (2004); Zhang and Burke (2004):

[TABLE]

where $\lambda>0$. For small values of $\lambda$, the system is weakly interacting, and exchange-type approximations are accurate.

The results are shown in Table 1. Because this is 1d, there are no degeneracies or multiplets. However, this model was purposely constructed to have near-degeneracies between the multiple and single excitations. With the harmonic confining potential, as $\lambda\to 0$, many levels approach one another. As shown, the double excitation of level 2 is very close to the single of level 3, and the double at level 4 is very close to the single at level 5.

In the 4th and 5th columns of the table, we report exact exchange results. The former is TDDFT, using the exact KS potential and the exact ground-state exchange in an adiabatic approximation. The latter is DEC/SEHX. We see that both are excellent approximations to the lowest excitations, and give almost identical results for the single excitations. This is because $\lambda=0.2$, ensuring that correlation effects are relatively weak. But, unlike adiabatic TDDFT, DEC/SEHX also yields predictions for the double excitations. Just like in the atoms, the errors are substantially larger for the doubles.

Because this model has only two electrons, we can calculate the exact DEC numerically with Eq. (8), by calculating the exact energies, densities and the xc potential of the model first. We then evaluate Eq. (8) numerically using these exact quantities (see supplemental material sup ). These DEC/exact results are in column 6, and agree within a mH with the exact results. This shows that exact DEC does handle doubles correctly, so that the failing in DEC/SEHX is due to the lack of correlation. The last column of the table shows results with the DSPA, a frequency-dependent model XC kernel designed for weakly-correlated systems with strong coupling between a single and double excitation, often called dressed TDDFT Maitra et al. (2004); Cave et al. (2004). This works extremely well here, as this system was designed to illustrate its accuracy. Study of the difference in the results between these two should provide a route to improving DEC approximations for double excitations.

A discerning reader might have noted that, throughout this work, we have avoided discussion of $N$- and $v$-representabilityParr and Yang (1989); Gilbert (1975); Levy (1982). These issues have been partially explored Oliveira et al. (1988); Hadjisavvas and Theophilou (1985) within EDFT in general, but not for this particular ensemble. But none of the calculations here ran into any representation difficulties, such as an inability to find a KS system with the required density. There is little reason to fear such problems in practice. Furthermore, as we use only DEC, any such difficulties in EDFT in general are likely to be least problematic for our applications.

There is obviously much work to be done to see if DEC can become competitive with standard TDDFT calculations. It should be applied to molecules with standard ground-state functionals, to see if the results are as accurate or if semilocal ground-state approximations destroy the accuracy found here within SEHX. Other challenges for TDDFT, such as charge transfer excitations, should be carefully tested. In such a case, we are less hopeful that DEC will provide accurate results as, like TDDFT, it also begins from (unrelaxed) KS transitions of the ground-state. Other ensembles might also yield direct ensemble corrections, or properties other than simple excitation energies might be accessible.

Several other EDFT-based methods for excitations were recently proposed, such as the linear interpolation method Senjean et al. (2015), the Helmholtz free-energy minimization method Pastorczak et al. (2013), and the ensemble-referenced Kohn-Sham method (REKS) Filatov (2015). The REKS method is a multi-reference extension to ground-state DFT and EDFT (see also Pastorczak and Pernal (2016)), while the others are within standard EDFT. Each has its own advantages, and the REKS method has been shown to work well in strongly-correlated systems Filatov et al. (2015). However, all these methods require extra self-consistent calculations aside from the ground-state one. The simplification achieved in this paper by changing the ensemble type suggests that similar simplifications may also be possible in these methods. Another route for future research would bypass the use of ensemble functionals altogether by developing approximate methods based on the DEC.

In summary, DEC (in Eq. (8)), is a formally exact approach to excitation energies from DFT, as illustrated by our model harmonic trap calculation. For example, where the fundamental and optical gaps match (insulating solids without excitons), DEC yields a new approach to the problem of finding accurate gaps within DFTPerdew et al. (2017), relating the derivative discontinuities with respect to particle numberPerdew et al. (1982) and those with respect to optical excitationLevy (1995). While DEC and TDDFT are both post-processing steps after a ground-state KS calculation, DEC is less expensive and applicable to traditionally difficult problems such as multiple excitations and spin-multiplets. Unlike TDDFT, EDFT is based on a variational principle Gross et al. (1988a), so the DEC derived in this work may be more reliable than TDDFT corrections, which are based on response theory. The calculations shown in this paper merely demonstrate the DEC method: SEHX yields better accuracy than TDDFT/ALDA for single excitations in atoms, and approximates doubles (albeit less accurately than singles). Simpler approximations, avoiding solution of OEP-type equations, might produce usefully accurate results for valence excitations in molecules. Thus DEC represents an exciting alternative to TDDFT.

The authors thank Cyrus Umrigar for providing us the exact KS potentials of the He and Be atoms. Z.-H.Y. is currently supported by Science Challenge Project No. TZ2016003 (China). Z.-H.Y. and C.A.U. were supported by NSF grant DMR-1408904. A.P.J. was supported by the University of California President’s Postdoctoral Fellowship. Part of this work was performed under the auspices of the U.S. Department of Energy by Lawrence Livermore National Laboratory under Contract DE-AC52-07NA27344. K.B. was supported by DOE grant DE-FG02-08ER46496.