Article

pubs.acs.org/JCTC

Simulating One-Photon Absorption and Resonance Raman

Scattering Spectra Using Analytical Excited State Energy Gradients
within Time-Dependent Density Functional Theory
Daniel W. Silverstein,† Niranjan Govind,*,‡ Hubertus J. J. van Dam,‡ and Lasse Jensen*,†
†
Department of Chemistry, The Pennsylvania State University, 104 Chemistry Building, University Park, Pennsylvania 16802, United
States
‡
William R. Wiley Environmental Molecular Sciences Laboratory, Pacific Northwest National Laboratory, 902 Battelle Boulevard,
Richland, Washington 99352, United States
*
S Supporting Information

ABSTRACT: A parallel implementation of analytical time-

dependent density functional theory gradients is presented for
the quantum chemistry program NWChem. The implementa-
tion is based on the Lagrangian approach developed by Furche
and Ahlrichs. To validate our implementation, we first
calculate the Stokes shifts for a range of organic dye molecules
using a diverse set of exchange-correlation functionals
(traditional density functionals, global hybrids, and range-
separated hybrids) followed by simulations of the one-photon
absorption and resonance Raman scattering spectrum of the
phenoxyl radical, the well-studied dye molecule rhodamine 6G, and a molecular host−guest complex (TTF⊂CBPQT4+). The
study of organic dye molecules illustrates that B3LYP and CAM-B3LYP generally give the best agreement with experimentally
determined Stokes shifts unless the excited state is a charge transfer state. Absorption, resonance Raman, and fluorescence
simulations for the phenoxyl radical indicate that explicit solvation may be required for accurate characterization. For the host−
guest complex and rhodamine 6G, it is demonstrated that absorption spectra can be simulated in good agreement with
experimental data for most exchange-correlation functionals. However, because one-photon absorption spectra generally lack
well-resolved vibrational features, resonance Raman simulations are necessary to evaluate the accuracy of the exchange-
correlation functional for describing a potential energy surface.

1. INTRODUCTION properties at costs comparable to the ground state gradient.

The excited state properties of molecules lead to diverse Hutter et al. employed this framework for TDDFT gradients
applications,1,2 including solution phase detection of molecular within the plane wave basis formalism.38,39 Seth et al. extended
binding events,3−5 photodynamic therapy,6,7 and biological the Lagrangian formulation for Slater-type orbitals and spin-flip
imaging.8−10 A detailed understanding of molecule-radiation TDDFT,40 and Petrenko et al. implemented the expressions for
interactions often requires analysis of experimental data and including frozen core orbitals and chain of spheres exchange.41
applications of theory, which has been demonstrated by several Extensions of the Lagrangian formulation of TDDFT gradients
groups.11−15 Since its introduction16 and implementation for have also been devised to include range separated XC
molecular systems,17,18 time-dependent density functional functionals,42,43 the Tamm-Dancoff approximation,44 molec-
theory (TDDFT) has become a powerful approach to study ular-orbital free implementations,45 and environmental effects
low-lying valence excited states19−25 as well as core-level including continuum46,47 and explicit solvent models.48,49
excitations26−31 efficiently and accurately. Although the Lagrangian approach is widely used, TDDFT
Excited-state gradients, which are needed for spectroscopy gradients have also been developed employing implicit
simulations and the dynamics of photochemical processes can differentiation,50 Green’s function methods,51 and real-space
also be evaluated analytically at a computationally reasonable TDDFT.52
cost within the framework of TDDFT. Pioneering efforts by Analytical excited state gradients are also useful for describing
Amos et al. demonstrated the accuracy of TDDFT energy vibronic features for one-photon absorption (OPA) and
gradients for determining excited state geometries and dipoles resonance Raman scattering spectra.53−57 Raman scattering is
using time-dependent Hartree−Fock (TDHF) and a variety of an inelastic light scattering process, where radiation with
exchange-correlation (XC) functionals.32−35 Furche and co-
workers36,37 introduced the Lagrangian formulation of the Received: September 1, 2013
excitation energy, enabling efficient calculation of excited state Published: November 6, 2013

© 2013 American Chemical Society 5490 dx.doi.org/10.1021/ct4007772 | J. Chem. Theory Comput. 2013, 9, 5490−5503
Journal of Chemical Theory and Computation Article

frequency ωL is absorbed, and then radiation is scattered with δ 2E xc

frequency ωs = ωL − ωv, with ωv corresponding to a vibrational f στxc =
δρσ (r′) δρτ (r)
frequency of the molecule. When ωL is tuned to an electronic (4)
excitation energy, resonance enhancement is observed, and the For all quantities, indices i, j, k, l, ... indicate occupied orbitals; a,
process is called resonance Raman scattering (RRS). Vibra- b, c, d, ... are virtual (unoccupied) orbitals; σ, τ, and υ are spin
tional information that is coupled to an electronic transition is indices; and p, q, r, s, ... are arbitrary orbitals. The molecular
therefore obtained by this method. RRS has been demonstrated orbital (MO) energies are εpσ, the percentage of the Hartree−
as a valuable technique in characterizing fluorescent dyes,58−61 Fock (HF) exchange is given by cx, and the two-electron
biomolecules,62−64 and photon-induced electron transfer integrals are written in Mulliken notation. Upon solving the
processes.11,65−68 TDDFT equations, the excitation energy ω and excitation
In this work, we present our parallel implementation of vectors (sometimes called configuration interaction (CI)
analytical TDDFT excitation energy gradients based on the vectors) X and Y are obtained.
Lagrangian formalism of Furche and Ahlrichs36 in the The derivation of TDDFT gradients depends on the
NWChem quantum chemistry program.69 We also present evaluation of the Lagrangian functional given by
the utility of TDDFT excitation energy gradients to determine
the Stokes shifts for a set of organic dyes, and by combining the L[X, Y, ω , C , Z , W]
excitation energy gradients with a time-dependent wavepacket 1
formalism,70−75 we demonstrate how OPA and RRS can be = ∑ [(X + Y )iaσ (A + B)iaσjbτ (X + Y )jbτ
2 iaσjbτ
simulated efficiently for fluorescent dyes and charge-transfer
complexes.53 The rest of the paper is organized as follows. In + (X − Y )iaσ (A − B)iaσjbτ (X − Y )jbτ ]
section 2, we outline the theory for the Lagrangian formalism
for excitation energy gradients and its application toward − ω[∑ (X + Y )iaσ (X − Y )iaσ − 1] + ∑ ZiaσFiaσ
obtaining vibronic coupling parameters. The computational iaσ iaσ
details are presented in section 3. Finally, in section 4 we
present and discuss features of our implementation based
− ∑ Wpqσ(Spqσ − δpq)
pqσ , p ≤ q
simulations of optical properties for organic dyes, radicals, and
molecules with charge-transfer states using several XC = G[X , Y , ω] + ∑ ZiaσFiaσ − ∑ Wpqσ(Spqσ − δpq)
functionals. iaσ pqσ , p ≤ q
(5)
2. THEORY The Lagrangian functional L depends on the excitation energy
2.1. Analytical Excitation Energy Gradients with and excitation vectors, the MO coefficients C, and Lagrange
TDDFT. In this section, we present the key expressions in multipliers Z and W. In eq 5, Fpqσ and Spqσ are the Fock and
the derivation of gradients of the excitation energy using overlap matrices, respectively. The Fock matrix is defined as
TDDFT based on the original derivation in ref 36. The
TDDFT equations can be written as follows: Fpqσ = hpqσ + ∑ [(pqσ |iiτ) − cxδστ(piσ |iqσ )] + Vpqxcσ
iτ
1 (6)
G[X, Y, ω] = [(X + Y)T (A + B)(X + Y)
2 In eq 6, h is the core Hamiltonian, the second term is the
+ (X − Y)T (A − B)(X − Y)] Coulomb integral, the third term is the HF exchange
contribution, and the last term, the XC potential, is given by
− ω[(X + Y)T (X − Y) − 1]
δE xc
=
1
∑ [(X + Y )iaσ (A + B)iaσjbτ (X + Y )jbτ
xc
Vpq xc
σ = ⟨pσ |Vσ |qσ ⟩ = ∫ dr ϕpσ (r ) ϕ (r )
δρσ (r ) qσ
2 iaσjbτ (7)
where ϕpσ is a MO. The Lagrange multipliers enforce two
+ (X − Y )iaσ (A − B)iaσjbτ (X − Y )jbτ ]
constraints:
− ω[∑ (X + Y )iaσ (X − Y )iaσ − 1] 1. Z enforces that the ground state Kohn−Sham (KS)
iaσ (1) equation has been solved.
2. W enforces that the MOs are orthonormal.
In eq 1, the matrices (A + B) and (A − B), which are
sometimes called orbital rotation Hessians, have the following The Lagrange multipliers ensure that the Lagrangian
matrix representation functional is forced to obey the variational principle. In other
words, the Lagrangian functional is stationary with respect to
(A + B)iaσjbτ = δστδabδij(εaσ − εiσ ) + 2(iaσ |jbτ ) + 2fiaxcσjbτ changes in the MO coefficients:
∂L ∂G ∂Fiaτ ∂Srsτ
− cxδστ[(abσ |ijσ ) + (jaσ |ibσ )] (2) =0= + ∑ Ziaτ − ∑ Wrsτ
∂Cμpσ ∂Cμpσ iaτ
∂Cμpσ rsτ , r ≤ s
∂Cμpσ
and (8)

