Neutral excitation density-functional theory: an efficient and variational first-principles method for simulating neutral excitations in molecules

We introduce neutral excitation density-functional theory (XDFT), a computationally light, generally applicable, first-principles technique for calculating neutral electronic excitations. The concept is to generalise constrained density functional theory to free it from any assumptions about the spatial confinement of electrons and holes, but to maintain all the advantages of a variational method. The task of calculating the lowest excited state of a given symmetry is thereby simplified to one of performing a simple, low-cost sequence of coupled DFT calculations. We demonstrate the efficacy of the method by calculating the lowest single-particle singlet and triplet excitation energies in the well-known Thiel molecular test set, with results which are in good agreement with linear-response time-dependent density functional theory (LR-TDDFT). Furthermore, we show that XDFT can successfully capture two-electron excitations, in principle, offering a flexible approach to target specific effects beyond state-of-the-art adiabatic-kernel LR-TDDFT. Overall the method makes optical gaps and electron-hole binding energies readily accessible at a computational cost and scaling comparable to that of standard density functional theory. Owing to its multiple qualities beneficial to high-throughput studies where the optical gap is of particular interest; namely broad applicability, low computational demand, and ease of implementation and automation, XDFT presents as a viable candidate for research within materials discovery and informatics frameworks.

Over the years, several such first-principles schemes 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 , Δ self-consistent field (ΔSCF) -DFT 28,29 , constrained DFT defined using virtual Kohn-Sham states 30 , and the maximum overlap method 31,32 . These are underpinned by the existence of a variational DFT, with an equivalent non-interacting Kohn-Sham (KS) state, for an individual excited state of interacting electrons 6,7,28 . Each method, including the one introduced here which builds upon the foundation that others have provided, has its relative strengths and weaknesses in terms of computational cost and ease of both implementation and convergence. We refer the reader to ref. 33 for a recent review of TDDFT, and to ref. 26 for a foundational comparison between TDDFT and DFT-based variational approaches.
In this article, we introduce neutral excitation DFT (XDFT), an inexpensive, fully first-principles extension of constrained DFT (cDFT) [34][35][36] for calculating neutral excitation energies in finite systems such as molecules and clusters. XDFT simulates an excitation in the KS system by reducing the electronic population of the ground state valence subspace by one electron, while keeping the total number of electrons unchanged. XDFT is quite unlike conventional cDFT in that no prior assumptions are made as to the spatial form of source or drain regions for charge constraint. It captures screening effects at all orders, unlike LR-TDDFT, but retains a single Slater determinant of Kohn-Sham orbitals and so it is readily portable to many standard DFT codes (this transpires to be drawback, as we will see, in the study of singlet excitations, motivating future developments). XDFT scales with the atom count, N, as per ground-state DFT, namely as O N ( ) 3 when no quantum nearsightedness is exploited (we have implemented it within an O N ( ) code). This contrasts with methods like TDDFT, which typically scales as O N ( ) 4 37 and the Bethe-Salpeter equation (BSE), which goes as O N ( ) 6 38 . In addition, unlike LR-TDDFT and BSE, which can be highly memory intensive when unoccupied states are treated explicitly, XDFT has a memory overhead comparable to that of standard DFT. Crucially, it avoids the calculation of unoccupied ground-state KS orbitals entirely, and we never calculate them here in practice. Therefore, in terms of computational efficiency, XDFT offers significant advantages over the aforementioned existing methods, as long as only the lowest-energy excited state of a given symmetry is of particular interest. XDFT offers ready compatibility with high-throughput frameworks 39 , the study of large-systems, and variational KS methods beyond DFT.
Applying this technique, we calculate the lowest singlet and triplet single excitation energies of a representative set of organic molecules. Here we find a surprisingly good agreement with LR-TDDFT values, notwithstanding the extreme simplicity of our method. We also move straightforwardly beyond single-particle excitations, while still only using a small number of coupled DFT calculations. As we will see, XDFT can accurately reproduce the energies of excited states with a predominantly double-excitation character in the canonical test-bed atom beryllium, which are inaccessible [16][17][18] to adiabatic LR-TDDFT.
In the next section, we describe the XDFT formalism. This is followed by a section that explores the connection between XDFT and exact theorems of excited state DFT and outlines the relevant approximations. The next section outlines certain important details of the calculation. Results, including single and double excitation energies and difference density are presented in the penultimate section. Finally, we present our conclusions.

