Excitation energies from thermally assisted-occupation density functional theory: Theory and computational implementation

assisted-occupation density functional theory: Theory and computational

implementation

Cite as: J. Chem. Phys. 153, 084120 (2020); https://doi.org/10.1063/1.5140243

Submitted: 26 November 2019 . Accepted: 12 August 2020 . Published Online: 31 August 2020 Shu-Hao Yeh, Aaditya Manjanath , Yuan-Chung Cheng , Jeng-Da Chai , and Chao-Ping Hsu

Excitation energies from thermally

assisted-occupation density functional theory:

Theory and computational implementation

Cite as: J. Chem. Phys. 153, 084120 (2020);doi: 10.1063/1.5140243 Submitted: 26 November 2019 • Accepted: 12 August 2020 • Published Online: 31 August 2020

Shu-Hao Yeh,^1,2 Aaditya Manjanath,¹ Yuan-Chung Cheng,² Jeng-Da Chai,^3,4,a) and Chao-Ping Hsu^1,a) AFFILIATIONS

1Institute of Chemistry, Academia Sinica, Taipei 11529, Taiwan

2Department of Chemistry, National Taiwan University, Taipei 10617, Taiwan

3Department of Physics, National Taiwan University, Taipei 10617, Taiwan

4Center for Theoretical Physics and Center for Quantum Science and Engineering, National Taiwan University, Taipei 10617, Taiwan

a)Authors to whom correspondence should be addressed:[email protected]and[email protected]

ABSTRACT

The time-dependent density functional theory (TDDFT) has been broadly used to investigate the excited-state properties of various molecular systems. However, the current TDDFT heavily relies on outcomes from the corresponding ground-state DFT cal- culations, which may be prone to errors due to the lack of proper treatment in the non-dynamical correlation effects. Recently, thermally assisted-occupation DFT (TAO-DFT) [J.-D. Chai, J. Chem. Phys. 136, 154104 (2012)], a DFT with fractional orbital occupations, was proposed, explicitly incorporating the non-dynamical correlation effects in the ground-state calculations with low computational complexity. In this work, we develop TDTAO-DFT, which is a TD, linear-response theory for excited states within the framework of TAO-DFT. With tests on the excited states of H2, the first triplet excited state (1³Σ⁺u) was described well, with non- imaginary excitation energies. TDTAO-DFT also yields zero singlet–triplet gap in the dissociation limit for the ground singlet (1¹Σ⁺g) and the first triplet state (1³Σ⁺u). In addition, as compared to traditional TDDFT, the overall excited-state potential energy surfaces obtained from TDTAO-DFT are generally improved and better agree with results from the equation-of-motion coupled-cluster singles and doubles.

Published under license by AIP Publishing.https://doi.org/10.1063/1.5140243., s

I. INTRODUCTION

Over the past decades, Kohn–Sham density functional the- ory (KS-DFT)^1,2 has been extensively used in the prediction of various ground-state (GS) properties of solids as well as finite- sized molecules.^3–5Its time-dependent (TD) extension, known as time-dependent density functional theory (TDDFT),^6–8has been a popular approach for computing excited-state properties, includ- ing the absorption and emission spectra,⁹ photochemical reac- tions,¹⁰ dynamics,¹¹ and energy and electron transfer,¹² due to its low computational cost and the availability of a plethora of computer codes in this area. The one-to-one correspondence between the TD density and the TD external potential was

rigorously demonstrated by Runge and Gross in 1984 in their the- orem.⁶ The linear-response framework was further introduced,^7,8 which brought forth a paradigm shift in the simulation of exci- tations of quantum systems from a density-functional perspec- tive^13–15 and is the main reason behind the popularity of this method.

However, conventional TDDFT is derived from ground-state (GS) KS-DFT, which is a single-determinant-based method. As a result, it can fail to describe the excited-state phenomena heav- ily governed by non-dynamical (or static) correlation, such as photochemistry processes involving photoinduced bond breaking, and problems associated with conical intersection.^7,9,16,17 A pro- totypical example is the bond dissociation process of the H2

molecule. It is known that the excitation energy of the low- est triplet state of H2, computed using conventional TDDFT,⁹ would become imaginary beyond a H–H bond distance of 1.75 Å, a phenomenon arising from a spin symmetry-breaking solution in the ground state,^18,19 a typical characteristic of nondynami- cal correlation effects. In contrast, in wavefunction-based meth- ods, the (nearly) degenerate determinants are considered on an equal footing when performing a self-consistent field (SCF) calcu- lation, and this is the basis of multi-configuration (MC) SCF or complete active space (CAS) SCF-based methodologies. However, these methods can be prohibitively expensive for large systems as their computational cost scales factorially with the size of active space.

KS-DFT with proper exchange energy functionals may rea- sonably model systems with non-dynamical correlation, albeit at the expense of enormous computation efforts. For example, the works by Becke^20,21 and the works by Kong and co-workers^22–24 demonstrated parametric functionals that need to be solved self- consistently within the single-determinant framework. Although these works significantly improved the bond dissociation trends of simple diatomic molecules, compared to the Hartree–Fock the- ory, they still deviate appreciably at the bond dissociation limit compared to a full configuration interaction (FCI) calculation.^22,23 Moreover, the SCF associated with these functionals adds to the computational effort that can scale dramatically with the size of molecules.

On the other hand, various approaches have been developed to cope with the non-dynamical correlation effects without the high computational cost of an exact exchange functional. The CAS- DFT model is one such method,²⁵wherein some amount of cor- relation has been accounted for by a density functional calcula- tion. As a result, the dynamical correlation associated with the MC representation of the system might be “double counted.”^26,27 To mitigate this issue, the multi-configuration pair-density func- tional theory^26,27and multi-configuration range-separated DFT^28,29 were developed. While the former utilizes the so-called on-top pair-density functional, the latter separates the electron interac- tion operator into short- and long-range parts, which are treated with DFT and wavefunction theory, respectively. Although the idea of using such a “hybrid” scheme seems to be an attrac- tive prospect,^26–29 they can be computationally demanding for increasing system sizes because of the initial generation of MC wavefunctions.