(A − B)iaσjbτ = δστδabδij(εaσ − εiσ ) + cxδστ[(jaσ |ibσ ) The most significant consequence of eq 8 is that derivatives of
the Lagrangian functional with respect to geometric perturba-
− (abσ |ijσ )] (3) tions or external fields (for polarizabilities) can be evaluated by
keeping the MO coefficients fixed at the ground state optimized
In eqs 2 and 3, fxc
στ is the XC kernel values. This enables gradients of the excitation energy to be
5491 dx.doi.org/10.1021/ct4007772 | J. Chem. Theory Comput. 2013, 9, 5490−5503
Journal of Chemical Theory and Computation Article

calculated at a cost comparable to a ground state energy Q iaσ = Hia+σ(T) + ∑ {(X + Y )ibσ Hab+ σ[(X + Y)]
gradient. Evaluation of eq 8 is required to determine the b
Lagrange multipliers and makes frequent use of the trans- −
+ (X − Y )ibσ Habσ [(X − Y)]}
formation from atomic orbital (AO) indices (μ, ν, λ, κ, ...) to
MO indices. For an arbitrary quantity Vpqσ, this transformation +2 ∑ [giaxcσkcτldυ](X + Y )kcτ (X + Y )ldυ
kcτldυ (14)
can be written as
Q aiσ = ∑ {(X + Y )jaσ H ji+σ[(X + Y )]
Vpqσ = ∑ CμpσVμνσCνqσ j
μν (9)
+ (X − Y)jaσ H −jiσ[(X − Y)]} (15)
The derivatives on the right-hand side of eq 8 can be
evaluated by employing the definitions of G[X,Y,ω], the Fock Q abσ = ω ∑ {(X + Y )iaσ (X − Y )ibσ + (X − Y )iaσ
matrix, and overlap matrix. For the Fock and overlap matrices, i

the expressions are given by × (X + Y )ibσ } + ∑ εbσ {(X + Y )iaσ (X + Y )ibσ

i
∂Fiaτ
∑ Ziaτ ∑ Cμqσ + (X − Y )iaσ (X − Y )ibσ } (16)
iaτ μ ∂Cμpσ
The linear transformations are
= ∑ (δστδipδqa + δστδapδiq)εqσ Ziaτ + ∑ δjpH pq+ σ[Z]
iaτ j +
H pqσ [V] = ∑ [2(rsτ|qpσ ) + 2f pqxcσrsτ − cxδσυ[(rqυ|psυ)
(10) rsτ

+ (rpυ|qsυ)]]Vrsτ (17)
∂S
∑ ∑ Wrsτ rsτ Cμqσ = (1 + δpq)Wpqσ −
H pqσ [V ] = ∑ {[cxδστ[(rqσ |psσ ) − (qsσ |prσ )]]}Vrsτ
μ rsτ , r ≤ s
∂Cμpσ (11) rsτ (18)

For the functional G[X,Y,ω], the derivatives of the orbital where the vector T is the unrelaxed one-particle difference
rotation Hessians with respect to the MO coefficients are density matrix, defined as
required. 1
Tabτ = ∑ {(X + Y )iaτ (X + Y )ibτ
2 i
∂G[X, Y, ω]
Q pqσ = ∑ Cμqσ + (X − Y )iaτ (X − Y )ibτ } (19)
μ ∂Cμpσ
⎧ 1
1 ⎪ ∂ Tijτ = − ∑ {(X + Y )iaτ (X + Y )jaτ
= ∑ (X + Y )kcτ ⎨∑ [(A + B)kcτldυ ] 2 a
⎪
2 kcτldυ ⎩ μ ∂Cμpσ
+ (X − Y )iaτ (X − Y )jaτ } (20)
⎫
⎪ 1
× Cμqσ ⎬(X + Y )ldυ + ∑ (X − Y )kcτ Tiaτ = Taiτ = 0 (21)
⎪
⎭ 2 kcτldυ
and the third derivative of the XC functional is
⎧ ⎫
⎪ ∂ ⎪
× ⎨∑ [(A − B)kcτldυ ]Cμqσ ⎬(X − Y )ldυ xc δ 3E xc
⎪
⎩ μ ∂Cμpσ
⎪
⎭ gστυ =
δρσ (r ) δρτ (r′) δρυ (r″) (22)
(12)
By combining eqs 13 through 16 with eqs 10 and 11, four
These derivatives lead to lengthy expressions but may be linear equations are obtained that enable the Lagrange
simplified by looking at different blocks of the matrix Qpqσ. For multipliers Z and W to be determined.
the occupied−occupied block (p = i ≤ q = j), occupied−virtual Q ijσ + Hij+σ[Z] = (1 + δij)Wijσ (23)
block (p = i and q = a), virtual−occupied block (p = a and q =
i), and virtual−virtual block (p = a ≤ q = b), the expressions Q iaσ + εaσ Ziaσ + Hia+σ[Z] = Wiaσ (24)
become
Q aiσ + εiσ Ziaσ = Wiaσ (25)
Q ijσ = ω ∑ {(X + Y )iaσ (X − Y )jaσ + (X − Y )iaσ
a Q abσ = (1 + δab)Wabσ (26)
× (X + Y )jaσ } − ∑ εaσ {(X + Y )iaσ (X + Y )jaσ Taking the difference between eqs 24 and 25, it is possible to
a evaluate Z
+ (X − Y )iaσ (X − Y )jaσ } + Hij+σ[T] ∑ (A + B)iaσjbτ Zjbτ = −R iaσ
+2 ∑ [gijxcσkcτldυ](X + Y )kcτ (X + Y )ldυ jbτ (27)
kcτldυ (13) where the right-hand side vector is
5492 dx.doi.org/10.1021/ct4007772 | J. Chem. Theory Comput. 2013, 9, 5490−5503
Journal of Chemical Theory and Computation Article

R iaσ = Hia+σ[T] + ∑ {(X + Y )ibσ Hab+ σ[(X + Y)] Table 1. Definitions of the XC Functionalsa
b
functional definition
−
+ (X − Y )ibσ Habσ [(X − Y)]} LDA Slater + VWN 5
BP86 Becke 88 + Perdew 86
+2 ∑ [giaxcσkcτldυ](X + Y )kcτ (X + Y )ldυ
PBE PBEx + PBEc
kcτldυ
BLYP Becke 88 + LYP
− ∑ {(X + Y )jaσ H ji+σ[(X − Y)] B3LYP Slater (80%) + Becke 88 (nonlocal, 72%) + HF Exchange (20%)
j + LYP (81%) + VWN 1 RPA (19%)
PBE0 PBEx (75%) + HF Exchange (25%) + PBEc
+ (X + Y )jaσ H ji−σ[(X − Y)]} (28) BHLYP Slater (50%) + Becke 88 (nonlocal, 50%) + HF Exchange (50%)
+ LYP
The Lagrange multiplier W can be evaluated after we obtain Z CAM- CAM-Becke 88 (α = 0.19, β = 0.46, ω = 0.33) + LYP
B3LYP
1 LC-PBE CAM-PBE (α = 0.00, β = 1.00, ω = 0.30) + PBEc
Wijσ = (Q + Hij+σ[Z])
1 + δij ijσ LC-PBE0 CAM-PBE (α = 0.25, β = 0.75, ω = 0.30) + PBEc
(29)
BNL Baer−Neuhauser−Livshits (α = 0.00, β = 1.00, ω = 0.33) + LYP
LC- LC-ωPBE (α = 0.00, β = 1.00, ω = 0.30) + PBEc
1 ωPBE
Wabσ = Q
1 + δab abσ (30) LC- LC-ωPBE (α = 0.20, β = 0.80, ω = 0.20) + PBEc
ωPBEh
Wiaσ = Q aiσ + εiσ Ziaσ LC-BLYP CAM-Becke 88 (α = 0.00, β = 1.00, ω = 0.33) + LYP
(31) a
The table is separated between the local and semilocal functionals,
With Z and W, the gradient of the Lagrangian can be evaluated global hybrids, and long-range corrected hybrids. For the latter, we
to obtain the gradient of the excitation energy. The expression report the α, β, and ω parameters.
for the gradient of the excitation energy is
(TDA) can be used to evaluate gradients for restricted singlet,
ωξ = ∑ hμνξPμνσ − ∑ SμνξWμνσ restricted triplet, and unrestricted calculations.
μνσ μνσ 2.2. Calculating Vibronic Coupling Parameters. The
xc(ξ) vibronic coupling parameters (dimensionless displacements)
+ ∑ V μνσ Pμνσ + ∑ (μν|κλ)ξ Γμνσκλτ are related to the gradient of the excited state energy (E0k,
μνσ μνκλστ
which is equivalent to ω from section 2.1) along normal mode
xc(ξ)
+ ∑ f μνσκλτ (X + Y )μνσ (X + Y )κλτ Qa. This is shown in Figure 1. In this section, the normal mode
μνκλστ (32)

The superscript, ξ, indicates the derivative of a quantity with

