Chem. Rev.

2005, 105, 4009−4037 4009

Single-Reference ab Initio Methods for the Calculation of Excited States of

Large Molecules
Andreas Dreuw*
Institut für Physikalische und Theoretische Chemie, Johann Wolfgang Goethe-Universität, Marie Curie-Strasse 11,
60439 Frankfurt am Main, Germany

Martin Head-Gordon
Department of Chemistry, University of California, and Chemical Sciences Division, Lawrence Berkeley National Laboratory,
Berkeley, California 94720-1470

Received April 8, 2005

Contents solved exactly for real molecular systems due to the

electron-electron interaction and the resulting high
Downloaded by NOTTINGHAM TRENT UNIV at 14:54:45:996 on June 28, 2019

1. Introduction 4009 dimensionality. Consequently, approximations had to

2. Wave-Function-Based Methods 4012 be introduced, which gave birth to the research field
2.1. Configuration Interaction Singles 4012 of quantum chemistry. The most prominent ap-
2.1.1. Derivation of the CIS Equations 4012 proximation is the Hartree-Fock approach,1,2 which
2.1.2. Properties and Limitations 4013 treats each electron independently moving in an
from https://pubs.acs.org/doi/10.1021/cr0505627.

2.2. Time-Dependent Hartree−Fock 4014 average field of all other electrons and the nuclei.
2.2.1. Concepts and Derivation of TDHF 4014 This leads to an uncoupling of the many-body SE to
2.2.2. Properties and Limitations 4015 many single-particle equations, the so-called Har-
3. Time-Dependent Density Functional Theory 4016 tree-Fock (HF) equations, and concomitantly to the
3.1. Formal Foundations 4016 familiar and in chemistry widely used single-particle
picture of molecular orbitals. However, an inherent
3.1.1. The Runge−Gross Theorem 4016
approximation in the HF method is the neglect of
3.1.2. The Action Integral 4018 electron correlation, that is, the explicit electron-
3.1.3. The Time-Dependent Kohn−Sham 4018 electron interactions. Much effort has been under-
Equation taken to recover the missing electron correlation and
3.2. Derivation of the Linear-Response TDDFT 4020 as a consequence a plethora of quantum chemical ab
Equation initio methods have emerged. Examples of such
3.3. Relation of TDDFT and TDHF, TDDFT/TDA, 4022 wave-function-based methods are Møller-Plesset
and CIS perturbation theory (MP),3,4 configuration interaction
3.4. Properties and Limitations 4023 (CI),5,6 and coupled-cluster approaches (CC).7,8
3.5. Charge-Transfer Excited States in TDDFT 4024 A conceptually different approach to include elec-
4. Analysis of Electronic Transitions 4026 tron correlation is represented by density functional
4.1. Molecular Orbitals 4027 theory (DFT),9-11 which relies on the electron density
4.2. Transition Density 4027 as a fundamental quantity. In DFT, exchange and
4.3. Difference Density 4028 correlation effects are gathered in the so-called
4.4. Attachment/Detachment Density Plots 4029 exchange-correlation (xc) functional. Since the exact
5. Illustrative Examples 4029 xc functional is unknown, it is fitted empirically to a
5.1. The Initial Steps of the Ultrafast 4029 set of experimental data or it is modeled on the basis
Photodissociation of CO-Ligated Heme of model systems such as the uniform electron gas
5.2. Charge-Transfer Excited States in 4031 and other known properties. Once determined, it is
Zincbacteriochlorin−Bacteriochlorin then employed as a universal xc functional. Depend-
Complexes ing on the choice of the ansatz, many different
6. Brief Summary and Outlook 4034 approximate xc functionals are available today. In
7. Acknowledgments 4035 general, xc functionals can be divided into three
8. References 4035 different classes: local functionals, gradient-corrected
functionals, and hybrid functionals. In fact, the art
of performing a DFT calculation is closely related to
1. Introduction choosing the appropriate xc functional for the system
Since the electronic Schrödinger equation (SE) was under investigation. However, DFT is a formally
first written down, it was clear that it cannot be exact theory, that is, if the exact xc functional were
used, exact results would be achieved. Density func-
* E-mail: andreas.dreuw@theochem.uni-frankfurt.de. tional theory is built on the famous Hohenberg-Kohn
10.1021/cr0505627 CCC: $53.50 © 2005 American Chemical Society
Published on Web 10/06/2005
4010 Chemical Reviews, 2005, Vol. 105, No. 11 Dreuw and Head-Gordon

ingredients for the formulation of a many-body theory

based on the electron density alone. Present-day DFT
calculations are almost exclusively done within the
so-called Kohn-Sham formalism,13 which corre-
sponds to an exact dressed single-particle theory. In
analogy to HF theory, the electrons are treated as
independent particles moving in the average field of
all others but now with correlation included by virtue
of the xc functional. This again gives rise to single-
electron molecular orbitals and orbital energies. All
the methods described so far aim at precise calcula-
tion of the electronic ground state of molecular
systems, and in principle, they can be divided into
wave-function-based methods and density-based meth-
Andreas Dreuw was born in 1972 in Neuss, Germany. He received his ods. Especially Hartree-Fock and DFT form the
undergraduate education in chemistry at the universities of Düsseldorf basis on which excited-state calculations are per-
and Heidelberg. In 1997, he joined the theoretical chemistry group of formed, as we will show in the following paragraphs.
Prof. Lorenz Cederbaum in Heidelberg and received his Ph.D. in 2001. Knowledge about the energetic position of elec-
Then he moved to the University of California at Berkeley to join the
research group of Prof. Martin Head-Gordon as an “Emmy Noether” tronically excited states relative to the ground state,
postdoctoral fellow of the Deutsche Forschungsgemeinschaft. Since 2003, as well as information about geometric and electronic
he has been at the University of Frankfurt where he is leading an properties of excited states, is necessary for the
independent research group funded within the “Emmy Noether” program. explanation and interpretation of electronic spectra
His research interests comprise the development of theoretical methods of molecular systems. Furthermore, optically forbid-
for the calculation of excited states of large molecules as well as their
application to biological problems such as electron and energy transfer
den (dark) states, which are thus experimentally not
processes in phostosynthesis. or only very poorly accessible, very often play impor-
tant roles in determining the dynamics of electroni-
cally excited systems. Quantum chemical calculations
of such excited states can, for instance, provide useful
information and do indeed often contribute to the
fundamental understanding of excited-state dynam-
ics.
Sometimes a ground-state method can be tweaked
to calculate a particular electronically excited state,
if the ground-state formalism can be forced to con-
verge onto an energetically higher solution by intro-
duction of constraints. These constraints may be
given, for example, as a different spin multiplicity;
hence one can use a ground state method to calculate
the ground state of every possible spin multiplicity.
Equally well, the excited state of interest can belong
Martin Head-Gordon was born in Canberra, Australia (1962), and grew to a different irreducible representation of the spatial
up in Halifax, Canada, and Melbourne, Australia. He received B.Sc. (Hons) symmetry group than the ground state. Then, it is
and M.Sc. degrees in Chemistry from Monash University, Melbourne, possible to converge the ground-state calculation for
where he worked in the laboratory of Prof. Ron Brown. In 1989, he the lowest solution of each irreducible representation
obtained his Ph.D. in Chemistry from Carnegie-Mellon University, working
with Prof. John Pople. After postdoctoral research with Dr. John Tully at of the point group of the given molecular system.
Bell Laboratories, Head-Gordon became an Assistant Professor of Following this procedure, excited-state properties
Chemistry at the University of California, Berkeley, in 1992. He was such as equilibrium geometries and static electric
promoted to Full Professor in 2000, and holds a joint appointment in the moments are accessible. The excitation energy is
Chemical Sciences Division at Lawrence Berkeley National Laboratory. given as the difference of the total energies of the
He was a NSF Young Investigator (1993), a Packard Fellow (1995), and
a Sloan Fellow (1995), received the Medal of the International Society of electronic ground state and the excited state obtained
Quantum Molecular Sciences (1998), and was a Miller Research Professor in two independent calculations. This method is also
(2001−2002). His research centers on the development and application commonly referred to as ∆-method. Logically, this
of electronic structure methods and their realization as practical computer methodology breaks down immediately if one is
algorithms, with particular interests in linear scaling, excited states, and interested in excited states of the same spin multi-
electron correlation. plicity and same irreducible representation of the
theorems, HK I and HK II.12 The first ensures a one- spatial symmetry group as the one of the electronic
to-one mapping between the electron density and the ground state. This is particularly often the case if one
external potential containing the electron-nuclei seeks to compute optical spectra of large molecular
attraction and any additional magnetic or electric systems, which mostly do not exhibit any spatial
field, while HK II guarantees the existence of a symmetry. It is also commonly the case that excited-
variational principle for electron densities analogous state configurations are not well represented by the
to the famous Raleigh-Ritz principle for wave func- model commonly used for ground-state wave func-
tions. Together, HK I and HK II, are the necessary tions, for example open-shell singlets.
Calculation of Excited States of Large Molecules Chemical Reviews, 2005, Vol. 105, No. 11 4011

Today, several quantum chemical approaches for to large molecules and do not explicitly include
the calculation of excited states are available that do correlation through the ground-state wave function.
not require any a priori constraints and that yield In quantum chemical calculations, the size of a
energies and oscillator strengths of several excited molecule is defined by the number of basis functions
states in one single calculation. In analogy to ground- rather than by the number of atoms, and we are
state methods, excited-state methods can also be interested in methods that can handle ca. 5000 basis
divided into wave-function-based methods and elec- functions in a computation time of maximum a few
tron-density-based methods. Typical wave-function- days on a standard personal computer (PC). If one
based methods are CI,6,5 multireference CI (MR- assumes a typical basis set size of about 15 basis
CI)14,15 or multireference MP approaches,16,17 multi- functions per second row atom, one can then inves-
configurational self-consistent field (MCSCF) meth- tigate molecules with up to 300 such atoms. These
ods,2 for example, complete active space SCF (CASS- restrictions limit the applicable ab initio methods
CF),18 and complete active space perturbation theory essentially to the wave-function-based methods CIS35
of second order (CASPT2).19 These approaches are in and time-dependent Hartree-Fock (TDHF), as well
principle based on the explicit inclusion of excited as to the density-based approach time-dependent
states in the many-body wave function as additional DFT (TDDFT), which builds upon the electron den-
so-called “excited” Slater determinants constructed sity obtained from ground-state DFT. CIS and TDHF
from the HF ground state by swapping occupied with are fairly old methods that have been widely used
virtual orbitals. The expansion coefficients of the not only in quantum chemistry but also in nuclear
Slater determinants are then calculated via the physics, where TDHF is better known as the random
Raleigh-Ritz variation principle, which in the case phase approximation (RPA) and CIS as the Tamm-
of CI corresponds to the diagonalization of the Dancoff approximation (TDA) to it.38 TDDFT, how-
Hamiltonian matrix in the basis of the excited ever, is a very modern method that was developed
determinants.1 In multiconfigurational SCF ap- about 20 years ago39-42 and today has become one of
proaches also the expansion coefficients of the mo- the most prominent methods for the calculation of
lecular orbitals setting up the Slater determinants excited states of medium-sized to large molecules.
are reoptimized1 making these calculations prohibi- Several reviews are available in the literature that
tively computationally expensive for large molecules. focus on the theoretical derivation and development
Other prominent wave-function-based approaches of the method,41-44 but so far, no review focuses on
are the equation-of-motion20-22 and linear-response23-26 the applicability of TDDFT to various systems, its
coupled cluster theories (EOM-CC and LR-CC, re- practical benefits, and its limitations. Also the dif-
spectively), which depending on the level of trunca- ferent available techniques for the analysis of com-
tion in the CC expansion can yield very accurate plicated electronic transitions has never been re-
results. Closely related to these coupled-cluster theo- viewed. The aim of this review is thus threefold.
ries is the symmetry-adapted cluster configuration First, CIS, TDHF, and TDDFT are rigorously intro-
interaction (SAC-CI) approach.27 Propagator theories duced by outlining their derivations and theoretical
emerging from the Green’s function formalism, such footing, where special emphasis is put on their
as the algebraic diagrammatic construction (ADC) relations to each other. Second, different methods for
scheme,28-30 also provide an elegant route to the the analysis of complicated electronically excited
calculation of excited state properties. However, all states are reviewed and comparatively discussed.
these wave-function-based methods are limited to Third, we focus on the applicability of the presented
fairly small molecules due to their high computa- methods and show their limitations and point out
tional costs. their strengths and weaknesses.
The cheapest excited-state methods that include The review is organized as follows. In the next
correlation via the wave function available today are section (section 2), the wave-function-based methods
the CIS(D) approach31,32 and the approximate coupled- CIS and TDHF are presented. The derivations and
cluster scheme of second order (CC2).33,34 The CIS- fundamental theoretical concepts are outlined in
(D) approach is a perturbative correction to config- sections 2.1.1 and 2.2.1 for CIS and TDHF, respec-
uration interaction singles (CIS)35 that approximately tively, and their properties and limitations are dis-
introduces effects of double excitations for the excited cussed in sections 2.1.2 and 2.2.2. Then we turn to
states in a noniterative scheme very similar to MP2, time-dependent DFT (section 3), where we first
in which doubly excited states are coupled to the review the theoretical foundations of the theory
ground state. It can also be compared to the pertur- (section 3.1) and then rederive the basic equations
bative triples correction (T) to the CCSD scheme.36 of linear response TDDFT in a density-matrix for-
Similar in spirit is the CC2 method, which is an mulation (section 3.2). The relations between the
approximation to CCSD and of which the linear introduced methods CIS, TDHF, TDDFT/TDA, and
response functions yield excited states and oscillator TDDFT (section 3.3) and properties and limitations
strengths of about MP2 quality. CC2 is also closely of TDDFT are pointed out (section 3.4). In last
related to the ADC(2) scheme mentioned before.37 subsection 3.5 of this section, we focus on one of the
Modern implementations of these wave-function- most severe failures of TDDFT for charge-transfer
based approaches allow for the treatment of molec- excited states. In section 4, modern schemes for the
ular systems of up to about 50 atoms.34 analysis of electronic transitions are reviewed com-
In this review, we want to focus on single-reference prising the analyses of molecular orbitals (section
ab initio excited state methods, which are applicable 4.1), the transition density (section 4.2), the difference
4012 Chemical Reviews, 2005, Vol. 105, No. 11 Dreuw and Head-Gordon

density matrix (section 4.3), and attachment/detach- constructed as a linear combination of the ground
ment density plots (section 4.4). In the last section state Slater determinant and so-called “excited”
5, two illustrative examples of typical theoretical determinants, which are obtained by replacing oc-
studies of large molecules employing TDDFT are cupied orbitals of the ground state with virtual ones.
given, where special emphasis is put onto the ap- If one replaces only one occupied orbital i by one
plicability, advantages, and limitations of TDDFT. virtual orbital a and one includes only these “singly
excited” Slater determinants, Φai (r), in the CI wave
2. Wave-Function-Based Methods function expansion, one obtains the CIS wave func-
tion, ΨCIS, which thus reads
2.1. Configuration Interaction Singles
2.1.1. Derivation of the CIS Equations ΨCIS ) ∑cai Φai (r) (7)
ia
Configuration interaction singles (CIS) is the com-
putationally as well as conceptually simplest wave- The summation runs over index pairs ia and has the
function-based ab initio method for the calculation dimension n × v. This ansatz for the many-body wave
of electronic excitation energies and excited-state function is substituted into the exact time-indepen-
properties. The starting point of the derivation of the dent electronic Schrödinger equation,
CIS equations is the Hartree-Fock (HF) ground
state, Φ0(r), which corresponds to the best single Ĥ(r)Ψ(r) ) [T̂(r) + V̂el-nuc(r) + V̂el-el(r)]Ψ(r) )
Slater determinant describing the electronic ground
state of the system. It reads EΨ(r) (8)

Φ0(r) ) |φ1(r)φ2(r)...φn(r)| (1) where T̂ has the usual meaning of the kinetic energy
operator
For simplicity, we assume a closed-shell ground-state
1
T̂(r) ) -∑ ∇2i
electronic configuration, and thus, the φi(r) cor-
respond to doubly occupied spatial orbitals and n )
(9)
i 2
N/2 (N is the number of electrons). Φ0(r) is obtained
by solving the time-independent Hartree-Fock equa-
and V̂el-nuc corresponds to the electron-nuclei attrac-
tion, which is given by
tion
F̂(r)Φ0(r) ) E0Φ0(r) (2)
ZK
with V̂el-nuc(r) ) -∑∑ (10)
i K |ri - rK|
n
F̂(r) ) ∑f̂i(r) (3) where the index i runs over all electrons and K over
i all nuclei. ZK is the charge of nucleus K. The
electron-electron interaction, V̂el-el(r) is given as
and
N N
1
f̂i(r) ) ĥi(r) + Ĵi(r) - K̂i(r) (4) Vel-el ) ∑∑ (11)
i j>i |ri - rj|
In this equation, ĥ(r) contains the kinetic energy of
the ith electron and its electron-nuclei attraction, Projection onto the space of singly excited determi-
while the Coulomb operator Ĵi(r) and exchange nants, that is, multiplication of eq 8 from the left with
operator K̂i(r) describe the averaged electron-elec-
tron interactions. They are defined as 〈Φbj |, yields

Ĵi(r)φi(r) ) [∑∫
N

j
dr′
|r - r′| ]
φ/j (r′)φj(r′)
φi(r) (5)
∑〈Φbj |Ĥ|Φai 〉cai ) ECIS∑cai δijδab
ia ia
(12)

[∑∫ ]
and with
N φ/j (r′)φi(r′)
K̂i(r)φi(r) ) dr′ φj(r) (6) 〈Φbj |Ĥ|Φai 〉 ) (E0 + a - i)δijδab + (ia||jb) (13)
j |r - r′|

Solution of the Hartree-Fock equations for the one readily obtains an expression for the excitation
ground-state Slater determinant (eq 1) within a given energies ωCIS ) ECIS - E0
basis set of size K yields n occupied molecular
orbitals, φi(r), and v ) K - n virtual orbitals, φa(r). ∑{(a - i)δijδab + (ia||jb)}cai ) ωCIS∑cai δijδab (14)
Here and in the following, we use the indices i, j, k, ia ia
etc. for occupied orbitals, a, b, c, etc. for virtual ones,
and p, q, r, etc. for general orbitals. In configuration a and i are the orbital energies of the single-electron
interaction, the electronic wave function is then orbitals φa and φi, respectively, and (ia||jb) corre-
Calculation of Excited States of Large Molecules Chemical Reviews, 2005, Vol. 105, No. 11 4013