Another category of computational methods exists, which can cope with non-dynamical correlation with the additional advan- tage that they are low-cost methods. They include the spin–

flip, ionization-potential, and electron-affinity based approaches, which are aimed to start with a high-spin, with 1-less or 1- more electron single-determinant references such that the non- dynamical correlation problem is minimal.^30–32These approaches require a well-balanced treatment of the orbitals in the refer- ence, and they can offer high-quality solutions in many cases.

However, the requirement of balanced treatment of orbitals in the reference is not always feasible, and thus, applications are limited.

In this regard, the thermally assisted-occupation density functional theory (TAO-DFT)³³was developed by Chai in 2012 to alleviate the formidable challenge of balancing the computational

cost and simultaneously incorporating the non-dynamical correla- tion effects with reasonable accuracy. In contrast to traditional KS- DFT, the underlying principle of TAO-DFT is in the usage of frac- tional orbital occupations according to a given fictitious temperature (θ) to effectively incorporate the different electronic configurations of a system. This approach ensures that some “excitations” in the form of fractional populations of electrons in the low-lying virtual orbitals are considered along with the GS of the system, similar to a multi-determinant expansion of the wavefunction. The inclusion of fractional occupancies is a computationally cheaper alternative to a multi-determinant expansion for accounting non-dynamical corre- lation effects. As a result, TAO-DFT has a computational cost similar to that of KS-DFT, which is O(N^3–4). In TAO-DFT, the entropy contribution [e.g., see Eq. (26) of Ref.33] can reasonably capture the non-dynamical correlation energy of a system, which was dis- cussed and numerically investigated in Ref.33, even when the sim- plest local density approximation (LDA) exchange-correlation (XC) energy functional is used. The XC energy functionals at the higher rungs of Jacob’s ladder, such as the generalized-gradient approxima- tion (GGA),³⁴global hybrid,³⁵and range-separated hybrid^35,36XC energy functionals, can also be employed in TAO-DFT. Moreover, a self-consistent scheme that determines the fictitious temperature in TAO-DFT has been recently proposed to improve the perfor- mance of TAO-DFT for a wide range of applications.³⁷Since TAO- DFT is similar to KS-DFT in computational efficiency, TAO-DFT has been recently adopted for the study of the electronic proper- ties of various nanosystems with pronounced radical nature.^36,38–45 In particular, the electronic properties (e.g., singlet–triplet energy gaps, vertical ionization potentials, vertical electron affinities, fun- damental gaps, and active orbital occupation numbers) of linear acenes and zigzag graphene nanoribbons (i.e., systems with polyrad- ical character) obtained from TAO-DFT^33–35,38have been shown to be in reasonably good agreement with those obtained from other accurate electronic structure methods, such as the particle–particle random-phase approximation (pp-RPA)⁴⁶ XC energy functional in KS-DFT, the density matrix renormalization group (DMRG) algorithm,^47,48the variational two-electron reduced density matrix (2-RDM) method,^49,50and other high-level methods.^51–54

II. GROUND-STATE REFERENCE: TAO-DFT

In TAO-DFT,³³the electron density is represented by the ther- mal equilibrium density of an auxiliary system ofNenon-interacting electrons at a fictitious temperature θ (in energy units),

ρ(r) = ∑

i

fiϕ^∗i(r)ϕi(r). (1)

Here,fi(a value between 0 and 1) is the fractional occupation num- ber of theith orbital ϕiand is given by the Fermi–Dirac distribution function,

fi= {1 + exp[(εi−μ)/θ]}⁻¹, (2) where μ is the chemical potential for electrons and is determined by

∑ifi=Nefor a given θ, orbital energies {εi}, and total electron num- berNe. This choice for the fractional occupation function and the corresponding one-particle density matrix has been extensively used in other methods such as finite-temperature DFT (FT-DFT)⁵⁵and

floating occupation molecular orbital-complete active space config- uration interaction (FOMO-CASCI).⁵⁶ With thisassisted occupa- tion number and generalized density expression, the total ground- state energy functional can be written as

EG[ρ] = T_TAO[{f_i,ϕi}]+Vext[ρ] + E^KSHxc+Eθ[ρ], (3) whereTTAOis the kinetic free energy functional of non-interacting electrons [equivalent toA^θsas defined in Eq. (24) of Ref.33],Vext[ρ]

is the energy functional of the external potential (or nuclei potential), E^KSHxcis the sum of Hartree and XC energy functionals in KS-DFT, andE_θis the θ-dependent energy functional.³³Alternatively (to the original derivation³³), from Eq.(3), upon performing the functional derivatives with respect to the orbitals (ϕi), we can also obtain the SCF equations in TAO-DFT,

[−1

2∇²_r+ vext(r)+ vHxc^KS(r)+ vθ(r)]ϕi(r) =εiϕi(r), (4) where vext, v_Hxc^KS, and vθare the potentials (or functional derivatives) of corresponding energy functionals (i.e.,Vext[ρ], E^KSHxc, andEθ[ρ], respectively) in Eq.(3), and {ϕi} and {εi} are the TAO orbitals and orbital energies, respectively, which can be solved self-consistently through SCF. The algorithm is similar to KS-DFT, with the only differences being the vθ(r) term in the Hamiltonian and the determi- nation of chemical potential μ, making this approach attractive and easy in implementation. We have provided a variational perspective of TAO-DFT inAppendix A, which complements the derivation in Ref.33.

III. EXCITED-STATE THEORY: TDTAO-DFT A. Mathematical formalism

In the present work, we propose TDTAO-DFT, which is a time-dependent linear-response theory for TAO-DFT, allowing excitation energy calculation using Casida’s formulation,⁸ within the framework of TAO-DFT. In TDTAO-DFT, the TD density is given by

ρ(r, t) = ∑

p

fpϕ^∗p(r,t)ϕp(r,t), (5)

