XDFT : an efficient first-principles method for neutral excitations in molecules

State-of-the-art methods for calculating neutral excitation energies are typically demanding and limited to single electron-hole pairs and their composite plasmons. Here we introduce excitonic density-functional theory (XDFT) a computationally light, generally applicable, first-principles technique for calculating neutral excitations based on generalized constrained DFT. In order to simulate an M -particle excited state of an N -electron system, XDFT automatically optimizes a constraining potential to confine N −M electrons within the ground-state Kohn-Sham valence subspace. We demonstrate the efficacy of XDFT by calculating the lowest single-particle singlet and triplet excitation energies of the well-known Thiel molecular test set, with results which are in excellent agreement with time-dependent DFT. Furthermore, going beyond the capability of adiabatic time-dependent DFT, we show that XDFT can successfully capture double excitations. Overall our method makes optical gaps, excition bindings and oscillator strengths readily accessible at a computational cost comparable to that of standard DFT. As such, XDFT appears as an ideal candidate to work within high-throughput discovery frameworks and within linear-scaling methods for large systems.

The first principles calculation of excited-state energies of quantum systems holds crucial importance for the study of solar cells 1 , organic light emitting diodes 2 , and chromophores in biological systems 3 , to name but a few. With some exceptions density-functional theory (DFT), which is the primary ab initio workhorse for computing ground state properties 4,5 , typically falls short on such tasks, although efforts are underway to extend the foundation of DFT to excited states [6][7][8][9][10][11] . The most commonly used first-principles method for calculating excitation energies, at least of finite systems, is perturbative time-dependent density functional theory (TDDFT) [12][13][14] . However, TDDFT has two significant limitations: 1) its considerable computational costs, which severely limit the size of the systems that it can investigate 15 , and 2) its inability to treat double (twoelectron) or higher-order excitations within adiabatic approximations to the exchange-correlation (XC) kernel [16][17][18] . Over the years several first-principles schemes based on time-independent DFT have been developed for calculating neutral excitation energies, such as ensemble DFT [19][20][21] , restricted open-shell Kohn-Sham DFT [22][23][24][25] , constricted variational DFT 26,27 , ∆SCF-DFT 28,29 and the maximum overlap method 30,31 . All of these approaches have their strengths and weaknesses in terms of computational costs and ease of both implementation and convergence. XDFT, like some of the latter methods, is motivated by the existence of a variational DFT, with a minimum principle and an equivalent noninteracting Kohn-Sham (KS) state, for an individual excited state of interacting electrons 6,7,28 . We refer the reader to Ref. (32) for a recent review of TDDFT, and to Ref. (26) for a foundational comparison between DFTbased variational approaches and TDDFT.
A neutral excitation, within the quasiparticle picture, is the promotion of one or more electrons from occupied levels to empty ones, resulting in the creation of bound electron-hole pairs, or excitons, and consequent energy-level relaxation. In this Letter, we introduce excitonic DFT (XDFT), an inexpensive, fully first-principles method based on constrained DFT (cDFT) [33][34][35] for calculating neutral excitation energies in finite systems such as molecules and clusters. XDFT scales with the atom count N as per ground-state DFT, namely as O(N 3 ). This contrasts with methods like TDDFT, which typically scales as O(N 4 ) 36 and the Bethe-Salpeter equation (BSE), which goes as O(N 6 ) 37 . In addition, unlike TDDFT and BSE, which are highly memory intensive, XDFT has a memory overhead comparable to that of standard DFT. Crucially, it avoids the calculation of unoccupied KS orbitals entirely. Therefore, in terms of computational efficiency, XDFT offers significant advantages over other methods and appears to be readily compatible with high-throughput frameworks, the study of large-systems, and KS methods beyond DFT.
Much like cDFT 34,35 , XDFT searches for the groundstate energy of a system subject to confining a given number of electrons with spin σ, N σ c , to a desired subspace. Such a constraining condition may be written as where 'Tr' denotes the trace,ρ σ is the spin-dependent Fermionic density operator andP is a projection operator onto the desired subspace. Then, the ground state of the system subject to the constraint is found at the stationary point of the functional where V c is a Lagrange multiplier. For a fixed V c , the second term on the right-hand side serves to modify the ground state potential by adding the term V cP . One then minimizes W [ρ, V c ] with respect toρ, just as E[ρ] in regular DFT. At the V c -dependent minima W [ρ, V c ] can be thought of as a function 35 , W (V c ), of V c . The maxima of W (V c ) occur at stable states of the constrained system 38 , at which the value of W is the total energy of interest. In conventional cDFT, the subspace spanned byP is a spatial region. IfP spans two spatial regions with opposite weighting, then one can enforce a chargeseparated density configuration for the simulation of charge-transfer excitations [39][40][41][42][43] . In this work, in order to access excitations beyond charge-separated states, we lift the restriction of the spatial confinement by defining the cDFT subspace in terms of the KS eigenstates. For a neutral N -electron system XDFT locates the energy of the lowest M -electron excited state by confining N − M electrons within the valence KS subspace of the unconstrained DFT ground-state. This circumvents the need for any prior, empirical specification of subspaces as in conventional cDFT. The projector that defines XDFT iŝ whereρ 0 is the ground state density operator, |ψ i is the i th KS orbital and f i is its occupation number. Once the energy of the first excited state, W , is determined by optimizing V c , then the lowest excitation energy, E * , can be evaluated as a total energy difference from the ground state DFT energy, E 0 , namely as XDFT is formally an orbital-dependent DFT, and its energy is separately invariant under arbitrary unitary transformations among the occupied Kohn-Sham orbitals of the ground-state and of the constrained ground-state. The KS wave-function of the excited state obtained from XDFT is orthonormal to that of the ground state. This is because it is a Slater determinant (SD) composed of KS orbitals the highest of which is, as a result of the constraint definition, orthonormal to each of the KS orbitals comprising the ground state KS wave-function.
XDFT can be used to simulate combinations of charge and spin excitations. Given a closed-shell ground state, triplet single-electron excitations and singlet double excitations are straightforward to access with a single constraint. These both incur the cost of just two DFT calculations -the ground-state one and the constrained one. In both cases, the electron-promotion constraint can be applied to the sum of the density operators for each spin, and the triplet state can be selected by setting m s = 1. Slightly more work is required to access singlet single and triplet double excitations, as we now discuss.
Given a closed-shell ground state, the final state of a singlet single excitation can not be represented by a single SD, and so the corresponding excitation energy cannot be obtained straightforwardly from a single cDFT calculation. Fortunately, for the non-interacting KS system both the closed-shell singlet excited state S=0 Ψ KS ms=0 (with a non-interacting energy S=0 E KS ms=0 ) and the open- ) can be written out as a linear combination of the same pair of SDs, within a frozen-orbital treatment. These two SDs are then degenerate, with a non-interacting energy SD E KS ms=0 . Invoking the multiplet sum method 32,44 , we can thereby express the non-interacting energy of a closed-shell singlet state approximately as In order to access one of these degenerate SDs that make up the singlet in practice, we apply the XDFT constraint to one spin channel only, using M = 1 and m s = 0. At this point, we make a final assumption that Eq. (3) may be used to approximately evaluate the energy of the interacting system. Keeping in mind that the three triplet states for m s = −1, 0, 1 are degenerate, we arrive at The advantage of Eq. (4) is that it involves only the energies of two single-SD states that are available using XDFT. This approximation has the same regime of validity as the frozen-orbital approximation for the KS orbitals. Each term in Eq. (4) derives from an interacting system that is obtained from an equivalent unrestricted KS system having the same density and spin density. Ultimately, the task of determining the triplet and singlet single excitation energies reduces to running the following three first-principles calculations (see Fig. 1): 1. A DFT calculation to determine the ground-state energy E 0 and density operatorρ 0 .
2. An XDFT calculation with m s = 1, confining N −1 electrons to the total valence subspace of the DFT run. This gives the energy of the lowest-lying interacting triplet state, S=1 E ms=1 .
3. An XDFT calculation with m s = 0, confining (N/2) − 1 electrons to the spin-up valence subspace of the DFT run, to find the energy SD E ms=0 .  [54] and [55]). The PBE 50 XC functional has been used in both cases. The dark and the light dots denote singlet and triplet gaps, respectively. The diagonal line indicates perfect agreement between XDFT and TDDFT.