sponds to the antisymmetrized two-electron integrals, triplet states for closed-shell molecules by allowing
which are defined as positive and negative combination of R and â excita-
tions from one doubly occupied orbital. Due to its
(ia||jb) ) ∫∫dr dr′ conceptually simple ansatz (eq 7) and the listed

[ ]
properties, the CIS excited-state wave functions are
φi(r)φa(r)φj(r′)φb(r′) - φi(r)φj(r)φa(r′)φb(r′) well defined. Hence, their wave functions and corre-
(15)
|r - r′| sponding energies are directly comparable, which is
a particularly necessary prerequisite if one is inter-
Equation 14 can be nicely written in matrix nota- ested in transitions between excited states.
tion as an eigenvalue equation An analytic expression for the total energy of
excited states can be obtained from eq 14 by adding
AX ) ωX (16)
E0 and multiplying from the left with the correspond-
in which we use the unusual symbol A for the matrix ing CIS vector. It reads
representation of the Hamiltonian in the space of the
singly excited determinants to make the connection ECIS ) EHF + ∑(cai )2(a - i) + ∑ cai cbj (ia||jb) (20)
to later occurring equations more clear. ω is the ia ia,jb
diagonal matrix of the excitation energies, and X is
the matrix of the CIS expansion coefficients. The As a consequence, ECIS is analytically differentiable
matrix elements of A are given as with respect to external parameters, for example,
nuclear displacements and external fields, which
Aia,jb ) (a - i)δijδab + (ia||jb) (17) makes the application of analytic gradient techniques
for the calculation of excited-state properties such as
The excitation energies are finally obtained by solv- equilibrium geometries and vibrational frequencies
ing the following secular equation possible. Analytic first derivatives for CIS excited
states have been published,45,46 and also second
(A - ω)X ) 0 (18) derivatives are available.46 They have been imple-
mented in various computer codes, for instance,
that is, by diagonalization of the matrix A. The Q-Chem, CADPAC, or TURBOMOLE.47-49
obtained eigenvalues correspond to the excitation In general, excitation energies computed with the
energies of the excited electronic states, and its CIS method are usually overestimated, that is, they
eigenvectors to the expansion coefficients according are usually too large by about 0.5-2 eV compared
to eq 7. with their experimental values (see, for instance, refs
45, 50, and 51). This is on one hand because the
2.1.2. Properties and Limitations
“singly excited” determinants derived from the Har-
In the previous section, the derivation of the CIS tree-Fock ground state are only very poor first-order
working equations has been outlined, and it has been estimates for the true excitation energies, since the
shown how excitation energies and excited-state virtual orbital energies a are calculated for the (N
wave functions are obtained. Here, we want to + 1)-electron system instead of for the N-electron
compile some useful properties of CIS: system.1 Consequently, the orbital energy difference
(1) Since the CIS wave function is determined by (a - i), which is the leading term in eq 17, is not
variation of the expansion coefficients of the ansatz related to an excitation energy, if the canonical
(eq 7), its total energy corresponds to an upper bound Hartree-Fock orbitals are used as reference. In other
of the true ground-state energy by virtue of the words, the canonical HF orbitals are not a particu-
Raleigh-Ritz principle. Also all excited-state total larly good basis for the expansion of the correlated
energies are true upper bounds to their exact values. wave function, which then needs a high flexibility to
(2) Owing to Brillouin’s theorem, which states that compensate for this disadvantage, that is, doubly and
“singly excited” determinants Φai (r) do not couple to higher excited determinants are demanded in the
the ground state, that is, wave function. On the other hand, since electron
correlation is generally neglected within the CIS
〈Φai (r)|Ĥ|Φ0(r)〉 ) 0 (19) method, the error will be the differential correlation
energy, which must be at least on the order of the
the excited state wave functions are Hamiltonian correlation energy of one and sometimes several
orthogonal to the ground state. (3) In contrast to electron pairs. Such energies are typically on the
other truncated configuration interaction methods, order of 1 eV per electron pair, and hence, one should
for instance, CI with single and double excitations expect errors of this magnitude.
(CISD), CIS is size-consistent, that is, the total Furthermore, CIS does not obey the Thomas-
ground-state energy of two noninteracting systems Reiche-Kuhn dipole sum rule, which states that the
is independent of whether they are computed to- sum of transition dipole moments must be equal to
gether in one calculation or are computed indepen- the number of electrons.52-54 Thus, transition mo-
dently. This is immediately plausible having Bril- ments cannot be expected to be more than qualita-
louin’s theorem in mind, owing to which the ground- tively accurate. While a full presentation of the
state energy within CIS is nothing else but the computational algorithms used to evaluate the CIS
Hartree-Fock energy. And since Hartree-Fock is energy is beyond our present scope, it is still impor-
size-consistent,1 so is CIS. (4) Another useful property tant to understand the dependence of computational
of CIS is that it is possible to obtain pure singlet and cost on molecule size. The computational cost for CIS
4014 Chemical Reviews, 2005, Vol. 105, No. 11 Dreuw and Head-Gordon

calculations on large molecules is dominated by the equation, that is, the linear response. In the follow-
evaluation and processing of two-electron integrals ing, we will follow this common convention. In most
in the atomic orbital basis.45 Per state, this cost scales quantum chemical applications, the TDHF approach
with the square of the molecule size for sufficiently is used to calculate electronic excitation spectra and
large systems, because the number of significant two- frequency-dependent polarizabilities of molecular
electron integrals also grows quadratically for large systems (for the latter see, for instance, ref 61).
enough molecules (for very small molecules the Furthermore, the linear response TDHF equations
growth is fourth-order). This is under the assumption are also well-known in physics under the name
that an atomic orbital basis of fixed size (per atom) random phase approximation (RPA) and have been
is chosen. Instead, if the atomic orbital basis on a applied in various fields.38,59,60
given atom is increased in size for a given molecule,
then the computational cost grows as O(n4), making 2.2.1. Concepts and Derivation of TDHF
large basis set calculations very expensive. This can The starting point of the derivation of the TDHF
be reduced to O(n3) by employing auxiliary basis equation is the general time-dependent electronic
expansions, as is discussed in section 3.4 for TDDFT. Schrödinger equation for molecular systems
Additionally there are linear algebra operations the
cost of which scales with the cube of basis set size ∂
that are nonetheless relatively small in comparison ĤΨ(r,t) ) i Ψ(r,t) (21)
∂t
to two-electron matrix element contractions.
The small cost of the linear algebra is a result of where Ĥ is the time-dependent Hamiltonian
the use of efficient Davidson-type algorithms55 to
iteratively obtain a small number of excited-state Ĥ(r,t) ) Ĥ(r) + V̂(r,t) (22)
eigenvalues and eigenvectors, in comparison to the
where Ĥ(r) is defined as in eq 8 and V̂(r,t) corresponds
very large rank of the CIS matrix, n(occ) × n(virt).
to an arbitrary single-particle time-dependent opera-
Direct diagonalization would scale as the cube of this
tor, for example, time-dependent electric field
rank, n(occ)3 × n(virt)3, which would be with the sixth
power of molecule size. Thus iterative diagonalization N
is critical to the applicability of CIS (and the related V̂(r,t) ) ∑v̂i(r,t) (23)
TDHF and TDDFT methods) to large molecules. i
Additionally, iterative diagonalization allows memory
requirements to be kept very modest, scaling with With the approximation that Ψ(r,t) can be written
the square of system size. This is because the full as a single Slater determinant (which is the well-
CIS matrix, A, which grows as the square of the rank, known Hartree-Fock assumption) of the form
n(occ)2 × n(virt)2, or the fourth power of molecule size,
is never constructed directly. The iterative diagonal- Φ(r,t) ) |φ1(r,t)φ2(r,t)...φN(r,t)| (24)
ization consists of the repeated contraction of the CIS
matrix A against a trial vector x for each state to a time-dependent variant of the Hartree-Fock equa-
produce a residual vector r. This step in the AO basis tion is obtained56 that reads
becomes the contraction of two-electron integrals
with a density-like matrix representing x. On today’s ∂
F̂(r,t)Φ(r,t) ) i Φ(r,t) (25)
standard computers, this allows for the treatment of ∂t
fairly large molecules of about 300 first row atoms
or 5000 or so basis functions. In addition to the definition of the time-indepen-
dent Fock operator (eq 3), the operator F̂(r,t) contains
2.2. Time-Dependent Hartree−Fock the time-dependent single-particle potential V̂(r,t).
Furthermore, the Coulomb and exchange operators
The time-dependent Hartree-Fock equations were are analogously defined as in eqs 5 and 6, respec-
written down for the first time by Dirac as early as tively, but acquire a time dependence since the
1930 following a density matrix and equation-of- single-particle wave functions φi(r,t) are now time-
motion formalism.56 Since then several different dependent.
derivations have been given (see, for instance, refs Let us assume that at t ) 0 a molecular system is
38 and 57-60), most notably the one by Frenkel in a stationary state given by a single Slater deter-
using a time-dependent variation principle.57 The minant Φ0(r) that obeys the time-independent
time-dependent Hartree-Fock equations constitute Hartree-Fock equation (eq 2). Now a small time-
an approximation to the exact time-dependent Schrö- dependent perturbation is applied, and the un-
dinger equation making use of the assumption that perturbed orbitals of the Slater determinant will
the system can at all times be represented by a single respond to this perturbation but change only slightly,
Slater determinant composed of time-dependent since the perturbation is weak. The TDHF equations
single-particle wave functions. These equations, how- are obtained via time-dependent perturbation theory
ever, which will be given later in detail, do not to first order, that is, the linear response of the
correspond to the scheme that is generally referred orbitals and of the time-dependent Fock operator are
to when today’s quantum chemists speak of time- taken into account. The latter comprises the time-
dependent Hartree-Fock (TDHF). What is meant, dependent perturbing potential itself as well as the
though, are the equations that are obtained in first- response of the Coulomb and exchange operators due
order time-dependent perturbation theory from Dirac’s to the change in the orbitals. One possible derivation
Calculation of Excited States of Large Molecules Chemical Reviews, 2005, Vol. 105, No. 11 4015

of the TDHF equations is via a density matrix ground state by virtue of some classes of “doubly
formulation, which is outlined later in section 3.2 for excited“ Slater determinants, and in this context, one
TDDFT, and as mentioned earlier, several other can indeed speak about “de-excited” states.59 In
routes can be found in the literature (see, for ex- practice, however, one uses the Hartree-Fock ground
ample, refs 2, 38, 59, and 60). Therefore we postpone state as reference state, and all necessary expectation
a detailed derivation until section 3.2. After some values are evaluated with respect to it. This ap-
algebra and Fourier transformation from the time to proximation is also known as the quasi-boson ap-
the energy domain, one arrives at a non-Hermitian proximation, which in general is reasonable if corre-
eigenvalue equation, which yields the excitation lation effects are only small in the corresponding
energies and transition amplitudes as eigenvalues ground state. This again is indicated by the magni-
and corresponding eigenvectors. It can be written tude of the Y amplitudes, which are a measure of the
conveniently in matrix notation as ground-state correlation and which, as a conse-
quence, should be small compared to the X ampli-

[ A B X
B* A* Y][ ] [
)ω
1 0 X
0 -1 Y ][ ] (26)
tudes.
The time-dependent HF method exhibits similar
properties as the CIS scheme. It is a size-consistent
where the matrix elements are defined as follows method, and one can obtain pure singlet and triplet
states for closed-shell molecules. However, TDHF
Aia,jb ) δijδab(a - i) + (ia||jb) encounters problems with triplet states, and in
general, triplet spectra are only very poorly predicted
Bia,jb ) (ia||bj) (27)
by TDHF calculations. This is because the HF ground
state is used as the reference, which in many cases
The leading term on the diagonal of the A matrix is even leads to triplet instabilities.62 In contrast to CIS,
the difference of the energies of the orbitals i and a, TDHF obeys the Thomas-Reiche-Kuhn sum rule of
which are the ones from which and to which the the oscillator strengths,63 and thus, one should expect
electron is excited, respectively. The second term of improved transition moments compared to CIS.
the A matrix and the elements of the B matrix, the It is worthwhile to note that time-dependent HF
antisymmetrized two-electron integrals (eq 15), stem has a close connection to ground-state HF stability
from the linear response of the Coulomb and ex- theory, in which it is analyzed whether a converged
change operators to the first-order changes in the HF solution is stable in the sense that it corresponds
single-particle orbitals. It is worthwhile to note, and to a minimum in parameter space.2,62,64 The structure
it will indeed be important later on, that the response of the equations to be solved for a stability analysis
of the exchange operator corresponds to a Coulomb- are very similar to the ones for TDHF, and they also
like term and vice versa. contain the first (gradient) and second derivatives
Furthermore, comparison of the TDHF eq 26 with (Hessian) of the energy with respect to variational
the CIS eigenvalue eq 16 reveals that the latter is parameters.
contained in the first one: when the B matrix of the Since the energy of the excited states in TDHF are
TDHF equation is set to zero, it reduces to the CIS given by an analytical expression that is similar to
scheme. In physics, this approximation is well-known the one for CIS (eq 20), analytical first derivatives
as the Tamm-Dancoff approximation,38,59 which will are accessible and have indeed been derived and
also be discussed later in section 3.2 in the context implemented in connection with TDDFT.65,66 As we
of TDDFT. Here, we instantly realize that there exist will see later, the analytical derivatives of TDHF are
two different derivation routes to CIS. One is via necessary ingredients for analytical gradients in
projection of the Hamiltonian operator Ĥ onto the TDDFT when hybrid functionals are employed (sec-
space of “singly excited” Slater determinants, a tion 3.4).
procedure that we will term CI formalism in the In all applications where the orbitals do not exhibit
following, and the other is by time-dependent re- triplet instabilities, the difference matrix (A - B, eq
sponse theory to first order and subsequent neglect 26) becomes positive definite, which allows reduction
of the B matrix. of the non-Hermitian TDHF eq 26 into a Hermitian
2.2.2. Properties and Limitations eigenvalue equation of half the dimension

In the previous subsection, the basic concepts and (A - B)1/2(A + B)(A - B)1/2Z ) ω2Z (28)
the working equation of linear response TDHF has
been derived, and it has been shown that it yields with
excitation energies and transition vectors. In com-
parison with CIS, TDHF is an extension since it Z ) (A - B)-1/2(X + Y) (29)
contains not only “singly excited” states but also
“singly de-excited” states, which are constructed by From a technical point of view, the TDHF eq 26 can
interchanging the orbital indices i and a. This state- then in analogy to CIS again be solved using the
ment ought to reflect the attempt to give a math- Davidson procedure,55 but the computational cost is
ematical procedure a physical meaning, but of course, roughly twice the one of a CIS calculation since in
the “de-excitations” are nonphysical since one cannot addition to the CIS calculation one has to form trial
de-excite the Hartree-Fock ground state. Historically vector-matrix products also for the B matrix.
one can justify this statement, since TDHF (RPA) has Although TDHF is in some sense an extension of
been constructed to include correlation effects in the CIS or equally well CIS is an approximation to
4016 Chemical Reviews, 2005, Vol. 105, No. 11 Dreuw and Head-Gordon

TDHF, TDHF has not been very successful in the of the existence of a noninteracting reference system,
quantum chemistry community, that is, it has not the ground state of which is a single Slater determi-
been applied very often. Probably, this is because nant and the electron density of which is by con-
excitation energies calculated with TDHF are slightly struction exactly equal to the electron density of the
smaller than the ones obtained with CIS, but they interacting real system, leads to the derivation of the
are still overestimations. The effect of the additional well-known Kohn-Sham (KS) equations.
B matrix must be only small, since its elements are However, traditional KS-DFT is limited to time-
supposed to be small. Otherwise the underlying independent systems, that is, ground states, and if
“quasi-boson approximation” is a bad one. But if the one wants to establish an analogous time-dependent
elements of the B matrix are small, one can equally theory, time-dependent versions of the first and
well use the computationally cheaper CIS scheme. second HK theorems must be formulated and proven
More seriously, the problem of the above-mentioned and a time-dependent KS equation must be derived.
poor treatment of triplet states does not exist in CIS. In the following three subsections, we present the
In summary, TDHF does usually not constitute a Runge-Gross theorem, which is a time-dependent
significant improvement over CIS that would justify analogue to HK I, we analyze the role of the action
its increased computational cost. integral in a time-dependent variational principle,
and we will outline the derivation of the time-
3. Time-Dependent Density Functional Theory dependent Kohn-Sham equations. The key steps of
the mathematical derivations, as well as the associ-
Twenty years after the formulation of the Runge- ated physical concepts will be outlined.
Gross theorem,39 which laid the theoretical founda-
tion for time-dependent density functional theory 3.1.1. The Runge−Gross Theorem
(TDDFT), it has become one of the most prominent The Runge-Gross theorem can be seen as the time-
and most widely used approaches for the calculation dependent analogue of the first Hohenberg-Kohn
of excited-state properties of medium to large molec- theorem and constitutes the cornerstone of the formal
ular systems, for example, excitation energies, oscil- foundation of the time-dependent Kohn-Sham for-
lator strengths, excited-state geometries, etc. Every malism.39 It states that the exact time-dependent
week a large number of publications containing electron density, F(r,t), determines the time-depend-
successful applications of TDDFT appear in the ent external potential, V(r,t), up to a spatially con-
literature. In this section, first the formal foundations stant, time-dependent function C(t) and thus the
of TDDFT are reviewed, which comprise the Runge- time-dependent wave function, Ψ(r,t), up to a time-
Gross theorem, the role of the action integral, and dependent phase factor. The wave function is thus a
the time-dependent Kohn-Sham equation (section functional of the electron density
3.1). They represent the necessary ingredients to
derive the TDDFT equations in the linear response Ψ(r,t) ) Ψ[F(t)](t) e-iR(t) (30)
formulation (section 3.2). Afterward, relations be-
tween CIS, TDHF, TDDFT/TDA, and TDDFT are with (d/dt)R(t) ) C(t). The density, as well as the
discussed in section 3.3, and properties and limita- potential, has to fulfill certain requirements, which
tions of TDDFT are outlined (section 3.4), where we will encounter in detail as we proceed with the
special emphasis is put on the failure of TDDFT for proof of the Runge-Gross theorem. The proof starts
CT excited states, which is one of the most severe from the general time-dependent Schrödinger equa-
problems of TDDFT (section 3.5). tion

3.1. Formal Foundations ∂

