2420
J. Chem. Theory Comput. 2009, 5, 2420–2435
Extensive TDDFT Benchmark: SingletExcited States of Organic Molecules
Denis Jacquemin,*,† Valerie Wathelet,† Eric A. Perpete,† and Carlo Adamo*,‡ ´ ` Groupe de ChimiePhysique Theorique et Structurale, Facultes UniVersitaires ´ ´ NotreDame de la Paix, rue de Bruxelles, 61, B5000 Namur, Belgium, and Ecole Nationale Superieure de Chimie de Paris, Laboratoire Electrochimie et Chimie ´ Analytique, UMR CNRSENSCP no. 7575, 11, rue Pierre et Marie
Extensive TDDFT Benchmark: SingletExcited States ofOrganic Molecules
Denis Jacquemin,*
,†
Vale´rie Wathelet,
†
Eric A. Perpe`te,
†
and Carlo Adamo*
,‡
Groupe de ChimiePhysique The´orique et Structurale, Faculte´s Uni
V
ersitaires NotreDame de la Paix, rue de Bruxelles, 61, B5000 Namur, Belgium, and Ecole Nationale Supe´rieure de Chimie de Paris, Laboratoire Electrochimie et Chimie Analytique, UMR CNRSENSCP no. 7575, 11, rue Pierre et Marie Curie,F75321 Paris Cedex 05, France
Received June 10, 2009
Abstract:
Extensive TimeDependent Density Functional Theory (TDDFT) calculations havebeen carried out in order to obtain a statistically meaningful analysis of the merits of a largenumber of functionals. To reach this goal, a very extended set of molecules (
∼
500 compounds,
>
700 excited states) covering a broad range of (bio)organic molecules and dyes have beeninvestigated. Likewise, 29 functionals including LDA, GGA,
meta
GGA, global hybrids, and longrangecorrected hybrids have been considered. Comparisons with both theoretical referencesand experimental measurements have been carried out. On average, the functionals providingthe best match with reference data are, one the one hand, global hybrids containing between22% and 25% of exact exchange (X3LYP, B98, PBE0, and mPW1PW91) and, on the otherhand, a longrangecorrected hybrid with a lessrapidly increasing HF ratio, namely LC
ω
PBE(20).Pure functionals tend to be less consistent, whereas functionals incorporating a larger fractionof exact exchange tend to underestimate signiﬁcantly the transition energies. For most treatedcases, the M05 and CAMB3LYP schemes deliver fairly small deviations but do not outperformstandard hybrids such as X3LYP or PBE0, at least within the vertical approximation. With theoptimal functionals, one obtains mean absolute deviations smaller than 0.25 eV, though theerrors signiﬁcantly depend on the subset of molecules or states considered. As an illustration,PBE0 and LC
ω
PBE(20) provide a mean absolute error of only 0.14 eV for the 228 states relatedto neutral organic dyes but are completely off target for cyaninelike derivatives. On the basisof comparisons with theoretical estimates, it also turned out that CC2 and TDDFT errors are ofthe same order of magnitude, once the abovementioned hybrids are selected.
1. Introduction
Developing methodological approaches able to accuratelydeliver the transition energies corresponding to electronicallyexcitedstates remains a major challenge for theoreticalchemists. Historically, the ﬁrst computational schemesdeveloped relied on semiempirical theories.
1
The mostsuccessful model, namely ZINDO,
2
was purposeddesignedto allow quick estimates of the main features of UV/visiblespectra and remains popular today. However, the quantitativeaspect of the obtained results (absorption wavelengths andtransition probabilities) was found to be highly systemdependent, a problematic feature.
3