where ϕp(r, t) are the TD orbitals (for the fictitious particles) and fp are the corresponding fractional occupation numbers, which are assumed to be time-independent, and their values are taken from those obtained from the corresponding ground-state TAO-DFT calculation [Eq. (1)]. In order to facilitate the map- ping between the original interacting system of electrons mov- ing under the influence of a TD external potential and the auxiliary system of non-interacting particles, an action varia- tional principle in TAO-DFT should be established. Following the variational principle, the TD effective potential for the non- interacting TAO system can be partitioned into the following parts:

v_eff^TAO(r,t) = vext(r,t) + v^TAOHxcθ[ρ](r, t), (6)

where v^TAO_Hxcθ is the functional derivative of the Hxcθ-action, which contains the time-dependent Hartree potential, exchange- correlation potential, and θ potentials for the fractional occupation.

Further details are included in Appendix B1 accompanying this work.

Similar to conventional TDDFT, with the equality connecting the effective potential and the functional derivative of TD action, an equation of motion for TDTAO-DFT can be expressed as

i∂

∂tϕp(r,t) = [−1

2∇²_r+ vext(r,t) + v^TAO_Hxcθ[ρ](r, t)]ϕp(r,t)

= ˆF(t) ϕp(r,t). (7)

We note that v^TAO_Hxcθ[ρ](r, t) is also a TD generalization of the potential associated with the Hartree, exchange, correlation, and θ- functionals in GS TAO-DFT. The equation of motion is reformu- lated in terms of the one-particle density matrix P(t),⁹

i∂

∂tP(t) = [F(t), P(t)], (8)

where F(t), the time-dependent “Fock matrix,” is the matrix rep- resentation of the one-particle operator (ˆF) in Eq.(7). The general time-evolution of the state of a system is given by

P(t) = P^○+ δP(t) (9)

and

F(t) = F^○+ δVext(t) + δFHxcθ[P](t), (10) where P^○ and F^○ denote the initial conditions for solving Eq.(8) and δP(t), δVext(t), and F_Hxcθ[P](t) are the time-dependent changes in the matrices of density, external field, and Hartree-exchange- correlation potential in matrix representation interaction, respec- tively, in the system. The initial state (att = t0) is commonly con- sidered to be the unperturbed GS of the system for convenience. In terms of the GS TAO orbitals,

P^○_pq=δpq⋅fp,F^○_pq=δpq⋅εp. (11) If the electronic eigenspectrum of a system is desired, the ampli- tude of the change in the external field |δVext(t)| is assumed to be infinitesimally small.^6–9It is therefore suitable to consider a linear- response relation between δF_Hxcθ[P](t) and δP(t). Using the GS TAO orbital basis, this can be obtained as

δF^Hxcθrs (t) = ∑

pq

∫ dτ(δFrs^Hxcθ(t)

δPpq(τ) )δPpq(τ). (12) Employing the time-domain Fourier transformations

δPqr(ω) =∫ dt e^−iωt[δPqr(t)], (13)

δVqr(ω) =∫ dt e^−iωt[δVqr(t)], (14)

δF^Hxcθrs

δPpq

(ω) =∫ dt e^{−iω(t−τ)}[δF^Hxcθ_rs (t)

δPpq(τ) ], (15) one could recast Eq.(8)into

∑

q

⎡⎢

⎢

⎢⎢

⎣

F^○pq⋅δPqr(ω) − δPpq(ω) ⋅ F^○qr

+⎛

⎝

δVpq(ω) + ∑

st

⎛

⎝ δF^Hxcθpq

δPst

(ω)⎞

⎠

⋅δPst(ω)⎞

⎠ P^○_qr

−Ppq^○

⎛

⎝

δVqr(ω) + ∑

st

⎛

⎝ δF^Hxcθqr

δPst

(ω)⎞

⎠

⋅δPst(ω)⎞

⎠

⎤⎥

⎥⎥

⎥⎦

=ω ⋅ δPpr(ω) (16)

by neglecting all second-order (or higher) terms. Upon invok- ing the GS definitions in Eq. (11) and assuming all δVpq(ω) to be infinitesimally small, the corresponding working equation becomes

(εp−εr)δPpr(ω) − ( fp−fr)

⎡⎢

⎢⎢

⎢⎣

∑

st

⎛

⎝ δF_pr^Hxcθ

δPst (ω)⎞

⎠ δPst(ω)

⎤⎥

⎥⎥

⎥⎦

=ω ⋅ δPpr(ω). (17)

A conventional linear-response relation [which is the inverse of Eq.(12)]^7,57gives the TD density–density response function. The details of this derivation are provided inAppendix B2.

Similar to conventional TDDFT, we apply the adiabatic approximation to the xcθ-kernel (i.e., the xcθ-kernel is assumed to be frequency-independent),^8,9,58

δF^Hxcθpr

δPst

(ω) ≈ δF^Hxcθpr

δPst

RR RR RR RR RR RP^○

≈ ∫ d³r d³r^′ϕ^∗r(r)ϕp(r)𝕗Hxcθ(r, r^′)ϕt(r^′)ϕ^∗s(r^′)

= (rp∣𝕗Hxcθ∣ts). (18)

The working equation would be reduced to an eigenvalue equation,

∑

st

[(εp−εr) ⋅δps,st− (fp−fr)(rp∣𝕗Hxcθ∣ts)] ⋅ Ω^Rk,st=ωk⋅Ω^Rk,pr, (19)

where Ω^Rpr = δPpr and k denotes the kth eigenvalue. This can be represented in the matrix form as Casida’s equation,⁸

⎛

⎝ Aˆ Bˆ Bˆ^∗ Aˆ^∗

⎞

⎠ (

X Y) =ωk(

ˆI ˆ0 ˆ0 −ˆI)(

Xk

Yk

), (20)

whereXk,pr=Ω^R_k,p>r,Yk,rp =Ω^R_k,p<rdenotes upward and downward transitions, respectively. The coupling matrices are defined as