i Ψ(r,t) ) Ĥ(r,t)Ψ(r,t) (31)
∂t
Traditional ground-state density functional theory
in the Kohn-Sham formulation (KS-DFT) relies on where
the Hohenberg-Kohn (HK) theorems and the exist-
ence of a noninteracting reference system, the elec- Ĥ(r,t) ) T̂(r) + V̂el-el(r) + V̂el-nuc(r) + V̂(t) (32)
tron density of which equals the electron density of
the real system (for reviews in the field of ground- T̂(r), V̂el-nuc, and V̂el-el(r) correspond to the kinetic
state DFT, the reader is referred to refs 9-11, 67, energy operator, the electron-nuclei attraction, and
and 68). The first HK theorem (HK I) establishes a electron-electron repulsion as defined in eqs 9-11,
one-to-one mapping between the exact electron den- respectively. V̂(t) is a time-dependent external po-
sity, F(r), and the exact external potential, Vext(r), and tential and is given as a sum of one-particle poten-
since Vext(r) determines the exact ground-state wave tials
function Ψ(r), the exact ground-state wave function
N
is a functional of the electron density, Ψ[F](r). The
second HK theorem (HK II) ensures the existence of V̂(t) ) ∑v̂(ri,t) (33)
a variational principle such that the electronic energy i
of a system calculated with a trial density is always
higher than the total energy obtained with the exact N is the number of electrons and is constant with
density. Both theorems together allow the construc- time. The electron density is given as
tion of an exact many-body theory using the electron
density as the fundamental quantity. The assumption F(r,t) ) ∫|Ψ(r1,r2,r3,...rN,t)|2 dr2 dr3 ... drN (34)
Calculation of Excited States of Large Molecules Chemical Reviews, 2005, Vol. 105, No. 11 4017

In the following, spin variables will be omitted for Subtraction of these equations yields
clarity. To prove the Runge-Gross theorem, it must
be demonstrated that two densities FA(r,t) and FB(r,t) ∂ A
[j (r,t) - jB(r,t)] )
evolving from a common initial state Ψ0 under the ∂t
influence of two different potentials vA(r,t) and vB- -i〈Ψ(r,t)|[ĵ(r),{ĤA(r,t) - ĤB(r,t)}]|Ψ(r,t)〉 (42)
(r,t) are always different if the two potentials differ
by more than a purely time-dependent function, that and evaluation of this expression at t ) t0 gives an
is expression that relates the time evolution of the
different current densities to the external potentials.
vA(r,t) * vB(r,t) + C(t) (35) It reads
The first assumption to be made is that the potentials ∂ A
can be expanded in a Taylor series in time around t0 [j (r,t) - jB(r,t)]t)t0 ) -F0∇(vA(r,t) - vB(r,t)) (43)
∂t
according to
k Consequently, if the potentials vA(r,t) and vB(r,t) differ
∞
1 ∂ v(r,t) at t ) t0, the right-hand side of eq 43 cannot vanish
v(r,t) ) ∑ |t0(t - t0)k identically, and hence, the current densities jA(r,t)
k)0k! ∂tk
and jB(r,t) will be different infinitesimally later than
∞
1 t0. Thus we have established a one-to-one mapping
) ∑ k!vk(r)(t - t0)k (36) between time-dependent potentials and current den-
sities. It is worthwhile to note that eq 43 holds in
k)0
this form only if eq 36 is satisfied for k ) 0. If the
Since vA(r,t) and vB(r,t) differ by more than a time- integer k for which eq 36 holds is greater than zero,
dependent function, some of the expansion coef- the (k + 1)th time derivative of (jA(r,t) - jB(r,t)) must
ficients, vAk (r) ) (∂kvA(r,t))/∂tk|t0 and vBk (r) ) (∂kvB(r,t))/ be evaluated at t ) t0 and the equation of motion is
∂tk|t0, must differ by more than a constant. Hence, to be applied (k + 1) times. Corresponding math-
there exists one smallest integer k such that ematical expressions can be found, for example, in
ref 41, which have a very similar form as eq 43 and
vAk (r) - vBk (r) * const (37) lead to the same conclusion. To be able to make this
argument, the potential v(r,t) must be expandable in
From here, the proof proceeds in two steps. First it a Taylor expansion in time according to eq 36.
will be shown that the current densities, jA(r,t) and Since now a one-to-one mapping between time-
jB(r,t), corresponding to vA(r,t) and vB(r,t) are always dependent external potentials and time-dependent
different, and in a second step, it will be derived that current densities is established, it remains to be
different current densities require different electron proven that different current densities require dif-
densities. In general, the current density is defined ferent electron densities. For this objective, the
as continuity equation (eq 39) is applied to FA(r,t) and
FB(r,t), and subtraction of the resulting equations and
1 differentiation with respect to time yields
j(r,t) ) [Ψ*(r,t)∇Ψ(r,t) - ∇Ψ*(r,t)Ψ(r,t)] (38)
2i
∂2 A ∂
By definition, the current density j(r,t) and the
2
[F (r,t) - FB(r,t)] ) -∇ [jA(r,t) - jB(r,t)] (44)
electron density F(r,t) obey the so-called continuity ∂t ∂t
equation69,70
Insertion of eq 43 yields
∂
F(r,t) ) -∇j(r,t) (39)
∂t ∂2 A
2
[F (r,t) - FB(r,t)] ) ∇[F0∇(vA(r,t) - vB(r,t))] (45)
which states that the temporal change of the electron ∂t
density in a certain volume is equal to the flux of
which corresponds to the desired relation between
current density through the surface of that volume.
time-dependent electron densities and time-depend-
Here, at initial time t ) t0, the current densities and
ent external potentials. If one can show that the
the electron densities are given as
right-hand side of this equation cannot vanish identi-
cally, it would be proven that FA(r,t) and FB(r,t) are
jA(r,t0) ) jB(r,t0) ) j0(r) different if the corresponding external potentials are
different. This proof is done by reductio ad absurdum
FA(r,t0) ) FB(r,t0) ) F0(r) (40)
assuming that the right-hand side does vanish.
According to Gauss’ theorem, the following equation
and the time-evolution of the current densities jA(r,t)
is valid
and jB(r,t) is given by their equation of motion
∂ A
j (r,t) ) -i〈Ψ(r,t)|[ĵ(r),ĤA(r,t)]|Ψ(r,t)〉 ∫F0(∇(vA - vB))2 d3r )
-∫(vA - vB)∇(F0∇(vA - vB)) d3r +
∂t
∂ B
j (r,t) ) -i〈Ψ(r,t)|[ĵ(r),ĤB(r,t)]|Ψ(r,t)〉 (41) I(vA - vB)F0∇(vA - vB) dS (46)
∂t
4018 Chemical Reviews, 2005, Vol. 105, No. 11 Dreuw and Head-Gordon

For real, experimentally realizable potentials, the Ψ(r,t) is a solution of the time-dependent Schrödinger
surface integral vanishes, because the fall off of the eq 31 with the initial condition
asymptote is at least as 1/r, and the second term of
the right-hand side vanishes by assumption. Because Ψ(r,t0) ) Ψ0(t) (51)
the integrand is nonnegative, it follows that for all r
then the wave function corresponds to a stationary
F0[∇(vA(r,t) - vB(r,t))]2 ) 0 (47) point of the quantum mechanical action integral.

and since F0 is greater than zero

t1
0 〈
A ) ∫t dt Ψ(r,t) i | ∂t∂ - Ĥ(r,t)|Ψ(r,t)〉 (52)

∇(vA(r,t) - vB(r,t)) ) 0 (48) which is a functional of F(r,t) owing to the above

proven Runge-Gross theorem, that is,
and thus

vA(r,t) ) vB(r,t) + const (49)

t1
0 〈 ∂t |
A[F] ) ∫t dt Ψ[F](r,t) i - Ĥ(r,t) Ψ[F](r,t)
∂
| 〉 (53)

which is in contradiction to the assumption (eq 35). Consequently, the exact electron density F(r,t) can be
Consequently, the right-hand side of eq 45 cannot obtained from the Euler equation
vanish identically, and for different time-dependent
external potentials at t ) t0, one obtains different δA[F]
time-dependent electron densities infinitesimally )0 (54)
δF(r,t)
later than t0. With this, the one-to-one mapping
between time-dependent densities and time-depend- when appropriate boundary conditions are applied.
ent potentials is established, and thus, the potential Furthermore from eq 32, the action integral can be
and the wave function are functionals of the density. split into two parts, one that is universal (for a given
number of electrons) and the other dependent on the
F(r,t) T v[F](r,t) + C(t) T Ψ[F](r,t) e-iR(t) applied potential v(r,t) ) Vel-nuc(r) + V(r,t)

A[F] ) B[F] - ∫t dt ∫d3r F(r,t)v(r,t)

t1
Recently, van Leeuwen presented a generalization of (55)
0
the Runge-Gross theorem and proved that a time-
dependent density F(r,t) obtained from a many- The universal functional B[F] is independent of the
particle system can under mild restriction on the potential v(r,t) and is given as
initial state always be reproduced by an external
potential in a many-particle system with different
two-particle interaction. For two states with equiva- 〈
B[F] ) ∫t dt Ψ[F](r,t) i
t1
0 | ∂t∂ - T̂(r) -
| 〉
lent initial state and the same two-particle interac-
tion, van Leeuwen’s theorem reduces to the Runge- V̂el-el(r) Ψ[F](r,t) (56)
Gross theorem.71
Furthermore, the expectation value of any quan- In summary, the variation of the action integral with
tum mechanical operator is a unique functional of respect to the density according to eq 54 is a
the density because the phase factor in the wave prescription of how the exact density can be obtained.
function cancels out according to In the next section, this stationary action principle
will be applied to derive a time-dependent Kohn-
O(t) ) 〈Ψ[F](r,t)|Ô(t)|Ψ[F](r,t)〉 ) O[F](t) (50) Sham equation in analogy to the time-independent
counterpart.
Strictly speaking, the expectation value implicitly
depends also on the initial state, Ψ0, that is, it is a 3.1.3. The Time-Dependent Kohn−Sham Equation
functional of F and Ψ0. For most cases, however, In analogy to the derivation of the time-indepen-
when Ψ0 is a nondegenerate ground state, O[F](t) is dent Kohn-Sham equations, it is assumed that a
a functional of the density alone, because Ψ0 is a time-dependent noninteracting reference system ex-
unique functional of its density F0 by virtue of the ists with external one-particle potential vS(r,t) of
traditional first Hohenberg-Kohn theorem. which the electron density FS(r,t) is equal to the exact
3.1.2. The Action Integral electron density F(r,t) of the real interacting system.
According to the generalization of the Runge-Gross
In the previous section, the one-to-one mapping theorem by van Leeuwen,71 the existence of a time-
between time-dependent potentials and time-depend- dependent noninteracting reference system is usually
ent functionals has been established, which repre- ensured. The noninteracting system is then repre-
sents the first step in the development of a time- sented by a single Slater determinant Φ(r,t) consist-
dependent many-body theory using the density as ing of the single-electron orbitals φi(r,t); thus, its
a fundamental quantity. A second requirement is density is given by
the existence of a variational principle in analogy to
the time-independent case, in which it is given by N
the above-described second Hohenberg-Kohn theo- F(r,t) ) FS(r,t) ) ∑|φi(r,t)|2 (57)
rem. In general, if the time-dependent wave function i
Calculation of Excited States of Large Molecules Chemical Reviews, 2005, Vol. 105, No. 11 4019

Provided that the one-particle potential vS(r,t) exists, This equation is only solved by the exact interacting
the single-electron orbitals are then given as the density. Comparison with eq 62 gives an expression
solution of the time-dependent one-particle Schröd- for the time-dependent single-particle potential
inger equation
F(r′,t) δAxc[F]
∂ 1 vS(r,t) ) v(r,t) + ∫ d3r′ + (66)
(
i φi(r,t) ) - ∇2i + vS(r,t) φi(r,t)
∂t 2 ) (58) |r - r′| δF(r,t)
Inserting this equation into the time-dependent
On the other hand, the noninteracting density, single-particle Schrödinger eq 58 yields the time-
which is by assumption equal to the exact density, dependent Kohn-Sham equations
is also determined by the Euler eq 54, in which the
action integral is varied with respect to the density. ∂
For the noninteracting system (V̂el-el ) 0 in eq 56), i φ (r,t) )
∂t i

( )
the action functional takes on the following appear-
F(r′,t) δAxc[F]
- ∇2i + v(r,t) + ∫ d3r′
1
ance: + φ (r,t)
2 |r - r′| δF(r,t) i
AS[F] ) BS[F] - ∫t dt ∫ d r F(r,t)vS(r,t) (59)
t1 3
∂
0 i φ (r,t) ) F̂KSφi(r,t) (67)
∂t i
where
in which the density is given according to eq 57. The
BS[F] ) ∫t
t1
0 〈 ∂
| |
dt Ψ[F](r,t) i - T̂(r) Ψ[F](r,t) (60)
∂t 〉
time-dependent Kohn-Sham equations are, in anal-
ogy to the time-independent case, single-particle
equations in which each electron is treated individu-
and v(r,t) corresponds as usual to the time-dependent ally in the field of all others. The kinetic energy of
external potential as defined in eq 33. Applying the the electrons is represented by -1/2∇2i ; the external
stationary action principle (eq 54) yields time-dependent potential v(r,t) and the Coulomb
interaction between the charge distribution of all
δAS[F] δBS[F] other electrons with the electron under consideration
)0) - vS(r,t) (61) are explicitly contained. In analogy to the traditional
δF(r,t) δF(r,t)
time-independent Kohn-Sham scheme all exchange
If now a time-dependent single-particle potential vS- and correlation effects (explicit Coulomb interaction
(r,t) exists that allows for the construction of the time- between the electrons) are collected in (δAxc[F])/(δF-
dependent single-particle Schrödinger eq 58, this (r,t)). To this end, no approximation has been intro-
potential is a unique functional of the density by duced and consequently the time-dependent Kohn-
virtue of the Runge-Gross theorem and can accord- Sham theory is a formally exact many-body theory.
ing to eq 61 be expressed as However, the exact time-dependent “exchange-cor-
relation” action functional (also called the xc kernel)
δBS[F] is not known, and approximations to this functional
vS(r,t) ) | (62) have to be introduced. The first approximation gen-
δê(r,t) ê(r,t))F(r,t) erally made is the so-called adiabatic local density
approximation (ALDA) in which the originally non-
evaluated at the exact interacting density F(r,t). To local (in time) time-dependent xc kernel is replaced
obtain more information about the properties of BS- with a time-independent local one based on the
[F], one considers the action functional (eq 55) of the assumption that the density varies only slowly with
interacting system and as a first step rewrites it in time. This approximation allows the use of a standard
the following form local ground-state xc potential in the TDDFT frame-
work. The available approximate xc functionals will
A[F] ) BS[F] - ∫t dt ∫d3r F(r,t)v(r,t) -
t1
0
be discussed in section 3.4.
The time-dependent Kohn-Sham equations can be
F(r,t)F(r′,t)
2
∫
1 t1
t0
dt ∫d3r ∫d3r′
|r - r′|
- Axc[F] (63) conveniently expressed in matrix notation in a basis
of, say, M time-independent single-particle wave
functions {øi(r)} such that
Axc is the so-called “exchange-correlation” part of
the action integral and is defined as M
φp(r,t) ) ∑cpj(t)øj(r) (68)
F(r,t)F(r′,t)
Axc[F] ) BS[F] - ∫t dt ∫ d3r ∫ d3r′
1 t1 j
-
2 0 |r - r′|
Then, the time-dependent KS equation reads
B[F] (64)
∂
If eq 63 is inserted into the stationary action principle i C ) FKSC (69)
corresponding to the Euler eq 54, one obtains ∂t
Here, the ith column of the matrix C contains the
δBS[F] F(r′,t) δAxc[F]
) v(r,t) + ∫ d r′ + 3
(65) time-dependent expansion coefficients of φi(r,t) and
δF(r,t) |r - r′| δF(r,t) FKS is the matrix representation of the time-depend-
4020 Chemical Reviews, 2005, Vol. 105, No. 11 Dreuw and Head-Gordon

ent Kohn-Sham operator in the given basis. Multi- independent Kohn-Sham Hamiltonian matrix are
plication of eq 69 from the right with C† and then given as

{
subtraction from the resultant equation of its Her-
mitian transpose leads to the Dirac form of the time- 1 2 M ZK
dependent Kohn-Sham equation in density matrix F(0)
pq ) ∫ d 3
r φ /
p (r) - ∇ - ∑ +
form. This equation reads 2 K)1|r - RK|

∑{FpqPqr - PpqFqr} ) i ∂tPpr

q
∂
(70) ∫ d3r′
F(r′)
|r - r′|
+
δExc

δF(r) }
φq(r) (74)

In the basis of the orthonormal unperturbed single-

in which the density matrix Ppr is in general related particle orbitals of the ground state, these matrices
to the electron density via are simply given as

M
F(0)
pq ) δpqp (75)
F(r,t) ) ∑cp(t)c/q(t)øp(r)ø/q(r)
p,q
and
M
) ∑Ppqøp(r)ø/q(r) (71) P(0)
ij ) δij
p,q
P(0) (0) (0)
ia ) Pai ) Pab ) 0 (76)
To obtain excitation energies and oscillator strengths
employing the time-dependent KS approach, two Again, we follow the convention that indices i, j, etc.
different strategies can be followed. One possibility correspond to occupied orbitals, a, b, etc. correspond
is to propagate the time-dependent KS wave function to virtual orbitals and p, q, r, etc. refer to general
in time, which is referred to as real-time TDDFT.72,73 orbitals, and p is the orbital energy of the one-
This technique still has the status of an expert’s electron orbital p.
method but is beginning to be used in chemistry and Now, an oscillatory time-dependent external field
biophysics, and some successful applications have is applied, and the first-order (linear) response to this
been reported recently in the literature.74,75 In the perturbation is analyzed. In general perturbation
following section, however, we want to focus on the theory, the wave function or in this case the density
analysis of the linear response of the time-dependent matrix is assumed to be the sum of the unperturbed
KS equation. This leads to the linear-response TD- ground state and its first-order time-dependent
DFT equations, which correspond to the widely used change,
TDDFT scheme implemented in most standard quan-
tum chemistry codes. Ppq ) P(0) (1)
pq + Ppq (77)
3.2. Derivation of the Linear-Response TDDFT The same holds for the time-dependent Kohn-
Equation Sham Hamiltonian, which to first order is given as
In this section, the derivation of the linear response the sum of the ground-state KS Hamiltonian and the
TDDFT equation is presented. Using a density first-order change
matrix formalism, it is shown how the excitation
energies are obtained from the linear time-dependent Fpq ) F(0) (1)
pq + Fpq (78)
response of the time-independent ground-state elec-
tron density to a time-dependent external electric Substituting eqs 77 and 78 into the time-dependent
field. Before the time-dependent electric field is Kohn-Sham eq 70 and collecting all terms of first
applied, the system is assumed to be in its electronic order yield
ground state, which is determined by the standard
time-independent Kohn-Sham equation, which in
the density matrix formulation takes on the following ∑[F(0) (1) (1) (0) (1) (0) (0) (1)
pq Pqr - Ppq Fqr + Fpq Pqr - Ppq Fqr ] )
q
appearance: ∂
i P(1)
pr (79)
∑ {F(0)
pq P(0)
qr - P(0)
pq qr }
F(0) )0 (72) ∂t
q
The first-order change of the Kohn-Sham Hamilto-
nian consists of two terms. The first contribution
with the idempotency condition
corresponds to the applied perturbation, the time-
dependent electric field itself, and it has been shown
∑P(0) (0) (0)
pq Pqr ) Ppr (73) that it is sufficient to consider only a single Fourier
q component of the perturbation,2 which is given in
matrix notation as
F(0) (0)
pq and Ppq correspond to the Kohn-Sham Hamil-
tonian and density matrix of the unperturbed 1
gpq ) [fpq e-iωt + f/qp eiωt] (80)
ground state, respectively. The elements of the time- 2
Calculation of Excited States of Large Molecules Chemical Reviews, 2005, Vol. 105, No. 11 4021