The XDFT Formalism in Brief
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 and with consequent energy relaxation due to screening. To simulate this excited state, XDFT searches for the ground-state energy of the same system, but now subject to the extra condition of a given number of electrons with spin σ, σ N c being confined to a pre-defined subspace. This constraining condition may be written as c where 'Tr' denotes the trace, ρ σ is the spin-dependent fermionic density operator and  is a projection operator onto the desired subspace. In conventional cDFT, the subspace spanned by  is a spatial region defined at the researcher's discretion. If  spans two spatial regions with opposite weighting, for example, then one can enforce a charge-separated density configuration for the simulation of charge-transfer excitations [40][41][42][43][44] .
Here we come to the key development of XDFT. In order to access excitations beyond charge-separated states of obvious spatial character, we free cDFT of all human assumptions by defining the subspace in terms of the ground-state KS eigenstates only. In doing so we retain the physical information encoded in the KS eigensystem, which is a key ingredient in LR-TDDFT but discarded in normal cDFT. More technically speaking, in a neutral N-electron system XDFT locates the energy of the lowest M-electron excited state (M may even be non-integer in principle) 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 and restores first-principles status. The projector that defines XDFT is where ρˆ0 is the ground state density operator, ψ i is the i th KS orbital and f i is its occupation number and the sum is over all KS levels. Note that, at zero temperature, = f 1(0) i for occupied (unoccupied) levels. As in conventional cDFT, the ground state of the system subject to the "exciting" constraint in XDFT is found at the stationary point of the functional www.nature.com/scientificreports www.nature.com/scientificreports/ 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 c . One then minimises ρŴ V [ , ] c with respect to ρˆ, just as ρÊ[ ] is minimized in regular DFT. At the V c -dependent minima ρŴ V [ , ] c can be thought of as a function 36 , W(V c ), of V c . The maxima of W(V c ) occur at stable states of the constrained system 45 , at which the value of W is the excited-state total energy of interest. Once the energy of the first excited state, W, is determined by optimising V c , then the lowest excitation (photon absorption) energy, ⁎ E , can be evaluated as a total energy difference from the ground state DFT energy, E 0 , namely as = − ⁎ E W E 0 . We note that excitations beyond the lowest-energy one of a given spin symmetry can be simulated in XDFT by employing multiple constraints. For example, if the Kohn-Sham valence subspaces of the ground state and the first XDFT excited state are projected onto by  0 and  1 , respectively, then the energy of the second excited state can be found by confining − N 1 electrons within the subspace of  0 using a Lagrange multiplier V c 1 and, separately, confining − N 1 electrons within the subspace of  1 using a multiplier V c 2 . 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 XDFT can be used to simulate combinations of charge and spin excitations. Given a closed-shell ground state, the double ( = M 2) excitations and the triplet single excitation 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 1 s . We have found that, for a common test set with different types of chemical bonds, XDFT allows us to directly access the lowest-lying excitations using coupled ground-state calculations. Non-linear response effects such as many-particle excitations are treated on the same footing as linear-response effects. In essence, while it has long been known that ground-state DFT contains the information needed to calculate neutral excitations, and relatively complex methods have been developed to explore it, we show here that constrained Kohn-Sham DFT at a level easily implementable in any normal DFT code provides a variational approach to exploit this. We note in passing that, since the excited state KS wavefunctions are accessible through XDFT, one can potentially use them to calculate approximate oscillator strengths. This is an interesting topic for future investigations.