Apr,st= (εp−εr)δpsδrt+Bpr,ts, (21) Bpr,st= −(fp−fr)(rp∣𝕗Hxcθ∣st). (22) These matrices are similar in form to those derived from con- ventional Casida’s equation, which most TDDFT works are based on.^9,59However, we consider the fractional occupation number dif- ference (Δf ) pre-factor in Eq.(20), which is equivalent to the original Casida’s equation in Ref.8. It is to be noted that the occupation num- bers are explicitly sourced from GS TAO. In Eq.(19), the superscript R in Ω^Rimplies that the eigenvectors obtained are the right eigenvec- tors. Using the density–density response function (Appendix B2), an eigenvalue-like equation that is complementary to that in Eq.(19) can be derived. The details are included inAppendix B3.

B. Idempotency in TDTAO-DFT

In KS theory, an idempotent one-electron density matrix (PP

= P)⁹ is derived from the single-determinant ansatz of the wave- function, so for anyfirst-order changes in the one-electron density matrix,

P^○⋅δP + δP ⋅ P^○−δP = 0, (23) which when represented in terms of KS orbitals becomes

(np+nq−1) ⋅ δPpq=0, (24) where {np} are the integer occupation numbers (either 0 or 1).

Within this particular condition, the conventional Casida’s scheme allows transitions between only occupied (ni= 1) and virtual (na= 0) orbitals. On the other hand, due to fractional occupation numbers, the one-electron density matrix in TAO-DFT violates this idempo- tency condition for nonvanishing θ. Therefore, a relaxed condition in terms of TAO orbitals is proposed as

(fp+fq−1) ⋅ δPpq∝θ, (25) where the KS limit of TDTAO-DFT is recovered for θ → 0. This condition implies that transitions withfp+fr tending to 1 would be dominant. These transitions require one of thep and q orbitals to bestrongly occupied, 1/2 ≤ fr ≤1, with the otherweakly occu- pied, 0 ≤ fr<1/2. More details on therelaxed idempotency condition for TDTAO-DFT can be found inAppendix Caccompanying this work.

IV. COMPUTATIONAL DETAILS

We implement this formalism in the development version of Q-Chem 5.2.⁶⁰ All numerical results are calculated with the cc- pVDZ basis set, which was determined by performing a comprehen- sive convergence test of different sets. The two-electron integrals are evaluated with the standard quadrature Euler–Maclaurin–Lebedev grid (50 194),⁶¹consisting of 50 Euler–Maclaurin⁶²radial grid points and 194 Lebedev⁶³angular grid points.

V. H2BOND DISSOCIATION USING TDTAO-DFT

We demonstrate how some of the challenges plaguing TDDFT are rectified with our method through the GS bond dissociation

process of the H2 molecule. This system has been studied exten- sively for many years using a plethora of methods. Successfully capturing the mechanism of bond dissociation within the frame- work of DFT has been elusive owing to the lack of incorporation of non-dynamical correlation effects. Within TAO-DFT, however, this challenge was resolved by choosing an appropriate θ of 40 mhartree.^33,34 It was further shown that, at the bond dissociation limit, the multi-reference character was more pronounced.^33,34

In TDDFT, one encounters the challenge of imaginary frequen- cies (i.e., excitation energies) for the triplet states that occurs in most of the results obtained from adiabatic local density approx- imation (ALDA) functionals (kernels).^18,19,64 This issue is related to the symmetry breaking where the difference in spin densities (i.e., ρα−ρβ) is not equal to zero for a large interatomic distance.

In other words, the unrestricted (asymmetric) solution obtained using KS-DFT becomes lower in total energy than the restricted (symmetric) one, as demonstrated by Casida et al. using a two- level model.¹⁹TAO-DFT significantly rectifies this issue for a large enough θ value.³³

Figure 1 shows the potential energy surface (PES) of the first triplet excited state (1³Σ⁺u) for H2 bond dissoci- ation using TDTAO-DFT and TDDFT (θ = 0 mhartree).

The TDDFT results show imaginary frequencies beyond the H–H bond distance of ∼1.5 Å. This is attributed to a poor ground- state reference, as mentioned previously, due to the lack of incor- poration of the non-dynamical correlation effects beyond this bond distance. In addition, this phenomenon is observed in TDTAO-DFT simulations for θ = 0 mhartree, 10 mhartree, 20 mhartree, and 30 mhartree. However, for θ ≥ 40 mhartree, the imaginary-frequency issue is resolved.

We also note here that the requirement for a real-value 1³Σ⁺u

excitation energy mandates a higher threshold value for θ than

FIG. 1. Potential energy surface of the first triplet excited state (1³Σ⁺_u) computed using TDDFT (θ = 0 mhartree) and TDTAO-DFT with the Perdew–Burke–Ernzerhof (PBE) XC-functional, cc-pVDZ basis set, and gradient expansion approximation (GEA) version of the E_θ functional³⁴for TAO calculations. The inset shows a zoomed-in view for the large bond-distance regime.

that obtained through a self-consistent scheme,³⁷which is around 15.5 mhartree. While a lower θ value is needed to describe the ground-state bond dissociation curves, our observation indicates that a higher θ value is needed for excitation properties and an opti- mal determination scheme for θ remains to be developed. One such direction is to include the excited-state information in the post-SCF variational scheme similar to that outlined in Eq. (9) of Ref.56.

Another advantageous aspect of TDTAO-DFT is that the energy of the first triplet excited (1³Σ⁺u) state in the dissociation limit correctly approaches the GS singlet energy.Figure 2shows the singlet–triplet (1¹Σ⁺g–1³Σ⁺u) vertical gap as a function of H–H bond dissociation computed using ground-state TAO-DFT, coupled–