Triplet
Next, we use Eq. (4) to approximate the energy S=0 E ms=0 of the excited singlet. Finally, we calculate the triplet and singlet neutral gaps, respectively, as We have implemented the XDFT formalism in the linear-scaling first-principles code onetep 45 , which variationally optimizes a minimal set of localized, nonorthogonal generalized Wannier Functions (NGWF), expanded in terms of psinc functions 46,47 , to minimize the total energy. onetep is equipped with an automated conjugate-gradients method for updating the cDFT (or XDFT) Lagrange multiplier 38,48,49 . We have used this, together with the Perdew-Burke-Ernzerhof (PBE) XC functional 50 to calculate the lowest singlet excitation energies of the 28 closed-shell organic molecules contained in the well-known Thiel set 51 . Our calculations are performed using scalar relativistic norm-conserving pseudopotentials, a plane-wave cutoff energy of 1500 eV and a radius of 14.0 a 0 for the NGWFs. The Martyna-Tuckerman periodic boundary correction scheme 52 was used with a parameter of 7.0 a 0 . The constrained KS system was found to contain symmetry-protected partialfilling of degenerate highest occupied state in certain molecules, and so we used finite-temperature ensemble DFT as implemented in onetep 53 in all cases.
In Fig. 2, we show a scatter plot of the singlet and triplet excitation energies calculated with XDFT against those obtained with linear-response TDDFT and adiabatic PBE in Ref. [54] (singlets) and Ref. [55] (triplets). The TDDFT results are generally in agreement with experimental values (see the supporting information in Ref. [56]). The triplet energies show an excellent agree- 3. Exciton charge density for the propanamide molecule, calculated with an isosurface value of ±0.05 eÅ −3 . Panel (a) shows the charge-density difference between the groundstate KS lumo and homo orbitals, while (b) and (c) show, respectively, the singlet and triplet exciton densities generated as a difference between the XDFT and DFT total densities. ment with TDDFT because, being SDs, the KS triplet excited states are directly accessible using a single XDFT calculation. The figure also demonstrates that, in spite of the multiplet sum approximation, XDFT calculates yields singlet energies with a remarkably good accuracy.