In this equation, the matrix fpq is a one-electron the electronic transitions occur for an infinitesimal
operator and describes the details of the applied perturbation, and making use of the fact that in the
perturbation. Furthermore, the two-electron part of basis of the canonical orbitals F(0) (0)
pp ) p and Pii ) 1
the Kohn-Sham Hamiltonian reacts on the changes (eqs 75 and 76), one obtains a non-Hermitian eigen-
in the density matrix, on which it explicitly depends. value equation, the TDDFT equation,
The changes in the KS Hamiltonian due to the
change of the density are given to first order as

∂F(0)
[ A B X
B* A* Y ][ ] [
)ω
1 0 X
0 -1 Y ][ ] (88)

pq ) ∑
pq
∆F(0) P(1) (81) the structure of which is equivalent to the TDHF eq
st
st ∂Pst 26 introduced in section 2.2. Here, the elements of
the matrices A and B are given as
such that the first-order change in the KS Hamilto-
nian is altogether given as Aia,jb ) δijδab(a - i) + (ia|jb) + (ia|fxc|jb)

F(1) (0) Bia,jb ) (ia|bj) + (ia|fxc|bj) (89)

pq ) gpq + ∆Fpq (82)

Turning to the time-dependent change of the where the two-electron integrals are again given in
density matrix induced by the perturbation of the KS Mulliken notation. In comparison with the TDHF eq
Hamiltonian, this is to first order given as 26, the definitions of the matrix elements differ only
in their last terms. While in TDHF the last terms
1 -iωt correspond to the response of the nonlocal HF ex-
P(1)
pq ) [dpq e + d/qp eiωt] (83) change potential, which yields a Coulomb-like term,
2
they correspond in TDDFT to the response of the
where dpq represent perturbation densities. Inserting chosen xc potential, which replaces the HF exchange
eqs 80-83 into eq 79 and collecting the terms that potential in KS-DFT. In the ALDA approximation
are multiplied by e-iωt yield the following expression (see section 3.1.3), the response of the xc potential

[ ( )
corresponds to the second functional derivative of the
∂F(0) exchange-correlation energy, which is also called the
∑ + fpq + ∑
pq
F(0)
pq dqr - dpqF(0)
qr dst P(0)
qr - xc kernel, and is given as

( )]
q st ∂Pst
∂F(0) (ia|fxc|jb) )
fqr + ∑
qr
P(0) dst ) ωdpr (84) δ2Exc
pq
st ∂Pst ∫ d r d r′
3 3
φ/i (r)φa(r)
δF(r)δF(r′)
φ/b(r′)φj(r′) (90)

The terms multiplied by eiωt lead to the complex

Explicit expressions for the xc kernel are given, for
conjugate of the above equation. The idempotency
example, in ref 76.
condition (eq 73) gives an expression for the first-
An alternative elegant route to the derivation of
order change of the density matrix of the form
the linear response expressions for TDDFT (eq 88)
∑{P(0)
pq Pqr + Ppq Pqr } ) Ppr
(1) (1) (0) (1)
(85)
via the energy-dependent density-density response
function ø(r,r′,ω) of the interacting system, which
q
contains all physical information about how the exact
which restricts the form of the matrix dpq in eq 84 density F(r,ω) changes upon small changes in the
such that occupied-occupied and virtual-virtual external potential vext(r,ω), has been presented by
blocks dii and daa are zero, and only the occupied- Marques and Gross.44 The quantities are energy-
virtual and virtual-occupied blocks, dia and dai, dependent since they correspond to the Fourier
respectively, contribute and are taken into account. transforms of the corresponding time-dependent ones.
Remembering the diagonal nature of the unperturbed The change in the density can equally well be
KS Hamiltonian and density matrixes, one obtains calculated using the response of the noninteracting
the following pair of equations: Kohn-Sham system, øKS(r,r′,ω) and is given as

F(0) (0)
aa xai - xaiFii + δF(r,ω) ) ∫ d3r′ øKS(r,r′,ω)δvS(r′,ω) (91)

( {
fai + ∑
∂Fai

bj ∂Pbj
xbj +
∂Fai
∂Pjb
ybj P(0)
ii ) ωxai (86) }) From eq 91 a formally exact expression for the exact
density response function of the interacting system
can be derived that reads
F(0) (0)
ii yai - yaiFaa -

( { })
ø(r,r′,ω) ) øKS(r,r′,ω) + ∫ d3 r′′ ∫d3 r′′′ ø(r,r′′,ω)
∂Fia ∂Fia
ii fia + ∑
P(0) xbj + ) ωyai (87) 1
bj ∂Pbj ∂Pjb
ybj
[
|r′′ - r′′′| ]
+ fxc(r′′,r′′′,ω) øKS(r′′′,r′,ω) (92)

where we have set xai ) dai and yai ) dia to follow Knowing that the exact density response function
conventional nomenclature. In the zero-frequency possesses poles at the exact excitation energies of the
limit (fai ) fia ) 0), that is, under the assumption that system,44 one can starting from eq 92 through a series
4022 Chemical Reviews, 2005, Vol. 105, No. 11 Dreuw and Head-Gordon

of algebraic manipulations arrive at a pseudo-eigen-

value equation similar to eq 88, which yields exact
excitation energies. If furthermore the single-pole
approximation (SPA)77 and the ALDA approximation
as described in section 3.1.3 are employed, one
obtains exactly the linear-response TDDFT, eq 88.
In analogy to TDHF and CIS (section 2.2), the
Tamm-Dancoff approximation (TDA) to TDDFT has
recently been introduced.78 Again, it corresponds to Figure 1. Schematic sketch of the relation between
Hartree-Fock (HF) and density functional theory (DFT),
neglecting the matrix B in eq 88, that is, only the time-dependent Hartree-Fock and time-dependent DFT
occupied-virtual block of the initial d matrix (eq 84) (TDDFT), and configuration interaction singles (CIS) and
is taken into account. This leads to a Hermitian the Tamm-Dancoff approximation to TDDFT (TDA).
eigenvalue equation
containing the response of the Hartree-Fock ex-
AX ) ωX (93) change potential, as well as the one of the chosen xc
potential at a rate determined by the factor cHF
where the definition of the matrix elements of A is determined in the hybrid xc functional. Comparing
still the same as in eq 89. It is worthwhile to note the definitions of the matrix elements of TDHF (eqs
that TDA/TDDFT is usually a very good approxima- 27) and TDDFT (eqs 89) with eq 96, it becomes
tion to TDDFT.78,79 A possible reason may be that in apparent that the latter equation contains TDHF and
DFT correlation is already included in the ground TDDFT as limiting cases if cHF ) 1 or cHF ) 0,
state by virtue of the xc functional, which is not the respectively. In other words, eq 96 corresponds to a
case in HF theory. Since the magnitude of the Y linear combination of eqs 27 and 89 thereby combin-
amplitudes and the elements of the B matrix are a ing TDHF and TDDFT in one hybrid scheme.
measure for missing correlation in the ground state,
they should be even smaller in TDDFT than in TDHF
and, thus, be less important. TDDFT is also more
3.3. Relation of TDDFT and TDHF, TDDFT/TDA,
resistant to triplet instabilities than TDHF. and CIS
Equations 88 and 89 represent the TDDFT formal- In the previous section, we have seen that TDDFT
ism, which is solved to obtain excitation energies ω and TDHF are closely related by a hybrid scheme
and transition vectors |XY〉 when the unperturbed KS that emerges when the unperturbed Kohn-Sham
Hamiltonian (eq 74), from which the response is Hamiltonian contains both a nonlocal Hartree-Fock
derived, contains a so-called pure DFT xc potential exchange potential and a local xc potential. In Figure
and not also parts of Hartree-Fock exchange. How- 1, the relations between the introduced methods CIS
ever, today it is very common to include parts of and TDHF as well as TDDFT/TDA and TDDFT are
Hartree-Fock exchange in the xc potential, which schematically shown. Starting from the usual time-
corresponds to using so-called hybrid functionals. For independent Hartree-Fock scheme, the ground-state
instance, the widely used B3LYP functional contains Kohn-Sham equations are in principle obtained by
20% HF exchange (cHF ) 0.2 in eq 94).80 In this case, exchanging the nonlocal HF exchange potential, vHF x
the unperturbed KS Hamiltonian acquires an ad- with the local Kohn-Sham xc potential vxc. Analyses
ditional term and takes on the following appearance:
of the linear response of the ground-state density

{
calculated either with HF or with DFT to an external
1 2 M -ZK
F(0)
pq ) ∫ d 3
r φ /
p (r) - ∇ + ∑ +
time-dependent perturbation lead to the TDHF or
TDDFT schemes, respectively, as is indicated by the
2 K)1|r - RK|
vertical arrows in Figure 1.
F(r′) F(r,r′) In analogy to the ground-state methods, TDHF and
∫ d3r′ - cHF∫ d3 r′ + TDDFT are similarly related since the TDHF equa-
|r - r′| |r - r′|

}
tions can be converted into the TDDFT equations by
δExc simply replacing the response of the HF exchange
(1 - cHF) φq(r) (94) potential, ∂vHF
x /∂F1, by the response of the xc poten-
δF(r) tial from DFT, ∂vxc/∂F1. The Tamm-Dancoff ap-
Using this more general time-independent KS Hamil- proximation, that is, the neglect of the B matrix in
tonian and following the same derivation route as the TDHF or TDDFT equations, leads in the TDHF
above, one arrives at more general expressions for case to CIS, while in the TDDFT case one obtains
the TDDFT equations, which allow for the use of the TDDFT/TDA equations. Of course, CIS and
hybrid functionals. Although the non-Hermitian eigen- TDDFT/TDA are equally related as TDHF and TD-
value eq 88 remains the same, the elements of its DFT are. However, in the case of CIS there exists a
matrices A and B are now given as derivation route other than the Tamm-Dancoff ap-
proximation to TDHF, namely, via a CI formalism,
Aia,jb ) δijδab(a - i) + (ia|jb) - cHF(ij|ab) + that is, direct projection of the molecular Hamiltonian
onto the “singly excited” Slater determinants. This
(1 - cHF)(ia|fxc|jb) (95)
is only possible because the projection of the Coulomb
operator ∑1/(r1 - r2) contained in the molecular
Bia,jb ) (ia|bj) - cHF(ib|aj) + (1 - cHF)(ia|fxc|bj) (96)
Hamiltonian yields terms equivalent to the response
Calculation of Excited States of Large Molecules Chemical Reviews, 2005, Vol. 105, No. 11 4023

of the Hartree-Fock exchange and Coulomb poten- eigenvalue problem (eq 88) can be converted into a
tial. Obviously, this is not the case for TDDFT/ Hermitian one, since the orbitals are usually real. It
TDA, since the response of the xc potential contains reads
the second derivative of the approximate local xc
functional. Projection of the KS determinant would (A - B)1/2(A + B)(A - B)1/2Z ) ω2Z (97)
result in equations analogous to the CIS equation
but the matrix elements evaluated for the KS orbit- with
als.

3.4. Properties and Limitations Z ) (A - B)-1/2(X + Y) (98)

In section 3.2, it was presented how the linear The solution of the hybrid TDDFT equations scales
response formulation of time-dependent density func- very similarly to the solution of the CIS and TDHF
tional theory allows for the calculation of excitation problems, which were already discussed in sections
energies and transition vectors in the DFT frame- 2.1.2 and 2.2.2, respectively. However, hybrid TDDFT
work. At present, TDDFT represents one of the most is slightly more expensive than TDHF since in
prominent approaches for this task, especially when addition to the response of the HF exchange, the
excited states of medium-sized or large molecular response of the xc potential also needs to be evalu-
systems are under investigation. ated. In common with ground-state DFT calculations,
Since the exact local xc potential, which is the key this is usually done numerically on an atom-centered
ingredient in DFT-based methods, is not known, an three-dimensional quadrature grid.89 With thresh-
approximate xc functional has to be chosen in any olding, this step can eventually scale linearly with
practical calculation. Many different flavors of xc the size of the molecule for sufficiently large systems,
functionals are available today, for example, local and thus eventually it becomes insignificant. The
functionals such as Slater-Vosko-Wilk-Nussair prefactor, however, is quite large, particularly for
(SVWN),56,81 gradient-corrected ones (GGAs) such as finer grids, and thus for medium-sized molecules, this
Becke-Lee-Yang-Parr (BLYP),82,83 Perdew-Burke- step can dominate the computation.
Enzerhof (PBE),84,85 or Becke-Perdew 1986 (BP86),82,86 When in addition only pure xc functionals are
and hybrid functionals such as Becke3-Lee-Yang- employed that do not include parts of HF exchange,
Parr (B3LYP).80 Moreover, the development of im- the matrix (A - B) becomes diagonal (eq 89) and one
proved xc functionals is still a very active field of can fully exploit this to avoid its multiplication with
research. At present, the functionals B3LYP and PBE the trial vectors within the iterative Davidson scheme.
are the most widely used xc functionals in standard In this case, the Tamm-Dancoff approximation does
ground-state DFT applications, and although all not imply a decrease in computational cost. On the
these functionals have been developed with respect contrary, as soon as hybrid functionals are used and
to the electronic ground state, they are also employed (A - B) is thus not diagonal, the matrix multiplica-
in TDDFT calculations, which is a consequence of the tion must be performed and then TDDFT/TDA does
ALDA approximation (section 3.1.3). In many cases, save computation time by a factor of approximately
results obtained with TDDFT are quite sensitive to two. Another algorithm to solve the TDDFT or TDHF
the choice of the xc functionals, in particular, when response equations has been provided by Olsen et
local or GGA functionals are compared with hybrid al.,90 which exploits the specific paired structure of
functionals (see, for instance, refs 87 and 88). This the non-Hermitian eigenvalue eq 90. This algorithm
aspect is also discussed in section 5.1. Therefore, the does not require the involved matrices to be explicitly
reliability of TDDFT calculations should always be available and thus allows for accurate calculations
checked by comparison with either wave-function- with high dimensions.
based benchmark calculations or experimental data, Additional computation time can be saved by a
as well as by the sensitivity of the results to the factor of 3-8 depending on the system of interest if
choices of xc functional. the resolution-of-the-identity (RI) approximation91 is
Although approximate xc functionals are employed used to evaluate the Coulomb-like terms of eq 96.92,93
in TDDFT, it has been proven to yield accurate The RI approximation involves computational effort
results for valence-excited states the excitation ener- that scales with the cube of system size for Coulomb
gies of which lie well below the ionization potential. interactions (and the fourth power for HF exchange)
For such states, the typical error of TDDFT lies for fixed basis set size. It is particularly valuable for
within the range of 0.1-0.5 eV, which is almost medium and large atomic orbital basis sets because
comparable with the error of high-level correlated the cost for fixed molecule size grows only as O(n3),
approaches such as EOM-CCSD or CASPT2. How- while conventional AO integral processing grows as
ever, to reach such high accuracy within the linear O(n4). The additional error introduced in the excita-
response formulation of TDDFT, one needs to include tion energies by the RI approximation is negligible
large sets of virtual orbitals. Still, in TDDFT, the since the larger errors in the total energies of the
accuracy is reached at very favorable computational ground and excited states are subtracted out, making
cost making TDDFT applicable to fairly large mol- TDDFT applicable to very large molecular systems.
ecules. On a modern computer, molecular systems with up
In general, the hybrid TDDFT scheme is solved to 3000 basis functions (approximately 200 first row
analogously to the TDHF approach as it has been atoms) can be treated by TDDFT in the standard
outlined in section 2.2.2. Again, the non-Hermitian implementation and systems with up to 4500 or so
4024 Chemical Reviews, 2005, Vol. 105, No. 11 Dreuw and Head-Gordon

basis functions (300 atoms) can be reached by TDDFT