cluster singles and doubles (CCSD), and TDTAO-DFT. To compute the 1¹Σ⁺g–1³Σ⁺u gap at the ground-state level (in order to mitigate the problem of imaginary frequencies in TDDFT), it was recom- mended to use the unrestricted ground-state SCF formalism for H2and other small molecules.^22,23,64However, this does not guar- antee the convergence of the energy of the 1³Σ⁺u state to that of the 1¹Σ⁺g state at the bond dissociation limit for H2for TAO-DFT (Fig. 2). This gap may violate thecovalent nature of the³Σ⁺u state, where the energies of covalent states 1³Σ⁺u and GS (1¹Σ⁺g) should be the same at the bond dissociation limit.¹⁸ In other words, at this limit, the electrons are located in the 1s orbitals of the corre- sponding atoms and are, therefore, isolated enough with respect to one another. This gap increases with θ due to the increase in the energy of 1³Σ⁺_uand a simultaneous decrease in the energy of 1¹Σ⁺_g

FIG. 2. The energy gap as a function of the H–H bond distance (in Å) between the singlet ground state (1¹Σ⁺_g) and the first excited triplet state (1³Σ⁺_u) calculated using TDTAO-DFT and unrestricted TAO-DFT with the PBE XC-functional and GEAθ- functional. The equation-of-motion–coupled cluster singles doubles (EOM-CCSD) results are presented as a benchmark. The cc-pVDZ basis set was employed for all calculations. The inset shows a zoomed-in view for the large bond-distance regime.

FIG. 3. Potential energies of (a) singlet and (b) triplet excited states, computed using TDTAO-DFT (with θ = 40 mhartree and GEA θ-functional), EOM-CCSD, and conventional TDDFT. (c) shows only the¹Σ⁺_gstates. (d) is a zoomed-in region showing the avoided crossing between two EOM-CCSD states that are not completely captured by either TDTAO-DFT or conventional TDDFT. The orange shaded regions in (a), (c), and (d) indicate portions of the EOM-CCSD curves that have double excitation character. DE in (a) signifies that the CCSD state is double excitation in nature. PBE is selected as the XC-functional for all DFT calculations, and cc-pVDZ is selected as the basis set for all calculations.

(this θ-dependent decrease is also observed for the total energy of 1³Σ⁺u calculated with TDTAO-DFT inFig. 1). On the other hand, the trend obtained for TDTAO-DFT (Fig. 2) is in excellent agree- ment with that obtained using the equation-of-motion coupled- cluster singles and doubles (EOM-CCSD) method or observed in experiments.⁶⁵ EOM-CCSD is used here as a benchmark method since it is equivalent to FCI for a two-electron system such as H2.

For the sake of completeness, we also computed the PESs of other excited states for H2. The lowest six singlet and triplet excited states in TDTAO-DFT and TDDFT are demonstrated with low- lying PESs from EOM-CCSD inFigs. 3(a) and 3(b). The overall feature of singlet and triplet states from TDTAO-DFT is in excel- lent agreement with the EOM-CCSD results, except for the charge- transfer state (1¹Σ⁺u) and the missing states with double excitation character [purple curve with unfilled squares and curves highlighted in orange with unfilled diamonds inFigs. 3(a)and3(b)]. We specu- late that the problem with the 1¹Σ⁺ustate could be due to the usage of

the simple adiabatic approximation to the xcθ-kernel^6,19,66,67as well as the time-independent occupation numbers in our formalism.^68–70 The missing CCSD double excited states also indicate the inability of TDTAO-DFT to capture the avoided crossing between the first two¹Σ⁺gexcited states [orange shaded regions as shown inFigs. 3(a), 3(c), and3(d)]. A more detailed investigation is certainly required for resolving these challenges.

VI. RELATIONSHIP BETWEEN θ AND IMAGINARY FREQUENCIES: A QUALITATIVE DESCRIPTION

We perform a detailed analysis of the PESs with different θ values to acquire more insight about the qualitative relationship between θ and the imaginary roots. Two molecular systems were chosen for this analysis, H2and N2, and their S–T gaps are shown in Fig. 4. The problem of imaginary frequencies is fixed with TDTAO- DFT for a suitable choice of θ, irrespective of the system under con- sideration, thereby indicating its versatility. However, we note that

FIG. 4. S–T gap of (a) H₂and (b) N₂with the bond distance and differentθ (in mhartree) values, calculated using TDTAO-DFT with the PBE XC-functional, cc-pVDZ basis set, and GEA version of the Eθfunctional.³⁴ The filled symbols indicate potential energy surfaces without any imaginary frequencies.

θ is a system-dependent quantity and a robust algorithm is needed to ascertain it. Based on the optimal choice of θ, we observe that the S–T gap vanishes at the bond dissociation limit for N2[Fig. 4(b)], similar to that in H2 [Fig. 4(a)]. This is also in agreement with experiments.⁷¹

VII. CONCLUDING REMARKS

In summary, a time-dependent linear-response theory for pre- dicting excited-state properties based on the TAO-DFT framework, TDTAO-DFT, is proposed. This theory takes advantage of TAO- DFT, where the spin-symmetry-breaking problem of orbitals in ground-state SCF is resolved. As a result, TDTAO-DFT provides a correct description of low-lying triplet excited states, without imag- inary energies, at the bond dissociation limit for a molecule. This was demonstrated through the dissociation curve of the hydrogen molecule, in which a reasonable lowest triplet state (1³Σ⁺u) is cap- tured by TDTAO-DFT, but is not so, for TDDFT. Additionally,

TAO-DFT (with a large fictitious temperature θ) may produce an erratic gap between the 1³Σ⁺u and ground states at the dissocia- tion limit, which is resolved by TDTAO-DFT. The PESs for higher excited states of stretched H2 are also improved significantly as compared to TDDFT.

SUPPLEMENTARY MATERIAL

Thesupplementary materialincludes additional results and the numerical data presented in this work.

AUTHORS’ CONTRIBUTIONS

S.-H.Y. and A.M. contributed equally to this work.

ACKNOWLEDGMENTS

C.-P.H. acknowledges support from Academia Sinica and the Investigator Award (Grant No. AS-IA-106-M01) and the Ministry of Science and Technology of Taiwan (Project No. 105-2113-M- 001-009-MY4). J.-D.C. acknowledges support from the Ministry of Science and Technology of Taiwan (Grant No. MOST107-2628-M- 002-005-MY3) and National Taiwan University (Grant No. NTU- CDP-105R7818). A.M. acknowledges additional financial support from the Academia Sinica Distinguished Postdoctoral Fellowship.