Exciton density
A plot of the difference in charge density between the excited state and the ground state provides a good approximation to the exciton charge density. In Fig. 3 we show such plots for a representative molecule of the Thiel set, propanamide. Fig. 3(a) shows an approximation to the exciton density based on the ground-state KS orbitals, which neglects the orbital relaxation and exciton binding. Since it captures these effects, the singlet (b) and triplet (c) isosurfaces generated using XDFT (and, in the case of the singlet, the multiplet sum method applied to the total electron densities) reflect a greater degree of exciton density localization than (a). Due to Pauli exclusion, furthermore, the singlet (b) exciton density attains a greater spatial localization than the triplet one (c).
In tests on excitations for which non-adiabatic linearresponse TDDFT is known to perform well, XDFT with m s = 1 yields the precisely the same total energy as an unconstrained DFT calculation with m s = 1 (i.e., DFT with a fixed spin-moment of 1 µ B ). Similarly, it does not affect the total energy if we apply an additional constraint to the spin population within the ground-state KS manifold. Based on this evidence, it appears that fixing the cDFT subspace to the ground-state KS valence manifold does not introduce any appreciable approximation to the definition of an excited state, at least in this 'adiabatic' regime. This is notwithstanding the fact that KS wave-function orthonormality does not imply the orthonormality of interacting wave-functions.
Before concluding our discussion on single-electron excitations, we note that higher-energy excitations can be simulated in XDFT by employing multiple constraints.
For example, if the valence subspaces of the ground state and the first XDFT excited state are projected onto bŷ P 0 andP 1 , respectively, then the energy of the second excited state can be found by confining N − 1 electrons within the subspace ofP 0 using a Lagrange multiplier V 1 c and, separately, confining N − 1 electrons within the subspace ofP 1 using a multiplier V 2 c . In general, the total-energy of the I th excited state system of a given spin symmetry will be found at the stationary point of Finally, we explore the ability of XDFT to calculate energies of excitations with strong double (two-electron) character, which are considered to be inaccessible to adiabatic TDDFT [16][17][18] . This is the case by construction within the linear response regime, but not necessarily so within full-response TDDFT. The XDFT method is nonperturbative, in that the Hartree and XC potentials are calculated self-consistently with the density in the excited state. Thus, XDFT is not limited to single excitations. In the benchmark case of atomic beryllium, the first double excitation promotes two electrons from the 2s to the 2p orbitals 57 . We calculated the singlet and triplet double excitations of Be by means of (the 1s electrons were pseudized, rendering multiplet summation unnecessary): 1. A ground-state DFT calculation.
2. An XDFT calculation with m s = 0, confining 0 electrons to the total valence subspace of the DFT run. This gives the energy S=0 E 2 ms=0 of the lowestlying interacting doubly-excited singlet state.
3. An XDFT calculation with m s = 1, confining 0 electrons to the total valence subspace of the DFT run. This yields the energy S=1 E 2 ms=1 .
The singlet and triplet double excitation energies were obtained as S=0 E 2 * = S=0 E 2 ms=0 − E 0 and S=1 E 2 * = S=1 E 2 ms=1 − E 0 , respectively. In Fig. 4 we plot the single and double excitation energies of Be calculated with semi-local and hybrid XC-functionals. These agree well with those calculated with ensemble DFT in Ref. [20], for all four excitation types. The singlet single-electron PBE excitation energy is also in very close agreement with our own linear-response TDDFT(PBE) result, indicating that the multiplet sum approximation is successful in this system. We note, however, that while our singlet double excitation energies agree well with experimental values, this is much less the case for our triplet double energies. Experimentally, the singlet 2s 2 → 2p 2 excitation is lower in energy than the triplet one, and this has been explained as resulting from a mixing of the singlet double with higher singlet single excitations 65 . Our results would support the opposite conclusion, however, since it is the triplet state which is poorly described. XDFT is capable of accessing excitations of non-integer electron character (e.g, mixed single and double excitations) with the aid of ensemble DFT 53,66 , in principle, and this is a promising avenue for future investigation.
In summary, we introduce the XDFT method for calculating the excited-state energies of finite systems by means of a small number of coupled DFT calculations. XDFT generalizes constrained DFT, in essence, by removing the necessity for users to pre-define the targeted subspaces. Unlike linear-response TDDFT or BSE, no reference is made to unoccupied orbitals. XDFT closely reproduces the TDDFT values for triplet and also, with the help of an additional approximation, in most cases the singlet excitation energies of the Thiel molecular test set. Unlike adiabatic TDDFT, however, XDFT can readily access the energies of double excitations, effectively circumventing the requirement for non-adiabaticity in TDDFT. We demonstrate this in a successful application to the well-known test case of the beryllium atom.
This work is supported by the European Research Council project quest. We acknowledge G. Teobaldi and N. D. M. Hine for their implementation of cDFT in FIG. 4. The excitation energies of atomic Be from experiment, linear-response TDDFT (the ONETEP linear-scaling implementation 58,59 ), and XDFT. The top panels show singleelectron excitation energies and the bottom panels show double (two-electron) excitation energies. The functionals tested were the generalized gradient approximation parameterizations PBE 50 and RPBE 60 , and the hybrid functionals PBE0 61 , B3LYP 62 , and B1PW91 63 . For the single excitations, we have included the results of linear-response TDDFT calculations within adiabatic PBE, where double excitations are inaccessible. TDA here refers to the same calculations but within the Tamm-Dancoff approximation 64 . We also provide experimental values taken from Ref. [65].