respect to a nuclear coordinate. When written as (ξ), this
indicates that the derivative is taken with respect to a nuclear
coordinate with the MO coefficients held fixed at their ground
state values. This is a consequence of the density functional
theory (DFT) contributions to the gradients implicitly
depending on the MO coefficients. Here, two new quantities
are defined, the relaxed one-particle difference density matrix
(P = T + Z) and the effective two-particle difference density
matrix Γμνσκλτ, defined as Figure 1. A depiction of the independent-mode displaced harmonic
oscillator model, assuming a resonance Raman scattering process.
1 Here, ωL is the incident frequency of light, ωs is the scattered
Γμνσκλτ = {2(X + Y )μνσ (X + Y )κλτ + 2PμνσDκλτ
2 frequency of light, and Δ0k
a is the dimensionless displacement of
normal mode qa.
− cxδστ[PμκσDνλτ + PμλσDνκτ + (X + Y )κνσ
× (X + Y )μλτ + (X + Y )κμσ (X + Y )λντ in mass-weighted units is written as Qa, while it is written as qa
for normal (dimensionless) coordinates. On the basis of
− (X − Y )κνσ (X − Y )μλτ + (X − Y )κμσ approximations where the ground and excited state potential
energy surfaces are harmonic and that there is no change in the
× (X − Y )λντ ]} (33) normal-mode frequencies when the ground and excited states
The expressions in this section were implemented into the are compared (i.e., no Duschinsky rotations76), the dimension-
NWChem quantum chemistry program.69 Evaluation of eqs 28 less displacements are calculated as
and 32 require third order functional derivatives, which were ⎛ ∂E 0k ⎞
derived, implemented, and tested using numerical differ- Vq = ⎜⎜ ⎟⎟ = −ωa0Δa0k
entiation of second order derivatives. Several types of XC a
⎝ a ⎠q = 0
∂q
a (34)
functionals are included in the implementation: the local
density approximation (LDA), generalized gradient approx- where ωa0 is the normal-mode frequency and is the Δ0k
a
imation (GGA), global hybrids, and range-separated hybrids. A dimensionless displacement of mode qa for the nth excited
complete list of the functionals is included in Table 1. Both full state. In general, it is simpler to evaluate the gradient with
linear response TDDFT and the Tamm-Dancoff approximation respect to Qa more directly. This requires an additional step,77
5493 dx.doi.org/10.1021/ct4007772 | J. Chem. Theory Comput. 2013, 9, 5490−5503
Journal of Chemical Theory and Computation Article

∂E 0k ⎛ dQ ⎞ ∂E 0k ℏ ∂E 0k ℏ energy, which come from performing an analytical excited state

= ⎜⎜ a ⎟⎟ = = VQ geometry optimization, are used to determine the dimension-
∂qa ⎝ dqa ⎠ ∂Q a 2πcνa0 ∂Q a 2πcνa0 a less displacements and only require one calculation. A side-by-
(35) side comparison of the numerical and analytical methods shows
that if you account for the redundant calculations (the ground
The only variability in this expression is how the excited state
state geometry optimization and normal-mode analysis), the
gradient in terms of mass-weighted coordinates VQa is evaluated. analytical approach results in a significant reduction in both
It has been demonstrated previously that numerical gradients CPU and wall time. This speedup results from analytical
can be used for vibronic simulations of spectroscopy.15,78−83 gradients requiring roughly the same computational effort as a
Other research groups have employed analytical gradients of ground state gradient calculation,36 which can decrease the
the excited state.53,54 Numerical gradients of the excitation computation time by 2 orders of magnitude as we have
energy require finite differentiation, where excitation energies shown.85
are calculated by stepping in the positive and negative The analytical method involves two different techniques:
directions along the normal modes. analytical excited state geometry optimization (method I) and
2.2.1. Numerical Excited State Energy Gradients. The analytical excited state gradients at the ground state equilibrium
numerical evaluation of excited state energy gradients revolves position (method II). We will focus on method II here. For
around calculating the linear response of a molecule up to a method II, the dimensionless displacement is calculated by
particular excitation. For any molecule, there are 3N − 5 or 3N projecting the gradient of the excited state energy at the ground
− 6 (N is the number of atoms) normal modes. Performing a state equilibrium position onto the normal modes. This method
numerical evaluation of the excited state energy gradient is written as
requires either 6N − 10 or 6N − 12 calculations, because this
1
requires calculations of the excitation energy in the positive and VQ a = ∑ VX , iLia
negative directions along every normal mode. This can be i
mi (38)
reduced significantly by focusing on a portion of the normal
modes in a wavenumber region. where mi is the mass of the atom corresponding to component i
For numerical evaluation of the excited state gradients, and VX,i is the Cartesian excited state energy gradient at the
formulas analogous to those discussed by Reiher et al.84 for the ground state equilibrium position. This is obtained from the
derivatives of the polarizability along normal mode Qa are used, first cycle of the excited state geometry optimization, assuming
substituting the excitation energy in place of the polarizability. that the ground state equilibrium geometry is the starting point
These derivatives are expressed using a three-point central of that optimization. From the gradient VQa, the dimensionless
differences formula as displacement is calculated using eq 34.
⎛ ∂E 0n ⎞ E 0k (R eq + sR R a) − E 0k (R eq − sR R a) 3. COMPUTATIONAL DETAILS
⎜⎜ ⎟⎟ =
⎝ ∂Q a ⎠Q 2sQ a|Q anorm| In Table 1, we present the definitions of all the XC functionals
a= 0 that are available with our analytic TDDFT gradients
(36) implementation and the ones used in this work. For the
In eq 36, Req is a vector of the optimized coordinates of the long-range (LC) corrected XC functionals, we use the
molecule’s ground state, Ra is the Cartesian normal mode following definition:
vector, sR is the Cartesian step size (typically 0.01 is used), and | 1 1 − [α + β erf(ωr12)] α + β erf(ωr12)
Qnorm
a | is the norm of the normalized mass-weighted normal = +
r12 r12 r12 (39)
mode vector. The latter quantity is used to give the equation
appropriate dimensions but is clearly equal to one. In this where r12 is the electronic separation, ω is the attenuation
expression, sQa is parameter, α is a parameter allowing the incorporation of HF
exchange, and β is a parameter allowing for the incorporation of
⎡ 3N (Q norm)2 ⎤−1/2 DFT exchange, with the requirements 0 ≤ α ≤ 1, 0 ≤ β ≤ 1,
sR
= sR ⎢∑ ⎥
a,i
sQ a = and 0 ≤ α + β ≤ 1 being satisfied.
|R a| ⎣⎢ i = 1 mi ⎥⎦
(37) As mentioned earlier, geometry optimizations, normal-mode
analysis, and excited state energy gradient calculations were
Either the expression in the middle or the right side may be performed in a development version of the NWChem quantum
used to find sQa depending on whether the normalized mass- chemistry program,69 while the OPA and RRS simulations were
weighted normal mode vector Qa or the unnormalized performed using a locally developed program, TDSPEC.86
Cartesian normal mode vector Ra is available. If the normalized For the simulations presented in section 4.1, the organic dye
mass-weighted normal mode vector is used, each squared molecules were optimized in the ground and excited states
component is divided by the mass of the ith atom (mi). using the BLYP, LC-BLYP, B3LYP, and CAM-B3LYP XC
2.2.2. Analytical Excited State Energy Gradients. If functionals with the 6-31G* basis set. In section 4.3, the host−
analytical geometry optimizations are available in a quantum guest complex was optimized in its ground state, and vibrational
chemistry program, a significant advantage can be gained in frequencies were determined using B3LYP with dispersion
calculating the dimensionless displacements. In the numerical corrections through the DFT-D2 method87 and 6-311G* basis
method described above, the majority of the time spent set. Excitation energy gradients were then determined using the
involves the 6N − 10 or 6N − 12 calculations to obtain all of B3LYP and LC-ωPBE functionals to simulate OPA and RRS.
the simulation parameters. In general, for all but the smallest Results in section 4.4 for rhodamine 6G are based on the
molecules, collection of that data requires substantial amounts ground state optimized geometry and vibrational frequencies
of computer time. Analytical gradients of the excited state using B3LYP/6-311G*. All normal-mode frequencies are scaled
5494 dx.doi.org/10.1021/ct4007772 | J. Chem. Theory Comput. 2013, 9, 5490−5503
Journal of Chemical Theory and Computation Article

Figure 2. Structures of organic dye molecules used in this work. The numbered dyes correspond to those presented in Table 2.