This work also benefited from discussions facilitated through the National Center for Theoretical Sciences, Taiwan.

APPENDIX A: A VARIATIONAL PERSPECTIVE OF TAO-DFT

In this the section, we briefly present the derivation of the TAO- DFT KS-like equations based on an alternative, variational principle.

The same variational approach is also employed in the derivation of the linear-response theory, which will be presented inAppendix B of this work.

According to the partition of energy functional,^33,34the func- tional derivative of the total energy functional can be expressed as

δE[ρ]

δϕi(r)=δTs^TAO

δϕi(r)+ δVext

δϕi(r)+δEHxcθ[ρ]

δϕi(r) , (A1) whereT^TAOs is the kinetic (free) energy functional, andVext+EHxcθ

is the energy associated with the effective potential. The explicit derivative of the kinetic (free) energy functional would be

δT_s^TAO δϕ^∗_j(r^′)

= δ

δϕ^∗_j(r^′) (∑

i

fi∫ dr ϕ^∗i(r) ˆt ϕi(r)

+ θ{∑

i

filnfi+ (1 −fi)ln(1 −fi)})

=fj⋅ ˆt ϕj(r^′)+ ∑

i

⎡⎢

⎢⎢

⎣ δfi

δϕ^∗_j(r^′)⋅ ∫ dr ϕ^∗i(r) ˆt ϕi(r)

⎤⎥

⎥⎥

⎦ + θ ∑

i

⎡⎢

⎢⎢

⎣

{lnfi−ln(1 −fi)} ⋅ δfi

δϕ^∗_j(r^′)

⎤⎥

⎥⎥

⎦

, (A2)

where ˆt = −∇²/2 and δϕj(r^′)/δϕ^∗_j(r^′) =0. Similarly, the derivatives of the energy term associated with external potential as well as the Hxcθ energy term are, respectively,

δVext

δϕ^∗_j(r^′) = ∫ dr^′′δρ(r^′′) δϕ^∗_j(r^′)

⋅ δ

δρ(r^′′)[∫ dr v^ext(r)ρ(r)]

= ∫ dr v^ext(r) δρ(r) δϕ^∗_j(r^′)

=fj⋅vext(r^′)ϕj(r^′) + ∑

i

⎛

⎝ δfi

δϕ^∗_j(r^′)⋅ ∫ dr vext(r)ϕ^∗i(r)ϕi(r)⎞

⎠ (A3)

and δEHxcθ[ρ]

δϕ^∗_j(r^′) = ∫ dr δρ(r) δϕ^∗_j(r^′)

⋅δEHxcθ[ρ]

δρ(r)

= ∫ dr δρ(r) δϕ^∗_j(r^′)

⋅v_Hxcθ[ρ](r)

=fj⋅vHxcθ(r^′)ϕj(r^′) + ∑

i

⎡⎢

⎢⎢

⎣ δfi

δϕ^∗_j(r^′)∫ drvHxcθ[ρ](r)ϕ^∗i(r)ϕi(r)

⎤⎥

⎥⎥

⎦

. (A4) Combining the three terms above, an explicit expression of the total energy functional is derived,

δE[ρ]

δϕ^∗_j(r^′)

= δTs^TAO

δϕ^∗_j(r^′)+ δVext

δϕ^∗_j(r^′)+δEHxcθ[ρ]

δϕ^∗_j(r^′)

=fj⋅ ˆt ϕj(r^′)+ ∑

i

⎡⎢

⎢⎢

⎣ δfi

δϕ^∗_j(r^′)⋅ ∫ dr ϕ^∗i(r) ˆt ϕi(r)

⎤⎥

⎥⎥

⎦ + θ ∑

i

⎡⎢

⎢⎢

⎣

(lnfi−ln(1 −fi)) ⋅ δfi

δϕ^∗_j(r^′)

⎤⎥

⎥⎥

⎦

+fj⋅vext(r^′)ϕj(r^′)

+ ∑

i

⎡⎢

⎢

⎢⎣ δfi

δϕ^∗_j(r^′)⋅ ∫ dr v^ext(r)ϕ^∗i(r)ϕi(r)

⎤⎥

⎥

⎥⎦

+fj⋅v_Hxcθ(r^′)ϕj(r^′)+ ∑

i

⎡⎢

⎢

⎢⎣ δfi

δϕ^∗_j(r^′)⋅ ∫ dr vHxcθ[ρ](r)ϕ^∗i(r)ϕi(r)

⎤⎥

⎥

⎥⎦

=fj⋅ [ˆt + vext(r^′)+ vHxcθ[ρ](r^′)]ϕj(r^′)+ ∑

i

δfi

δϕ^∗_j(r^′)

⋅ {θ ln( fi

1 −fi

)+∫ dr ϕ^∗_i(r)[ˆt + vext(r)+ vHxcθ[ρ](r)]ϕi(r)}. (A5) Enforcing the normalization conditions for both density and orbital functions, a Lagrangian is introduced,

L[ρ] = E[ρ] − ∑

ij

[λij∫ dr ϕ^∗i(r)ϕj(r) −δij] −μ[∫ dr ρ^TAO(r) −Ne], (A6) where {λij} and μ are Lagrange multipliers. Considering the functional derivative with respect to orbital functions

δL[ρ]

δϕ^∗_j(r^′)

=fj⋅ [ˆt + vext(r^′)+ vHxcθ[ρ](r^′)]ϕj(r^′)+ ∑

i

δfi

δϕ^∗_j(r^′)⋅ {θ ln( fi

1 −fi

)

+∫ dr ϕ^∗i(r)[ˆt + vext(r)+ vHxcθ[ρ](r)]ϕi(r)} − ∑

i

λji⋅ϕj(r^′) −μ ∑

i

δfi

δϕ^∗_j(r^′)

−μfj⋅ϕj(r^′)