5
More recently, calculations carried out for organic dyes have indicated that PM5
6
could be a promising approach,
7
but such a claim remainsto be tested on a broader set of transitions and molecules.At the other extreme of the theoretical palette, one ﬁndshighly correlated ab initio approaches such as SACCI,EOMCC, MRCI, or CASPT2 that allow very accurate
* Corresponding author email: denis.jacquemin@fundp.ac.be(D.J.); carloadamo@enscp.fr (C.A.).
†
FUNDP, Namur.
‡
ENSCP, Paris.
J. Chem. Theory Comput.
2009,
5,
2420–2435
2420
10.1021/ct900298e CCC: $40.75
2009 American Chemical Society
Published on Web 08/11/2009
estimates but are limited to rather small systems due to theirextreme computational cost. Two of the most extensiveinvestigations performed with such highaccuracy approaches(CC3 and CCSDR) have been published very recently byThiel and its coworkers.
8,9
Although these contributionscertainly represent a very large computational effort, theyhave been “limited” to molecules of about 15 atoms(naphthalene was the largest compound treated with CC3)and relied on a diffusefree polarized triple
ζ
basis set. Whileit is true that CASPT2 and CC2, the two “lighter” wavefunction approaches, could possibly be applied to moleculescontaining about 40 atoms,
10

18
the current implementationsof these
ab intio
theories often do not permit a systematicinclusion of medium effects. This is a problematic drawback,as it is wellknown that excitedstate properties tend to bemore solventsensitive than their groundstate counterparts.
19
Clearly, the difﬁculties to apply the highly correlatedapproaches to a broad set of molecules in a reallifeenvironment have yet not been completely solved. In termsof computational cost, one ﬁnds an intermediate betweensemiempirical theories and wave function approaches, namelytimedependent density functional theory (TDDFT).
20

23
TDDFT is the most widely applied ab initio tool formodeling the electronic spectra of organic and inorganicmolecules
24,25
and can be extended to incorporate environmental effects either through a modeling of the bulk environment
19,26

29
or through a variety of QM/MMapproaches.
30

35
Despite its successes and versatility, TDDFT is limited and suffers an important drawback: the qualityof the obtained results is profoundly functionaldependent.Indeed, the appropriate selection of the exchangecorrelationform is often crucial to grasp chemically sound conclusions.For most excitedstates, hybrid functionals that incorporatea fraction of exact exchange (EE) tend to provide moreaccurate estimates than pure functionals. Anyway, transitionwavelengths to excitedstates presenting a doubly excitedcharacter or a signiﬁcant chargetransfer nature are traditionally poorly estimated, as are the electronic spectrum of molecules having a strong multideterminantal nature. For surethese deﬁciencies are related to the approximate nature of today’s implementation, as illustrated by the recently developed longrangecorrected hybrids (LCH),
36

44
that appearto correctly appraise the chargetransfer properties. Contraryto the global hybrids (GH), LCH presents an EE percentagedepending on the interelectronic distance, allowing a physically correct asymptotic behavior when the two electronsare far apart.It is quite astonishing that only a limited number of contributions collated the
pros
and
cons
of functionals inthe TDDFT framework, for a signiﬁcant set of molecules.In Table 1, we summarize the selected methodologies andtraining sets for twelve investigations tackling this question.One can certainly ﬁnd many other TDDFT studies usingGH or LCH but often speciﬁc to a speciﬁc class of molecules.
5,45