when the RI approximation is used.
The reason for the accuracy of TDDFT excitation
energies is that the difference of the Kohn-Sham
orbital energies, which are the leading term of the
diagonal elements of the A matrix in eq 96, are
usually excellent approximations for excitation ener-
gies.94,95 This is because the virtual KS orbital ener-
gies are evaluated for the N-electron system and,
thus, correspond more to the single-particle energy
of an excited electron than to the energy of an
additional electron as in Hartree-Fock theory, where
the virtual orbital energies are evaluated for the N
+ 1 electron system. Consequently, orbital energy
differences are a much better estimate for valence-
excited states in KS-DFT than in HF theory. Al-
though TDDFT performs usually very well for valence-
excited states, it is now well-known that TDDFT has
severe problems with the correct description of Ry- Figure 2. Schematic sketch of a typical valence-excited
dberg states, valence states of molecules exhibiting state, in which the transition occurs on one of the indi-
extended π-systems,96,97 doubly excited states,98,99 and vidual molecules, that is, the orbitals i, j and a, b are
charge-transfer excited states.100-103 For such states, located on the same molecule in contrast to a charge-
the errors in the excitation energies can be as large transfer excited state in which an electron is transferred
as a few electronvolts, and the potential energy from orbital i on molecule A to orbital a on molecule B.
surfaces can exhibit incorrect curvature. The prob- When the molecules A and B are spatially separated from
each other the orbitals i and j do not overlap with a and b.
lems with Rydberg states and extended π-systems
can be attributed to the wrong long-range behavior
of TDDFT (eqs 88 and 96) are available,65,66,110 which
of current standard xc functionals, since they decay
allow for the efficient calculation of first-order prop-
faster than 1/r, where r is the electron-nucleus
erties of excited states, for instance, their equilibrium
distance. For example, asymptotically corrected func-
geometries and dipole moments. The equilibrium
tionals such as van Leeuwen-Baerends 1994 (LB94)104
structures, dipole moments, and harmonic frequen-
or statistical averaging of orbital potentials (SAOP)105
cies of excited states calculated with TDDFT are of
or local exact exchange potentials106,107 yield substan-
generally high quality, which is comparable to that
tially improved Rydberg state excitation energies.108
of ground-state KS-DFT calculations.66 So far, the
States with substantial double excitation character
analytic gradients have been implemented in CAD-
cannot be treated within linear response theory in
PAC and TURBOMOLE.48,49
the usual ALDA approximation, since only singly
excited states are contained in the linear response
formalism.98,99 They can however be recovered when 3.5. Charge-Transfer Excited States in TDDFT
the xc kernel is allowed to be frequency/energy As already mentioned in section 3.4, TDDFT yields
dependent, since the xc kernel is in fact strongly substantial errors for charge-transfer excited
frequency-dependent close to a double excitation.98,99 states.87,100,101,103,111 The excitation energies for such
In the case of excited charge-transfer (CT) states, states are usually drastically underestimated, and
which are the topic of section 3.5, the excitation moreover, the potential energy curves of CT states
energies are much too low and the potential energy do not exhibit the correct 1/R dependence along a
curves do not exhibit the correct 1/R asymptote when charge-separation coordinate R,101,102,112 for in reality
R corresponds to a distance coordinate between the the positive and negative charges in a CT state
positive and negative charges of the CT state. electrostatically attract each other and, thus, separa-
Since TDDFT and TDHF are very closely related tion of these charges must result in an attractive 1/R
linear response theories, it is not surprising that dependence. The wrong shape of the potential energy
TDDFT possesses all the properties of TDHF (section curves excludes the CT states from the reliable
2.2.2). The linear response equations can be derived calculation of all molecular properties that involve
variationally,66 but due to the approximate nature their geometric derivative employing TDDFT or
of the employed xc functionals, comparison with the TDDFT/TDA.
exact total energy is not possible. One can say, These failures of TDDFT in the calculation of CT
though, that TDDFT is variational within the “model excited states can be understood by analysis of the
chemistry”109 defined by the approximate xc func- basic eqs 88 and 96 for the case of charge-transfer
tional. In analogy to TDHF, TDDFT is size-consis- excited states. In contrast to a valence-excited state,
tent, and it can yield pure spin singlet and triplet in a long-range, say for example, an intermolecular
states for closed-shell molecules. Furthermore, TD- charge-transfer state, an electron is transferred from
DFT obeys the Thomas-Reiche-Kuhn sum rule of an occupied orbital i of molecule A to a virtual orbital
the oscillator strengths.42 a of another molecule B (Figure 2). For simplicity,
Since 2002, analytic first geometric derivatives of let us assume that the overlap between orbitals on
the excitation energies given by the hybrid expression molecule A and orbitals on molecule B is zero.
Calculation of Excited States of Large Molecules Chemical Reviews, 2005, Vol. 105, No. 11 4025

Consequently, for such a state, all terms of eq 96

containing products of occupied and virtual orbitals
vanish. This comprises the second and fourth term
of the definition of the A matrix (eq 96) being the
response of the Coulomb potential and the xc poten-
tial of the KS operator, respectively. Only the first
term, which is the difference of the one-particle
energies of the donor orbital i on A and the acceptor
orbital a on B, and the third term originating from
the nonlocal HF exchange part of the Kohn-Sham
operator contribute to the A matrix of eq 86. This
term contributes to the matrix elements of A, since
orbitals i and j are both on A and the orbitals a and
b are on B. In fact, this term is a Coulomb-like term,
since the created holes (orbitals i and j corresponding
to the positive charge in the CT state) interact with
the electrons (orbitals a and b reflecting the negative
charge), which relates to the electrostatic attraction
within the CT state. Consequently, this term is Figure 3. Comparison of the long-range behavior of the
lowest CT excited state of an ethylene-tetrafluoroethylene
essential for the correct 1/R dependence of the dimer along the intermolecular separation coordinate
potential energy curves of CT states along the computed with TDDFT employing the SVWN (black), LB94
intermolecular separation coordinate. The same ar- (red), B3LYP (green), and “half-and-half” (blue) functionals
guments are valid for the elements of the matrix B with the curve obtained at the CIS (brown) level. The
(eq 96), and all terms are zero, that is, this matrix excitation energy of the lowest CT state at 4 Å is set to
does not contribute to CT states at all. For CT states, zero.
eq 96 reduces to the following simple expressions
contained in the second term of eq 100, which
Aiaσ,jbτ ) δστδijδab(aσ - iτ) - δστcHF(iσjσ|aτbτ) (99) corresponds to the linear response of HF exchange.
Therefore, the correct 1/R long-range behavior of the
Biaσ,jbτ ) 0 (100) potential energy surfaces can in principle be recov-
ered by inclusion of nonlocal HF exchange in the xc
It is now obvious from the definition of the A matrix potential, which will improve the asymptotic behavior
in eq 100 that the excitation energy of a CT state in according to the factor cHF of the exchange functional
TDDFT is simply given by the difference of the orbital (eq 100). This explains the observed asymptotic
energies of the electron-accepting and electron-donat- behavior of potential energy curves of CT states of a
ing molecular orbitals, a and i, respectively, when tetrafluoroethylene-ethylene complex calculated with
a pure local xc functional (for instance, SVWN, the SVWN, LB94, B3LYP, and “half-and-half” func-
BLYP,82 or LB94) is employed, that is, cHF ) 0. tional compared to configuration interaction singles
Within HF theory this is already a rough estimate (CIS) (Figure 3).101
for the energy of the CT state at large distances, since The 1/R failure of TDDFT employing pure standard
Koopman’s theorem states that -i and -a cor- xc functionals has in fact been understood as an
respond to the ionization potential of molecule A and electron-transfer self-interaction error. To clarify this,
to the electron affinity of molecule B, respectively. let us first inspect the case of TDHF (cHF ) 1 in eq
This is because the occupied orbitals are calculated 90). The excitation energy of a long-range CT state,
for the N-electron system, while the virtual orbitals where an electron is excited from orbital i on molecule
are formally evaluated for the (N + 1)-electron A into orbital a on molecule B, is dominated by the
system. This is not the case in density functional orbital energy difference a - i. In general, a
theory following the Kohn-Sham formalism (DFT), contains the Coulomb repulsion of orbital a with all
since the same potential is used to calculate the occupied orbitals of the ground state including the
occupied and virtual orbitals. As a consequence, while orbital i, which is no longer occupied in the CT state.
the HOMO still corresponds to the IP, the LUMO is In other words, the electrostatic repulsion between
generally more strongly bound in DFT than in HF orbitals a and i, the integral (aa|ii), is contained in
theory and cannot be related to the EA. Since the the orbital energy difference although orbital i is
negative of the LUMO energy is therefore much empty in the CT state. That means that the trans-
larger than the true EA, the orbital energy difference ferred electron in orbital a experiences the electro-
corresponding to a CT state is usually a drastic static repulsion with itself still being in orbital i, that
underestimation of its correct excitation energy. is, it experiences molecule A as neutral. This electron-
Furthermore, since the excitation energy of a CT transfer self-interaction effect is canceled in TDHF
state is simply given by the constant difference of the by the response of the HF exchange term, the second
corresponding orbital energies, the potential energy term of the A matrix (eq 100), being (ii|aa), which
curves of CT states do not exhibit the correct 1/R gives rise to the particle-hole attraction. In pure
shape along a distance coordinate but are constant. TDDFT employing approximate xc functionals, HF
As already mentioned, the electrostatic attraction exchange is not present and the electron-transfer self-
between the positive charge (the holes i, j on A) and interaction effect is not exactly canceled leading to
the negative charge (the electrons a, b on B) is an incorrect long-range behavior of their potential
4026 Chemical Reviews, 2005, Vol. 105, No. 11 Dreuw and Head-Gordon

energy curves. 1 1 - erf(µr12) erf(µr12)

It is also worthwhile to note another peculiarity of ) + (101)
r12 r12 r12
CT states, namely, that the B matrix (eq 100)
vanishes for a long-range CT state in the TDDFT,
as well as TDHF, case. Since the neglect of the B where r12 ) |r1 - r2|. The first term of the right-hand
matrix is equivalent to applying the Tamm-Dancoff side (rhs) corresponds to the short-range part and is
approximation, TDHF and CIS, as well as TDDFT evaluated using the xc potential from DFT, while the
and TDDFT/TDA, yield identical results for the second term, the long-range part, is calculated with
excitation energies of long-range CT states, which is exact Hartree-Fock exchange. This idea is fairly old
the difference of the corresponding orbital energies. and had originally been suggested by Stoll and Savin
This fact can be exploited in TDDFT as a first already in 1985.117-120 The scheme (eq 101) has been
diagnostic for whether one deals with a problematic applied in combination with various xc functionals
CT state. yielding, for instance, LC-BLYP, which indeed cor-
rects the failures of TDDFT for CT excited states.114
However, if the exact local Kohn-Sham xc poten- A major drawback of this approach, however, is that
tial would be used, which unfortunately is not known, the standard xc functionals require a refitting pro-
the correct 1/R long-range behavior would be found. cedure. Similar in spirit is the approach of Baer and
This is due to the derivative discontinuities of the Neuhauser, who also include long-range Hartree-
exact exchange-correlation energy with respect to Fock exchange but who employ a different partition
particle number.112,113 In the case of an electron function for the Coulomb operator.116 An extension
transfer from orbital i on molecule A to orbital a on of this approach has been presented by Yanai et al.,
molecule B, the xc potential jumps discontinuously who combine B3LYP80 at short range with an in-
by the constant IPA - EAB, leading to a singularity creasing amount of exact HF exchange at long range
in the derivative of the xc potential with respect to resulting in a functional called CAM-B3LYP,115 which
the density, that is, the xc kernel fxc used in the gives excellent CT excitation energies in comparison
TDDFT calculation (eq 96). This singularity then with benchmark calculations. But since they use at
compensates the vanishing overlap between the long range at most 60% HF exchange, the long-range
orbitals i and a when the molecules A and B are being asymptotic behavior of the CT states is not fully
separated, and thus, the fourth term of eq 96 does in corrected.115 A slightly different route is taken by
fact contribute to a CT state when the exact xc Gritsenko and Baerends who suggest a new long-
potential would be employed, and the correct 1/R range corrected xc kernel that shifts the orbital
asymptote along the separation coordinate would energy of the acceptor orbital to a value related to
then be obtained for a CT state. Furthermore, this the electron affinity. The wrong asymptotic behavior
clearly explains why standard approximate xc func- of the CT states is corrected by a distance-dependent
tionals in TDDFT fail in describing long-range CT Coulomb term correcting for electron-transfer self-
states correctly since they do not contain the deriva- interaction.95
tive discontinuities.
A completely different ansatz to the solution of the
At present, several different pathways are starting CT problem may be represented by time-dependent
to emerge to address this substantial failure of current density functional theory (TDCDFT), which
TDDFT for CT states and to correct for it. The most has recently been implemented in the Amsterdam
obvious way is to improve the xc functional by density functional (ADF) package.121,122 We have
including exact exchange in the unperturbed Hamil- previously seen that a correct description of charge-
tonian, either in the form of nonlocal Hartree-Fock transfer excited states requires the nonlocal HF
exchange or of the exact local Kohn-Sham exchange exchange potential, and current functionals are non-
potential. The latter are known as local Hartree- local. However, at present no test calculation has yet
Fock (LHF) or also optimized effective potentials been performed to show whether TDCDFT would in
(OEPs), which have recently been introduced by principle yield correct CT excited states. Also, due
Görling106 and Ivanov et al.107 These potentials may to its high computational cost the question remains
be able to yield the correct asymptotic 1/R depen- whether this method will be applicable to large
dence of excited CT states, since they possess singu- molecules soon.
larities as soon as the overlap between the electron-
donating and electron-accepting orbitals i and a,
respectively, approaches zero. This, however, remains 4. Analysis of Electronic Transitions
to be explored in detail, and it seems numerically When a molecule is excited from the electronic
difficult to compensate the vanishing overlap with ground state to an energetically higher lying excited
the divergence of the xc kernel fxc such that the fourth state, the electronic many-body wave function changes.
term of the A matrix (eq 96) cancels the electron- To gain insight in the nature of the corresponding
transfer self-interaction error correctly and recovers electronic transition, this change in the wave function
the correct 1/R asymptote. Also, satisfactory match- needs to be analyzed. In principle, this is possible by
ing correlation functionals remain to be developed. direct analysis of the excited state wave function or
Inclusion of nonlocal Hartree-Fock exchange has its electron density. For this objective, the techniques
been realized in a few schemes so far.114-116 In all developed for the analysis of the wave function or
these schemes, the Coulomb operator of the Hamil- electron density of the ground state are simply
tonian is split into two parts, a short-range and a applied to the excited state of interest. However, here
long-range part, as, for example, in ref 114, we want to focus on approaches that are dedicated
Calculation of Excited States of Large Molecules Chemical Reviews, 2005, Vol. 105, No. 11 4027

to addressing differences in the ground state and 4.2. Transition Density

excited state directly, avoiding tedious analysis of the
excited-state wave function. These approaches com- Another possibility to obtain information about an
prise inspection of the molecular orbitals involved in electronic transition is to analyze the one-electron
the transition (section 4.1), the analysis of the transi- transition density. The transition density, T(r), couples
tion density (section 4.2) or the difference density the electronic ground state with the excited state of
(section 4.3), and the so-called attachment/detach- interest and is in general given as
ment density plots (section 4.4). The advantages and
disadvantages of these analysis tools will be il- T(r) ) N∫|Ψex(r1,r2,...rn)〉〈Ψ0(r1,r2,...rn)| dr2...drn
luminated in the following subsections. (102)
4.1. Molecular Orbitals If the electronic ground state is given as a single
Slater determinant and the electronically excited
Usually, the electronic ground state of a molecular state of interest as a linear combination of determi-
system can be represented by a single Slater deter- nants, then T(r) is given as a linear combination of
minant composed of single-electron wave functions single-determinant transition densities weighted by
describing the movement of the individual electrons their expansion coefficients
in the molecule. These single-electron wave functions
are the familiar molecular orbitals (MOs) well-known T(r) )
to chemists and often used to understand molecular
processes. It is thus natural to analyze the changes ∑cia∫|Φia(r1,r2,...rn)〉〈Φ0(r1,r2,...rn)| dr2...drn (103)
in the many-body wave function, that is, the Slater ia
determinant, in terms of the MOs that constitute it.
As we have already seen in section 2.1, eq 7, an In principle, one can analyze the transition density
ansatz for an excited state wave function can be directly to obtain valuable information about the
easily constructed by exchange of an occupied MO symmetry of the transition and about the way in
with a virtual one in the ground-state Slater deter- which a one-electron operator couples two different
minant. Usually, an excited state is a linear combi- states. It is however useful to analyze the transition
nation of several such “excited” determinants, and density matrix, which is given in the molecular
the expansion coefficients are obtained by the corre- orbital basis as
sponding calculation, in case of eq 7 by a CIS
calculation. However, sometimes only one or two Tia ) 〈φi|T̂(r)|φa〉 (104)
excited determinants contribute significantly to the
excited-state wave function (i.e., cia ) 0.8-1.0), and via so-called “natural transition orbitals” (NTOs)
in such cases, it is easily possible to analyze the analogous to the well-known natural orbitals ob-
nature of the corresponding transition in terms of the tained by diagonalization of the ground-state single-
occupied and virtual MOs that have been inter- electron density.123 Since the transition density ma-
changed in the excited state determinants compared trix is a rectangular nocc × nvirt matrix, it cannot
to the ground state. Such a case is, for instance, the simply be diagonalized. Instead, the “corresponding
charge-transfer excited state of the zincbacteriochlo- orbital transformation” by Amos and Hall124 can be
rin-bacteriochlorin complex investigated in section applied, which is based on a singular value decom-
5.2. The wave function of this excited state is position of the transition density matrix
composed of essentially one Slater determinant in
which the highest occupied MO (HOMO) of the T ) USV† (105)
ground state is replaced with the lowest unoccupied
MO (LUMO), since it has an expansion coefficient of where U and V are nocc × nocc and nvirt × nvirt unitary
0.99 at the theoretical levels of TDDFT as well as matrices, respectively, and S is a singular matrix
CIS. Hence, the excited state corresponds to a single- containing the singular values of T. In general, one
electron transition from the HOMO to the LUMO (see can write
also Figure 9 in section 5.2).
Although the analysis of an electronic transition TT† ) USV†VS†U†
via MOs is computationally easy and thus convenient
and straightforward for states with only one or two ) USS†U†
significantly contributing Slater determinants, it can
become very tedious for states that are represented ) US2U† (106)
by several determinants with expansion coefficients
of similar size. Examples for this scenario are the from which follows that
wave functions of the repulsive excited states of CO-
ligated heme, the investigation of which is described U†TT†U ) S2 (107)
in section 5.1. These states are represented by a
linear combination of four “singly excited” determi- that is, the unitary transformation U diagonalizes
nants with coefficients between 0.3 and 0.7, thus the matrix TT† and consequently contains its eigen-
making it almost impossible to obtain a clear picture vectors as columns. The singular values of S are
of the nature of these electronic transitions based on given as the square root of the eigenvalues of TT†.
MOs. Therefore, one can determine the unitary transfor-
4028 Chemical Reviews, 2005, Vol. 105, No. 11 Dreuw and Head-Gordon

mation matrix U and the singular matrix S by

solving the eigenvalue equation

TT†b
ui ) λib
ui i ) 1, 2, ..., nocc (108)

Following an argument for the matrix T†T analogous

to eq 106 one finds that the unitary matrix V can be
determined by solving a corresponding eigenvalue
equation

T†Tυ
bi ) λibυi i ) 1, 2, ..., nvirt (109)