Table 2. Absorption Maxima, Emission Maxima, and Stokes Shifts (S.S.) for a Test Set of Dye Molecules Using BLYP, LC-
BLYP, B3LYP, and CAM-B3LYP Using the 6-31G* Basis Seta
BLYP LC-BLYP B3LYP CAM-B3LYP exp.
dye Abs. Em. S.S. Abs. Em. S.S. Abs. Em. S.S. Abs. Em. S.S. S.S.
1 3.09 2.12 0.97 3.72 3.04 0.68 3.55 2.72 0.83 3.79 3.07 0.72 0.86b
2 3.68 3.62 0.06 5.00 4.79 0.21 4.16 4.06 0.10 4.73 4.56 0.17 0.23c
3 2.84 2.16 0.68 4.00 3.54 0.46 3.36 2.85 0.51 3.85 3.48 0.37 0.44d
4 2.63 2.49 0.14 3.71 3.34 0.37 3.11 2.90 0.21 3.58 3.26 0.32 0.19e
5 1.94 1.87 0.07 1.97 1.77 0.20 2.09 2.01 0.08 2.08 1.94 0.14 0.02f
6 1.99 1.26 0.73 2.70 2.45 0.25 2.31 2.09 0.22 2.64 2.41 0.23 0.10e
7 3.16 (I) − 3.97 3.31 0.66 3.54 2.93 0.61 3.88 3.29 0.59 0.68g
8 2.14 1.90 0.14 3.15 2.60 0.55 2.49 2.16 0.33 2.78 2.45 0.33 0.36h
9 1.94 1.89 0.05 1.96 1.83 0.13 2.09 2.02 0.07 2.07 1.97 0.10 0.02f
10 1.47 1.29 0.18 1.97 1.94 0.03 1.94 1.76 0.18 2.13 2.11 0.02 −
ME 0.07 0.08 0.01 0.02
MAE 0.18 0.10 0.06 0.08
RMSD 0.26 0.12 0.07 0.09
a
Values are given in eV. Instabilities in the optimized excited state are labeled (I). Mean errors (ME), mean absolute errors (MAE), and root mean
square deviations (RMSD) are reported in eV in comparison to experiment. bReference 93. See text for discussion. cReference 94. dReference 95.
e
Reference 96. fReference 100. gReferences 97 and 98. hReference 99.

by 0.98 to account for anharmonicity in the normal modes and functionals to determining accurate Stokes shifts. Sections
limitations of the basis set.88−91 4.2−4.4 use the excited state energy gradients to evaluate
vibronic coupling parameters, illustrating the capabilities of the
4. RESULTS AND DISCUSSION
code for studying radicals, molecules with charge-transfer
The results presented in this section show applications of the
analytical excitation energy gradient implementation to a set of excitations, and fluorescent dyes. In section 4.5, the timings of
selected systems. Section 4.1 highlights the available excited the current development version of the TDDFT gradient code
state geometry optimization code for several types of XC are shown.
5495 dx.doi.org/10.1021/ct4007772 | J. Chem. Theory Comput. 2013, 9, 5490−5503
Journal of Chemical Theory and Computation Article

Figure 3. (a) Absorption (green and violet lines) and fluorescence (red and blue lines) spectra from 2,3,5-trifluorobenzaldehyde in the syn and anti
conformations. (b) Simulated fluorescence spectra of the syn and anti conformations compared with the experimental fluorescence spectrum.93
Results are calculated with B3LYP/aug-cc-pVTZ.