53
As can be seen, not only the training setbut also the details of the methodologies selected forbenchmarking (including the size of the basis set and thepossible modeling of solvent effects) differ signiﬁcantly fromone work to the other. We believe it is especially strikingthat most studies include only a very small number of functionals (typically three) and that only four works usedmore than 100 excitations to obtain statistically meaningfulconclusions. Considering the different training sets andprocedures, it is to be expected that the conclusions of theseinvestigations are not perfectly uniform. While the obtainedmean absolute deviation (MAE) for the “best” functional istypically close to 0.25 eV, the actual ﬁndings are in factpartly antagonistic, making it difﬁcult to appreciate the“general” functional performance in the TDDFT framework:1. Tozer and coworkers concluded that CAMB3LYP
40
leads to much smaller deviations than B3LYP
54
for a varietyof transitions of mediumsize chromogens.
55
The averageB3LYP error being completely unacceptable (
>
1.0 eV) forboth Rydberg and chargetransfer states.
55,56
For valencetransitions, all tested functionals (PBE, B3LYP, and CAMB3LYP) provided similar MAE (0.27, 0.26, and 0.27 eV,respectively).
55
2. Rohrdanz and Herbert found that an accurate descriptionof both the groundstate and excitedstate properties of largemolecules was uneasy with common LCH functionals
57
andsubsequently design a LCH functional working for bothground and excitedstates.
58
This new LCH functionalprovided a MAE of about 0.3 eV.
58
3. For the
λ
max
related to
π
f
π
f
transitions in 100 organicdyes, we found, within the vertical approximation, that PBE0outperforms LCH and provides a MAE close to 0.15 eV,
59
the errors being of the same order of magnitude for
n
f
π
f
transitions.
60,61
For the same set of
π
f
π
f
transitions,CAMB3LYP provided signiﬁcantly larger deviations (0.26eV).
59
4. Thiel’s group used BP86, B3LYP, and BHHLYP andthey obtained MAE of 0.52 eV, 0.27, and 0.50 eV,respectively, for more than 100 transitions in small molecules,
62
using their own “best theoretical estimates”
8
asreference values.5. Dierksen and Grimme concluded from an extensivevibronic investigation of (mainly) hydrocarbons that theoptimal global hybrid should contain between 30% and 40%of EE.
63
Comparing their vertical (0

0) TDDFT data totheir solventcorrected experimental references, we calculatedMAE of 0.43 (0.57) eV, 0.21 (0.34) eV, and 0.31 (0.18) eVfor BP86, B3LYP, and BHHLYP, respectively.6. Very recently,
64
Goerigk et al. used the CASPT2 resultof ref 8 to benchmark doublehybrid functionals
65
and founda MAE of 0.22 eV for B2PLYP and B2GPPLYP, signiﬁcantly smaller than with B3LYP (0.30 eV) and conﬁrmedthis ﬁnding on a set of ﬁve large chromophores.Consequently, given an arbitrary molecule, it remainsdifﬁcult to know without testing what is (are) “reasonably”the most adequate functional(s) to evaluate the electronicspectra. Should one choose a GH or a LCH? Would the errorbe much larger with a GGA than with a GH? What is the“expected” accuracy with today’s computational procedure?Are
ab intio
functionals outperforming (or not) parametrizedfunctionals? Should the chosen functional vary for moleculesof different size? Of course, all these questions have beentackled in part in the abovementioned works, but with nogeneric answer embraced by a large community. Here, we
Extensive TDDFT Benchmark
J. Chem. Theory Comput., Vol. 5, No. 9, 2009
2421
have performed benchmarks that are more complete than anypreviously published data, both from the point of view of the number of molecules considered and of the set of pureand hybrid functionals incorporated.
2. Methodology
2.1. Strategy.
As can be seen in Table 1, two philosophiescan be used to benchmark TDDFT functionals: versusexperiment (VE) or versus theory (VT). Both approacheshave advantages and disadvantages. Trying to closely matchexperiment (VE) is generally desired in most practicalapplications and allows to include in the training set a widerange of molecules and compounds. On the other hand, onewould normally need to compute the full vibronic spectra(and not “simply” vertical transitions) and to perfectly modelthe experimental setup (pressure, temperature, full environmental effects, ...), both tasks being impossible for a largeset of solvated molecules. Additionally, it is not alwaysstraightforward to pinpoint the theoretical transition actuallycorresponding to the experimental measures, especially forhighly excited states. Comparisons with accurate wavefunction estimates (VT) allows straightforward and physicallymeaningful comparisons (same conditions, same transitions)but is obviously limited by the availability of theoretical data,i.e. only small molecules can be included. In many cases,CC2 results have been used as reference values for mediumsize molecules, a strategy that we think unsatisfying. Indeed,we computed a MAE of 0.27 eV (0.30 eV) between the CC2/ TZVP and the CASPT2/TZVP (“best estimates”) values forthe 103 singletexcited excitedstates of ref 8.
66
Even forlowlying excitedstates, CC2 is often off the theoretical limitby 0.1 eV,
8
a value equal to onehalf or onethird of thetypical TDDFT error.In the following, we will use both philosophies so to beas general as possible. In what concerns the versus theoryscheme, we have selected Thiel’s set (VT set in thefollowing) and mimic exactly the computational procedure(basis set and geometry). For the VE set, we have used acomputational strategy that is at the limit of today’s possibilities for such a set of molecules, trying to circumventthe possible limitations of our computational procedure. Forthe sake of consistency, we have chosen to use a uniformmethodology (basis set, solvent effects, ...) for all VEmolecules.
2.2. General Computational Procedure.
All calculationshave been performed with the Gaussian suite of programs,using both the commercial and development versions
67,68
with a tight selfconsistent ﬁeld convergence threshold (10

8
to 10

10
au). For the VE set, we have followed a wellestablished threestep approach:
25
i) the groundstate geometry of each compound has been optimized until the residualmean force is smaller than 1.0
×
10