=fj⋅ [ˆt + vext(r^′)+ vHxcθ[ρ](r^′)]ϕj(r^′)+ ∑

i

{εi−θ[1

θ(εi−μ)]} ⋅ δfi

δϕ^∗_j(r^′)− ∑

i

λji⋅ϕi(r^′) −μ ∑

i

δfi

δϕ^∗_j(r^′)

−μfj⋅ϕj(r^′)

=fj⋅ [ˆt + vext(r^′)+ vHxcθ[ρ](r^′)]ϕj(r^′)+ μ ∑

i

δfi

δϕ^∗_j(r^′)− ∑

i

λji⋅ϕi(r^′) −μ ∑

i

δfi

δϕ^∗_j(r^′)

−μfj⋅ϕj(r^′)

=fj⋅ [ˆt + vext(r^′)+ vHxcθ[ρ](r^′)]ϕj(r^′) − (λjj+ μfj) ⋅ϕj(r^′)+

i≠j

∑

i

λji⋅ϕi(r^′). (A7)

and using the variational condition δL[ρ]/δϕ^∗_j(r^′) =0, one obtains

[ˆt + vext(r^′)+ vHxcθ[ρ](r^′)]ϕj(r^′) = (fj⁻¹λjj+ μ) ⋅ ϕj(r^′)+fj⁻¹ i≠j

∑

i

λji⋅ϕi(r^′). (A8)

Note that the second term in Eq.(A7)indicates that it is necessary to introduce the entropy term θ[∑ifilnfi+ (1 −fi)ln(1 −fi)] to the kinetic functional in order to preserve the correct variational prop- erty such that the derivative terms arising from vextand vHxcθ[last terms in Eqs.(A3)and(A4)] are compensated.

With a canonical orbital assumption (because the orbitals are orthonormal to one another), the equation can be recast into an eigenvalue equation, similar to a KS-like equation,

ˆh^TAO[ρ](r)ϕi(r) =εi⋅ϕi(r), (A9)

where ˆh^TAO= ˆt + vext+ vHxcθand εi=f_i⁻¹λii+ μ.

APPENDIX B: DETAILED DERIVATION OF LR-TDTAO-DFT

1. Variational principle for TAO action functional and TD effective potential

Starting from the action variational principle⁷²and its modified form,⁷³we have the general definitions of action functionals for a physical system,

A[ρ] =∫

τ

0 dt ⟨Ψ(t)∣(∂

∂t− ˆH)∣Ψ(t)⟩, (B1)

B[ρ] = A[ρ] +∫ dt∫ dr vext(r,t)ρ(r, t), (B2)

δB[ρ]

δρ(r, t) =vext(r,t) + i⟨Ψ[ρ](τ)∣δΨ[ρ; τ]

δρ(r, t) ⟩, (B3) where Ψ[ρ; τ] represents the wavefunction at time t and τ denotes the upper bound of the time integral.

For a TDTAO system, the definition of universal action func- tionals can be written similarly, following that of the conventional TDDFT scheme,⁷³

δATAO[ρ]

δρ(r, t) =i⟨ΨTAO[ρ](τ)∣δΨTAO[ρ; τ]

δρ(r, t) ⟩, (B4)

BTAO[ρ] = ATAO[ρ] +∫

τ

0 dtdr veff(r,t)ρ(r, t). (B5) The TD effective potential for TAO can be expressed as

v_eff(r,t) =δBTAO[ρ]

δρ(r, t) +i⟨ΨTAO[ρ](τ)∣δΨTAO[ρ; τ]

δρ(r, t) ⟩. (B6) One can define the difference between the two functionals as

AHxcθ[ρ] = BTAO[ρ] − B[ρ], (B7) which is the TAO extension of Hartree-exchange-correlation func- tionals. Summarizing the equations above, similar to TDDFT, one can recast the effective potential in TAO as

v_Hxcθ(r,t) =δAHxcθ[ρ]

δρ(r, t)

=veff(r,t) + i⟨ΨTAO[ρ](τ)∣δΨTAO[ρ; τ]

δρ(r, t) ⟩

−vext(r,t) − i⟨Ψ[ρ](τ)∣δΨ[ρ; τ]

δρ(r, t) ⟩. (B8)

2. Density–density response function

Here, we show that the linear-response equation can also be constructedinversely,

δρ(r t) =∫ dr^′dt^′χ^TAO_s (rt, r^′t^′)δv^TAO_eff (r^′t^′), (B9) where

χs^TAO(rt, r^′t^′) ≡ δρ^TAO(r,t) δv_eff^TAO(r^′,t^′)

(B10) is the density–density response function for a non-interacting TAO system. With the density expression in terms of TD orbitals, one obtains

χ^TAOs (rt, r^′t^′) = δρ^TAO(r,t) δv^TAO_eff (r^′,t^′)

=

δ[∑pfpϕ^∗p(r,t)ϕp(r,t)]

δv_eff^TAO(r^′,t^′)

= ∑

p

fp

δ[ϕ^∗p(r,t)ϕp(r,t)]

δv^TAO_eff (r^′,t^′)

= ∑

p

fp[

δϕ^∗p(r,t)

δv_eff^TAO(r^′,t^′)ϕ^○p(r,t) + ϕ^○∗p (r,t) δϕp(r,t)

δv_eff^TAO(r^′,t^′)

], (B11)

where ϕ^○p(r,t) and its complex conjugate represent the evolution of the TD orbitals in the absence of any TD perturbation (i.e., the TD external field). Applying the first-order perturbation theory, the TD orbital functions in a TD external field can be described by the equation

ϕp(r,t) = [1 − i∫

t 0 dt^′

r≠p

∑

rs

ϕ^○r(r)e^{−iεr(t−t′)}

× [∫ dr^′ϕ^○∗r (r^′)δv^TAOeff (r^′,t^′)ϕ^○s(r^′)]e^−iεst′

× ∫ dr ϕ^○∗s (r)]ϕ^○_p(r)

=ϕ^○p(r) −i e^−iεrt∫

t 0dt^′ ∑

r, r≠p

ϕ^○r(r)