4.1. Absorption and Emission by Organic Dyes. In this hydroxybenzylidene-2,3-dimethylimidazolinone) are more sen-
section, we investigate the dependence of Stokes shifts on the sitive to the amount of HF exchange, and the most accurate
XC functional for 10 organic chromophores. The test set functional is the range separated hybrid LC-BLYP. For dyes 3
includes the chromophores used in ref 92, which are illustrated and 7, this is a consequence of the S1 ← S0 transition involving
in Figure 2. Dye molecules 1−6, 8, and 9 have experimentally an intramolecular charge transfer event. It was not possible to
measured absorption and fluorescence spectra in the gas phase, optimize the bright excited state of dye 7 using BLYP due to
allowing a direct comparison of the Stokes shifts with the the presence of a nearby spurious charge transfer transition. On
calculated results. Dye 7 (4′-hydroxybenzylidene-2,3-dimethy- the other hand, dye 2 involves excitations of electrons in a
limidazoline) is a model of the chromophore in green highly delocalized π system, where traditional XC functionals
fluorescent protein and only has available absorption and perform poorly in general. Asymptotic corrections to the XC
fluorescent data in the protein environment.97,98 While the functional are known to be important for describing charge
simulations presented here neglect environmental effects, such transfer excitations42,103−113 and potential energy surfa-
as implicit or explicit solvation, this only has a small effect on ces.114,115
the results presented in this section. It should be stressed that The measured absorption spectrum93 of dye 1 shows a
in situations where simulated results are compared with weakly allowed S1 ← S0 transition and strongly allowed S2 ← S0
experiments measured in polar solvents or where environ- transition. Fluorescence for dye 1 was measured by exciting to
mental interactions are important, solvation effects generally the S2 state, but the observed spectrum is largely obscured by
must be included for accurate results.46,101,102 Table 2 presents the overlapping phosphorescence spectrum. On the basis of the
the calculated absorption and emission maxima along with low intensity and proximity to the absorption feature of the S1
Stokes shifts in comparison to the measured results. All excited state, it is likely that the emission comes from the S1 state rather
state geometry optimizations were performed for the lowest than S2, so the excited state optimization in the TDDFT
energy strongly allowed excitation. This usually corresponded calculations was performed for the S1 state. The intense
to the S1 state of the chromophore. phosphorescence and overlapping onset of the S2 state
The experimentally measured Stokes shifts span values absorption with the S1 state absorption make it difficult to
ranging between nearly coincident absorption and emission determine the Stokes shift experimentally. We therefore
maxima (free pthalocyanine, dye 5) to large shifts of ∼0.80 eV simulated the absorption and fluorescence including vibronic
(2,3,5-trifluorobenzaldehyde, dye 1). In general, DFT captures effects to estimate the Stokes shift.
the correct magnitudes of the Stokes shifts but is sensitive to Figure 3a shows the simulated absorption and fluorescence
the choice of XC functional. The BLYP functional, on average, spectra of dye 1. Simulations for both the syn and anti
yields the worst agreement with experimental measurements conformations of the aldehyde were used, although the anti
with a mean absolute error (MAE) of 0.18 eV and root-mean- conformation is ∼0.10 eV (774 cm−1) more stable. The spectra
square deviation (RMSD) of 0.26 eV. Inclusion of 20% HF were simulated with the 0−0 peak of the vibronic progression
with the B3LYP functional improves the results substantially for absorption and fluorescence coincident with the exper-
with a MAE of 0.07 eV and RMSD of 0.08 eV. The range imental measurements (26 040 and 26 015 cm−1, respectively)
separated hybrid CAM-B3LYP yields comparable performance for both the syn and anti conformations.93 For the anti
with B3LYP, while LC-BLYP (which includes 100% HF conformation, the absorption and fluorescence spectra show a
exchange at large interelectronic separation) slightly worsens progression in the carbonyl stretching frequency (1730 cm−1).
the agreement with experimental results (MAE = 0.10 eV and The peak of the absorption occurs between the third and fourth
RMSD = 0.12 eV). It is suggested through these results that HF peaks, while the peak of fluorescence occurs at the third peak of
exchange is required to improve the calculation of excited state the progression. This yields an estimated Stokes shift between
properties in these systems. 0.86 and 1.07 eV, in agreement with the values determined
In general, the Stokes shift is best described by the global using the TDDFT optimization.
hybrid B3LYP for the test set presented here. However, dye 2 Figure 3b compares the simulated and measured fluorescence
(coronene), dye 3 (coumarin 153), and dye 7 (4′- spectra of dye 1. The analysis presented in ref 93 determined
5496 dx.doi.org/10.1021/ct4007772 | J. Chem. Theory Comput. 2013, 9, 5490−5503
Journal of Chemical Theory and Computation Article

Figure 4. (a) Absorption spectrum for the D3 ← D0 transition of the phenoxyl radical complex simulated with B3LYP/aug-cc-pVDZ (blue line) in
comparison to the experimental spectrum (black line) adapted from ref 124. (b) Simulated resonance Raman spectrum using B3LYP compared with
the experimental spectrum124 at 399 nm.

that the fluorescence almost entirely results from the anti In general, the magnitudes of the differences between simulated
conformation based on DFT calculated frequencies and the and experimental Stokes shifts are comparable to those of
small amount of syn conformers at room temperature. Among section 4.1. However, it is found that the BLYP functional gives
the features marked by arrows in the measured fluorescence the best agreement with experimental data. The overestimation
spectrum, the first, third, fifth, and sixth peaks can be explained of the Stokes shift by the different functionals may be attributed
by emission by the anti conformer. However, the fourth to the neglect of solvent effects on the absorption and
vibronic peak corresponds strongly with the syn conformer. emission.49,125,126
This indicates that emission from both conformers is observed Figure 4 presents the simulated OPA and RRS spectra for the
experimentally, so an isomerization process between the more D3 ← D0 transition of the phenoxyl radical compared with the
stable anti conformer to the syn conformer likely occurs in the experimental spectra of ref 124. The simulated OPA spectrum
excited state. yields many of the observed features in the vibrational
4.2. Optical Processes of the Phenoxyl Radical. The progression, although the features are blueshifted from the
phenoxyl radical is a model system for investigating the tyrosyl positions observed in the experiment. This behavior was also
radical,116 which is important in biological processes such as the observed in a recent study by Santoro et al. regardless of
formation of deoxyribonucleotides,117 water oxidation during whether Duschinsky rotations were included in the simu-
photosynthesis,118 and an important intermediate in the lation.122 Tentatively, we attribute this disagreement to
oxidation of benzene.119 Characterization of the phenoxyl interactions between the solvent and phenoxyl radical that are
radical often proceeds using spectroscopy due to its short not included in the present model. Such effects, including both
lifetime. Simulated and experimental OPA120−122 and explicit and bulk interactions, which have been included in
RRS54,122−124 have facilitated the understanding of the several previous studies of absorption and emission processes,
electronic and structural properties of the phenoxyl radical have been demonstrated to be important.46,49,101,102 The RRS
and enabled characterization of the tyrosyl radical also. spectrum (Figure 4b) shows good agreement with the
Tripathi and Schuler characterized the phenoxyl radical and experimental spectrum. The four most prominent features
its isotopomers using OPA and RRS in aqueous solution (see (528, 990, 1398, and 1505 cm−1) are vibronically coupled to
Figure 4).123,124 The RRS measured on resonance with the 2B1 the electronic transition, and this is captured in the simulation.
← 2A2 transition possesses a broad background between 1300 These results are confirmation that scattering from the
and 1900 cm−1 (between 420 and 440 nm) and is argued to phenoxyl radical results from an A term (Franck−Condon)
result from weak fluorescence of the phenoxyl radical. On the mechanism. Weaker features labeled in the experimental
basis of the experimental OPA and RRS spectra, the Stokes spectrum are absent in the simulations and have been
shift of the phenoxyl radical is 0.21 eV. Table 3 presents demonstrated to arise from Herzberg−Teller coupling terms
simulated results compared with the experimental Stokes shift. that are not included here.122 The peak at 1585 cm−1 in the
simulation is not well resolved in the experiment due to the
broad fluorescence background in that region of the RRS
Table 3. Absorption Maxima, Emission Maxima, and Stokes
spectrum.
Shift (S.S.) for the Phenoxyl Radical Resulting from the D3
4.3. Charge Transfer in a Molecular Host−Guest
← D0 (2B1 ← 2A2) Transition Using BLYP, LC-BLYP,
Complex. The host−guest complex of tetrathiafulvalene
B3LYP, and CAM-B3LYP Using the 6-31G* Basis Seta
(TTF) with cyclobis(paraquat-p-phenylene) (CBPQT4+) is a
method Abs. Em. S.S. model system for studying analyte binding processes. This
BLYP 3.37 3.07 0.30 complex is shown in Figure 5. It was demonstrated previously
B3LYP 3.59 3.27 0.32 that resonance Raman scattering is a powerful analytical
CAM-B3LYP 3.81 3.44 0.37 technique for determining binding energies while also obtaining
LC-PBE 3.85 3.45 0.40 structural information about complexation.127 The formation of
Exp. 3.11 2.90 0.21 the TTF⊂CBPQT4+ complex results in an excited state located
at 865 nm, which is caused by an intermolecular charge transfer
a
Values are given in eV. The experimental data is from ref 124. (CT) state between TTF and CBPQT4+. The simulated
5497 dx.doi.org/10.1021/ct4007772 | J. Chem. Theory Comput. 2013, 9, 5490−5503
Journal of Chemical Theory and Computation Article

4.4. Simulated Spectra for Rhodamine 6G. The

fluorescent dye rhodamine 6G (R6G, shown in Figure 7) has

Figure 5. Structure of the TTF⊂CBPQT4+ host−guest complex.

Figure 7. Structure of the fluorescent dye rhodamine 6G.

resonance Raman spectrum of the complex using BP86 was been the center of many recent Raman scattering studies.
shown to result in intense features at approximately 500 cm−1, TDDFT calculations of static and frequency-dependent
attributed to the poor description of CT states with traditional polarizabilities were used to simulate nonresonant Raman and
DFT XC functionals. Here, we simulate resonance Raman RRS, allowing the resonance enhancement (∼105) to be
spectra of the complex using the global hybrid B3LYP and quantified.128 It has been demonstrated that vibronic effects are
range separated hybrid LC-ωPBE, which we have previously
important for capturing the correct relative peak intensities with
used to study photon driven CT events.115
RRS for R6G also.79−81 TDDFT simulations with vibronic
Figure 6a shows the simulated absorption spectra of
coupling have also been instrumental in understanding the
TTF⊂CBPQT4+. While B3LYP calculations find the absorption
complicated nature of resonance hyper-Raman scattering and
maximum to be 1.13 eV, LC-ωPBE determines the absorption
maximum to be 2.09 eV. The experimental absorption surface-enhanced hyper-Raman scattering in R6G.15,83,129
maximum is 1.43 eV, so both calculated results are shifted We examined the strongly allowed S1 ← S0 transition with
such that the absorption maximum coincides with the different XC functionals. Because each functional yields a
experimental absorption maximum. Simulations with both XC different excitation energy, the OPA spectra are shifted so that
functionals are capable of reproducing the experimental the absorption maximum is aligned with the experimental
spectrum. However, the blue edge of the spectrum is spectrum. Figure 8 shows the simulated OPA spectrum with
underestimated. each functional. The functionals are split into four panels and
Figure 6b illustrates the simulated resonance Raman spectra show very different behaviors. For the local and semilocal
of the complex in comparison with the experimental spectrum. functionals (LDA, BP86, and PBE), the spectra show a main
Both B3LYP and LC-ωPBE yield spectra comparable to peak with a vibronic shoulder at approximately 490 nm. Each
experimental data, in contrast with the simulated BP86 functional in this case overestimates the intensity of the
spectrum in ref 127. HF exchange is therefore important for vibronic shoulder compared to the experimental spectrum.
describing relative peak intensities in resonance Raman spectra. The global hybrid (B3LYP and PBE0) functionals show very
Strong features are observed at 471 cm−1 (TTF S atom ring similar behavior. Both functionals significantly improve the
breathing), 701 cm −1 (TTF C−S stretch), 1477 cm −1 description of the vibronic feature in the OPA spectrum.
(CBPQT4+ CH2 scissoring), and 1637 cm−1 (CBPQT4+ pyridyl Interestingly, the long-range corrected hybrid (CAM-B3LYP,
C−C stretching). For the simulated B3LYP spectrum, the BNL, LC-PBE, LC-PBE0, LC-ωPBE, LC-ωPBEh) functionals
mode at 1637 cm−1 is overestimated significantly. The same are less consistent in their description of the OPA spectrum.
mode is underestimated slightly with LC-ωPBE. In general, the BNL and LC-PBE0 yield poor spectra in comparison with
relative peak intensities using LC-ωPBE are in better agreement experimental results, where the vibronic structure is more
with experimental results, illustrating the importance of range intense than the main peak leading to a broad feature. This is
separated hybrids for predicting vibronic features for CT states. contrasted by the other long-range corrected hybrids, which

Figure 6. (a) Absorption spectra of the TTF⊂CBPQT4+ complex simulated with B3LYP and LC-ωPBE in comparison to the experimental spectrum
(black line). (b) Simulated resonance Raman spectra using B3LYP and LC-ωPBE compared with the experimental spectrum127 at 785 nm.
Acetonitrile solvent lines in the resonance Raman spectrum are designated with MeCN.

5498 dx.doi.org/10.1021/ct4007772 | J. Chem. Theory Comput. 2013, 9, 5490−5503

Journal of Chemical Theory and Computation Article

Figure 8. Absorption spectra of rhodamine 6G using different XC functionals in comparison with experimental data.130 Spectra are shifted so that
the absorption maximum coincides with the experimental absorption maximum.

Figure 9. Resonance Raman spectra of rhodamine 6G using different XC functionals. Simulations are performed assuming a wavelength of 514.5 nm.

have spectra comparable to the global hybrid functionals. It was these methods. Surprisingly, however, CAM-B3LYP, LC-PBE,
thought, because long-range corrected hybrids generally use LC-ωPBE, and LC-ωPBEh do not possess problems in the
100% HF exchange at large interelectronic separations, that a OPA spectrum, so further investigation requires the vibrational
pure HF calculation might elaborate on deficiencies of some of information provided by RRS.
these functionals. The HF calculation gives similar results to Figure 9 shows the RRS spectra of rhodamine 6G using the
BNL and LC-PBE0, demonstrating the connection between same functionals. For comparison, the experimental surface-
5499 dx.doi.org/10.1021/ct4007772 | J. Chem. Theory Comput. 2013, 9, 5490−5503
Journal of Chemical Theory and Computation Article

enhanced Raman scattering (SERS) spectrum is shown, which

is measured at 514.5 nm. Although many peaks are present in
the SERS spectrum, we focus on the modes at 616, 772, 1356,
1515, 1579, and 1656 cm−1. Aside from the mode at 1663
cm−1, the local and semilocal functionals underestimate the
intensity of every feature in the SERS spectrum. The global
hybrid, B3LYP, which we have previously used in other studies
to obtain good agreement with experimental results,15,83
compares very well with the experimental spectrum. However,
the modes at 1373 and 1518 cm−1 are underestimated.
Although the global hybrids PBE0 and B3LYP use similar
amounts of HF exchange (25% and 20%, respectively), the
simulated spectra from both functionals vary substantially
below 1000 cm−1. Many features that are not observed in the Figure 10. Resonance Raman spectra of rhodamine 6G calculated with
SERS spectrum appear in the PBE0 spectrum. Examining the B3LYP. The figures are (a) analytical gradient, (b) analytical gradient
semilocal functionals BP86 and PBE, we observe that the with displacement of the 1656 mode scaled by 70%, (c) numerical
spectra are comparable. It may be surprising then that B3LYP gradient, and (d) simulated RRS with parameters from ref 59.
and PBE0 give very different looking spectra, since B3LYP uses Simulations were performed at 532 nm.
the same exchange functional as BP86 and PBE0 uses the same
exchange functional as PBE. Tentatively, we attribute these found that numerical evaluation of the vibronic coupling
differences to sensitivity of the different XC functionals to the parameters using B3LYP/6-311G* gives good agreement with
structure and normal modes used in the calculations as well as experimental spectra. It is expected that analytical gradients are
to the basis sets. For different functionals, the ground state more accurate than numerical gradients, so the better
optimized structure and normal modes do not need to be agreement using numerical gradients with experimental results
identical, and this may impact the quality of the response is likely fortuitous.
Table 4 shows a comparison of the dimensionless displace-
properties we calculate. Further work is being done to quantify
ments from the numerical and analytical methods. In general,
the effect of the geometry and basis set choice for large
molecules.
RRS spectra using the long-range corrected hybrids and HF Table 4. Dimensionless Displacements for Important Modes
are split into two groups. The functionals CAM-B3LYP, LC- in Rhodamine 6G Using Numerical and Analytical Gradients
PBE, LC-ωPBE, and LC-ωPBEh give similar spectra. Recall Determined Using B3LYPa
that these functionals also yielded reasonable OPA spectra in frequency (cm−1) numerical Δ analytical Δ absolute difference (%)
comparison to the experimental data (Figure 8). Although the 616 −0.24 −0.24 0.00 (0.0)
ratios of the peaks in the RRS spectra are sensitive to the 772 −0.16 −0.12 0.04 (33.3)
functional used, it appears that the observed modes are 1129 −0.12 −0.10 0.02 (20.0)
identical for these functionals. The other functionals, BNL, LC- 1182 0.09 0.07 0.02 (28.6)
PBE0, and HF, were characterized by poor agreement with the 1195 −0.10 −0.09 0.01 (10.0)
experimental OPA spectrum. From the RRS spectra, these 1298 −0.14 −0.12 0.02 (16.7)
functionals each have a large peak at 855 and 1730 cm−1. The 1356 0.18 0.17 0.01 (6.0)
broad vibronic features in the OPA spectra for BNL, LC-PBE0, 1515 −0.14 −0.13 0.01 (7.7)
and HF are likely associated with the 1730 cm−1 mode, since 1561 0.14 0.19 0.05 (26.3)
that feature differs from the other functionals. 1579 0.12 0.11 0.01 (9.1)
To better understand the analytical gradient results with 1656 −0.20 −0.31 0.11 (35.5)
B3LYP, we compare them to a RRS spectrum simulated with a
Percent differences are determined relative to the analytical gradient
numerical gradients. Figure 10a and c show the comparison of result.
the analytical and numerical gradient simulations with a
simulation using parameters from ref 59 (Figure 10d). We
see that the same modes are predicted to couple to the S1 ← S0 the parameters are similar. The modes at 1561 and 1656 cm−1
transition in the analytical and numerical methods, although the have the largest absolute deviations in the two methods. These
relative intensities differ. The analytical gradient method yields differences are comparable to those calculated between the
a spectrum comparable to the RRS spectrum reported in a analytical and numerical gradient methods by Petrenko and
previous study by our group.79 In comparison to the Neese for trans-hexatriene with the same functional.53,78
experimental RRS spectrum, the data generated from numerical To obtain better agreement between the two methods, we
gradients are in better agreement (the peak at ∼375 cm−1 in the tried scaling the intensity of the 1656 mode (Figure 10b). By
experimental spectrum was not included in the simulations). A reducing the intensity of that mode, the analytical gradient
possible explanation for this behavior is that the analytical spectrum becomes comparable to the numerical one. Although
gradient depends on two perturbations (the external electric the spectra still do not overlap perfectly, the agreement
field and change in geometry), meaning that it is a quadratic between both methods is improved by scaling the 1656 mode
response property. Previous work has demonstrated that by 70%. Small differences in the other relative peak intensities
quadratic response properties are sensitive to the size of the are caused by the variations in the dimensionless displacements
basis set,131,132 so the gradients and vibronic coupling (Table 4).
parameters may change significantly depending on the quality 4.5. Timings. In Table 5, CPU times for a TDDFT energy
of the basis set. Guthmuller and Champagne80,81 similarly + gradient and TDDFT gradient are presented for the P2TA
5500 dx.doi.org/10.1021/ct4007772 | J. Chem. Theory Comput. 2013, 9, 5490−5503
Journal of Chemical Theory and Computation Article

Table 5. CPU Timings (s) for a TDDFT Energy + Gradient multiprogram national laboratory operated for DOE by the
of P2TA (Dye 10 in Figure 2) and TDDFT Gradient of the Battelle Memorial Institute under Contract DE-AC06-76RLO-
First Excited Statea 1830. D.W.S. and L.J. acknowledge support from the U.S.
National Science Foundation under Grant Number CHE-
NCPUs TDDFT energy + gradient TDDFT gradient
0955689. N.G. acknowledges support from the U.S. Depart-
32 3260 1477 ment of Energy, Office of Basic Energy Sciences, under Grant
64 1714 742 No. DESC0008666 of the SciDAC program.

■
128 969 377
256 687 224
REFERENCES
a
The TDDFT energy was calculated for the lowest five states.
Calculations were performed using B3LYP/6-31G* (864 basis (1) Morton, S. M.; Silverstein, D. W.; Jensen, L. Chem. Rev. 2011,
functions). 111, 3962−3994.
(2) Helgaker, T.; Coriani, S.; Jørgensen, P.; Kristensen, K.; Olsen, J.;
Ruud, K. Chem. Rev. 2012, 112, 543−631.
(3) Champagne, B.; Plaquet, A.; Pozzo, J.-L.; Rodriguez, V.; Castet, F.
molecule (dye 10 in Figure 2) using the B3LYP functional and J. Am. Chem. Soc. 2012, 134, 8101−8103.
6-31G* basis set (864 basis functions). Although the (4) Kulesza, A.; Mitrić, R.; Bonačić-Koutecký, V. J. Phys. Chem. A
implementation described here is still a development version, 2009, 113, 3783−3788.
good scaling is observed up to 128 processors for this molecule. (5) Nielsen, K. A.; Sarova, G. H.; Martín-Gomis, L.; Stein, P. C.;
These timing tests were performed on the Chinook computer Sanguinet, L.; Levillain, E.; Sessler, J. L.; Guldi, D. M.; Sastre-Santos,
system at EMSL, PNNL.133 Á .; Jeppesen, J. O. J. Am. Chem. Soc. 2008, 130, 460−462.
(6) Yuan, H.; Chong, H.; Wang, B.; Zhu, C.; Liu, L.; Yang, Q.; Lv, F.;
5. CONCLUSION Wang, S. J. Am. Chem. Soc. 2012, 134, 13184−13187.
(7) Zhang, J.; Wong, K.-L.; Wong, W.-K.; Mak, N.-K.; Kwong, D. W.
In this work, we have presented an implementation of analytical
J.; Tam, H.-L. Org. Biomol. Chem. 2011, 9, 6004−6010.
excitation energy gradients within the TDDFT framework. As (8) Uppal, R.; Ciesienski, K. L.; Chonde, D. B.; Loving, G. S.;
part of our validation, we calculated the Stokes shifts for a range Caravan, P. J. Am. Chem. Soc. 2012, 134, 10799−10802.
of organic dye molecules using a diverse set of XC functionals (9) Liu, Z.; Zhang, C.; He, W.; Qian, F.; Yang, X.; Gao, X.; Guo, Z.
including traditional DFT functionals, global hybrids, and New J. Chem. 2010, 34, 656.
range-separated hybrids. We also simulated the OPA and RRS (10) Picot, A.; D’Aléo, A.; Baldeck, P. L.; Grichine, A.; Duperray, A.;
of the well-studied dye molecule rhodamine 6G, phenoxyl Andraud, C.; Maury, O. J. Am. Chem. Soc. 2008, 130, 1532−1533.
radical, and a molecular host−guest complex (11) Egolf, D. S.; Waterland, M. R.; Myers Kelley, A. J. Phys. Chem. B
(TTF⊂CBPQT4+). Our calculations demonstrate the impor- 2000, 104, 10727−10737.
tance of analytic TDDFT gradients for efficient computation of (12) Kleinman, S. L.; Ringe, E.; Valley, N.; Wustholz, K. L.; Phillips,
E.; Scheidt, K. A.; Schatz, G. C.; Van Duyne, R. P. J. Am. Chem. Soc.
vibronic coupling parameters. Our results show that the RRS
2011, 133, 4115−4122.
depends strongly on the XC functional used and thus serves as (13) Kamarchik, E.; Krylov, A. I. J. Phys. Chem. Lett. 2011, 2, 488−
a stringent benchmark of the accuracy of the XC functionals. 492.
For the molecule, rhodamine 6G, we observe that the best (14) Petrenko, T.; Krylova, O.; Neese, F.; Sokolowski, M. New J.
description of the potential energy surface is obtained using the Phys. 2009, 11, 015001.
global hybrid B3LYP functional. Surprisingly, we find that the (15) Milojevich, C. B.; Silverstein, D. W.; Jensen, L.; Camden, J. P. J.
long-range corrected functionals provide a rather poor Am. Chem. Soc. 2011, 133, 14590−14592.
description of the RRS of rhodamine 6G. (16) Runge, E.; Gross, E. K. U. Phys. Rev. Lett. 1984, 52, 997−1000.

■
(17) Casida, M. E. In Recent Advances in Density Functional Methods,
ASSOCIATED CONTENT Part I; Chong, D. P., Ed.; World Scientific: Singapore, 1995; Recent
Advances in Computational Chemistry, Vol. 1, Chapter 5, pp 155−
* Supporting Information
S
192.
Optimized structures of various dye molecules in the ground (18) Casida, M. E. In Recent Developments and Applications of Modern
state using the B3LYP functional. This material is available free Density Functional Theory; Seminario, J. M., Ed.; Elsevier: Amsterdam,
of charge via the Internet at http://pubs.acs.org/.

■
1996; Theoretical and Computational Chemistry, Vol. 4, pp 391−439.
(19) Hirata, S.; Head-Gordon, M. Chem. Phys. Lett. 1999, 314, 291.
AUTHOR INFORMATION (20) Stratmann, R. E.; Scuseria, G. E.; Frisch, M. J. J. Chem. Phys.
Corresponding Authors 1998, 109, 8218−8224.
*E-mail: niri.govind@pnl.gov. (21) Masunov, A.; Tretiak, S. J. Phys. Chem. B 2004, 108, 899−907.
*E-mail: jensen@chem.psu.edu. (22) Zhao, L. L.; Jensen, L.; Schatz, G. C. Nano Lett. 2006, 6, 1229−
1234.
Notes (23) Aikens, C. M. J. Phys. Chem. A 2009, 113, 10811−10817.
The authors declare no competing financial interest.

■
(24) Marques, M. A. L.; Gross, E. K. U. Annu. Rev. Phys. Chem. 2004,
55, 427−455.
ACKNOWLEDGMENTS (25) Casida, M.; Huix-Rotllant, M. Annu. Rev. Phys. Chem. 2012, 63,
The authors thank Jon Camden (University of Tennessee, 287−323.
Knoxville) for the SERS spectra of rhodamine 6G and Amar (26) Ray, K.; DeBeer George, S.; Solomon, E. I.; Wieghardt, K.;
Neese, F. Chem.Eur. J. 2007, 13, 2783−2797.
Flood (Indiana University) for the OPA and RRS spectra of the (27) Stener, M.; Fronzoni, G.; de Simone, M. Chem. Phys. Lett. 2003,
TTF⊂CBPQT4+ host−guest complex. A portion of the work 373, 115−123.
was performed using EMSL, a national scientific user facility (28) Besley, N. A.; Asmuruf, F. A. Phys. Chem. Chem. Phys. 2010, 12,
sponsored by the U.S. Department of Energy, Office of 12024−12039.
Biological and Environmental Research, and located at the (29) Liang, W.; Fischer, S. A.; Frisch, M. J.; Li, X. J. Chem. Theory
Pacific Northwest National Laboratory (PNNL). PNNL is a Comput. 2011, 7, 3540−3547.

5501 dx.doi.org/10.1021/ct4007772 | J. Chem. Theory Comput. 2013, 9, 5490−5503

Journal of Chemical Theory and Computation Article

(30) Lopata, K.; Van Kuiken, B. E.; Khalil, M.; Govind, N. J. Chem. (68) Zong, Y.; McHale, J. L. J. Chem. Phys. 1997, 107, 2920−2929.
Theory Comput. 2012, 8, 3284−3292. (69) Valiev, M.; Bylaska, E.; Govind, N.; Kowalski, K.; Straatsma, T.;
(31) Van Kuiken, B. E.; Valiev, M.; Daifuku, S. L.; Bannan, C.; Van Dam, H.; Wang, D.; Nieplocha, J.; Apra, E.; Windus, T.; de Jong,
Strader, M. L.; Cho, H.; Huse, N.; Schoenlein, R. W.; Govind, N.; W. Comput. Phys. Commun. 2010, 181, 1477−1489.
Khalil, M. J. Phys. Chem. A 2013, 117, 4444−4454. (70) Heller, E. J. J. Chem. Phys. 1975, 62, 1544−1555.
(32) Van Caillie, C.; Amos, R. D. Chem. Phys. Lett. 1999, 308, 249− (71) Heller, E. J. J. Chem. Phys. 1976, 64, 63−73.
255. (72) Lee, S.-Y.; Heller, E. J. J. Chem. Phys. 1979, 71, 4777−4788.
(33) Van Caillie, C.; Amos, R. D. Chem. Phys. Lett. 2000, 317, 159− (73) Tannor, D. J.; Heller, E. J. J. Chem. Phys. 1982, 77, 202−218.
164. (74) Heller, E. J. Acc. Chem. Res. 1981, 14, 368−375.
(34) Amos, R. D. Chem. Phys. Lett. 2002, 364, 612−615. (75) Heller, E. J.; Sundberg, R.; Tannor, D. J. Phys. Chem. 1982, 86,
(35) Burcl, R.; Amos, R. D.; Handy, N. C. Chem. Phys. Lett. 2002, 1822−1833.
355, 8−18. (76) Duschinsky, F. Acta Phys. U.R.S.S. 1937, 7, 551−556.
(36) Furche, F.; Ahlrichs, R. J. Chem. Phys. 2002, 117, 7433−7447. (77) Neugebauer, J.; Hess, B. A. J. Chem. Phys. 2004, 120, 11564.
(37) Rappoport, D.; Furche, F. J. Chem. Phys. 2005, 122, 064105. (78) Petrenko, T.; Neese, F. J. Chem. Phys. 2007, 127, 164319.
(38) Hutter, J. J. Chem. Phys. 2003, 118, 3928−3934. (79) Kane, K. A.; Jensen, L. J. Phys. Chem. C 2010, 114, 5540−5546.
(39) Odelius, M.; Laikov, D.; Hutter, J. J. Mol. Struct.: THEOCHEM (80) Guthmuller, J.; Champagne, B. ChemPhysChem 2008, 9, 1667−
2003, 630, 163−175. 1669.
(40) Seth, M.; Mazur, G.; Ziegler, T. Theor. Chem. Acc. 2011, 129, (81) Guthmuller, J.; Champagne, B. J. Phys. Chem. A 2008, 112,
331−342. 3215−3223.
(41) Petrenko, T.; Kossmann, S.; Neese, F. J. Chem. Phys. 2011, 134, (82) Jarzȩki, A. A. J. Phys. Chem. A 2009, 113, 2926−2934.
054116. (83) Milojevich, C. B.; Silverstein, D. W.; Jensen, L.; Camden, J. P. J.
(42) Chiba, M.; Tsuneda, T.; Hirao, K. J. Chem. Phys. 2006, 124, Phys. Chem. C 2013, 117, 3046−3054.
144106. (84) Reiher, M.; Neugebauer, J.; Hess, B. A. Z. Phys. Chem. 2003,
(43) Nguyen, K. A.; Day, P. N.; Pachter, R. Int. J. Quantum Chem. 217, 91−103.
2010, 110, 2247−2255. (85) Silverstein, D. W. Applications of Time-dependent Quantum
(44) Liu, F.; Gan, Z.; Shao, Y.; Hsu, C.-P.; Dreuw, A.; Head-Gordon, Mechanics to Resonantly-enhanced Linear and Nonlinear Optical
M.; Miller, B. T.; Brooks, B. R.; Yu, J.-G.; Furlani, T. R.; Kong, J. Mol. Processes. Ph.D. thesis, The Pennsylvania State University, University
Phys. 2010, 108, 2791−2800. Park, PA, 2013.
(45) Liu, J.; Liang, W. Z. J. Chem. Phys. 2011, 134, 044114. (86) Silverstein, D. W.; Jensen, L. TDSPEC: Time-Dependent
(46) Scalmani, G.; Frisch, M. J.; Mennucci, B.; Tomasi, J.; Cammi, R.; Spectroscopic Simulations of Linear and Nonlinear Optical Processes.
Barone, V. J. Chem. Phys. 2006, 124, 094107. The Pennsylvania State University, University Park, PA, 2011.
(47) Marenich, A. V.; Cramer, C. J.; Truhlar, D. G.; Guido, C. A.; (87) Grimme, S. J. Comput. Chem. 2006, 27, 1787−1799.
Mennucci, B.; Scalmani, G.; Frisch, M. J. Chem. Sci. 2011, 2, 2143− (88) Besler, B. H.; Scuseria, G. E.; Scheiner, A. C.; Schaefer, H. F. J.
2161.
Chem. Phys. 1988, 89, 360−366.
(48) Si, D.; Li, H. J. Chem. Phys. 2010, 133, 144112.
(89) Stanton, J. F.; Watts, J. D.; Bartlett, R. J. J. Chem. Phys. 1991, 94,
(49) Minezawa, N.; Silva, N. D.; Zahariev, F.; Gordon, M. S. J. Chem.
404−413.
Phys. 2011, 134, 054111.
(90) Scott, A. P.; Radom, L. J. Phys. Chem. 1996, 100, 16502−16513.
(50) Doltsinis, N. L.; Kosov, D. S. J. Chem. Phys. 2005, 122, 144101.
(91) Alecu, I. M.; Zheng, J.; Zhao, Y.; Truhlar, D. G. J. Chem. Theory
(51) Ismail-Beigi, S.; Louie, S. G. Phys. Rev. Lett. 2003, 90, 076401.
(52) Sitt, A.; Kronik, L.; Ismail-Beigi, S.; Chelikowsky, J. R. Phys. Rev. Comput. 2010, 6, 2872−2887.
(92) Kowalczyk, T.; Tsuchimochi, T.; Chen, P.-T.; Top, L.; Voorhis,
A 2007, 76, 054501.
(53) Petrenko, T.; Neese, F. J. Chem. Phys. 2012, 137, 234107. T. V. J. Chem. Phys. 2013, 138, 164101.
(54) Ma, H.; Liu, J.; Liang, W. J. Chem. Theory Comput. 2012, 8, (93) Itoh, T. J. Mol. Struct. 2010, 261, 53−58.
4474−4482. (94) Itoh, T. J. Mol. Struct. 2008, 252, 115−120.
(55) Savarese, M.; Aliberti, A.; De Santo, I.; Battista, E.; Causa, F.; (95) Ernsting, N. P.; Asimov, M.; Schäfer, F. P. Chem. Phys. Lett.
Netti, P. A.; Rega, N. J. Phys. Chem. A 2012, 116, 7491−7497. 1982, 91, 231−235.
(56) Avila Ferrer, F. J.; Cerezo, J.; Stendardo, E.; Improta, R.; (96) Haucke, G.; Graness, G. Angew. Chem., Int. Ed. 1995, 34, 67−68.
Santoro, F. J. Chem. Theory Comput. 2013, 9, 2072−2082. (97) Creemers, T. M. H.; Lock, A. J.; Subramaniam, V.; Jovin, T. M.;
(57) Kupfer, S.; Guthmuller, J.; González, L. J. Chem. Theory Comput. Völker, S. Nat. Struct. Biol. 1999, 6, 557−560.
2012, 9, 543−554. (98) Ma, Y.; Rohlfing, M.; Molteni, C. J. Chem. Theory Comput. 2009,
(58) Poizat, O.; Sourisseau, C. J. Phys. Chem. 1984, 88, 3007−3014. 6, 257−265.
(59) Shim, S.; Stuart, C. M.; Mathies, R. A. ChemPhysChem 2008, 9, (99) Williams, R.; Goldsmith, G. J. J. Chem. Phys. 1963, 39, 2008−
697−699. 2011.
(60) Angeloni, L.; Smulevich, G.; Marzocchi, M. P. J. Raman (100) Eastwood, D.; Edwards, L.; Gouterman, M.; Steinfeld, J. J. Mol.
Spectrosc. 1979, 8, 305−310. Struct. 1966, 20, 381−390.
(61) Lueck, H. B.; McHale, J. L.; Edwards, W. D. J. Am. Chem. Soc. (101) Jacquemin, D.; Planchat, A.; Adamo, C.; Mennucci, B. J. Chem.
1992, 114, 2342−2348. Theory Comput. 2012, 8, 2359−2372.
(62) Ibrahim, M.; Xu, C.; Spiro, T. G. J. Am. Chem. Soc. 2006, 128, (102) Silva, D. L.; Murugan, N. A.; Kongsted, J.; Rinkevicius, Z.;
16834−16845. Canuto, S.; Ågren, H. J. Phys. Chem. B 2012, 116, 8169−8181.
(63) Ruban, A. V.; Berera, R.; Ilioaia, C.; van Stokkum, I. H. M.; (103) Lange, A. W.; Herbert, J. M. J. Am. Chem. Soc. 2009, 131,
Kennis, J. T. M.; Pascal, A. A.; van Amerongen, H.; Robert, B.; Horton, 3913−3922.
P.; van Grondelle, R. Nature 2007, 450, 575−578. (104) Silverstein, D. W.; Jensen, L. J. Chem. Phys. 2010, 132, 194302.
(64) Balakrishnan, G.; Jarzecki, A. A.; Wu, Q.; Kozlowski, P. M.; (105) Govind, N.; Valiev, M.; Jensen, L.; Kowalski, K. J. Phys. Chem.
Wang, D.; Spiro, T. G. J. Phys. Chem. B 2012, 116, 9387−9395. A 2009, 113, 6041−6043.
(65) Markel, F.; Ferris, N. S.; Gould, I. R.; Myers, A. B. J. Am. Chem. (106) Jensen, L.; Govind, N. J. Phys. Chem. A 2009, 113, 9761−9765.
Soc. 1992, 114, 6208−6219. (107) Andzelm, J.; Rinderspacher, B. C.; Rawlett, A.; Dougherty, J.;
(66) Lilichenko, M.; Tittelbach-Helmrich, D.; Verhoeven, J. W.; Baer, R.; Govind, N. J. Chem. Theory Comput. 2009, 5, 2835−2846.
Gould, I. R.; Myers, A. B. J. Chem. Phys. 1998, 109, 10958−10969. (108) Kowalski, K.; Krishnamoorthy, S.; Villa, O.; Hammond, J. R.;
(67) Zong, Y.; McHale, J. L. J. Chem. Phys. 1997, 106, 4963−4972. Govind, N. J. Chem. Phys. 2010, 132, 154103.

5502 dx.doi.org/10.1021/ct4007772 | J. Chem. Theory Comput. 2013, 9, 5490−5503

Journal of Chemical Theory and Computation Article

(109) Glaesemann, K. R.; Govind, N.; Krishnamoorthy, S.; Kowalski,

K. J. Phys. Chem. A 2010, 114, 8764−8771.
(110) Lopata, K.; Govind, N. J. Chem. Theory Comput. 2011, 7,
1344−1355.
(111) Lopata, K.; Reslan, R.; Kowalska, M.; Neuhauser, D.; Govind,
N.; Kowalski, K. J. Chem. Theory Comput. 2011, 7, 3686−3693.
(112) Refaely-Abramson, S.; Sharifzadeh, S.; Govind, N.; Autschbach,
J.; Neaton, J. B.; Baer, R.; Kronik, L. Phys. Rev. Lett. 2012, 109, 226405.
(113) Stein, T.; Autschbach, J.; Govind, N.; Kronik, L.; Baer, R. J.
Phys. Chem. Lett. 2012, 3, 3740−3744.
(114) Dreuw, A.; Head-Gordon, M. J. Am. Chem. Soc. 2003, 126,
4007−4016.
(115) Silverstein, D. W.; Jensen, L. J. Chem. Theory Comput. 2010, 6,
2845−2855.
(116) Sander, W.; Roy, S.; Polyak, I.; Ramirez-Anguita, J. M.;
Sanchez-Garcia, E. J. Am. Chem. Soc. 2012, 134, 8222−8230.
(117) Engström, M.; Himo, F.; Gräslund, A.; Minaev, B.; Vahtras, O.;
Ågren, H. J. Phys. Chem. A 2000, 104, 5149−5153.
(118) Nugent, J. H.; Ball, R. J.; Evans, M. C. Biochim. Biophys. Acta
2004, 1655, 217−221.
(119) Herbinet, O.; Husson, B.; Ferrari, M.; Glaude, P.-A.; Battin-
Leclerc, F. Proc. Combust. Inst. 2013, 34, 297−305.
(120) Radziszewski, J. G.; Gil, M.; Gorski, A.; Spanget-Larsen, J.;
Waluk, J.; Mróz, B. J. J. Chem. Phys. 2001, 115, 9733−9738.
(121) Dierksen, M.; Grimme, S. J. Chem. Phys. 2004, 120, 3544−
3554.
(122) Santoro, F.; Cappelli, C.; Barone, V. J. Chem. Theory Comput.
2011, 7, 1824−1839.
(123) Tripathi, G.; Schuler, R. H. Chem. Phys. Lett. 1983, 98, 594−
596.
(124) Tripathi, G. N. R.; Schuler, R. H. J. Chem. Phys. 1984, 81, 113.
(125) Doltsinis, N. L.; Markwick, P. R. L.; Nieber, H.; Langer, H. In
Challenges and Advances In Computational Chemistry and Physics;
Shukla, M. K., Leszczynski, J., Eds.; Springer: Netherlands, 2008; Vol.
5, pp 265−299.
(126) Gustavsson, T.; Bányász, Á .; Lazzarotto, E.; Markovitsi, D.;
Scalmani, G.; Frisch, M. J.; Barone, V.; Improta, R. J. Am. Chem. Soc.
2005, 128, 607−619.
(127) Witlicki, E. H.; Hansen, S. W.; Christensen, M.; Hansen, T. S.;
Nygaard, S. D.; Jeppesen, J. O.; Wong, E. W.; Jensen, L.; Flood, A. H.
J. Phys. Chem. A 2009, 113, 9450−9457.
(128) Jensen, L.; Schatz, G. C. J. Phys. Chem. A 2006, 110, 5973−
5977.
(129) Milojevich, C. B.; Silverstein, D. W.; Jensen, L.; Camden, J. P.
ChemPhysChem 2011, 12, 101−103.
(130) Dixon, J. M.; Taniguchi, M.; Lindsey, J. S. Photochem. Photobiol.
2005, 81, 212−213.
(131) Sałek, P.; Vahtras, O.; Helgaker, T.; Ågren, H. J. Chem. Phys.
2002, 117, 9630−9645.
(132) Paterson, M. J.; Christiansen, O.; Pawlowski, F.; Jørgensen, P.;
Hättig, C.; Helgaker, T.; Sałek, P. J. Chem. Phys. 2006, 124, 054322.
(133) Chinook computer system. See http://www.emsl.pnl.gov/
capabilities/computing/msc/supercomputing/supercomputer.jsp.

5503 dx.doi.org/10.1021/ct4007772 | J. Chem. Theory Comput. 2013, 9, 5490−5503