5
au (socalled
tight
threshold in Gaussian); ii) the vibrational spectrum isanalytically determined to conﬁrm that the structure is a trueminimum; and iii) the vertical transition energies to thevalence excited states are computed with TDDFT. For theVT set, the geometries have been taken from ref 8 and stepiii) directly performed.As the majority of experimental data are obtained incondensed phase, we have included bulk solvent effects inour VE model (all VT calculations are in gasphase). Thiswas performed at each stage, including geometry optimizations and Hessian calculations, using the wellknown Polarizable Continuum Model (PCM),
19
that is able to obtaina valid approximation of solvent effects as long as no speciﬁcinteractions link the solute and the solvent molecules.Typically solvent

solute hydrogen bonds tend to inﬂuencemore signiﬁcantly the
n
f
π
f
transitions than their
π
f
π
f
counterparts, and we have tried to select aprotic solvent forthe former, at least when different experimental values areavailable. The list of solvent selected is given in theSupporting Information. The default PCM Gaussian parameters have generally been used, though for a few calculationsif was necessary to change the atomic raddi (UAKS insteadof UA0) or to switch off the presence of smoothing sphere(NoAddSph) to converge the force minimizations. For therecords, note that some default PCM parameters might differbetween the two versions of the program used. All TDDFTcalculations have been performed within the nonequilibriumapproximation, valid for absorption spectra.
19
2.3. Functionals and Basis Sets.
As we want to assessthe
pros
and
cons
of a series of DFT approaches, a veryextended set of functionals has been used. Apart from theTimeDependent Hartree

Fock approach (TDHF, refereedto as HF in the following), the selected functionals can beclassiﬁed in ﬁve major categories: LDA, GGA,
meta
GGA,GH, and LCH. In the ﬁrst category, that is expected to bethe less efﬁcient we have selected only one functional,SVWN5.
69,70
We have chosen four GGAs, namely BP86,
71,72
BLYP,
71,73
OLYP
73,74
and PBE,
75
whereas we have pickedup three popular
meta
GGA: VSXC,
76
τ
HCTH
77
andTPSS.
78
Twelve global hybrids have been used: TPSSh(10%),
79
O3LYP (11.61%),
80
τ
HTCHh (15%),
77
B3LYP(20%),
54,81
X3LYP (21%),
82
B98 (21.98%),
83
mPW1PW91(25%),
84
PBE0 (25%),
85,86
M05 (28%),
87
BMK (42%),
88
BHHLYP (50%),
89
and M05

2X (56%).
90
The LCH constitute the last category and use a growing fraction of EEwhen the interelectronic distance increases. This is formallyperformed by partitioning the twoelectron operator as
36,40,91
The ﬁrst term of the rhs of this equation describes the socalled shortrange effect and is modeled through DFTexchange, whereas the second term corresponds to the longrange contribution calculated with the HF exchange formula.In eq 1
ω
is the range separation parameter, while
R
and
R +
β
deﬁne the EE percentage at
r
12
)
0 and
r
12
)
∞
,respectively. The LC model uses
R )
0.00,
β
)
1.00, and
ω
)
0.33 au

1
in eq 1
37,38
and has been applied to bothGGA and
meta
GGA to give LCBLYP, LCOLYP, LCPBE, LC
τ
HCTH, and LCTPSS. The approach designedby Vydrov and Scuseria,
42,43
namely LC
ω
PBE, with
ω
)
0.40 au

1
and
R )
0,
β
)
1, has been used as well. Notethat in LC
ω
PBE, the shortrange exchange functional canbe rigorously derived
41,92
by integration of the model1
r
12
)
1

[
R +
β
erf(
ω
r
12
)]
r
12
+R +
β
erf(
ω
r
12
)
r
12
(1)
Extensive TDDFT Benchmark
J. Chem. Theory Comput., Vol. 5, No. 9, 2009
2423