× [∫ dr^′ϕ^○∗r (r^′)δv^TAOeff (r^′,t^′)ϕ^○p(r^′)]e^{−i(εp−εr)t′}, (B12) and the corresponding orbital response functions are expressed explicitly in terms of initial orbitals (ground-state TAO orbitals) and orbital energies,

δϕp(r,t) δv^TAO_eff (r^′,t^′)

= −iΘ(t − t^′)e^−iεrt ∑

r, r≠p

ϕ^○r(r)ϕ^○∗r (r^′)ϕ^○p(r^′)e^{−i(εp−εr)t′},

δϕ^∗p(r,t) δv^TAO_eff (r^′,t^′)

=iΘ(t − t^′)e^iεrt ∑

r, r≠p

ϕ^○∗r (r)ϕ^○r(r^′)ϕ^○∗p (r^′)e^{i(εp−εr)t′}.

(B13)

Combining Eqs.(B11)and(B13), the time-domain non-interacting response function in TDTAO-DFT can be evaluated as follows:

χs^TAO(rt, r^′t^′) = ∑

p

fp

⎧⎪

⎪

⎨

⎪⎪

⎩

⎡⎢

⎢⎢

⎢⎣

iΘ(t − t^′)e^iεrt ∑

r, r≠p

ϕ^○∗r (r)ϕ^○r(r^′)ϕ^○∗p (r^′)e^{i(εp−εr)t′}

⎤⎥

⎥⎥

⎥⎦ ϕ^○p(r,t)

+ ϕ^○∗p (r,t)

⎡⎢

⎢⎢

⎢⎣

−iΘ(t − t^′)e^−iεrt ∑

r, r≠p

ϕ^○r(r)ϕ^○∗r (r^′)ϕ^○p(r^′)e^{−i(εp−εr)t′}

⎤⎥

⎥⎥

⎥⎦

⎫⎪

⎪

⎬

⎪⎪

⎭

=i

r≠p

∑

pr

fpΘ(t − t^′){ϕ^○p(r)ϕ^○∗r (r)ϕ^○r(r^′)ϕ^○∗p (r^′)e^{i(εr−εp)(t−t′)}−ϕ^○∗p (r)ϕ^○r(r)ϕ^○∗r (r^′)ϕ^○p(r^′)e^{−i(εr−εp)(t−t′)}}. (B14)

Performing a Fourier transformation, the corresponding frequency- domain expression becomes

χ^TAOs (r, r^′, ω) =

r≠p

∑

pr

fp{

ϕ^○p(r)ϕ^○∗r (r)ϕ^○r(r^′)ϕ^○∗p (r^′) ω − (εr−εp)+iη

−

ϕ^○∗p (r)ϕ^○r(r)ϕ^○∗r (r^′)ϕ^○p(r^′) ω + (εr−εp)+iη }

=

r≠p

∑

pr

(fp−fr)

ϕ^○∗p (r^′)ϕ^○r(r^′)ϕ^○∗r (r)ϕ^○p(r) ω − (εr−εp)+iη ,

where η → 0. (B15)

We note that there are noself-transition terms in both Eqs.(B15)and (B12)since every TD orbital is considered as an orthonormalized function at any given instant of time. As a result, an explicit response function for a non-interacting reference system (TAO system) is obtained, and the resulting expression is similar to the conventional TDDFT.⁷

3. Alternative path to Casida’s equation Recall the partition of effective potential⁵⁷

δv^TAOeff (r, ω) = δvext(r, ω) +∫ dr1δρ(r1, ω) ⋅ 𝕗Hxcθ(r, r1, ω)

=δvext(r, ω) + δvHxcθ[ρ](r, ω), (B16) where 𝕗Hxcθis theFock matrix defined in Eq.(17). Since an infinitesi- mal external field change is considered [δvext(r1, ω) → 0],^8,9Eq.(B9) can be recast into

δρ(r, ω) =∫ dr1∫ dr2χ^TAOs (r, r1, ω)δρ(r2, ω) ⋅ 𝕗Hxcθ(r1, r2, ω).

(B17)

If ∫ dr 𝕗Hxcθ(r^′, r, ω) is operated on both sides of the equation, one obtains an iterative formula

δvHxcθ(r, ω) =∫ dr1∫ dr2𝕗Hxcθ(r, r1, ω)

×χs^TAO(r1, r2, ω)δvHxcθ(r2, ω). (B18) Recalling the explicit expression of the non-interacting response function in Eq.(B15), Eq.(B18)can be reformulated into

δv^Hxcθrs (ω) =

q≠p

∑

pq

∫ dr∫ dr1ϕ^○∗r (r)ϕ^○s(r)ϕ^○∗q (r1)ϕ^○p(r1)

×𝕗Hxcθ(r, r1, ω)

⎡⎢

⎢⎢

⎢⎣

(fp−fq) ⋅δvpq^Hxcθ(ω) ω − (εq−εp)+iη

⎤⎥

⎥⎥

⎥⎦

, (B19)

where δv^Hxcθrs (ω) = ∫ dr ϕ^○∗r (r)ϕ^○s(r)δvHxcθ(r, ω) is the Hxcθ poten- tial projected on the single-particle basis set. Similar to the deriva- tion in the main manuscript, the two-electron integral is defined as follows:

(rs∣𝕗Hxcθ(ω)∣pq) ≡∫ dr ϕ^○∗r (r)ϕ^○s(r)

× ∫ dr¹ϕ^○∗q (r1)ϕ^○p(r1)𝕗_Hxcθ(r, r1, ω). (B20) With a rescaling factor [ω − (εs−εr) +iη], an iterative equation in a finite basis set is obtained,

[ω − (εs−εr)] ⋅Ω^L_rs(ω) =

q≠p

∑

pq

[(rs∣𝕗_Hxcθ(ω)∣pq)( fp−fq)] ⋅Ω^L_pq(ω), (B21) where

Ω^Lrs(ω) ≡ δv^Hxcθrs (ω)

ω − (εs−εr)+iη. (B22)