Excitation Energies from Thermally-Assisted-Occupation Density Functional Theory:
Theory and Computational Implementation
arXiv:2007.07590v2 [physics.chem-ph] 14 Aug 2020
Shu-Hao Yeh,1, 2, ∗ Aaditya Manjanath,1, ∗ Yuan-Chung Cheng,2 Jeng-Da Chai,3, 4, † and Chao-Ping Hsu1, ‡
1
Institute of Chemistry, Academia Sinica, Taipei 11529, Taiwan
Department of Chemistry, National Taiwan University, Taipei 10617, Taiwan
3
Department of Physics, National Taiwan University, Taipei 10617, Taiwan
4
Center for Theoretical Physics and Center for Quantum Science and Engineering,
National Taiwan University, Taipei 10617, Taiwan
(Dated: August 1, 2020)
2
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 density functional theory (DFT) calculations
which may be prone to errors due to the lack of proper treatment in the non-dynamical correlation
effects. Recently, thermally-assisted-occupation density functional theory (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 time-dependent (TD) TAO-DFT, which
is a time-dependent, 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 (13 Σ+
u ) was described well,
with non-imaginary excitation energies. TDTAO-DFT also yields zero singlet-triplet gap in the
3 +
dissociation limit, for the ground singlet (11 Σ+
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 (EOM-CCSD).
I.
INTRODUCTION
Over the past decades, Kohn-Sham density functional
theory (KS-DFT)1,2 has been extensively used in the prediction of various ground-state properties of solids as well
as finite-sized molecules.3–5 Its time-dependent (TD) extension, known as TDDFT6–8 , has been a popular approach for computing excited-state properties, including
the absorption and emission spectra9 , photochemical reactions10 , dynamics11 , energy and electron transfer12 , etc.,
due to its low computational cost and the availability of
a plethora of computer codes in this area. The one-toone correspondence between the TD density and the TD
external potential was rigorously demonstrated by Runge
and Gross in 1984 in their theorem6 . The linear-response
framework was further introduced 7,8 , which brought forth
a paradigm shift in the simulation of excitations of quantum systems from a density-functional perspective13–15
and is the main reason behind the popularity of this
method.
However, conventional TDDFT is derived from groundstate (GS) KS-DFT which is a single-determinant–based
method. As a result, it can fail to describe the excitedstate phenomena heavily governed by non-dynamical (or
static) correlation, such as photochemistry processes involving photoinduced bond breaking, and problems associated with conical intersection7,9,16,17 . A prototypical example is the bond dissociation process of the H2
molecule. It is known that the excitation energy of the
lowest triplet state of H2 , computed using conventional
TDDFT9 , would become imaginary beyond a H-H bond
distance of 1.75 Å, a phenomenon arising from a spin
symmetry-breaking solution in the ground state18,19 , a
typical characteristic of nondynamical correlation effects.
In contrast, in wavefunction-based methods, the (nearly)
degenerate determinants are considered on an equal footing when performing a self-consistent field (SCF) calculation, 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
reasonably model systems with non-dynamical correlation,
albeit at the expense of enormous computation efforts.
For example, the works by Becke20,21 and the works by
Kong and coworkers22–24 demonstrated parametric functionals which 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 theory, they still deviate appreciably at the bond dissociation
limit compared to a full configuration interaction (FCI)
calculation22,23 . Moreover, the SCF associated with these
functionals adds to the computational effort which can
scale dramatically with the size of molecules.
2
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 method25 ,
wherein some amount of correlation has been accounted
for, by a density functional calculation. 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 functional theory26,27 and multi-configuration range-separated
DFT28,29 were developed. While the former utilizes the socalled on-top pair-density functional, the latter separates
the electron interaction operator into short- and longrange parts which are treated with DFT and wavefunction theory, respectively. Although the idea of using such
a “hybrid” scheme seems to be an attractive prospect26–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 advantage 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 singledeterminant references such that the non-dynamical correlation problem is minimal30–32 . These approaches require a well-balanced treatment of the orbitals in the reference, 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)33 was developed by
Chai in 2012 to alleviate the formidable challenge of balancing the computational cost and simultaneously incorporating the non-dynamical correlation effects with reasonable accuracy. In contrast to traditional KS-DFT, the
underlying principle of TAO-DFT is in the usage of fractional 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 correlation effects. As a result, TAO-DFT has a computational cost
similar to that of KS-DFT, which is O(N 3-4 ). In TAODFT, the entropy contribution (e.g., see Eq. (26) of Ref.
33), can reasonably capture the non-dynamical correlation energy of a system, which was discussed and numer-
ically investigated in Ref. 33, even when the simplest local density approximation (LDA) XC energy functional is
used. The XC energy functionals at the higher rungs of
Jacob’s ladder, such as the generalized-gradient approximation (GGA)34 , global hybrid35 , and range-separated
hybrid35,36 XC 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 performance of TAO-DFT
for a wide range of applications37 . Since TAO-DFT is
similar to KS-DFT in computational efficiency, TAO-DFT
has been recently adopted for the study of the electronic
properties of various nanosystems with pronounced radical nature36,38–45 . In particular, the electronic properties
(e.g., singlet-triplet energy gaps, vertical ionization potentials, vertical electron affinities, fundamental gaps, and
active orbital occupation numbers) of linear acenes and
zigzag graphene nanoribbons (i.e., systems with polyradical character) obtained from TAO-DFT33–35,38 have 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)46 XC energy functional in KS-DFT, the density
matrix renormalization group (DMRG) algorithm47,48 , the
variational two-electron reduced density matrix (2-RDM)
method49,50 , and other high-level methods51–54 .
II.
GROUND-STATE REFERENCE: TAO-DFT
In TAO-DFT33 , the electron density is represented by
the thermal equilibrium density of an auxiliary system of
Ne non-interacting electrons at a fictitious temperature θ
(in energy units):
X
fi φ∗i (r)φi (r).
(1)
ρ(r) =
i
Here, fi (a value between 0 and 1) is the fractional occupation number of the ith orbital φi , and is given by the
Fermi-Dirac distribution function
n
o−1
fi = 1 + exp (εi − µ)/θ
,
(2)
where µ is thePchemical potential for electrons, and is determined by
i fi = Ne for a given θ, orbital energies
{εi }, and total electron number Ne . This choice for the
fractional occupation function and the corresponding oneparticle density matrix has been extensively used in other
methods such as finite-temperature DFT (FT-DFT)55
and floating occupation molecular orbital-complete active
space configuration interaction (FOMO-CASCI)56 . With
this assisted occupation number and generalized density
expression, the total ground-state energy functional can
be written as
KS
EG [ρ] = TTAO[{fi ,φi }] + Vext [ρ] + EHxc
+ Eθ [ρ],
(3)
3
where TTAO is the kinetic free energy functional of noninteracting electrons (equivalent to Aθs as defined in
Eq. (24) of Ref. 33), Vext [ρ] is the energy functional of the
KS
external potential (or nuclei potential), EHxc
is the sum
of Hartree and XC energy functionals in KS-DFT, and Eθ
is the θ-dependent energy functional33 . Alternatively (to
the original derivation33 ), 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
KS
− ∇r + vext (r) + vHxc (r) + vθ (r) φi (r) = εi φi (r), (4)
2
KS
where vext , vHxc
, and vθ are the potentials (or functional derivatives) of corresponding energy functionals
KS
(i.e., Vext [ρ], EHxc
, and Eθ [ρ], 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 determination of chemical potential µ, making this
approach attractive and easy in implementation. We have
provided a variational perspective of TAO-DFT in Appendix 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 formulation8 , within the framework of TAO-DFT. In TDTAODFT, the TD density is given by
X
ρ(r, t) =
fp φ∗p (r, t)φp (r, t),
(5)
p
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 mapping between the original interacting system of electrons moving under the influence
of a TD external potential and the auxiliary system of
non-interacting particles, an action variational principle
in TAO-DFT should be established. Following the variational principle, the TD effective potential for the noninteracting TAO system can be partitioned into the following parts:
TAO
TAO
veff
(r, t) = vext (r, t) + vHxcθ
[ρ](r, t)
(6)
TAO
where vHxcθ
is the functional derivative of the Hxcθaction, which contains the time-dependent Hartree potential, exchange-correlation potential, and the θ potentials
for the fractional occupation. Further details are included
in Appendix B 1 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
∂
1
TAO
i φp (r, t) = − ∇2r + vext (r, t) + vHxcθ
[ρ](r, t) φp (r, t)
∂t
2
= F̂ (t) φp (r, t).
(7)
TAO
vHxcθ
[ρ](r, t)
We note that
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 reformulated in terms of the one-particle density
matrix P(t)9 :
∂
i P(t) = F(t), P(t)
(8)
∂t
where F(t), the time-dependent “Fock matrix”, is the
matrix representation 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) and
F(t) = F◦ + δVext (t) + δFHxcθ [P](t),
(9)
(10)
where P◦ and F◦ denote the initial conditions for solving Eq. 8, δP(t), δVext (t), and δFe-e [P](t) are the timedependent changes in the matrices of density, external
field, and the electron-electron interaction, respectively,
in the system. The initial state (at t = t0 ) is commonly
considered to be the unperturbed GS of the system for
convenience. In terms of the GS TAO orbitals:
◦
◦
Ppq
= δpq · fp ; Fpq
= δpq · εp
(11)
If the electronic eigenspectrum of a system is desired,
the amplitude of the change in the external field |δVext (t)|
is assumed to be infinitesimally small6–9 . It is therefore
suitable to consider a linear response relation between
δFHxcθ [P](t) and δP(t). Using the GS TAO orbital basis this can be obtained as
!
Hxcθ
XZ
δFrs
(t)
Hxcθ
δPpq (τ ). (12)
δFrs (t) =
dτ
δPpq (τ )
pq
Employing the time-domain Fourier transformations
Z
(13)
δPqr (ω) = dt e−iωt δPqr (t)
Z
(14)
δVqr (ω) = dt e−iωt δVqr (t)
Z
Hxcθ
Hxcθ
δFrs
(t)
δFrs
(ω) = dt e−iω(t−τ )
(15)
δPpq
δPpq (τ )
4
one could recast Eq. 8 into
"
X
◦
◦
Fpq
· δPqr (ω) − δPpq (ω) · Fqr
q
!
Hxcθ
δFpq
◦
(ω) · δPst (ω) Pqr
δP
st
st
!
#
Hxcθ
X δFqr
◦
δVqr (ω) +
− Ppq
(ω) · δPst (ω)
δPst
st
X
+ δVpq (ω) +
= ω · δPpr (ω)
(16)
by neglecting all second-order (or higher) terms. Upon
invoking the GS definitions in Eq. 11 and assuming all
δVpq (ω) to be infinitesimally small, the corresponding
working equation becomes
"
#
Hxcθ
X δFpr
(εp − εr )δPpr (ω) − (fp − fr )
(ω) δPst (ω)
δPst
st
= ω · δPpr (ω).
(17)
A conventional linear-response relation (which is the inverse of Eq. 12)7,57 gives the TD density-density response
function. The details of this derivation are provided in
Appendix B 2.
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
Hxcθ
Hxcθ
δFpr
δFpr
(ω) ≈
δPst
δPst P◦
Z
≈ d3 r d3 r′ φ∗r (r)φp (r)fHxcθ (r, r′ )φt (r′ )φ∗s (r′ )
= rp|fHxcθ |ts .
(18)
The working equation would be reduced to an eigenvalue
equation
X
(εp − εr ) · δps,st − (fp − fr ) rp fHxcθ ts · ΩR
k,st
st
= ω k · ΩR
k,pr ,
(19)
where ΩR
pr = δPpr and k denotes the k-th eigenvalue. This
can be represented in the matrix form as Casida’s equation8 :
Xk
X
Î 0̂
 B̂