Although the matrices TT† and T†T have different Figure 4. Molecular structure (A) and NTOs of the hole
dimensions, namely, nocc × nocc and nvirt × nvirt, (B) and particle (C) densities of the lowest 3π-π* transition
respectively, their first nocc eigenvalues are identical, of [Ru(bpy)2dppz]2+.125
and the eigenvalues λnocc...λnvirt of the larger matrix
T†T are zero in the case of CIS or TDA ting the corresponding NTOs gives detailed insight
into the nature of the electronic transition.
1 g λi ) λ′i g 0 i ) 1, 2, ..., nocc An example for the successful application of NTOs
nocc
is the analysis of the nature of the energetically
∑ λi ) 1
lowest excited states of the ruthenium complex [Ru-
(bpy)2dppz]2+ (bpy ) 2,2′-bipyridine; dppz ) dipyrido-
i
[3,2-a:2′,3′-c]phenazine, which serves as very sensi-
λi ) 0 i ) (nocc + 1), ..., nvirt (110) tive luminescent reporter of DNA in aqueous solu-
tion.125 In Figure 4, the molecular structure of the
complex and the NTOs of the hole and particle of the
These additional zero eigenvalues arise from map- lowest 3π-π* transition are given. In this case, an
ping the original transition density matrix onto the analysis of the involved molecular orbitals is very
larger matrix T†T, and hence the additional eigen- tedious since the transition is represented by a linear
vectors and eigenvalues must belong to the kernel combination of several singly substituted determi-
of the map. In the case of TDHF or TDDFT, the sum nants, while on the contrary, the natural transtion
of the eigenvalues will not be identical to one but will orbitals capture 90% of the character of the transi-
deviate from that to the extent that the de-excitations tion. More examples for the application of analyses
are significant. Especially in TDDFT, these contribu- via NTOs can be found in refs 123, 125, and 126.
tions are usually small.
The occupied and virtual “natural transition orbit-
4.3. Difference Density
als”, φnto
i and φnto
a , respectively, are defined as
A complementary approach to the analysis of the
φnto
i ) φ iU i ) 1, 2, ..., nocc (111) transition density matrix T is the investigation of the
difference density matrix ∆, which is simply given
φnto
a ) φaV a ) 1, 2, ..., nocc (112) as the difference between the single electron density
matrices of the excited state Pex and the electronic
ground state P0
where in eq 112 the index a ends at nocc and not nvirt,
since nvirt - nocc virtual orbitals are mapped onto the ∆ ) Pex - P0 (113)
null vector. Following this procedure, a one-to-one
mapping between occupied and virtual orbitals is
established, because one occupied and one virtual Today, the analysis of an electronic transition by
NTO correspond to each eigenvalue λi. means of the difference density is frequently per-
In summary, the originally nocc × nvirt dimensional formed (see, for example, refs 127-130), and in
transition density matrix has been reduced to a nocc principle, valence and Rydberg excited state can be
× nocc dimensional matrix in the basis of so-called easily distinguished. Also the nature of the transition
“natural transition orbitals”, which correspond to nocc (n-π* or charge-transfer) is often readily apparent
particle-hole amplitudes. With each hole in the for simple molecular systems. However, the differ-
occupied space, one single corresponding particle in ence density is a complicated function with often
the virtual space can be associated by means of the intricate nodal surfaces, which makes its plotting and
associated eigenvalue λ. The importance of a par- analysis often very tedious especially for larger
ticular particle-hole excitation to an electronic ex- molecules. This is mostly because both the ground-
citation is reflected by the value of the corresponding state electron density that is removed upon excitation
eigenvalue. Usually, electronic transitions can be and the “new” electron density of the excited state
expressed by one single particle-hole pair in the are shown with different signs together in one
NTO basis with an associated eigenvalue of es- picture. It is also somewhat awkward that a density
sentially one, even such transitions that are highly acquires a negative sign since it is originally defined
mixed in the canonical molecular orbital basis. Plot- as the square of the wave function.
Calculation of Excited States of Large Molecules Chemical Reviews, 2005, Vol. 105, No. 11 4029

4.4. Attachment/Detachment Density Plots This density matrix A corresponds to the particle
levels occupied in the electronic transition. The
The physically most appealing and conceptually difference between the two new matrices A and D
easiest way to analyze the nature of a complicated corresponds to the original difference density matrix
electronic transition is via so-called attachment/ ∆:
detachment density plots.32,131 The basis of this
analysis is the diagonalization of the difference ∆)A-D (120)
density matrix ∆ given by eq 113 in section 4.3 via
In other words, the detachment density is that part
U†∆U ) δ (114) of the single-electron ground-state density that is
removed and rearranged as attachment density.
where U is a unitary transformation matrix contain- Together these densities characterize an electronic
ing the eigenvectors of the difference density matrix, transition as D f A, which permits the visualization
which again could be considered as “natural orbitals and analysis of electronic transitions more or less as
of the electronic transition” but which are generally if they correspond to just single-orbital replacements,
different from those in section 4.2. δ is the diagonal regardless of the extent of configuration mixing that
matrix containing the eigenvalues δp of ∆, which are occurs in the excited-state wave function and regard-
interpreted as occupation numbers of the eigenvec- less of how inappropriate the molecular orbitals are
tors. The difference density matrix is a square N × for the description of the transition. The analysis via
N dimensional matrix, where N is (nocc + nvirt). In attachment/detachment density plots can be applied
electronic transitions that do not involve ionization to any excited-state calculation that yields a one-
or electron attachment, the sum of all occupation particle difference density.
numbers must be zero; otherwise it corresponds to An application of attachment/detachment density
the net electron gain or loss of electrons n during the plots is given in section 5.1, where the excited states
transition, that is, of CO-ligated heme are analyzed (Figure 7). These
states are characterized by strong configuration
N
mixing of four excited Slater determinants, which
Tr(∆) ) ∑ δp ) n (115) makes an analysis in the molecular orbital basis
p)1
impossible. By means of attachment/detachment
In a next step, the diagonal difference density matrix plots, the nature of the repulsive states could be
δ is split into two matrices A and D. The matrix D, easily analyzed and understood in terms of an Fe-
the so-called detachment density, is defined as the CO back-bonding to anti-back-bonding transition.
sum of all eigenvectors of ∆ that possess negative Further applications of attachment/detachment plots
eigenvalues, weighted by their absolute value of their can be found, for example, in refs 132-135.
occupations. Consequently, if d is a diagonal matrix
with elements 5. Illustrative Examples
In this section, illustrative and educative examples
dp ) -min(δp,0) (116) for applications of CIS, TDHF, and TDDFT to large
molecular systems will be presented. First, an inves-
then tigation of the ultrafast photodissociation process of
CO-ligated heme employing TDDFT is described,
D ) UdU† (117) which is a well suited model study, since it demon-
The detachment density matrix D is a positive strates the necessity to perform a detailed investiga-
semidefinite N × N matrix, since the sign of the tion of the functional performance and it highlights
negative eigenvalues is changed, and thus, it has the usefulness of attachment/detachment density
positive entries everywhere, where ∆ has negative plots to analyze complicated electronic transitions.
ones. All other values of ∆ are set to zero. This Then, we will show TDDFT results for a zincbacte-
density corresponds to the electron density of the riochlorin-bacteriochlorin complex, where we are
ground state that is removed upon electron excitation interested in the spectral properties of the complex
and can thus be seen as a hole density, although not in comparison to the individual molecules. Here,
in general corresponding to removal of an integral severe problems with charge-transfer excited states
number of electrons. are encountered using TDDFT alone, and a practical
The attachment density matrix A is similarly hybrid approach combining the benefits of TDDFT
defined as the sum of all natural orbitals of the and CIS must be applied to obtain physically reason-
difference density matrix with positive occupation able results. During the whole section, emphasis is
numbers, weighted by the absolute value of their put on the theoretical approach and the choice of
occupation. Consequently, if a is a N × N diagonal method rather than on the chemically and physically
matrix with elements interesting results of the investigations, which can
be found in the original publications.
ap ) max(δp,0) (118)
5.1. The Initial Steps of the Ultrafast
then Photodissociation of CO-Ligated Heme
CO-ligated iron porphyrin (heme) is known to
A ) UaU† (119) undergo ultrafast dissociation of the CO ligand
4030 Chemical Reviews, 2005, Vol. 105, No. 11 Dreuw and Head-Gordon

Table 1. Calculated Geometrical Parameters (Å) of

the CO-Ligated Heme Model Complex in Comparison
with Experimentally Determined Values for the
Complete System
BLYP B3LYP expt
Fe-CO 1.81 1.80 1.73-1.93
Fe-NP 2.01 2.03 1.98-2.06
Fe-NIm 1.98 2.04 2.06-2.20
C-O 1.21 1.16 1.07-1.12

Table 2. Excitation Energies (eV) of the Q and B

States of the Model Complex Compared with the
Experimentally Observed Values
Figure 5. Model complex used for the theoretical inves- SVWN BLYP B3LYP
tigation of CO-ligated heme.
LANL2DZ 6-31G* LANL2DZ 6-31G* LANL2DZ 6-31G* expt

within 50 fs upon excitation of the system in the Qy 1.515 1.396 1.833 1.573 2.404 2.403 2.18
Qx 1.537 1.418 1.843 1.593 2.418 2.413 2.30
energetically lowest lying Q-states of the porphyrin Bx 3.009 2.909 3.049 2.927 3.347 3.314 2.96
ring.136 In a recent theoretical investigation employ- By 3.020 2.927 3.056 2.943 3.359 3.325 3.16
ing TDDFT, the electronically excited states of CO-
ligated heme were studied with the objective to
As second test of the functional performance, the
identify and characterize the electronic states that
excited states of the model complex were calculated
are involved in the ultrafast photodissociation pro-
using TDDFT employing the functionals from above,
cess.88,137 Within the calculation, a model complex
and the obtained values were compared with the
was used that reflects the structural characteristics
experimentally known values of the excitation ener-
of the heme group in the intact proteins hemoglobin
gies of the allowed transitions, the Q and B states.140
and myoglobin (Figure 5).
The Q and B states are known to be π-π* transitions
A crucial step in theoretical investigations employ- and can thus be expected to be reasonably well
ing TDDFT is the choice of an appropriate exchange- described by TDDFT. For these test calculations, the
correlation (xc) functional and a reasonable basis set. standard 6-31G* basis set and the LANL2DZ basis
Since one cannot expect that one approximate xc set were employed. The results of the test calcula-
functional describes all excited states of a molecule tions are summarized in Table 2.
equally well, it is important to test several xc func- In contrast to the geometrical parameters, which
tionals with respect to their performance for the varied only slightly when different xc functionals
system under investigation. In principle, calculated were tested, the excitation energies are sensitive to
properties are compared with their experimental the choice of the xc functional. While the local SVWN
values, for example, geometrical parameters or ex- and the gradient-corrected BLYP functional strongly
perimentally determined excitation energies of opti- underestimate the excitation energies of the Q states,
cally allowed transitions. In the presented investi- the B bands are in good agreement with the experi-
gation of the CO-ligated heme group, the geometrical mental values. The hybrid B3LYP functional, how-
parameters optimized with ground-state DFT, as well ever, overestimates the excitation energies of both
as the values of the experimentally observable elec- the Q states and the B states by about 0.2 eV. After
tronic transitions, the Q and B states, calculated with a closer look at the basis set dependence, the results
TDDFT, are compared with their experimental val- of the SVWN and BLYP calculations change strongly
ues. Three different xc functionals were tested: the when going from the LANL2DZ basis set to the
local Slater-Vosko-Wilk-Nusair (SVWN) func- 6-31G* basis set. This is not the case for the calcula-
tional,56,81 the gradient-corrected Becke-Lee-Yang- tions employing the B3LYP functional. Due to the
Parr (BLYP)82 functional, and the hybrid Becke3- consistent error of 0.2 eV for all tested states and the
Lee-Yang-Parr (B3LYP)80 functional. As first test, robustness of the B3LYP functional with respect to
the geometries of the complex were optimized using basis set change, the B3LYP functional has been
standard ground-state DFT with the three xc func- chosen to be the one used throughout the whole
tionals mentioned above and the LANL2DZ basis set, investigation of the photodissociation process of CO-
which makes use of the Los Alamos effective core ligated heme. This investigation clearly illustrates
potential for the inner electron shell of the iron atom the importance of a performance check of the xc
and uses the 6-31G basis set on all other lighter functionals prior to the actual investigation. As one
atoms. The calculated values for four selected pa- can see, it is not possible to pick one random
rameters (the Fe-CO bond, the bond between the functional and be guaranteed to obtain reliable
iron and the porphyrin nitrogens (Fe-NP), the bond results for excitation energies within the TDDFT
between the iron and the imidazole nitrogen (Fe- framework.
NIm), and the C-O bond of the CO ligand) are given To identify the relevant excited electronic states
in Table 1 and compared with their experimental for the ultrafast photodissociation process of CO-
values.138,139 It can be seen that the calculated values ligated heme, the lowest excited singlet states of the
depend only slightly on the quality of the xc func- system were calculated along the CO stretch coordi-
tional, that is, they are essentially independent of the nate employing TDDFT in combination with the
choice of the functional. B3LYP functional and the LANL2DZ basis set. The
Calculation of Excited States of Large Molecules Chemical Reviews, 2005, Vol. 105, No. 11 4031

Figure 6. Calculated potential energy curves along the