Formal Justification of the XDFT Method
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 state. The XDFT constraint, which, for simulation of the first excited state, amounts to ejecting one electron from the valence subspace, encodes a well-defined many-particle excitation of the non-interacting Kohn-Sham system. Now, while this certainly does not imply a well-defined excitation of the interacting system, it may provide a good approximation in cases where the non-interacting and interacting wave-functions are similar (modulo unitary transformations) in both the ground and excited states. This is expected to hold well for systems that do not exhibit strong static correlation, where furthermore we may expect a degree of cancellation of error when only looking at the differences of the ground and excited-state total energies. To investigate the validity of the XDFT formalism without having to assume equivalence between the Kohn-Sham and the interacting many-particle wave function, in the following, we seek to establish a rigorous connection between the formally exact theorems of excited state DFT and the XDFT method. We start with a very brief discussion of these exact theorems.
Exact theorems for excited state DFT. Ground-state density-functional theory involves finding the ground-state (GS) density of a system of interacting electrons through the number-conserving minimisation of the Levy constrained search 46 LL ext . Such minimisation is with respect to the density, ρ r ( ), and for a given external potential Here, T is the electron kinetic energy operator, V ee is the electron-electron interaction operator, and ρ F [ ] LL is the minimum of +T V ee provided that the state ψ produces the density ρ r ( ). Now, for some external potential v r ( ) ext , we may denote the k th excited state, energy and density by ψ k , E k and ρ r ( ) k , respectively. Furthermore, we may assume that, for the same or a different external potential ′ v r ( ) ext , there exists a stationary state ψ′ producing the same density ρ r ( ) k , but such that Then, necessarily, k k e e k LL www.nature.com/scientificreports www.nature.com/scientificreports/ and consequently, This leads us to an important result, shown by Perdew and Levy in ref. 6 , namely that a stationary excited state will correspond to a local density minimum of LL ext if and only if, for any external potential, there is no other stationary state that gives the same density and yields a lower value for +T V ee .
In an alternate approach to this 7,47 , the first excited-state energy of a system subject to an external potential v r ( ) ext can be found by minimising, with respect to density ρ r ( ), the functional 1 e xt ext 1 ext 1 e xt is a bifunctional defined as where the minimisation of +T V ee is to be performed over states ψ , which yield the density ρ r ( ) and are orthonormal to the ground state, ψ 0 , of the external potential, v r ( ) ext . Foundational justifications for numerous time-independent excited-state DFT schemes 48,49 including ∆SCF 28 , TOCIA 50 , and OCDFT 51 have been developed on the basis of this approach.
Remarkably, if the external potential is Coulombic, namely if, for a number of M nuclei, it has the form where the α th nucleus with charge α Z is located at α R , then two different stationary states in the presence of the same or different external potential cannot have the same density 52 . In other words, the stationary-state density uniquely specifies the stationary state and the external potential. In this situation, as shown in ref. 52 , one can, in principle, omit the v ext dependence of F 1 and find the first excited state density ρ r ( ) 1 by minimising 1 e xt ext 1 XDFT may be placed in a formal context using this equation, as we now explain.
The XDFT approximation. Here we show how the XDFT method follows from Eq. (12) with the use of certain well-defined approximations. Let v r ( ) KS be the KS potential and ψ 0 KS be the KS ground state corresponding to an interacting system with an external potential v r ( ) ext . Then, assuming representability where required, consider a non-interacting KS like system whose ground state density equals ρ 1 , which is that of the first excited state of the interacting system. This system is subject to a local potential, such that the exchange-correlation potential v r ( ) XC1 ensures that the correct density is recovered. Unfortunately, v r ( ) XC1 is not known to us. To facilitate the use of available approximations for exchange-correlation, therefore, let us consider a different non-interacting KS-like system such that its lowest-energy stationary state ψ 0 KS which satisfies the condition For such a non-interacting system, the XC-potential v r ( ) XC , which generates a KS potential , contrasts with the standard ground-state XC potential, v r ( ) XC , which is constructed in such way that the non-interacting ground-state density associated to coincides with the density of the interacting ground state for the external potential, v r ( ) ext . We note, using Eq. (12), that v r ( ) XC must be a unique functional of the density. We also note that, in contrast with the corresponding object in the treatment presented in ref. 53 , ψ 0 KS is not generally the first excited state of the original, unconstrained KS system since ψ 0 KS and ψ 0 KS are stationary states of KS Hamiltonians with different potentials v KS and v KS .
Using XDFT, we seek to obtain the non-interacting state ψ 0 KS that is orthonormal to the non-interacting ground state ψ 0 KS , by ejecting a single electron out of the valence subspace, i.e., by creating an electron-hole pair in the KS system. Formally speaking, therefore, the XDFT method effectively amounts to making the central approximation www.nature.com/scientificreports www.nature.com/scientificreports/ used in OCDFT 51 . We refer the reader in particular to the very informative Table 3 of ref. 51 , where several properties of various time-independent excited-state DFT approaches are carefully compared. Using the same notation as is used in that Table, the properties of the XDFT method are tabulated in Table 1.
Proof that XDFT encodes the orthonormality of KS Slater determinants. In the following, we prove that the the density generated by XDFT belongs to a non-interacting state that is orthonormal to ψ 0 KS . For an N electron system, let the excited-state valence orbitals obtained with XDFT be ψ , . Then, since a unitary transformation of orbitals preserves the density, it is sufficient to prove that there is at least one unitary transformation of ψ { } i that produces an "excited electron" orbital, i.e., an orbital that is orthonormal to all of the ground-state valence orbitals that generate ψ 0 KS . Let us define a candidate orbital ψ ⁎ , with a view to describing the excited-state component of the XDFT KS valence eigensystem, as Then, noting that, for zero-temperature fermionic systems the density matrix is idempotent, and so   Table 1. Using the same notation as is introduced in Table 3 of ref. 51 , the properties of XDFT from the leftmost column report (i) orthonormality with respect to the ground state; (ii) the number of electrons transferred from the occupied to the virtual orbitals; (iii) the orbital space in which the hole orbital belongs; (iv) the orbital space in which the occupied orbitals belong, and (v) whether the method may suffer from variational collapse to the ground state. (  and therefore ψ′ N is orthonormal to all of the ground-state valence orbitals. Hence, the Slater determinant constructed from the set of orbitals ψ′ { } i is necessarily orthonormal to the determinant ψ 0 KS .

Methodological Details
Implementation and parameters. The linear-scaling first-principles code ONETEP 54 , within which we have implemented the XDFT formalism, variationally optimizes a minimal set of localized, non-orthonormal generalized Wannier Functions (NGWF), expanded in terms of psinc functions 55,56 , to minimize the total energy. ONETEP is equipped with an automated conjugate-gradients method for optimizing the cDFT (or XDFT) Lagrange multiplier 45,57,58 . We have used this, together with the Perdew-Burke-Ernzerhof (PBE) XC functional 59 to calculate the lowest singlet excitation energies of the 28 closed-shell organic molecules contained in the wellknown Thiel set 60 . 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 61 was used with a parameter of 7.0 a 0 . The constrained KS system was found to contain symmetry-protected partial-filling of the degenerate highest occupied state in certain molecules, and so we used finite-temperature ensemble DFT as implemented in ONETEP 62 in all cases. In order to access one of these degenerate Sds that make up the singlet, in practice, we promote one electron by applying the XDFT constraint to one spin channel only, maintaining = m 0 s . We note that, for a spin-restricted treatment, since the two Sds are degenerate, each of them has the same energy as the singlet and the triplet state and Eq. (21) is trivially satisfied. However, at this point, we make a final assumption that Eq. (21) may be used to approximately evaluate the energy of the interacting system. Keeping in mind that the three triplet states for = − m 1,0,1 s are degenerate, the energy of the singlet first-excited state can then be approximated as  (22) is that it involves only the energies of two single-Sd states that are available using XDFT. Each term on the right hand side of Eq. (22) derives from an interacting system that is obtained from an equivalent unrestricted KS system having the same density and spin density. In passing, we note that the Sd state is sometimes referred to in the literature as a contaminated singlet. A formalism using restricted (spin-independent) KS orbitals might offer an energy E Sd that is more appropriate for use with Eq. 22.

Multiplet sum method.
Calculation flowchart. The XDFT flowchart involving coupled calculations for finding neutral gaps of finite systems is presented in Fig. 1. In practice, the task of determining the triplet and singlet single excitation energies boils down to calculating total energy differences.
1. First, a standard DFT run is performed to calculate the closed-shell ground-state energy E 0 and density operator ρˆ0. 2. ρˆ0 is used to run an XDFT calculation confining

Results and Discussion
Excitation energy of Thiel set molecules. We have used XDFT to calculate the lowest singlet excitation energies of the 28 closed-shell organic molecules contained in the well-known Thiel set 60 . In Fig. 2, we show a scatter plot of the singlet and triplet excitation energies calculated with XDFT against those obtained with LR-TDDFT and the adiabatic Perdew-Burke-Ernzerhof (PBE) XC functional 59 in ref. 65 (singlets) and ref. 66 (triplets). The LR-TDDFT results are broadly in agreement with experimental values (see the supporting information in ref. 67 ). The triplet energies, which do not rely on any additional approximation beyond XDFT such as multiplet sum, show almost perfect agreement with those calculated using ground-state DFT with = m 1 s . This may provide a practical way of validating the XDFT approximations when applied to a new system. The figure also demonstrates that, in spite of the multiplet sum approximation, XDFT yields singlet energies with a remarkably good accuracy, if adiabatic, semi-local LR-TDDFT is taken as a reasonable benchmark. In terms of computational efficiency, we note that, for the representative molecule p-Benzoquinone, an XDFT calculation run on 3  www.nature.com/scientificreports www.nature.com/scientificreports/ processors for the singlet excitation energy offers an approximate two-fold reduction in computation time compared to its LR-TDDFT counterpart. Unlike linear-scaling LR-TDDFT, XDFT does not require the prior optimization of a defined number of conduction-band orbitals, which is a process that can demand some trial-and-error before well-converged results are obtained. It is to be noted that, as a result of the charge-delocalization error of semi-local XC functionals, XDFT, in its currently-implemented form, is applicable only to finite systems.
Charge difference density. The difference between the local part of the constrained density operator, ρˆ, and that of the ground state density operator, ρˆ0, can be viewed as an approximation to the difference density, from which transition dipole moments for example can be calculated. In Fig. 3 we show such plots for a representative molecule of the Thiel set 60 , propanamide. Figure 3(a) shows an approximation to the difference density based on the ground-state KS orbitals, which neglects orbital relaxation and electron-hole 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 33,64 applied to the total electron densities) reflect a greater degree of difference density localisation than (a). Due to Pauli exclusion, furthermore, the singlet (b) difference density attains a greater spatial localisation than the triplet one (c). Such conclusions are also supported by quantitative analysis of the charge densities.
choice of Xc functional. The accuracy of results obtained with any DFT-based method is necessarily dependent on the approximate XC-functional used. We have found good agreement between XDFT and LR-TDDFT results for the Thiel set, using the semi-local PBE functional for both methods. This choice of functional is motivated by the fact that, for many DFT and adiabatic-kernel LR-TDDFT calculations, when compared to experimental results, PBE typically offers acceptable accuracy with relatively inexpensive calculations and is therefore a highly popular choice of functional.
However, non-local hybrid XC-functionals containing a fraction of KS exact-exchange typically improves the agreement with experimental results, at the expense of increasing the computational cost. This trend can be seen in Table 2, where we compare experimental singlet excitation energies with XDFT results evaluated using PBE and the B3LYP 74 hybrid functional for six Thiel set molecules. These molecules are those for which the agreement between the XDFT(PBE) and experimental results is particularly poor (giving a Mean Absolute Error or MAE of 0.85 eV, whereas the MAE is 0.50 eV for the entire Thiel set). The MAE between XDFT and experiment for the six molecules reduces to 0.33 eV using XDFT(B3LYP). It is more probable that this improvement is primarily due to an improved description of the fundamental gap and electron-hole binding, rather than an improvement in the performance of the XDFT or multiplet sum method approximations per se. It is to be noted that, although a hybrid functional improves the energies considerably, the discrepancies resulting from the multiplet sum approximation are present nonetheless. This can be seen by contrasting the MAE (w.r.t. experimental values) of the XDFT for the propanamide molecule, calculated with an isosurface value of ±0.05 eÅ −3 . Panel (a) shows the charge-density difference between the ground-state KS LUMO and HOMO orbitals, while (b,c) show, respectively, the singlet and triplet difference densities generated as the difference between the excited XDFT (ρˆ) and ground-state DFT (ρˆ0) densities. The change from (a) to (b,c) is due to electron-hole pair binding and screening, which are captured variationally and at all orders by XDFT. Figure  created using XCrysDen v1.6.2 68  www.nature.com/scientificreports www.nature.com/scientificreports/ results (0.33 eV) with the LR-TDDFT ones (0.08 eV). On the other hand, the fact that, for triplet excitations, the MAE of the XDFT energies is 0.08 eV indicates that the error in the singlet energies arises from the multiplet sum approximation and not from the XDFT approximation per se. It seems feasible that this agreement with experiment may be improved further within the XDFT framework, potentially using range-separated hybrid functionals, implicit dielectric screening, zero-point phonon corrections, and a more elaborate treatment of spin contamination that circumvents the multiplet sum method.
Simulation of a double excitation. Finally, we explore the ability of XDFT to calculate energies of excitations with strong double (two-electron) character, which are inaccessible by constuction to adiabatic-kernel LR-TDDFT [16][17][18] . The XDFT method is non-linear, unlike LR-TDDFT, in the sense that its Hartree and XC potentials are calculated self-consistently with the density in the excited state, i.e., not just corrected to first order using the interaction kernel f Hxc . Thus, XDFT is not limited to single excitations. As a proof of principle, we focus here on one particular excitation of a known, strong double (i.e. two-electron) character, nothing that a more comprehensive study of excitations of more mixed single-double character using XDFT would be necessary to fully establish the range of applicability of the method to such effects. We note, in passing, that XDFT is not restricted to exciting integer numbers of electrons, M, particularly when coupled with ensemble DFT. In the benchmark case of atomic beryllium, the first double excitation promotes two electrons from the 2s to the 2p orbitals 75 . For the Be atom, 1s electrons are described by a pseudopotential, rendering the multiplet sum method unnecessary. Consequently, using the ground state valence density operator ρˆ0 to confine zero electrons to the total valence subspace from a standard ground-state DFT calculation, we can directly access the energies of the lowest lying doubly excited singlet 0 , respectively. Our results 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 LR-TDDFT(PBE) result, indicating that the multiplet sum method is accurately applicable to this system. We note, however, that while our singlet double excitation energies agree surprisingly well with experimental values, this is much less the case for our triplet double energies. Experimentally, the singlet → s p 2 2 2 2 excitation is