= ωk
,
(20)
∗
∗
Y
Yk
0̂ −Î
B̂ Â
R
where Xk,pr = ΩR
k,p>r , Yk,rp = Ωk,p<r , denotes upward
and downward transitions, respectively. The coupling matrices are defined as
Apr,st = (εp − εr )δps δrt + Bpr,ts
Bpr,st = −(fp − fr ) rp|fHxcθ |st .
(21)
(22)
These matrices are similar in form to those derived from
conventional Casida’s equation which most TDDFT works
are based on9,59 . However, we consider the fractional occupation number difference (∆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 numbers are explicitly sourced from GS TAO. In Eq. 19, the superscript R
in ΩR implies that the eigenvectors obtained are the right
eigenvectors. Using the density-density response function
(Appendix B 2), an eigenvalue-like equation that is complementary to that in Eq. 19 can be derived. The details
are included in Appendix B 3.
B.
Idempotency in TDTAO-DFT
In KS theory, an idempotent one-electron density matrix (PP = P)9 is derived from the single-determinant
ansatz of the wavefunction, so for any first-order changes
of 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 idempotency
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 with fp +fr tending
to 1 would be dominant. These transitions require one of
the p and q orbitals to be strongly occupied, 1/2 ≤ fr ≤ 1,
with the other weakly occupied, 0 ≤ fr < 1/2. More details
on the relaxed idempotency condition for TDTAO-DFT
can be found in Appendix C accompanying this work.
IV.
COMPUTATIONAL DETAILS
We implement this formalism in the development version of Q-Chem 5.260 . All numerical results are calculated with cc-pVDZ basis set, which was determined by
performing a comprehensive convergence test of different
sets. The two-electron integrals are evaluated with the
standard quadrature grid EML(50,194)61, consisting of 50
Euler-Maclaurin62 radial grid points and 194 Lebedev63
angular grid points.
5
-0.3
H2 BOND 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 extensively for many years using a
plethora of methods. Successfully capturing the mechanism of bond dissociation within the framework of DFT
has been elusive owing to the lack of incorporation of nondynamical correlation effects. Within TAO-DFT, however, this challenge was resolved by choosing an appropriate θ of 40 mHartree33,34 . It was further shown that, at
the bond dissociation limit, the multi-reference character
was more pronounced33,34 .
In TDDFT, one encounters the challenge of imaginary
frequencies (i.e., excitation energies) for the triplet states
which occurs in most of the results obtained from 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 model19 . TAODFT significantly rectifies this issue, for a large enough θ
value33 . Fig. 1 shows the potential energy surface (PES) of
the first triplet excited state (13 Σ+
u ) for H2 bond dissociation 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 lack of incorporation of the non-dynamical correlation effects beyond this bond distance. In addition, this
phenomenon is observed in TDTAO-DFT simulations for
θ = 0, 10, 20, and 30 mHartree. However, for θ ≥ 40
mHartree, the imaginary-frequency issue is resolved.
We also note here that the requirement for a realvalue 13 Σ+
u excitation energy mandates a higher threshold
value for θ than that obtained through a self-consistent
scheme37 , 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 optimal
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 in Ref. 56.
Another advantageous aspect of TDTAO-DFT is that
the energy of the first triplet excited (13 Σ+
u ) state in the
dissociation limit correctly approaches the GS singlet en3 +
ergy. Fig. 2 shows the singlet-triplet (11 Σ+
g -1 Σu ) vertical
gap as a function of H-H bond dissociation computed us-
-1
-0.4
Total energy (Hartree)
V.
-1.01
-0.5
-1.02
-0.6
-1.03
-0.7
-1.04
-0.8
-1.05
3.25
3.75
4.25
4.75
5.25
5.75
6.25
-0.9
-1
-1.1
0
1
2
3
4
5
6
7
H-H bond distance (Å)
qTDTAO=0
qTDTAO=10
qTDTAO=20
qTDTAO=30
qTDTAO=40
qTDTAO=50
qTDTAO=60
(qTDTAO in mHartree)
qTDTAO=70
FIG. 1. Potential energy surface of the first triplet excited
state (13 Σ+
u ) computed using TDDFT (θ = 0 mHartree) and
TDTAO-DFT with the PBE XC-functional, cc-pVDZ basis
set and GEA version of Eθ functional34 for TAO calculations.
The inset shows a zoomed-in view for the large bond-distance
regime.
ing ground-state TAO-DFT, CCSD, and TDTAO-DFT.
3 +
To compute the 11 Σ+
g -1 Σu gap at the ground-state level
(in order to mitigate the problem of imaginary frequencies in TDDFT), it was recommended to use the unrestricted ground-state SCF formalism for H2 and other
small molecules22,23,64 . However, this does not guarantee
the convergence of the energy of the 13 Σ+
u state to that
of the 11 Σ+
g state at the bond dissociation limit for H2
for TAO-DFT (Fig. 2). This gap may violate the covalent
nature of the 3 Σ+
u state, where the energies of covalent
states 13 Σ+
and
GS (11 Σ+
u
g ) should be the same at the
bond dissociation limit18 . In other words, at this limit,
the electrons are located in the 1s orbitals of the corresponding atoms and are therefore, isolated enough with
respect to one another. This gap increases with θ due to
the increase in the energy of 13 Σ+
u and a simultaneous decrease in the energy of 11 Σ+
g (this θ-dependent decrease is
also observed for the total energy of 13 Σ+
u calculated with
TDTAO-DFT in Fig. 1). On the other hand, the trend obtained for TDTAO-DFT (Fig. 2) is in excellent agreement
with that obtained using the equation-of-motion coupledcluster singles and doubles (EOM-CCSD) method or observed in experiments65 . EOM-CCSD is used here as a
benchmark method since it is equivalent to FCI for a twoelectron system like 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
6
18
1
Sg+-3Su+ gap (eV)
16
VI. RELATIONSHIP BETWEEN θ AND
IMAGINARY FREQUENCIES: A QUALITATIVE
DESCRIPTION
3.25
2.75
14
2.25
12
1.75
1.25
10
0.75
8
0.25
6
-0.25
1.75
2.25
2.75
3.25
3.75
4.25
4.75
5.25
4
2
0
-2
0
1
2
3
4
5
6
H-H bond distance (Å)
CCSD
qTDTAO=40
qTDTAO=50
qTDTAO=60
qTDTAO=70
qGSTAO=40
qGSTAO=50
qGSTAO=60
qGSTAO=70
We perform a detailed analysis of the PESs with the
different θ values to acquire more insight about the qualitative relationship between θ and the imaginary roots.
Two molecular systems were chosen for this analysis, H2
and 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
consideration, thereby indicating its versatility. However,
we note that θ 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
experiments71 .
(qTDTAO/qGSTAO in mHartree)
FIG. 2. The energy gap as a function of H-H bond distance
(in Å) between the singlet ground state (11 Σ+
g ) and the first
excited triplet state (13 Σ+
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.
are demonstrated with low-lying PESs from EOM-CCSD
in Figs. 3 (a) and (b). The overall feature of singlet
and triplet states from TDTAO-DFT is in excellent agreement with the EOM-CCSD results, except for the chargetransfer state (11 Σ+
u ) and the missing states with double
excitation character (purple curve with unfilled squares
and golden yellow curves with unfilled diamonds in Figs.
3 (a) and (b)). We speculate that the problem with the
1 1 Σ+
u state could be due to the usage of the simple adiabatic approximation to the xcθ-kernel6,19,66,67 as well
as the time-independent occupation numbers in our formalism68–70 . The missing CCSD double excited states
also indicate the inability of TDTAO-DFT to capture the
avoided crossing between the first two 1 Σ+
g excited states
(orange shaded regions as seen in Figs. 3 (a), (c), and
(d)). A more detailed investigation is certainly required
for resolving these challenges.
VII.
CONCLUDING REMARKS
In summary, a time-dependent linear-response theory
for predicting excited-state properties based on the TAODFT 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 imaginary 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 (13 Σ+
u ) is captured by TDTAO-DFT, but is not so,
for TDDFT. Additionally, TAO-DFT (with a large fictitious temperature θ) may produce an erratic gap between the 13 Σ+
u and ground states at the dissociation 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
The Supplementary Material includes additional results
and the numerical data presented in this work.
ACKNOWLEDGMENTS
CPH acknowledges support from Academia Sinica and
the Investigator Award (AS-IA-106-M01) and the Ministry of Science and Technology of Taiwan (project
105-2113-M-001-009-MY4). JDC acknowledges support
7
2
2
1.5
CCSD
1 1Σu+
1 1Σu+
2 1Σg+
3 1Σg+
2 1Σu+
2 1Σu+
1, 2 1Πu+
1 1Πu+
3 1Σg+
6 1Σg+
4 1Σg+
2 1Σg+
5 1Σg+
2 1Σu+
3 1Σu+ (DE)
TDDFT
2 1Σg+
1 1Σu+
1, 2 1Πu+
Total energy (Hartree)
Total energy (Hartree)
TDTAO
3 1Σg+
1
0.5
0
−0.5
1.5
2 1Σg+
3 1Σg+
CCSD
3 1Σg+
2 1Σg+
6 1Σg+
4 1Σg+
5 1Σg+
1
1 1Σg+
TDDFT
3 1Σg+
0.5
0
−0.5
−1
−1
(a)
−1.5
0
1
2
3
4
5
6
−1.5
0
7
(c)
1
2
3
4
5
6
7
H-H bond distance (Å)
H-H bond distance (Å)
−0.25
1
TDTAO
1 3Σu+
CCSD
1 3Σg+
2 3Σu+
1, 2 3Πu+
2 3Σg+
1 3Σu+
1 3Σg+
2 3Σu+
1 3Πu+
3 3Σg+
TDDFT
1 3Σu+
2 3Σg+
CCSD
(double excitations)
1 3Σg+
3 3Σu+
2 3Σu+
1, 2 3Πu+
3 3Σg+
0.5
0
−0.5
Total energy (Hartree)
1.5
Total energy (Hartree)
TDTAO
−0.3
−0.35
−0.4
−0.45
−1
(b)
−1.5
0
1
2
3
4
5
H-H bond distance (Å)
6
7
−0.5
0.7
(d)
0.8
0.9
1.0
1.1
1.2
1.3
H-H bond distance (Å)
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 1 Σ+
g states. (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, that CCSD state is double excitation in nature. PBE is selected as the XC-functional for all
DFT calculations and cc-pVDZ as the basis set for all calculations.
8
(a)
qTDTAO=30
qTDTAO=31
qTDTAO=32
qTDTAO=33
qTDTAO=34
qTDTAO=35
qTDTAO=36
qTDTAO=37
qTDTAO=38
qTDTAO=39
qTDTAO=40
0.10
0.09
1
Sg+-3Su+ gap (eV)
0.08
0.07
0.06
0.05
0.04
0.03
0.02
0.01
0.00
1.5
1.7
1.9
2.1
2.3
2.5
2.7
2.9
3.1
3.3
3.5
H-H bond distance (Å)
(b)
0.10
qTDTAO=30
qTDTAO=31
qTDTAO=32
qTDTAO=33
qTDTAO=34
qTDTAO=35
qTDTAO=36
qTDTAO=37
qTDTAO=38
qTDTAO=39
qTDTAO=40
0.09
1
Sg+-3Su+ gap (eV)
0.08
0.07
0.06
0.05
0.04
0.03
0.02
0.01
0.00
1.5
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 in the latter section of this
work.
According to the partition of energy functional33,34 , the
functional derivative of the total energy functional can be
expressed as
δE[ρ]
δT TAO
δVext
δEHxcθ [ρ]
= s
+
+
,
δφi (r)
δφi (r)
δφi (r)
δφi (r)
where TsTAO is the kinetic (free) energy functional, and
Vext + EHxcθ is the energy associated with the effective
potential. The explicit derivative of the kinetic (free) energy functional would be
X Z
δ
δTsTAO
fi dr φ∗i (r) t̂ φi (r)
= ∗ ′
δφ∗j (r′ )
δφj (r )
i
X
fi ln fi + (1 − fi ) ln(1 − fi )
+θ
i
= fj · t̂ φj (r′ ) +
X
i
1.7
1.9
2.1
2.3
2.5
2.7
2.9
3.1
3.3
3.5
N-N bond distance (Å)
FIG. 4. S-T gap of (a) H2 and (b) N2 with the bond distance
and different θ (in mHartree) values, calculated using TDTAODFT with the PBE XC-functional, cc-pVDZ basis set and GEA
version of Eθ functional34 . The filled symbols indicate potential energy surfaces without any imaginary frequencies.
from the Ministry of Science and Technology of Taiwan (Grant No. MOST107-2628-M-002-005-MY3) and
National Taiwan University (Grant No. NTU-CDP105R7818). AM 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.
DATA AVAILABILITY STATEMENT
The data that supports the findings of this study are
available within the article and its supplementary material.
(A1)
+θ
X
δfi
·
δφ∗j (r′ )
Z
ln fi − ln(1 − fi ) ·
i
dr φ∗i (r) t̂ φi (r)
δfi
,
δφ∗j (r′ )
(A2)
where t̂ = −∇2 /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,
Z
Z
′′
δVext
δ
′′ δρ(r )
=
dr
·
dr
v
(r)ρ(r)
ext
δφ∗j (r′ )
δφ∗j (r′ ) δρ(r′′ )
Z
δρ(r)
= dr vext (r) ∗ ′
δφj (r )
= fj · vext (r′ )φj (r′ )
Z
X δfi
∗
· dr vext (r)φi (r)φi (r) ,
+
δφ∗j (r′ )
i
(A3)
and
Z
δEHxcθ [ρ]
δρ(r) δEHxcθ [ρ]
= dr ∗ ′ ·
δφ∗j (r′ )
δφj (r )
δρ(r)
Z
δρ(r)
= dr ∗ ′ · vHxcθ [ρ](r)
δφj (r )
= fj · vHxcθ (r′ )φj (r′ )
i
X h δfi Z
∗
+
drv
[ρ](r)φ
(r)φ
(r)
.
Hxcθ
i
i
δφ∗j (r′ )
i
(A4)
9
Combining the three terms above, an explicit expression
of the total energy functional is derived
δTsTAO
δVext
δEHxcθ [ρ]
δE[ρ]
=
+
+
δφ∗j (r′ )
δφ∗j (r′ ) δφ∗j (r′ )
δφ∗j (r′ )
Z
i
X h δfi
∗
′
·
dr
φ
(r)
t̂
φ
(r)
= fj · t̂ φj (r ) +
i
i
δφ∗j (r′ )
i
Xh
δfi i
ln fi − ln(1 − fi ) · ∗ ′
+θ
δφj (r )
i
the functional derivative with respect to orbital functions
δL[ρ]
= fj · t̂ + vext (r′ ) + vHxcθ [ρ](r′ ) φj (r′ )
∗
′
δφj (r )
(
X δfi
fi
·
θ
ln
+
δφ∗j (r′ )
1 − fi
i
)
Z
∗
+ dr φi (r) t̂ + vext (r) + vHxcθ [ρ](r) φi (r)
δfi
′
∗ (r′ ) − µfj · φj (r )
δφ
j
i
i
′
= fj · t̂ + vext (r ) + vHxcθ [ρ](r′ ) φj (r′ )
X
1
δfi
εi − θ (εi − µ) · ∗ ′
+
θ
δφj (r )
i
X δfi
X
′
λji · φi (r′ ) − µ
−
∗ (r′ ) − µfj · φj (r )
δφ
j
i
i
= fj · t̂ + vext (r′ ) + vHxcθ [ρ](r′ ) φj (r′ )
X
X δfi
λji · φi (r′ )
−
+µ
∗
′
(r
)
δφ
j
i
i
X δfi
′
−µ
∗ (r′ ) − µfj · φj (r )
δφ
j
i
= fj · t̂ + vext (r′ ) + vHxcθ [ρ](r′ ) φj (r′ )
−
+fj · vext (r′ )φj (r′ )
Z
i
X h δfi
∗
·
dr
v
(r)φ
(r)φ
(r)
+
ext
i
i
δφ∗j (r′ )
i
+fj · vHxcθ (r′ )φj (r′ )
Z
i
X h δfi
∗
+
·
dr
v
[ρ](r)φ
(r)φ
(r)
Hxcθ
i
i
δφ∗j (r′ )
i
= fj · t̂ + vext (r′ ) + vHxcθ [ρ](r′ ) φj (r′ )
X δfi
fi
· θ ln
+
δφ∗j (r′ )
1 − fi
i
Z
h
i
∗
+ dr φi (r) t̂ + vext (r) + vHxcθ [ρ](r) φi (r) .
X
λji · φj (r′ ) − µ
X
i6=j
X
λji · φi (r′ ).
− λjj + µfj · φj (r′ ) +
(A5)
(A7)
i
and using the variational condition δL[ρ]/δφ∗j (r′ ) = 0, one
obtains
Enforcing the normalization conditions for both density
and orbital functions, a Lagrangian is introduced
t̂ + vext (r′ ) + vHxcθ [ρ](r′ ) φj (r′ )
= (fj−1 λjj + µ) · φj (r′ ) + fj−1
i6=j
X
λji · φi (r′ ). (A8)
i
L[ρ] = E[ρ] −
i
Xh Z
λij dr φ∗i (r)φj (r) − δij
ij
−µ
hZ
i
dr ρTAO (r) − Ne ,
(A6)
Note that the second term in Eq. A7 indicates
that it is
P
necessary to introduce the entropy term θ[ i fi ln fi +(1−
fi ) ln(1 − fi )] to the kinetic functional, in order to preserve
the correct variational property such that the derivative
terms arising from vext and vHxcθ (last terms in Eqs. A3
and A4) are compensated.
With 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
ĥTAO [ρ](r)φi (r) = εi · φi (r),
where {λij } and µ are Lagrange multipliers. Considering
where ĥTAO = t̂ + vext + vHxcθ and εi = fi−1 λii + µ.
(A9)
10
2.
Appendix B: Detailed derivation of LR-TDTAO-DFT
1.
Variational principle for TAO action functional
and TD effective potential
Here we show that the linear response equation can also
be constructed inversely,
Starting from the action variational principle72 and its
modified form73 , we have the general definitions of action
functionals for a physical system,
δρ(r t) =
(B1)
TAO ′ ′
dr′ dt′ χTAO
(r t, r′ t′ ) δveff
(r t ), (B9)
s
χTAO
(r t, r′ t′ ) ≡
s
(B2)
(B3)
where Ψ[ρ; τ ] represents the wavefunction at time t and τ
denotes the upper bound of time-integral.
For a TDTAO system, the definition of universal action
functionals can be written similarly, following that of the
conventional TDDFT scheme73 ,
*
+
δBTAO [ρ]
δΨTAO [ρ; τ ]
.
+ i ΨTAO [ρ](τ )
δρ(r, t)
δρ(r, t)
(B6)
One can define the difference between the two functionals
as
AHxcθ [ρ] = BTAO [ρ] − B[ρ],
(B7)
+
δΨTAO [ρ; τ ]
δρ(r, t)
*
+
δΨ[ρ; τ ]
. (B8)
−vext (r, t) − i Ψ[ρ](τ )
δρ(r, t)
= veff (r, t) + i ΨTAO [ρ](τ )
hP
p
i
fp φ∗p (r, t)φp (r, t)
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 following
equation
Z t
r6=p
X
′
φp (r, t) = 1 − i
dt′
φ◦r (r) e−iεr (t−t )
0
rs
hZ
i
′
′
TAO ′ ′ ◦ ′
dr′ φ◦∗
(r
)δv
(r
,
t
)φ
(r
)
e−iεs t
r
eff
s
Z
◦∗
× dr φs (r) φ◦p (r)
×
=
*
(B10)
δρ
(r, t)
=
TAO (r′ , t′ )
TAO (r′ , t′ )
δveff
δveff
h
i
X δ φ∗p (r, t)φp (r, t)
=
fp
TAO (r′ , t′ )
δveff
p
X δφ∗p (r, t)
=
fp
φ◦ (r, t)
TAO (r′ , t′ ) p
δv
eff
p
δφp (r, t)
◦∗
+ φp (r, t) TAO ′ ′ , (B11)
δveff (r , t )
which is the TAO extension of Hartree-exchangecorrelation functionals. Summarizing the equations above,
similar to TDDFT, one can recast the effective potential
in TAO as
δAHxcθ [ρ]
vHxcθ (r, t) =
δρ(r, t)
δ
TAO
χTAO
(r t, r′ t′ ) =
s
0
The TD effective potential for TAO can be expressed as
δρTAO (r, t)
TAO (r′ , t′ )
δveff
is the density-density response function for a noninteracting TAO system. With the density expression in
terms of TD orbitals, one obtains
+
*
δΨTAO [ρ; τ ]
δATAO [ρ]
(B4)
= i ΨTAO [ρ](τ )
δρ(r, t)
δρ(r, t)
Z τ
BTAO [ρ] = ATAO [ρ] +
dtdr veff (r, t)ρ(r, t). (B5)
veff (r, t) =
Z
where
τ
D
E
∂
A[ρ] =
dt Ψ(t)
− Ĥ Ψ(t) ,
∂t
0
Z
Z
B[ρ] = A[ρ] + dt dr vext (r, t)ρ(r, t),
+
*
δΨ[ ρ; τ ]
δB[ρ]
,
= vext (r, t) + i Ψ[ρ](τ )
δρ(r, t)
δρ(r, t)
Z
Density-density response function
φ◦p (r)
×
−ie
hZ
−iεr t
Z
0
t
dt′
X
φ◦r (r)
r, r6=p
i
′
TAO ′ ′ ◦ ′
−i(εp −εr )t′
dr′ φ◦∗
,
r (r )δveff (r , t )φp (r ) e
(B12)
and the corresponding orbital response functions are expressed explicitly in terms of initial orbitals (ground-state
11
Combining Eq. B11 and Eq. B13, the time-domain noninteracting response function in TDTAO-DFT can be evaluated as follows
TAO orbitals) and orbital energies
δφp (r, t)
= −iΘ(t − t′ )e−iεr t
TAO (r′ , t′ )
δveff
X
′ ◦ ′ −i(εp −εr )t′
×
φ◦r (r)φ◦∗
r (r )φp (r )e
r, r6=p
δφ∗p (r, t)
TAO (r′ , t′ )
δveff
= iΘ(t − t′ )eiεr t
X
×
′
◦ ′ ◦∗ ′ i(εp −εr )t
φ◦∗
.
r (r)φr (r )φp (r )e
r, r6=p
χTAO
(r t, r′ t′ )
s
=
X
fp
p
(B13)
(
′
iΘ(t − t ) e
iεr t
X
◦ ′ ◦∗ ′ i(εp −εr )t′
φ◦∗
r (r)φr (r )φp (r )e
r, r6=p
φ◦p (r, t)
)
X
′
′
−iεr t
′ ◦ ′ −i(εp −εr )t
+ φ◦∗
φ◦r (r)φ◦∗
p (r, t) − iΘ(t − t ) e
r (r )φp (r )e
(B14)
r, r6=p
=i
r6=p
X
′
(
fp Θ(t − t )
pr
◦ ′ ◦∗ ′ i(εr −εp )(t−t′ )
φ◦p (r)φ◦∗
r (r)φr (r )φp (r )e
Performing a Fourier transformation, the corresponding
frequency-domain expression becomes
3.
−
◦
◦∗ ′ ◦ ′ −i(εr −εp )(t−t′ )
φ◦∗
p (r)φr (r)φr (r )φp (r )e
)
Alternative path to the Casida’s equation
Recall the partition of effective potential57
TAO
δveff
(r, ω)
χTAO
(r, r′ , ω)
s
=
r6=p
X
r6=p
′ ◦ ′ ◦∗
◦
X
φ◦∗
p (r )φr (r )φr (r)φp (r)
,
(fp − fr )
=
ω − (εr − εp ) + iη
pr
where η → 0.
dr1 δρ(r1 , ω) · fHxcθ (r, r1 , ω)
= δvext (r, ω) + δvHxcθ [ρ](r, ω).
(
◦ ′ ◦∗ ′
φ◦p (r)φ◦∗
r (r)φr (r )φp (r )
ω − (εr − εp ) + iη
pr
)
◦
◦∗ ′ ◦ ′
φ◦∗
p (r)φr (r)φr (r )φp (r )
−
ω + (εr − εp ) + iη
fp
= δvext (r, ω) +
Z
(B15)
(B16)
where fHxcθ is the Fock matrix defined in Eq. 17 in the
main manuscript. Since an infinitesimal external field
change is considered (δvext (r1 , ω) → 0)8,9 , Eq. B9 can be
recast into
Z
Z
δρ(r, ω) = dr1 dr2 χTAO
(r, r1 , ω)
s
×δρ(r2 , ω) · fHxcθ (r1 , r2 , ω).
(B17)
R
If dr fHxcθ (r′ , r, ω) is operated on both sides of the equation, one obtains an iterative formula
We note that there are no self-transition terms in both
Eq. B15 and Eq. 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 noninteracting reference system (TAO system) is obtained,
and the resulting expression is similar to the conventional
TDDFT7 .
δvHxcθ (r, ω) =
Z
dr1
Z
dr2 fHxcθ (r, r1 , ω)
×χTAO
(r1 , r2 , ω)δvHxcθ (r2 , ω) (B18)
s
Recalling the explicit expression of non-interacting response function in Eq. B15, Eq. B18 can be reformulated
12
into
Hxcθ
δvrs
(ω) =
q6=p Z
X
pq
dr
Z
◦
◦∗
◦
dr1 φ◦∗
r (r)φs (r)φq (r1 )φp (r1 )
Hxcθ
(fp − fq ) · δvpq
(ω)
,
× fHxcθ (r, r1 , ω)
ω − (εq − εp ) + iη
(B19)
R
Hxcθ
◦
where δvrs
(ω) =
dr φ◦∗
r (r)φs (r)δvHxcθ (r, ω) is the
Hxcθ potential projected on the single-particle basis set.
Similar to the derivation in main manuscript, the twoelectron integral is defined as follows:
rs fHxcθ (ω) pq
(B20)
Z
Z
◦
◦
≡ dr φ◦∗
dr1 φ◦∗
r (r)φs (r)
q (r1 )φp (r1 )fHxcθ (r, r1 , ω).
With a rescaling factor [ω − (εs − εr ) + iη], an iterative
equation in a finite basis set is obtained
h
i
ω − (εs − εr ) · ΩLrs (ω)
=
q6=p
X
pq
rs fHxcθ (ω) pq (fp − fq ) · ΩLpq (ω), (B21)
where
where Ppq is a matrix element of transition density matrices. This condition leads to the result that only transitions between occupied and virtual orbitals would contribute to a physical (single) excitation. On the other
hand, since the single-determinant assumption is removed
from TAO-DFT, we consider an alternative invariant assumption based on the recurrence relation of the derivative
of Fermi function
∂fp
= −(fp − fp2 )/θ,
∂εp
(C2)
or in matrix representation is
∂P0
= −(P0 − P20 )/θ,
∂F0
(C3)
where P0 −P20 on the left-hand-side implies a relaxed idempotency feature of TAO one-particle density matrix. In
other words, instead of equating to zero, P0 − P20 is associated with another constant, θ · ∂P0 /∂F0 . To employ the
relaxed condition in excited state TAO, we further assume
that the simple partial derivative form would be preserved
in the TD extension of ∂P0 /∂F0 . Recall the total functional derivative of the density matrix
δ(P − P2 ) = δP − P0 · δP − δP · P0 ,
(C4)
and combine it with Eq. C2
ΩLrs (ω) ≡
Hxcθ
δvrs
(ω)
ω − (εs − εr ) + iη
.
(B22)
Within ALDA, the corresponding eigenvalue equation
would be
X
(εq − εp )δqs δpr − rs fHxcθ pq (fq − fp ) ΩLk,pq
pq
= ωk · ΩLk,rs ,
(B23)
where k denotes the k-th eigenvalue. We note that this
eigenvalue equation is not exactly the same as Eq. 19 the
main text. However, because of the transpose relation
between the two matrices, they will generate the same
eigenspectra.
Appendix C: Relaxed idempotency condition
In conventional TDDFT, transitions between orbital are
pre-selected by the idempotency condition9 , which is derived from a single-determinant assumption, and can be
formulated as
(fp + fq − 1)δPpq = 0,
(C1)
δPpq − fp · δPpq − δPpq · fq = −(fp + fq − 1)δPpq
∂P
,
(C5)
= −θ · δ
∂F pq
where δ[∂P/∂F]pq is not a explicit derivative and is assumed as an infinitesimal constant. Therefore, the relaxed
condition is proposed as follows:
(fp + fq − 1)δPpq ∝ θ.
(C6)
Note that the original idempotency condition would be
preserved when the KS limit is considered (θ → 0). Based
on the relaxed condition, an excitation should be dominated by those p and q terms with (fp + fq ) tending to 1.
Therefore, to reduce the interference from spurious excitations70 , only transitions between strongly occupied orbitals
and strongly virtual orbitals, where (fp + fq − 1) is minimized, are considered in the current version of TDTAODFT. The criteria to classify orbitals are
1
or εp∈occ. ≤ µ,
2
1
≤ or εq∈vir. ≥ µ.
2
fp∈occ. ≥
fq∈vir.
(C7)
13
∗
†
‡
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
These authors contributed equally.
Corresponding author: jdchai@phys.ntu.edu.tw
Corresponding author: cherri@sinica.edu.tw
P. Hohenberg and W. Kohn, Phys. Rev. 136, B864 (1964).
W. Kohn and L. J. Sham, Phys. Rev. 140, A1133 (1965).
R. O. Jones, Rev. Mod. Phys. 87, 897 (2015).
A. D. Becke, J. Chem. Phys. 140, 18A301 (2014).
N.
Mardirossian
and
M.
Head-Gordon,
Mol. Phys. 115, 2315 (2017).
E. Runge and E. K. U. Gross, Phys. Rev. Lett. 52, 997
(1984).
M. Petersilka, U. J. Gossmann, and E. K. U. Gross, Phys.
Rev. Lett. 76, 1212 (1996).
M. E. Casida, “Time-Dependent Density Functional Response Theory for Molecules | Recent Advances in Density
Functional Methods,” (WORLD SCIENTIFIC, 1995) pp.
155–192.
A. Dreuw and M. Head-Gordon, Chem. Rev. 105, 4009
(2005).
M. P. Deskevich, M. Y. Hayes, K. Takahashi, R. T. Skodje,
and D. J. Nesbitt, J. Chem. Phys. 124, 224303 (2006).
I. Tavernelli, U. F. Röhrig, and U. Rothlisberger, Mol. Phys.
103, 963 (2005).
T. L. J. Toivonen, T. I. Hukka, O. Cramariuc, T. T. Rantala,
and H. Lemmetyinen, J. Phys. Chem. A 110, 12213 (2006).
S.
Hirata
and
M.
Head-Gordon,
Chem. Phys. Lett. 314, 291 (1999).
R. E. Stratmann, G. E. Scuseria, and M. J. Frisch,
J. Chem. Phys. 109, 8218 (1998).
C.-P. Hsu, S. Hirata,
and M. Head-Gordon,
J. Phys. Chem. A 105, 451 (2001).
B. G. Levine, C. Ko, J. Quenneville, and T. J. Martı́nez,
Mol. Phys. 104, 1039 (2006).
M. Filatov, J. Chem. Theory Comput. 9, 4526 (2013).
O. V. Gritsenko, S. J. A. van Gisbergen, A. Görling, and
E. J. Baerends, J. Chem. Phys. 113, 8478 (2000).
M. E. Casida, F. Gutierrez, J. Guan, F.-X. Gadea,
D. Salahub, and J.-P. Daudey, J. Chem. Phys. 113, 7062
(2000).
A. D. Becke, J. Chem. Phys. 122, 064101 (2005).
A. D. Becke, J. Chem. Phys. 138, 074109 (2013).
E. Proynov, Y. Shao, and J. Kong, Chem. Phys. Lett. 493,
381 (2010).
E. Proynov, F. Liu, Y. Shao, and J. Kong, J. Chem. Phys.
136, 034102 (2012).
J. Kong and E. Proynov, J. Chem. Theory Comput. 12, 133
(2016).
J. Gräfenstein and D. Cremer, Chem. Phys. Lett. 316, 569
(2000).
G. Li Manni, R. K. Carlson, S. Luo, D. Ma, J. Olsen, D. G.
Truhlar, and L. Gagliardi, J. Chem. Theory Comput. 10,
3669 (2014).
G. Li Manni, R. K. Carlson, S. Luo, D. Ma, J. Olsen, D. G.
Truhlar, and L. Gagliardi, J. Chem. Theory Comput. 12,
458 (2016).
E. Fromager, S. Knecht, and H. J. A. Jensen, J. Chem.
Phys. 138, 084101 (2013).
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
K. Sharkas, A. Savin, H. J. A. Jensen, and J. Toulouse, J.
Chem. Phys. 137, 044104 (2012).
J.
F.
Stanton
and
J.
Gauss,
J. Chem. Phys. 101, 8938 (1994).
M.
Nooijen
and
R.
J.
Bartlett,
J. Chem. Phys. 102, 3629 (1998).
Y. Shao, M. Head-Gordon,
and A. I. Krylov,
J. Chem. Phys. 118, 4807 (2003).
J.-D. Chai, J. Chem. Phys. 136, 154104 (2012).
J.-D. Chai, J. Chem. Phys. 140, 18A521 (2014).
J.-D. Chai, J. Chem. Phys. 146, 044102 (2017).
F. Xuan, J.-D. Chai, and H. Su, ACS Omega 4, 7675 (2019).
C.-Y. Lin, K. Hui, J.-H. Chung, and J.-D. Chai, RSC Adv.
7, 50496 (2017).
C.-S. Wu and J.-D. Chai, J. Chem. Theory Comput. 11,
2003 (2015).
C.-N. Yeh and J.-D. Chai, Sci. Rep. 6, 30562 (2016).
S. Seenithurai and J.-D. Chai, Sci. Rep. 6, 33081 (2016).
C.-S. Wu, P.-Y. Lee, and J.-D. Chai, Sci. Rep. 6, 37249
(2016).
S. Seenithurai and J.-D. Chai, Sci. Rep. 7, 4966 (2017).
S. Seenithurai and J.-D. Chai, Sci. Rep. 8, 13538 (2018).
C.-N. Yeh, C. Wu, H. Su, and J.-D. Chai, RSC Adv. 8,
34350 (2018).
J.-H. Chung and J.-D. Chai, Sci. Rep. 9, 2907 (2019).
Y. Yang, E. R. Davidson, and W. Yang, Proc. Natl Acad.
Sci. 113, E5098 (2016).
J. Hachmann, J. J. Dorando, M. Aviles, and G. K.-L. Chan,
J. Chem. Phys. 127, 134309 (2007).
W. Mizukami, Y. Kurashige, and T. Yanai, J. Chem. Theory Comput. 9, 401 (2013).
K. Pelzer, L. Greenman, G. Gidofalvi, and D. A. Mazziotti,
J. Phys. Chem. A 115, 5632 (2011).
J. Fosso-Tande, T.-S. Nguyen, G. Gidofalvi, and A. E. DePrince, J. Chem. Theory Comput. 12, 2260 (2016).
M. S. Deleuze, L. Claes, E. S. Kryachko, and J.-P. Franco̧is,
J. Chem. Phys. 119, 3106 (2003).
B. Hajgato, M. S. Deleuze, D. J. Tozer, and F. D. Proft, J.
Chem. Phys. 129, 084308 (2008).
B. Hajgato, D. Szieberth, P. Geerlings, F. D. Proft, and
M. S. Deleuze, J. Chem. Phys. 131, 224321 (2009).
B. Hajgato, M. Huzak, and M. S. Deleuze, J. Phys. Chem.
A 115, 9282 (2011).
N. D. Mermin, Phys. Rev. 137, A1441 (1965).
P. Slavı́ček and T. J. Martı́nez, J. Chem. Phys. 132, 234102
(2010).
C. Ullrich, Time-Dependent Density-Functional Theory:
Concepts and Applications, Oxford Graduate Texts (OUP
Oxford, 2012).
A. D. Becke, J. Chem. Phys. 104, 1040 (1996).
M. Casida and M. Huix-Rotllant, Annu. Rev. Phys. Chem.
63, 287 (2012).
Y. Shao et al., Mol. Phys. 113, 184 (2015).
P. M. Gill, B. G. Johnson, and J. A. Pople, Chem. Phys.
Lett. 209, 506 (1993).
C. W. Murray, N. C. Handy, and G. J. Laming, Mol. Phys.
78, 997 (1993).
14
63
64
65
66
67
68
V. I. Lebedev and D. Laikov, in Dokl. Math., Vol. 59 (1999)
pp. 477–481.
Z.-L.
Cai
and
J.
R.
Reimers,
J. Chem. Phys. 112, 527 (2000).
N. Kouchi, M. Ukai, and Y. Hatano, J. Phys. B 30, 2319
(1997).
D. J. Tozer, J. Chem. Phys. 119, 12697 (2003).
A. Dreuw, J. L. Weisman, and M. Head-Gordon, J. Chem.
Phys. 119, 2943 (2003).
K. J. H. Giesbertz, E. J. Baerends, and O. V. Gritsenko,
Phys. Rev. Lett. 101, 033004 (2008).
69
70
71
72
73
K. J. H. Giesbertz, K. Pernal, O. V. Gritsenko, and E. J.
Baerends, J. Chem. Phys. 130, 114104 (2009).
K. Giesbertz, O. Gritsenko, and E. Baerends, J. Chem.
Phys. 133, 174119 (2010).
I. N. Kadochnikov, B. I. Loukhovitski, and A. M. Starik,
Phys. Scr. 88, 058306 (2013).
E. Gross and W. Kohn, in Density Functional Theory of
Many-Fermion Systems, Advances in Quantum Chemistry,
Vol. 21, edited by P.-O. Löwdin (Academic Press, 1990) pp.
255–291.
G. Vignale, Phys. Rev. A 77, 062511 (2008).