Fe-CO stretch coordinate for CO-ligated heme. The data Figure 7. Attachment/detachment density plots of the
have been obtained using TDDFT/B3LYP/LANL2DZ. repulsive states calculated at the equilibrium Fe-CO
distance of 1.8 Å (top) and at an Fe-CO bond length of 2.5
geometry of the system was not optimized along the Å (bottom) using TDDFT/B3LYP/LANL2DZ. The isosur-
dissociation pathway, since geometry relaxation ef- faces shown are calculated for a 90% density enclosure.
fects can be expected to be small for this ultrafast
process. The calculated potential energy curves are back-bonding). Furthermore, it can be seen from the
shown in Figure 6. attachment/detachment density plots of the repulsive
states at 2.5 Å Fe-CO bond distance (lower panel of
Upon photoexcitation of CO-ligated heme into the
Figure 7) that the character of this state changes
Q states (1A′ and 1A′′ in Figure 6), the system decays
along the Fe-CO stretch coordinate and that it gets
nonradiatively into two quasi-degenerate repulsive
more iron d-orbital character.
states. These states are identified as 3A′ and 5A′′ at
In summary, application of TDDFT to the ultrafast
the equilibrium geometry of CO-ligated heme. Along
photodissociation of CO-ligated heme has led to a
the dissociation coordinate, a small energy barrier
detailed understanding of the initial steps of this
of 0.15 eV must be crossed, but this is probably easily
complicated process. Furthermore, the presented
accomplished by vibrational excitation going along
investigation can serve as a good example for a
with the photoexcitation. Calculation of the vibra-
typical TDDFT study, since the excited electronic
tional frequencies of the electronic ground state with
states under investigation are typical valence-excited
DFT/B3LYP/6-31G* corroborates this picture, since
states, which can be expected to be reasonably well
the period of the vibration associated with the Fe-
described with TDDFT. A detailed study of the
CO stretch coordinate is calculated to be about 70
dependence of the results on different xc functionals
fs, which strongly correlates with the time scale of
has been performed, which is a crucial step in any
the photodissociation process of 50 fs.
theoretical investigation employing TDDFT. The
Since the repulsive excited states involve linear usefulness of density-based analyses, here attach-
combinations of several molecular orbitals and are, ment/detachment density plots, has been demon-
thus, difficult to analyze in terms of orbitals, the strated.
nature of the repulsive states has been analyzed
using attachment/detachment density plots32 (see
also section 4.4) at the equilibrium geometry, as well
5.2. Charge-Transfer Excited States in
as at an Fe-CO bond length of 2.5 Å. The calculated
Zincbacteriochlorin−Bacteriochlorin Complexes
plots are displayed in Figure 7. Comparison of the Despite its failure for CT states as outlined in
detachment and attachment densities nicely explains section 3.5, TDDFT has recently been applied to
the repulsive character of the 5A′′ and 3A′ states at several molecular complexes in which photoinitiated
the equilibrium geometry (upper panel of Figure 7). electron transfer might occur and, thus, long-range
While the detachment density, which is removed CT excited states are relevant. In this spirit, phe-
from the ground-state density, is clearly dominated nylene-linked free-base porphyrin and zincporphyrin
by a bonding iron-carbon interaction, the attach- complexes, as well as their bacteriochlorin (BC)
ment density, which is added to the ground-state analogues, have been investigated employing TDDFT
density, has clearly antibonding character, which is in combination with the BLYP functional,141,142 and
seen as a node along the Fe-C bond. The bonding a plethora of spurious CT states was obtained in the
interaction between the iron atom and the CO ligand visible range of the electronic spectra of these com-
in the detachment density can be understood in pounds. In contrast, in a similar study of linked
chemical terms as back-bonding from an iron d- zincporphyrin dimers in which the symmetry-adapted
orbital into the π* orbital of CO. The antibonding cluster configuration interaction (SAC-CI) method27
interaction in the attachment density corresponds to has been employed, CT excited states were not found
the antibonding combination of these orbitals (anti- in the energy regime of the energetically lowest Q
4032 Chemical Reviews, 2005, Vol. 105, No. 11 Dreuw and Head-Gordon

Table 3. Comparison of the Energies of the Ten

Lowest Singlet Excited States of the Full
Phenylene-Linked ZnBC-BC Complex as Well as of
the Model Complex without the Phenylene Bridge
with the Individually Calculated and Experimentally
Determined Excitation Energies of the Q States of
the Monomersa
ZnBC-BC complex monomers
state full model transition calcd expt
1 1.33 (0.000) 1.32 (0.000) ZnBC f BC CT
2 1.46 (0.000) 1.47 (0.000) BC f ZnBC CT
3 1.86 (0.000) 1.90 (0.000) BC f ZnBC CT
4 1.94 (0.001) 1.96 (0.000) ZnBC f BC CT
5 2.05 (0.393) 2.07 (0.266) π-π* ZnBC (Qx) 2.07 (0.231) 1.65b
6 2.09 (0.131) 2.12 (0.170) π-π* BC (Qx) 2.10 (0.187) 1.6c
7 2.38 (0.059) 2.40 (0.038) π-π* BC (Qy) 2.39 (0.034) 2.3c
8 2.42 (0.019) 2.46 (0.018) π-π* ZnBC (Qy) 2.44 (0.026) 2.2b
9 2.43 (0.022) 2.42 (0.000) ZnBC f BC CT
10 2.58 (0.000) 2.66 (0.000) BC f ZnBC CT
a
Figure 8. Molecular structure of the (1,4)-phenylene- The oscillator strength of each transition is given in
linked zincbacteriochlorin-bacteriochlorin complex, as well brackets behind the corresponding energy. All calculations
as of the model complex used in some calculations. The have been performed at the level of TDDFT/BLYP/6-31G* and
distance coordinate R is here defined as the distance all energies are given in eV. b Experimental data for zinctetra-
between the formerly linked carbon atoms. phenylbacteriochlorin-pyridine were taken from ref 154.
c
Experimental data for bacteriopheophorbide were taken from
ref 155.
states of zincporphyrins,143 which is in accordance
with experimental findings (see ref 144 and refer-
ences herein). Typically, charge-transfer excited states
are experimentally identified via polarity dependent
solvatochromic shifts of their absorption band or by
means of ultrafast transient absorption spectros-
copy.145-147 In this section, we are going to outline a
recent reinvestigation of the excited states of a
phenylene-linked zincbacteriochlorin-bacteriochlorin
complex (Figure 8) employing TDDFT102 for two
reasons. On one hand, we want to demonstrate the
substantial failure of TDDFT employing standard xc
functionals for CT excited states, and on the other
hand, the usefulness of the practical scheme to
circumvent the problem introduced in here should be
pointed out. Figure 9. Highest occupied molecular orbital (HOMO,
bottom) and lowest unoccupied molecular orbital (LUMO,
As a first step, the excited states of the monomers top) of the ZnBC-BC complex. In the lowest excited CT
zincbacteriochlorin and bacteriochlorin have been state, an electron is transferred from the HOMO to the
calculated employing TDDFT with the BLYP82 xc LUMO, and thus, this transition corresponds to an electron
functional and the 6-31G* basis set at the separately transfer from ZnBC to BC.
optimized geometries. At this level of theory, ZnBC
possesses weakly allowed singlet excited states at energetically lowest state at 1.33 eV is a pure one-
2.07 and 2.44 eV, as well as strongly allowed transi- electron transition from the highest occupied molec-
tions at 3.62 and 3.88 eV, which are interpreted as ular orbital (HOMO) into the lowest unoccupied
the well-known Qx and Qy, as well as By and Bx, molecular orbital (LUMO), and analysis of these MOs
states, respectively. The corresponding states of (Figure 9) clearly shows that this state corresponds
bacteriochlorin are found at energies of 2.10 (Qx), 2.39 to an electron transfer from ZnBC to BC. In analogy,
(Qy), 3.67 (By), and 3.85 eV (Bx). Our calculated as states 2 and 3 at 1.46 and 1.86 eV represent BC-to-
well as experimental values for the excitation ener- ZnBC CT states, while state 4 at 1.94 eV is again a
gies of the Q states are given in Table 3. ZnBC-to-BC CT state. Of course, these states cannot
be present if the monomers are calculated individu-
All values are in reasonable agreement with other
ally but are a characteristic of complex formation.
calculations of the electronic absorption spectra of
ZnBC and BC.148-150 Using the same theoretical However, as pointed out in section 3.5, the calcu-
approach for the calculation of the 10 lowest singlet lated values of the CT states are given at too low
excited states of the full (1,4)-phenylene-linked ZnBC- energies, which can be easily confirmed by simple
BC complex (Figure 8), one obtains the values given electrostatic considerations. Assuming that the sepa-
in Table 3. As expected, the Q-states of the constitut- rated charges in the CT states could be treated as
ing monomers are found at almost exactly identical point charges, the distance-dependent excitation
energies in the linked complex. They are slightly red- energy of the energetically lowest ZnBC-to-BC CT
shifted by only 0.01-0.02 eV. In addition to the state ωCT(R) can simply be estimated via
monomer states, further energetically low-lying states
are found at 1.33, 1.46, 1.86, and 1.94 eV. The ωCT(R) > IPZnBC + EABC - 1/R (121)
Calculation of Excited States of Large Molecules Chemical Reviews, 2005, Vol. 105, No. 11 4033

where IPZnBC is the ionization potential of the ZnBC

monomer, EABC is the electron affinity of bacterio-
chlorin, and 1/R is the electrostatic attraction be-
tween them. Of course, R corresponds to the distance
between the charges, which in our approximation is
chosen to be the smallest distance between the
carbon atoms of the different tetrapyrrol rings in the
(1,4)-phenylene-linked ZnBC-BC complex, which has
a value of 5.84 Å (see also Figure 8). At the level of
DFT/BLYP/6-31G*, IPZnBC has a value of 5.57 eV and
EABC is given as -0.42 eV yielding an excitation
energy for the lowest ZnBC-to-BC CT state of about
2.69 eV at the equilibrium distance R of the linked
complex. In this equation, the cation and anion are
treated as point charges and the shortest possible
distance R is assumed, which of course leads to an
overestimation of the electrostatic attraction. As a
consequence, the estimated excitation energy ωCT is
a true lower bound to the correct value. Comparison
with the excitation energy calculated at the level of
TDDFT/BLYP/6-31G*, which has a value of 1.33 eV,
yields a minimum error of 1.36 eV for this CT state.
This demonstrates the tendency of TDDFT to under-
estimate the excitation energies of CT states drasti-
cally. It is worthwhile to note that the excitation
energies of the CT states computed with TDDFT
correspond in fact exactly to the difference of the
orbital energies involved in the excitations, as we
have already discussed in section 3.5. For instance,
the excitation energy of the lowest CT state (1.33 eV)
is solely given by the difference of the HOMO (-3.70
eV) and LUMO (-2.37 eV) energies of the electronic
ground state.
As the next step, the dependence of the excitation Figure 10. Potential energy curves of the lowest singlet
energies of the lowest CT states on the distance R excited states of the ZnBC-BC complex along the distance
between the separated charges will be investigated. coordinate R calculated with TDDFT/BLYP/6-31G* (top)
This requires the introduction of a suitable model and with the hybrid TDDFT/CIS approach (bottom, see text
complex, in which the distance between the positive for details). Short-dashed lines correspond to CT excited
states, solid lines are the valence-excited Q states, and the
and negative charges of the CT state can be easily long-dashed line represents the electronic ground state.
varied. For the (1,4)-phenylene-linked ZnBC-BC
complex, this is easily accomplished by neglecting the energies correspond to the orbital energy differences),
phenylene bridge, and choosing the distance between and they do not show the correct 1/R dependence.
the formerly connected carbon atoms as the distance This is due to the electron-transfer self-interaction,
coordinate R, which has a value of 5.84 Å in the which is not canceled in the employed pure BLYP
linked complex and was used previously to estimate functional as we have outlined in section 3.5. In the
the minimum excitation energy, ωCT, in eq 121. As lower part of Figure 10, potential energy curves are
can be seen in Table 3, the phenylene bridge has only shown that are obtained with a hybrid approach of
minor influence on the 10 lowest excited states of the TDDFT and CIS that will be described in detail
complex, which are slightly shifted by an average of within the next paragraph.
0.03 eV, while the largest shift of 0.08 eV occurred In contrast to TDDFT methods, CIS and TDHF
for state 10. As a consequence, the model complex yield the correct 1/R behavior of the potential energy
can be used in further calculations without losing curves of CT states with regard to a distance coor-
general validity of the obtained results. dinate, because of the full inclusion of HF exchange,
To investigate the asymptotic behavior of the which leads to the cancellation of electron-transfer
excited states and, in particular, of the excited CT self-interaction. On the other hand, the excitation
states of the ZnBC-BC complex, the 10 lowest energies calculated with CIS or TDHF are usually
excited states of the model complex (Figure 8) have much too large owing to the calculation of the
been calculated employing TDDFT with the BLYP energies of the virtual orbitals for the (N + 1)-electron
functional and the 6-31G* basis set along the dis- system instead of for the N-electron system. There-
tance coordinate R. The obtained potential energy fore, we suggest using a hybrid scheme that combines
curves (PECs) are displayed in the upper part of the benefits of both methods to obtain reasonable
Figure 10. As can be easily seen, the potential energy estimates for the energies and potential energy
curves of the CT states are constant along the surfaces of CT states relative to valence-excited
distance R between the monomers (the excitation states. This approach can only be seen as a practical
4034 Chemical Reviews, 2005, Vol. 105, No. 11 Dreuw and Head-Gordon

work-around to the CT problem in TDDFT until a lecular CT states compared to valence-excited states
more theoretically sound solution is found. However, and to obtain reasonable potential energy curves
it is a useful “diagnostic” for whether CT states are along an intermolecular separation coordinate. The
in the range of the valence excitations of interest. As hybrid approach has already been successfully ap-
the first step of this hybrid approach, a ground-state plied to CT states in large xanthophyll-chlorophyll
DFT calculation is performed for the energetically dimers.87,111
lowest CT state at a large intermolecular distance In Figure 10, the electron-transfer self-interaction-
(R0) by exchanging the orbital i located at molecule free PECs of the two energetically lowest CT states
A against the orbital a located at molecule B in the of the ZnBC-BC complex are shown that have been
â-part of the wave function. To converge the DFT calculated along the distance coordinate R with CIS
calculation onto the energetically lowest CT state, we and shifted by a ∆DFT offset. This offset has been
generally use the geometric direct minimization calculated for the lowest ZnBC-to-BC CT state to be
technique (GDM).151 At shorter distances when the 3.79 eV and for the lowest BC-to-ZnBC CT state to
orbitals start to overlap, GDM converges to the be 3.95 eV at the level of BLYP/6-31G*. These curves
closed-shell ground state of the dimer, presumably are then plotted together with the curves of the
because the CT solution is no longer a distinct wave Q-states given by the TDDFT/BLYP/6-31G* calcula-
function minimum. Otherwise this would provide a tion. The obtained picture is substantially different
convenient way to map out a complete potential from the one obtained by TDDFT alone. First of all,
energy surface of the CT state. A corrected excitation the two energetically lowest CT states are clearly
energy for the lowest CT state is then easily obtained above the valence-excited Q-states of the complex.
by subtraction of the total energies of the ground Second, the CT states do exhibit the correct asymp-
state and the CT state, which corresponds to the totic behavior such that the excitation energies do
known ∆DFT method.152 This excitation energy does increase with 1/R. Furthermore, from this picture one
not suffer from electron-transfer self-interaction. The can obtain a reasonable value for the lowest intramo-
self-interaction-free long-range value of the CT ex- lecular CT states of the 1,4-phenylene-linked ZnBC-
citation energy is then used as a constant offset (for BC complex. At the value of 5.84 Å, the value of R in
exchange-correlation effects) for the asymptotically the full complex, the lowest CT state possesses an
correct potential energy curve calculated at the CIS excitation energy of 3.75 eV and the second lowest
level: CT state has an excitation energy of 3.91 eV, which
agrees with the previously computed lower bound for
ωCT(R) ) ωCIS CIS the excitation energies of these states. Compared
CT (R) + (∆DFT(R0) - ωCT (R0)) (122)
with the TDDFT computed values of 1.33 and 1.46
eV for these states, this yields errors in the TDDFT
The quality of the total energy of the CT state is calculation of 2.42 and 2.43 eV, respectively. As a
roughly the same as the one of the ground state, since consequence, the spectrum as obtained by TDDFT
the energetically lowest CT state usually corresponds alone and published previously142 is an artifact of the
to a pure occupied-virtual single-electron transition approximate xc functionals employed in present-day
over the considered distance range, that is, the TDDFT.
excited state is well described by a single Slater
determinant. Furthermore, in the limit of the exact
xc functional, TDDFT and ∆DFT are equivalent. 6. Brief Summary and Outlook
Plotting the shifted CIS curve together with the At present, time-dependent density functional theory
curves of the valence-excited states calculated with (TDDFT) is the most prominent method for the
TDDFT yields a complete self-interaction-free picture calculation of excited states of medium-sized and
of all relevant excited states of the dimer. We have large molecules. Every week, a large number of
investigated the error that is introduced by using the publications appear that present successful applica-
CIS curve instead of the ∆DFT curve by studying the tions of the method in fields such as inorganic,
energetically lowest CT state of a tetrafluoroethyl- organic, and physical chemistry. With the enormous
ene-ethylene dimer. In this small system, it is advance in computer technology, systems up to 300
possible to converge the ground-state DFT calculation second-row atoms became tractable, which even
at even shorter distances onto the CT state with the allows for the treatment of excited-state problems in
help of a maximum-overlap method (MOM)153 and the fields of biochemistry, biophysics, and material
thereby to map out a complete potential energy sciences.
surface for the CT state. We found that by using the TDDFT is a formally exact theory that relies on
CIS curve with the self-interaction-free offset value the analysis of the time-dependent linear response
instead of the ∆SCF curve one introduces an ad- of the formally exact ground-state density to a time-
ditional error of at most 0.1 eV. Therefore, it can be dependent external perturbation, which after Fourier
expected that the accuracy of the hybrid approach transformation yields exact excited-state energies
for CT states is of about the same order as TDDFT and oscillator strengths. The derivation of the famous
is for valence-excited states. Runge-Gross theorem and the subsequent formula-
The proposed hybrid approach is of course not tion of a time-dependent Kohn-Sham equation were
useful for small systems, since more reliable wave- the cornerstones in the development of the TDDFT
function-based methods are available, but for large formalism. However, since the exact xc functional is
systems, this approach is at present a viable way to not known, approximate xc functionals need to be
gain insight into the energetic positions of intermo- employed in a practical calculation. Usually the
Calculation of Excited States of Large Molecules Chemical Reviews, 2005, Vol. 105, No. 11 4035

adiabatic local density approximation (ALDA) is ously provided by the National Energy Research
employed, and standard time-independent xc func- Scientific Computing Center (NERSC).
tionals derived for ground-state DFT are used. Con-
comitantly, errors in the excitation energies and 8. References
oscillator strengths are introduced. Still, for most
(1) Szabo, A.; Ostlund, N. S. Modern Quantum Chemistry; Mac-
valence-excited states, which lie well below the first Millan Publishing: New York, 1982.
ionization potential, TDDFT yields results with high (2) McWeeny, R.; Sutcliffe, B. T. Methods of Molecular Quantum
accuracy at relatively low computational cost in Mechanics; Academic Press: London, 1969.
(3) Møller, C.; Plesset, M. S. Phys. Rev. 1934, 46, 618.
comparison with highly accurate methods such as (4) Bartlett, R. J. Annu. Rev. Phys. Chem. 1981, 32, 359.
MRCI, CASPT2, or EOM-CCSD, which are applicable (5) Carsky, P. In Encyclopedia of Computational Chemistry; Schley-
only to small molecules up to 20 atoms. In fact, the er, P. v. R., Clark, N. L., Gasteiger, J., H. F. S., III, Schreiner,
P. R., Eds.; Wiley: Chichester, U.K., 1998; p 485.
quality of the excitation energies often lies within (6) Kutzelnigg, W. J. Mol. Struct. (THEOCHEM) 1988, 181, 33.
0.1-0.5 eV compared with experimental data. Nev- (7) Gauss, J. In Encyclopedia of Computational Chemistry; Schleyer,
P. v. R., Clark, N. L., Gasteiger, J., H. F. S., III, Schreiner, P.
ertheless, one has to be very careful using TDDFT R., Eds.; Wiley: Chichester, U.K., 1998; p 615.
with approximate xc functionals owing to its failures (8) Bartlett, R. J., Ed. Modern Ideas in Coupled-Cluster Methods;
for Rydberg states, systems with large π-systems, World Scientific: Singapore, 1997.
(9) Parr, R. G.; Yang, W. Density-Functional Theory of Atoms and
doubly excited states, and CT states. The latter Molecules; Oxford Science Publication: New York, 1989.
failure limits the applicability of TDDFT to large (10) Dreizler, R. M.; Gross, E. K. U. Density functional theory;
systems or small molecules in solution or protein Springer-Verlag: Heidelberg, Germany, 1995.
(11) Baerends, E. J.; Gritsenko, O. V. J. Phys. Chem. A 1997, 101,
environments dramatically, because erroneous inter- 5383.
and intramolecular CT excited states occur in the (12) Hohenberg, P.; Kohn, W. Phys. Rev. 1964, B136, 864.
(13) Kohn, W.; Sham, L. J. Phys. Rev. 1965, A140, 1133.
low-energy region of the electronic spectra. At present, (14) Buenker, R. J.; Peyerimhoff, S. D.; Butscher, W. Mol. Phys. 1978,
this prevents TDDFT from being a black-box method 35, 771.
for the calculation of excited-state properties in the (15) Buenker, R. J.; Peyerimhoff, S. D.; Bruna, P. J. In Computational
Theoretical Organic Chemistry; Csizmadia, I. G., Daudel, R.,
same fashion that time-independent DFT has become Eds.; Reidel: Dordrecht, The Netherlands, 1981; p 91.
a standard tool for the study of electronic ground (16) McDonall, J. J.; Peasley, K.; Robb, M. A. Chem. Phys. Lett. 1988,
states. 148, 183.
(17) Andersson, K.; Malmqvist, P.-Å.; Roos, B. O. J. Chem. Phys.
In view of these findings, much research is still 1992, 96, 1218.
(18) Roos, B. O. Adv. Chem. Phys. 1987, 69, 399.
dedicated to the improvement and development of (19) Andersson, K.; Roos, B. O. In Modern Electronic Structure
new xc functionals to eliminate the known failures Theory; Yarkony, D. R., Ed.; World Scientific: New York, 1995;
of TDDFT. However, a remaining open question is Vol. 1, p 55.
(20) Emrich, K. Nucl. Phys. A 1981, 351, 379.
whether there exists an approximate xc functional (21) Sekino, H.; Bartlett, R. J. Int. J. Quantum Chem. Symp. 1984,
that can describe all excited states equally well. Since 18, 255.
different excited states can possess very different (22) Geertsen, J.; Rittby, M.; Bartlett, R. J. Chem. Phys. Lett. 1989,
164, 57.
electronic structures, it seems unlikely that they can (23) Monkhorst, H. J. Int. J. Quantum Chem. Symp. 1977, 11, 421.
all be captured by one simple approximate xc func- (24) Dalgaard, E.; Monkhorst, H. J. Phys. Rev. A 1983, 28, 1217.
tional. An approximation may be an excellent one for (25) Koch, H.; Christiansen, O.; Jørgensen, P. Chem. Phys. Lett. 1995,
244, 75.
one class of excited states but at the same time a very (26) Christiansen, O.; Gauss, J.; Schimmelpfennig, B. Phys. Chem.
poor one for another equally important class of Chem. Phys. 2000, 2, 965.
(27) Nakatsuji, H.; Hirao, K. J. Chem. Phys. 1978, 68, 2053.
excited states. From this point of view, the discovery (28) Gromov, E. V.; Trofimov, A. B.; Vitkovskaja, N. M.; Schirmer,
of the “one” approximate xc functional that solves all J.; Köppel, H. J. Chem. Phys. 2003, 119, 737.
the problems seems unlikely; at least it poses a huge (29) Trofimov, A. B.; Stelter, G.; Schirmer, J. J. Chem. Phys. 2002,
117, 6402.
challenge to leading experts in the field. (30) Schirmer, J. Phys. Rev. A 1982, 26, 2395.
After its birth in 1984, TDDFT has now left its (31) Head-Gordon, M.; Rico, R. J.; Oumi, M.; Lee, T. J. Chem. Phys.
Lett. 1994, 219, 21.
childhood behind and advanced to a juvenile method (32) Head-Gordon, M.; Grana, A. M.; Maurice, D.; White, C. A. J.
that has proven to be a useful alternative approach Phys. Chem. 1995, 99, 14261.
(33) Christiansen, O.; Koch, H.; Jørgensen, P. Chem. Phys. Lett. 1995,
to standard wave-function-based methods for the 243, 409.
calculation of excited states. Whether its indicated (34) Haettig, C.; Weigend, F. J. Chem. Phys. 2000, 113, 5154.
failures are remainders of children’s diseases that can (35) del Bene, J.; Ditchfield, R.; Pople, J. A. J. Chem. Phys. 1971,
55, 2236.
be fully cured in the future and whether TDDFT (36) Raghavachari, K.; Trucks, G. W.; Pople, J. A.; Head-Gordon, M.
thereby possesses the potential to advance to a Chem. Phys. Lett. 1989, 157, 479.
standard black-box research tool for excited states (37) Schirmer, J.; Trofimov, A. B. J. Chem. Phys. 2004, 120, 11449.
(38) Fetter, A. L.; Walecka, J. D. Quantum Theory of Many-Particle
as its older cousin DFT has become for ground states, Systems; McGraw-Hill: New York, 1971.
the near future will show. (39) Runge, E.; Gross, E. K. U. Phys. Rev. Lett. 1984, 52, 997.
(40) Gross, E. K. U.; Kohn, W. Phys. Rev. Lett. 1985, 55, 2850.
(41) Gross, E. K. U.; Kohn, W. Adv. Quantum Chem. 1990, 21, 255.
7. Acknowledgments (42) Casida, M. E. In Recent Advances in Density Functional Methods,
Part I; Chong, D. P., Ed.; World Scientific: Singapore, 1995; p
155.
A. D. gratefully acknowledges financial support by (43) van Leeuwen, R. Int. J. Mod. Phys. B 2001, 15, 1969.
the Deutsche Forschungsgemeinschaft as an “Emmy (44) Marques, M. A. L.; Gross, E. K. U. Annu. Rev. Phys. Chem. 2004,
55, 427.
Noether” fellow. This work was also supported by the (45) Foresman, J. B.; Head-Gordon, M.; Pople, J. A.; Frisch, M. J. J.
Director, Office of Energy Research, Office of Basic Phys. Chem. 1992, 96, 135.
Energy Sciences, Chemical Science Division, of the (46) Maurice, D.; Head-Gordon, M. Mol. Phys. 1999, 96, 1533.
(47) Kong, J.; White, C. A.; Krylov, A. I.; Sherrill, D.; Adamson, R.
U.S. Department of Energy under Contract No DE- D.; Furlani, T. R.; Lee, M. S.; Lee, A. M.; Gwaltney, S. R.; Adams,
AC03-76SF00098. Computer time has been gener- T. R.; Ochsenfeld, C.; Gilbert, A. T. B.; Kedziora, G. S.; Rassolov,
4036 Chemical Reviews, 2005, Vol. 105, No. 11 Dreuw and Head-Gordon

V. A.; Maurice, D. R.; Nair, N.; Shao, Y.; Basley, N. A.; Maslen, (99) Maitra, N. T.; Zhang, F.; Cave, R. J.; Burke, K. J. Chem. Phys.
P. E.; Dombroski, J. P.; Daschel, H.; Zhang, W.; Korambath, P. 2004, 120, 5932.
P.; Baker, J.; Byrd, E. F. C.; Van Voorhis, T.; Oumi, M.; Hirata, (100) Tozer, D. J.; Amos, R. D.; Handy, N. C.; Roos, B. J.; Serrano-
S.; Hsu, C.-P.; Ishikawa, N.; Florian, J.; Warshel, A.; Johnson, Andres, L. Mol. Phys. 1999, 97, 859.
B. G.; Gill, P. M. W.; Head-Gordon, M.; Pople, J. A. J. Comput. (101) Dreuw, A.; Weisman, J. L.; Head-Gordon, M. J. Chem. Phys.
Chem. 2000, 21, 1532. 2003, 119, 2943.
(48) CADPAC6: The Cambridge Analytic Derivatives Package Issue (102) Dreuw, A.; Head-Gordon, M. J. Am. Chem. Soc. 2004, 126, 4007.
6, Cambridge, England, 1992. A suite of quantum chemistry (103) Sobolewski, A. L.; Domcke, W. Chem. Phys. 2003, 294, 73.
programs developed by R. D. Amos with contributions from (104) van Leeuwen, R.; Baerends, E. J. Phys. Rev. A 1994, 49, 2421.
Alberts, I. L.; Andrews, J. S.; Colwell, S. M.; Handy, N. C.; (105) Grüning, M.; Gritsenko, O. V.; van Gisbergen, S. J. A.; Baerends,
Jayatilaka, D.; Knowles, P. J.; Kobayashi, R.; Koga, N.; Laidig, E. J. J. Chem. Phys. 2001, 114, 652.
K. E.; Maslen, P. E.; Murray, C. W.; Rice, J. E.; Sanz, J.; (106) Görling, A. Phys. Rev. Lett. 1999, 83, 5459.
Simandiras, E. D.; Stone, A. J.; Su, M.-D. (107) Ivanov, S.; Hirata, S.; Bartlett, R. J. Phys. Rev. Lett. 1999, 83,
(49) Alrichs, R. TURBOMOLE, version 5.7.1; University of Karlsru- 5455.
he: Germany, 2004. (108) Della Sala, F.; Görling, A. Int. J. Quantum Chem. 2002, 91, 131.
(50) Hirata, S.; Head-Gordon, M.; Bartlett, R. J. J. Chem. Phys. 1999, (109) Pople, J. A.; Binkley, J. S.; Seeger, R. Int. J. Quantum Chem.
111, 10774. Symp. 1976, 10, 1.
(51) Stanton, J. F.; Gauss, J.; Ishikawa, N.; Head-Gordon, M. J. (110) Van Caillie, C.; Amos, R. D. Chem. Phys. Lett. 1999, 308, 249.
Chem. Phys. 1995, 103, 4160. (111) Dreuw, A.; Fleming, G. R.; Head-Gordon, M. J. Phys. Chem. B
(52) Thomas, W. Naturwissenschaften 1925, 13, 627. 2003, 107, 6500.
(53) Reiche, F.; Thomas, W. Z. Phys. 1925, 34, 510. (112) Tozer, D. J. J. Chem. Phys. 2005, 119, 12697.
(54) Kuhn, W. Z. Phys. 1925, 33, 408. (113) Perdew, J. P.; Parr, R. G.; Levy, M.; Balduz, J. L., Jr. Phys. Rev.
(55) Davidson, E. R. J. Comput. Phys. 1975, 17, 87. Lett. 1982, 49, 1691.
(56) Dirac, P. A. M. Proc. Cambridge Philos. Soc. 1930, 26, 376. (114) Tawada, Y.; Tsuneda, T.; Yanagisawa, S.; Yanai, T.; Hirao, K.
(57) Frenkel, J. Wave Mechanics, Advanced General Theory; Clar- J. Chem. Phys. 2004, 1210, 8425.
endon Press: Oxford, U.K., 1934. (115) Yanai, T.; Tew, D. P.; Handy, N. C. Chem. Phys. Lett. 2004, 393,
(58) Heinrichs, J. Chem. Phys. Lett. 1968, 2, 315. 51.
(59) Ring, P.; Schuck, P. The nuclear many-body problem; Springer- (116) Baer, R.; Neuhauser, D. Phys. Rev. Lett. 2005, 94, 043002.
Verlag: Heidelberg, Germany, 1980. (117) Stoll, H.; Savin, A. In Density Functional Methods in Physics;
(60) Thouless, D. J. The Quantum Mechanics of Many Body Systems; Dreizler, R. M., da Providencia, J., Eds.; Plenum: New York,
Academic Press: New York, 1972. 1985; p 177.
(61) Sekino, H.; Bartlett, R. J. J. Chem. Phys. 1986, 85, 976. (118) Savin, A. In Recent Developments and Applications of Modern
(62) Čižek, J.; Paldus, J. J. Chem. Phys. 1967, 47, 3976. Density Functional Theory; Seminario, J. M., Ed.; Elsevier:
(63) McLachlan, A. D.; Ball, M. A. Rev. Mod. Phys. 1964, 36, 844. Amsterdam, 1996.
(64) Seeger, R.; Pople, J. A. J. Chem. Phys. 1977, 66, 3045. (119) Leininger, T.; Stoll, H.; Werner, H.-J.; Savin, A. Chem. Phys.
(65) Van Caillie, C.; Amos, R. D. Chem. Phys. Lett. 2000, 317, 159. Lett. 1997, 275, 151.
(66) Furche, F.; Ahlrichs, R. J. Chem. Phys. 2002, 117, 7433. (120) Iikura, H.; Tsuneda, T.; Yanai, T.; Hirao, K. J. Chem. Phys. 2001,
(67) Eschrig, H. The Fundamentals of Density Functional Theory; 115, 5340.
Teubner Verlag: Stuttgart, Germany, 1996. (121) van Faassen, M.; de Boeij, P. L.; van Leeuwen, R.; Berger, J. A.;
(68) Koch, W.; Holthausen, M. C. A Chemist’s Guide to Density Snijders, J. G. J. Chem. Phys. 2003, 118, 1044.
Functional Theory; Wiley-VCH: Weinheim, Germany, 2000. (122) Baerends, E. J.; Autschbach, J.; Brces, A.; Bo, C.; Boerrigter, P.
(69) Baym, G., Ed. Lecture Notes on Quantum Mechanics; Perseus M.; Cavallo, L.; Chong, D. P.; Deng, L.; Dickson, R. M.; Ellis, D.
Publishing: Cambridge, Massachusetts 1974. E.; Fan, L.; Fischer, T. H.; Fonseca Guerra, C.; van Gisbergen,
(70) Atkins, P. W. Molecular Quantum Mechanics; Oxford University S. J. A.; Groeneveld, J. A.; Gritsenko, O. V.; Grning, M.; Harris,
Press: Oxford, U.K., 1983. F. E.; van den Hoek, P.; Jacobsen, H.; van Kessel, G.; Kootstra,
(71) van Leeuwen, R. Phys. Rev. Lett. 1999, 82, 3863. F.; van Lenthe, E.; McCormack, D. A.; Osinga, V. P.; Patchk-
(72) Yabana, K.; Bertsch, G. F. Phys. Rev. B 1996, 54, 4484. ovskii, S.; Philipsen, P. H. T.; Post, D.; Pye, C. C.; Ravenek, W.;
(73) Marques, M. A. L.; Castro, A.; Bertsch, G. F.; Rubio, A. Comput. Ros, P.; Schipper, P. R. T.; Schreckenbach, G.; Snijders, J. G.;
Phys. Commun. 2003, 151, 60. Sola, M.; Swart, M.; Swerhone, D.; te Velde, G.; Vernooijs, P.;
(74) Castro, A.; Marquez, M. A. L.; Alonso, J. A.; Bertsch, G. F.; Versluis, L.; Visser, O.; van Wezenbeek, E.; Wiesenekker, G.;
Yabana, K.; Rubio, A. J. Chem. Phys. 2002, 116, 1930. Wolff, S. K.; Woo, T. K.; Ziegler, T. ADF2004.01, SCM, Theoreti-
(75) Marques, M. A. L.; Lopez, X.; Varsano, D.; Castro, A.; Rubio, A. cal Chemistry; Vrije Universiteit: Amsterdam, The Netherlands,
Phys. Rev. Lett. 2003, 90, 258101. http://www.scm.com.
(76) Hirata, S.; Head-Gordon, M. Chem. Phys. Lett. 1999, 302, 375. (123) Martin, R. L. J. Chem. Phys. 2003, 118, 4775.
(77) Petersilka, M.; Grossmann, U. J.; Gross, E. K. U. Phys. Rev. Lett. (124) Amos, A. T.; Hall, G. G. Proc. R. Soc. A 1961, 263, 483.
1996, 76, 1212. (125) Batista, E. R.; Martin, R. L. J. Phys. Chem. A 2005, 109, 3128.
(78) Hirata, S.; Head-Gordon, M. Chem. Phys. Lett. 1999, 314, 291. (126) Clark, A. E.; Martin, R. L.; Hay, P. J.; Green, J. C.; Jantunen,
(79) Hsu, C.-P.; Hirata, S.; Head-Gordon, M. J. Phys. Chem. A 2001, K. C.; Kiplinger, J. L. J. Phys. Chem. A 2005, 109, 5481.
105, 451. (127) Climent, T.; Gonzalez-Luque, R.; Merchan, M. J. Phys. Chem.
(80) Becke, A. D. J. Chem. Phys. 1993, 98, 5648. A 2003, 107, 6995.
(81) Vosko, S. H.; Wilk, L.; Nusair, M. Can. J. Phys. 1980, 58, 1200. (128) Serrano-Andres, L.; Merchan, M.; Jablonski, M. J. Chem. Phys.
(82) Becke, A. D. Phys. Rev. A 1988, 38, 3098. 2003, 119, 4294.
(83) Lee, C.; Yang, W.; Parr, R. G. Phys. Rev. B 1988, 37, 785. (129) Molina, V.; Merchan, M. J. Phys. Chem. A 2001, 105, 3745.
(84) Perdew, J. P.; Burke, K.; Enzerhof, M. Phys. Rev. Lett. 1996, (130) Wiberg, K. B.; Hadad, C. M.; Foresman, J. B.; Chupka, W. A. J.
77, 3865. Phys. Chem. 1992, 96, 10756.
(85) Perdew, J. P.; Burke, K.; Enzerhof, M. Phys. Rev. Lett. 1997, (131) Grana, A. M.; Lee, T. J.; Head-Gordon, M. J. Phys. Chem. 1995,
78, 1396. 99, 3493.
(86) Perdew, J. P. Phys. Rev. B 1986, 33, 8822. (132) Halasinski, T. M.; Weisman, J. L.; Ruiterkamp, R.; Lee, T. J.;
(87) Dreuw, A.; Fleming, G. R.; Head-Gordon, M. Phys. Chem. Chem. Salama, F.; Head-Gordon, M. J. Phys. Chem. A 2003, 107, 3660.
Phys. 2003, 5, 3247. (133) Vaswani, H. M.; Hsu, C. P.; Head-Gordon, M.; Fleming, G. R. J.
(88) Dunietz, B. D.; Dreuw, A.; Head-Gordon, M. J. Phys. Chem. B Phys. Chem. B 2003, 107, 7940.
2003, 107, 5623. (134) Weisman, J. L.; Head-Gordon, M. J. Am. Chem. Soc. 2001, 123,
(89) Becke, A. D. J. Chem. Phys. 1988, 88, 2547. 11686.
(90) Olsen, J.; Jørgen, H.; Jensen, A.; Jørgensen, P. J. Comput. Phys. (135) Oumi, M.; Maurice, D.; Head-Gordon, M. Spectrochim. Acta A
1988, 74, 265. 1999, 55, 525.
(91) Eichkorn, K.; Treutler, O.; Öhm, H.; Häser, M.; R, A. Chem. Phys. (136) Franzen, S.; Kiger, L.; Poyart, C.; Martin, J.-L. Biophys. J. 2001,
Lett. 1995, 240, 283. 80, 2372.
(92) Jamorski, C.; Casida, M. E.; Salahub, D. R. J. Chem. Phys. 1996, (137) Dreuw, A.; Dunietz, B. D.; Head-Gordon, M. J. Am. Chem. Soc.
104, 5134. 2002, 124, 12070.
(93) Bauernschmitt, R.; Häser, M.; Treutler, O.; Ahlrichs, R. Chem. (138) Derewenda, Z.; Dodson, G.; Emsley, P.; Harris, D.; Nagai, K.;
Phys. Lett. 1997, 264, 573. Perutz, M.; Reynaud, J.-P. J. Mol. Biol. 1990, 211, 515.
(94) Savin, A.; Umrigar, C. J.; Gonze, X. Chem. Phys. Lett. 1998, 288, (139) Della Longa, S.; Arcovito, A.; Girasole, M.; Hazemann, J. L.;
391. Benfatto, M. Phys. Rev. Lett. 2001, 87, 155501.
(95) Gritsenko, O.; Baerends, E. J. J. Chem. Phys. 2004, 121, 655. (140) Makinen, M. W.; Eaton, W. Ann. N. Y. Acad. Sci. 1973, 206,
(96) Cai, Z.-L.; Sendt, K.; Reimers, J. R. J. Chem. Phys. 2002, 117, 210.
5543. (141) Yamaguchi, Y.; Yokomichi, Y.; Yokoyama, S.; Mashiko, S. Int.
(97) Grimme, S.; Parac, M. Chem. Phys. Chem. 2003, 3, 292. J. Quantum Chem. 2001, 84, 338.
(98) Cave, R. J.; Zhang, F.; Maitra, N. T.; Burke, K. Chem. Phys. (142) Yamaguchi, Y.; Yokoyama, S.; Mashiko, S. J. Chem. Phys. 2002,
Lett. 2004, 389, 39. 116, 6541.
Calculation of Excited States of Large Molecules Chemical Reviews, 2005, Vol. 105, No. 11 4037

(143) Miyahara, T.; Nakatsuji, H.; Hasegawa, J.; Osuka, A.; Aratani, (149) Hashimoto, T.; Choe, Y. K.; Nakano, H.; Hirao, K. J. Phys. Chem.
N.; Tsuda, A. J. Chem. Phys. 2002, 117, 11196. A 1999, 103, 1894.
(144) Cho, H. S.; Jeong, D. H.; Yoon, M.-C.; Kim, Y. H.; Kim, Y.-R.; (150) Parusel, A. B. J.; Grimme, S. J. Porphyrins Phthalocyanines
Kim, D.; Jeoung, S. C.; Kim, S. K.; Aratani, N.; Shinmori, H.; 2001, 5, 225.
Osuka, A. J. Phys. Chem. A 2001, 105, 4200. (151) Voorhis, T. V.; Head-Gordon, M. Mol. Phys. 2002, 100, 1713.
(145) Häberle, T.; Hirsch, J.; Pöllinger, F.; Heitele, H.; Michel-Beyerle, (152) Ziegler, T.; Rauk, A.; Baerends, E. J. Theor. Chim. Acta 1977,
M. E.; Anders, C.; Döhling, A.; Krieger, C.; Rückemann, A.; 43, 261.
Staab, H. A. J. Phys. Chem. 1996, 100, 18269. (153) Gilbert, A. T. B.; Gill, G. B. W.; Gill, P. M. W., manuscript in
(146) Polivka, T.; Zigmantas, D.; Herek, J. L.; He, Z.; Pascher, T.; preparation.
Pulleritis, T.; Cogdell, R. J.; Frank, H. A.; Sundström, V. J. Phys. (154) Vasudevan, J.; Stibrany, R. T.; Bumby, J.; Knapp, S.; Potenza,
Chem. B 2002, 106, 11016. J. A.; Emge, T. J.; Schugar, H. J. J. Am. Chem. Soc. 1996, 118,
(147) Thomson, M.; Novosel, M.; Roskos, H. G.; Müller, T.; Wagner, 11676.
M.; Fabrizi de Biani, F.; Zanello, P. J. Phys. Chem. A 2004, 108, (155) Scheer, H.; Inhoffen, H. H. In The Porphyrins; Dolphin, D., Ed.;
3281. Academic: New York, 1978; Vol. 2, p 45.
(148) Hasegawa, J.; Ozeki, Y.; Ohkawa, K.; Hada, M.; Nakatsuji, H.
J. Phys. Chem. B 1998, 102, 1320. CR0505627