Longest-Wavelength Electronic Excitations of Linear Cyanines: The

Sep 18, 2013 - Barry Moore , II , Haitao Sun , Niranjan Govind , Karol Kowalski , and Jochen ..... Z. C. Wong , W. Y. Fan , T. S. Chwee , Michael B. S...
0 downloads 0 Views 1MB Size
Article pubs.acs.org/JCTC

Longest-Wavelength Electronic Excitations of Linear Cyanines: The Role of Electron Delocalization and of Approximations in TimeDependent Density Functional Theory Barry Moore II and Jochen Autschbach* Department of Chemistry, University at Buffalo, State University of New York, Buffalo, New York 14260-3000, United States S Supporting Information *

ABSTRACT: The lowest-energy/longest-wavelength electronic singlet excitation energies of linear cyanine dyes are examined, using time-dependent density functional theory (TDDFT) and selected wave function methods in comparison with literature data. Variations of the bond-length alternation obtained with different optimized structures produce small differences of the excitation energy in the limit of an infinite chain. Hybrid functionals with range-separated exchange are optimally ‘tuned’, which is shown to minimize the delocalization error (DE) in the cyanine π systems. Much unlike the case of charge-transfer excitations, small DEs are not strongly correlated with better performance. A representative cyanine is analyzed in detail. Compared with accurate benchmark data, TDDFT with ‘pure’ local functionals gives too high singlet excitation energies for all systems, but DFT-based ΔSCF calculations with a local functional severely underestimates the energies. TDDFT strongly overestimates the difference between singlet and triplet excitation energies. An analysis points to systematically much too small magnitudes of integrals from the DFT components of the exchange-correlation response kernel as the likely culprit. The findings support previous suggestions that the differential correlation energy between the ground and excited state is not correctly produced by TDDFT with most functionals.

1. INTRODUCTION Electronic excitation energies, when calculated with timedependent (TD) density functional theory (DFT) linear response methods, are often said to have a tendency to be too low,1−3 in particular, with nonhybrid functionals. The reason lies in well publicized deficiencies inherent in TDDFT when used in conjunction with approximate functionals. However, there are some exceptions for various classes of excitations. Chiefly among those are the longest-wavelength/lowest-energy excitations of cyanines,4−9 for which TDDFT overestimates the singlet excitation energies. Formulas for simple linear cyanines are shown in Figure 1. Part of the cyanine problem can be attributed to the fact that calculated vertical excitation energies are not directly comparable with the energies corresponding to the experimental absorption band maxima, but the overestimation persists when comparing TDDFT with accurate theoretical vertical excitation energies.8 Likewise, the singlet excitations of related delocalized organic π chromophores such as rhodamine and rosamine dyes10 are overestimated with TDDFT and standard functionals. In both cases, the overestimation is least pronounced with nonhybrid functionals (around a few tenths of an electron volt) but can become larger when supposedly more accurate hybrid functionals such as B3LYP or PBE0 are used. However, trends among different rhodamines and rosamines were better reproduced with hybrid functionals.10 The lowestenergy electronic excitations of these systems, including the © XXXX American Chemical Society

Figure 1. Linear cyanine series CNn, with k = 1, 2, 3, 4 and n = 2k + 3 = 5, 7, 9, 11. Two equivalent resonance structures are shown. The number n is the sum of C and N atoms in the π backbone.

cyanine dyes, are assigned to essentially pure transitions from the highest occupied molecular orbital (HOMO) to the lowest unoccupied MO (LUMO) of the π system. Other cases where TDDFT has a tendency to overestimate excitation energies are ligand−field transitions for complexes of Co(III), Rh(III), and Received: July 23, 2013

A

dx.doi.org/10.1021/ct400649r | J. Chem. Theory Comput. XXXX, XXX, XXX−XXX

Journal of Chemical Theory and Computation

Article

similar systems,11,12 but here, the trends for nonhybrid versus hybrid functionals are different. There is no obvious red flag for the HOMO−LUMO (HL) transitions of cyanines, since such a valence excitation is not typically regarded as problematic for TDDFT. This is much unlike charge-transfer excitations whose energies can be grossly underestimated with popular standard functionals. Moreover, any CT-like character of a HL transition of an extended π system is expected to result in a calculated excitation energy that is too low13 rather than too high. Curiously, for cyanine dyes, the HL excitation energy is very well described by the simplest one-dimensional quantum model system, the particle in a box (PB), with a single empirical parameter to account for the spatial extent of the terminal groups.14,15 According to this model, and in reasonable agreement with experimental data for series of finite-length cyanines, the HL absorption wavelength increases linearly with the chain length and approaches infinity for the limit of an infinite chain (i.e., the excitation energy converges to zero). The linear relationship between absorption wavelength and chromophore chain length has been referred to as the vinyl shift.16 The PB model can be extended easily in order to account for bond-length alternation (BLA) in a π chromophore. BLA then creates a chainlength independent additional contribution to the HL gap and, consequently, a limiting finite excitation energy/finite wavelength for long π chromophore backbones. The orange to red colors of long-chain polyenes such as carotene and lycopene are examples for the presence of a limiting absorption wavelength for π-chromophores with sizable BLA. Cyanines cover a much larger spectral range, indicating that the infinite-chain excitation energy limit is small due to a small BLA. For some time, the longest-wavelength excitations of cyanines have been considered as a problem too difficult for TDDFT. Early work pointed out the difficulties of finding a rationale for this failure (e.g., the single-reference nature of TDDFT response calculations), even though complete active space (CAS) wave function calculations appeared to perform well in comparison with experiment.4 The authors of a more recent (2010) TDDFT study17 noted that for the tested cyanine series the dependence of the excitation energies on the functional (among those selected) was rather weak. TDDFT errors were still attributed to a possible multireference behavior of the cyanines. Masunov proposed a ΔSCF procedure (see Section 3.4) utilizing a density matrix obtained from TDDFT to improve the results.18 Other work deemed insufficient treatment of differential electron correlation between ground and excited states to be the main culprit.6 Meguellati et al.19 proposed a vibrational correction to improve the comparison between calculations and experimental band peak maxima. An extensive recent benchmark has deemed the active space used for the CAS calculations of ref 4 insufficient.7 In agreement with ref 17, it was noted that the inclusion of exact exchange in hybrid functionals had little effect on the excitation energies, which was taken as an indication that DFT selfinteraction errors are most likely small. A hybrid functional with range-separated exchange (RSE) was included in the calculation, albeit not one that was fully long-range corrected (LC). More recent computational work8 has indicated that TDDFT has no particularly severe problems with cyanine dyes if members of the Minnesota (M) family of functionals20 are used. Reliable gasphase vertical excitation energies for benchmark purposes were obtained from coupled-cluster (CC3) and diffusion Monte Carlo (DMC) calculations. The use of these reference energies showed that previous TDDFT results had afforded less of an error than

previously thought. Yet, the TDDFT singlet excitation energies remained above the reference values. At present is remains unclear how the description of the HL transition relates to other aspects of the functional in describing correctly the delocalized electronic structure of π chromophores. The present study has been motivated in particular by the question how the performance of TDDFT for cyanines relates to the DFT delocalization error (DE)21 (‘self-interaction errors’) and whether the DE severely affects the calculated spectroscopic properties. The spectral properties of the cyanine dyes are consistent with strongly delocalized electronic structures and small BLA. Based on the simple picture of the two resonance structures in Figure 1, one would expect that the delocalization and the BLA are interconnected.22 This leads to important questions regarding the performance of TDDFT. For instance, is there a quantifiable improvement of the Minnesota functionals over standard nonhybrid and hybrid functionals in terms of the DE that can be linked to the performance in the TDDFT excitation energy calculations? Do optimally tuned RSE functionals,23−25 which tend to afford very small DEs, perform comparably or better? A related question is the extent of BLA and its impact on the excitation energies. Recent work26 has pointed out that functionals with small DEs may not necessarily produce the most physically reasonable BLA. Further, we address the relation of the HL orbital energy gap to the calculated excitation energies in the context of optimal tuning of RSE functionals. For the CN9 system, a detailed analysis of various contributions to the TDDFT excitation energy is performed. We show that the DFT components of the exchange-correlation response kernel entering the TDDFT results, and corresponding terms in ΔSCFtype calculations, are likely much too small in magnitude. For the cyanine dyes, the resulting errors on the singlet excitation energies have the same direction as those in Hartree−Fock based response calculations, placing the TDDFT results above accurate benchmark excitation energies. Unlike systems that suffer from the TDDFT charge-transfer problem, minimization of the DE for the cyanines is not correlated with spectacular improvements. Section 2 provides computational details. In Section 3, we first analyze the impact of structural parameters on extrapolated limits of the HL excitation energy for long chains and compare with a particle in a box model. The extent of the DFT delocalization error is investigated next, followed by a comparison of calculated excitation energies. Finally, the CN9 system is analyzed in detail. Concluding remarks can be found in Section 4.

2. COMPUTATIONAL DETAILS Molecular structures for the linear cyanine series (CN5−CN11) were optimized in C2v symmetry with Gaussian 0927 using the ccpVQZ basis set.28 For the optimizations, Hartree−Fock (HF) theory, DFT, and a second-order Møller−Plesset (MP2) perturbative treatment of electron correlation were applied. ‘Optimal tuning’ of hybrid functionals with range-separated exchange (RSE) was carried out with a developer’s version of the Northwest Computational Chemistry (NWChem) package.29,30 DFT calculations with fractional electron numbers were also carried out with NWChem, using a code implemented previously by one of us.31 The RSE functionals used for this work employ a 3-parameter separation of the interelectronic distance r12 in the exchange functional.32 1 − [α + β erf(γr12)] α + β erf(γr12) 1 = + r12 r12 r12 B

(1)

dx.doi.org/10.1021/ct400649r | J. Chem. Theory Comput. XXXX, XXX, XXX−XXX

Journal of Chemical Theory and Computation

Article

Table 1. Bond Lengths and Bond Length Alterationsa, Compared to MP2, for the Cyanine Series CN5-CN11 Optimized with Various Electronic Structure Methods CN5 HF BP PBE PBE0 MP2 LC-PBE* LC-PBE0* CN7 HF BP PBE PBE0 MP2 LC-PBE* LC-PBE0* CN9 HF BP PBE PBE0 MP2 LC-PBE* LC-PBE0* CN11 HF BP PBE PBE0 MP2 LC-PBE* LC-PBE0*

N1−C1 1.303(−0.019) 1.325(0.003) 1.324(0.003) 1.312(−0.010) 1.314(−0.007) 1.308(−0.014) 1.304(−0.018) N1−C1 1.309(−0.012) 1.331(0.010) 1.330(0.009) 1.318(−0.004) 1.320(−0.002) 1.315(−0.007) 1.310(−0.012) N1−C1 1.315(−0.007) 1.335(0.014) 1.334(0.013) 1.322(0.001) 1.324(0.003) 1.321(−0.001) 1.315(−0.006) N1−C1 1.320(−0.002) 1.339(0.017) 1.338(0.016) 1.326(0.004) 1.328(0.006) 1.325(0.003) 1.319(−0.003)

C1−C2 1.379(−0.006) 1.390(0.005) 1.390(0.005) 1.381(−0.004) 1.382(−0.002) 1.377(−0.008) 1.375(−0.010) C1−C2 1.375(−0.010) 1.389(0.004) 1.388(0.003) 1.379(−0.005) 1.382(−0.003) 1.376(−0.008) 1.373(−0.012) C1−C2 1.369(−0.016) 1.386(0.001) 1.386(0.001) 1.376(−0.009) 1.379(−0.006) 1.374(−0.011) 1.370(−0.015) C1−C2 1.364(−0.021) 1.384(−0.001) 1.383(−0.001) 1.373(−0.011) 1.376(−0.008) 1.372(−0.013) 1.368(−0.017)

C2−C3 1.385(0.000) 1.395(0.010) 1.395(0.010) 1.386(0.001) 1.387(0.002) 1.383(−0.002) 1.380(−0.005) C2−C3 1.392(0.007) 1.399(0.014) 1.398(0.013) 1.390(0.005) 1.390(0.005) 1.388(0.004) 1.385(0.000) C2−C3 1.398(0.013) 1.401(0.017) 1.401(0.016) 1.393(0.009) 1.393(0.008) 1.393(0.008) 1.389(0.004)

C3−C4 1.380(−0.005) 1.394(0.009) 1.394(0.009) 1.384(−0.001) 1.386(0.001) 1.382(−0.003) 1.378(−0.007) C3−C4 1.374(−0.011) 1.392(0.007) 1.391(0.007) 1.381(−0.004) 1.384(−0.001) 1.379(−0.006) 1.375(−0.010)

C4−C5 1.387(0.002) 1.397(0.012) 1.397(0.012) 1.388(0.003) 1.389(0.004) 1.387(0.002) 1.383(0.002)

a

All distances in Å. Symmetry equivalent bonds are excluded for brevity. In parentheses: The difference of the bond length minus the average C−N or C−C MP2 value, respectively. A positive value in parentheses indicates a bond that is longer than the MP2 average. BLA is indicated by how much the values in parentheses differ for the bonds in the backbone, for a given molecule and a given theory level.

references, the procedure produces HOMO−LUMO energy gaps that are optimally close to the fundamental gap of the Nelectron system, IP − EA, (EA is the electron affinity), and it tends to produce system-dependent functionals that afford small DFT delocalization errors (see Section 3.2). The optimal tuning calculations and all fractional-electron DFT calculations utilized a triple-ζ valence polarized (TZVP) Gaussian-type basis.39 It is well-known that the use of these system-specific functionals is not suitable to determine energy differences between different systems.35,40,41 The main advantage, presently, is for response calculations, for applications where the HOMO−LUMO gap is supposed to represent IP − EA, and to help diagnose DFT related problems. Time-Dependent Density Functional Theory (TDDFT) based calculations of excitation energies were performed with NWChem, and for a selection of Minnesota functionals (M06HF, M06, and M06-2X) with QChem 4.1.0.42 QChem was was also used to compute configuration interaction singles (CIS) and CIS with perturbative treatment of doubles (CIS(D)) calculations.43 In all cases, the TZVP basis set was used. Basis set convergence for the first excitation was examined based on additional calculations with the aug-cc-pVTZ basis.28 See the Supporting Information (SI, Table S1).

The second term on the right-hand side is used for the long-range exact-exchange component of the functional. The range separation parameter γ is the inverse of a distance at which the functional switches from predominantly DFT-like to predominantly HF-like. We use α + β = 1, providing a full long-range correction (LC). Two hybrid variants of the Perdew-BurkeErnzerhof exchange correlation (XC) functional33 were used: LC-PBE with α = 0, and LC-PBE0 with α = 0.25, both with values of γ = 0.3 atomic units (au). The optimally tuned versions are labeled as LC-PBE* and LC-PBE0*. For the tuning, γ was varied in order to minimize 1

J2 =

∑ [εH(N + i) + IP(N + i)]2 i=0

(2)

Here, N is the electron number of the system, εH is the energy of the HOMO, and IP is the ionization potential. The optimal γ were determined by cubic spline interpolation, with an estimated uncertainty of 0.003. The procedure and a variety of sample applications are discussed in detail, for example, in refs 13, 23−25, 34−37. The purpose is to establish that the HOMO energy is as close as possible to the negative IP, a condition known to hold for exact Kohn−Sham (KS) and generalized Kohn−Sham (GKS) theory,23,38 both for the N-electron species and the system with an extra electron. As discussed in the cited C

dx.doi.org/10.1021/ct400649r | J. Chem. Theory Comput. XXXX, XXX, XXX−XXX

Journal of Chemical Theory and Computation

Article

Table 2. Excitation Energies ΔEn for CNn Cyanines with n = 5, 7, 9, 11, Calculated with LC-PBE/TZVP, for the Geometries of Table 1a geometry

ΔE5b

ΔE7

ΔE9

ΔE11

ac

bc

Bd

HF BP PBE PBE0 LC-PBE* LC-PBE0* MP2

5.31 5.16 5.16 5.26 5.30 5.33 5.26

4.13 4.01 4.01 4.09 4.12 4.15 4.09

3.43 3.32 3.32 3.39 3.42 3.44 3.39

2.96 2.85 2.86 2.92 2.94 2.96 2.92

2.341(20) 2.322(26) 2.321(26) 2.310(26) 2.300(21) 2.291(21) 2.288(25)

5.177(4) 5.223(5) 5.219(5) 5.182(5) 5.161(4) 5.152(4) 5.182(5)

0.653(27) 0.534(36) 0.533(36) 0.567(37) 0.570(30) 0.576(30) 0.538(36)

a

Fit parameters for a PB model with bond-length alternation (BLA). bLowest-energy singlet excitation energies, in eV. cÅ units. deV units. Fitting parameter standard deviations in the last or last two decimals are given in parentheses.

long bond distances, these functionals end up overestimating the chain lengths. A set of calculations of the longest-wavelength (lowest-energy, ‘first’) singlet excitation was performed with the various optimized structures, using the same electronic structure method in each case. For this task, we chose the RSE functional LC-PBE in a standard parametrization, as it is one of the better performers in terms of DE (see SI Figures S1 and S2), among the functionals tested. The results are collected in Table 2. The same table also lists parameters obtained from a fit of the excitation energies for the CNn series, for a given set of structure optimizations, to a PB model with a BLA perturbation. The fit is (atomic units)

For the CN9 analysis, TDDFT and DFT calculations were performed with a locally modified version of the Amsterdam Density Functional (ADF)44 package, using the PBE functional and a TZ2P Slater-type orbital basis. Various integrals for the HOMO and LUMO orbitals were calculated with an ADF addon previously developed for a study of TDDFT problems for 3d metal ligand field transitions, ref 11. In this code, integrals involving the XC response kernel were calculated with the Vosko−Wilk−Nusair (VWN) local density approximation.45

3. RESULTS AND DISCUSSION This section is divided into three parts: The first part is concerned with bond length alternation and the dependence of the first excitation energy on the chain length and chain structure. The second part considers tuning of the range-separation parameter and delocalization error associated with various electronic structure methods. The third part examines the excitation energy for one of the cyanines in more detail in order to understand the impact of various approximations better. 3.1. Cyanine Geometries, Bond Length Alternation. As pointed out in the Introduction, prior work on cyanine dyes, as well as a PB model with perturbing potentials to model the structure of the chain, points at an influence of BLA on the excitation energies and the infinite-chain limit (CN∞). In turn, a small BLA can be associated with a particularly strong delocalization of the π system along the cyanine backbone. As the DFT delocalization error (DE) is of some interest in this work, it is worthwhile to study the BLA in some detail and relate it to the predictions of a suitable PB model. It is noted up-front (for numerical data see Section 3.2) that ‘pure’ (nonhybrid) density functionals delocalize too much, HF theory delocalizes too little, and optimally tuned RSE functionals typically afford small DEs (as far as the frontier orbitals are concerned). Table 1 collects structural and BLA data for the set of cyanines CN5 to CN11. Based on previous work on cyanines,1 the MP2 structures are considered to be sufficiently accurate (bond lengths agree to within about 1 pm with CCSD(T)/6-31G(d)46) and used as references for the present work. As one might expect from the delocalization behavior of pure DFT vs HF, the pure functionals BP and PBE produce small BLA in the optimized structures, whereas HF produces large BLA. Further, we see the well-known performance of GGA functionals reflected in the data, viz. a global slight overestimation of bond distances, while HF produces an expected underestimation. The hybrid functionals yield some improvements over HF, but the BLA is generally larger than desired. In terms of BLA, BP and PBE perform best in comparison to MP2. While these functionals are shown below as delocalizing too much, the effect on the cyanine structural parameters is desired. However, due to the overall too

ΔEk =

π 2(2k + 5) +B 2(ak + b)2

(3)

with k = (n − 3)/2 and n equal to 5, 7, 9, 11. The parameter a accounts for the average linear extension of the repeat unit in the backbone, b accounts for the spatial extension of the end groups, and B represents the extrapolated infinite-chain limit due to BLA (see ref 15 for details). The key point is that the fit to the model correctly predicts that with increasing BLA a chain-length independent additional energy gap opens up. The trends for B are perfectly in line with the structural data, but overall minor. MP2, BP, and PBE, which produce the smallest structural BLA, give the smallest B-parameters in the fits. The fit for the HF geometry produces a 20% larger B, and the hybrid functionals are in between. The other two parameters are very robust, as they should be, varying only by 2% and 1% for a and b, respectively. Notwithstanding the limitations of 3parameter fits to four data points, these results show that the quantum chemical calculations of the cyanine excitation energies indeed exhibit trends that can be rationalized with a simple perturbation theory model based on a one-dimensional PB. 3.2. Delocalization Errors and Optimal Tuning. The MP2/cc-pVQZ optimized structures were used for studies of the DFT delocalization error (DE) and an ‘optimal tuning’ of the range-separation parameter γ in the RSE hybrid functionals. Similar calculations performed with structures optimized with the BP functions produced very similar results regarding the optimal values for γ and numerical data related to the DE and are therefore not discussed further. For brevity, detailed results are presented for the CN9 molecule, as a representative of the whole CN set. While numerical values for the other molecules differ, the overall trends for the set are similar. A full account of the calculations can be found in the SI. Numerical data related to the DE can be obtained by considering the curvature of the energy E(N) as a function of fractional electron number N. As explained elsewhere,21 the correct behavior of E(N) would afford straight-line segments D

dx.doi.org/10.1021/ct400649r | J. Chem. Theory Comput. XXXX, XXX, XXX−XXX

Journal of Chemical Theory and Computation

Article

better, as indicated by the smaller curvatures. Selected Minnesota functionals for which we were able to perform fractional electron number calculations are seen to improve over PBE and PBE0 as far as the DE is concerned, but still afford significant E(N) curvature. M06-HF exhibits negative curvature consistently, roughly of the same order as LC-PBE0. An optimal tuning of the RSE functionals was performed as detailed in Section 2. The dependence of J2 of eq 2 on the rangeseparation parameter γ of eq 1 is shown in Figure 3 for the CN9

between integers, with slopes corresponding to the negative IP of the species with next higher integer N. Curvature indicates presence of the DE. As discussed in ref 25, there is a close relation of the lack (or presence) of curvature in E(N) to the presence (or lack) of an integer discontinuity of the KS potential. Approximate local functionals such as PBE, BP, have been developed to perform well at integer N, but in order to accomplish this, curvature must be present.25 Figure 2 shows a plot of E(N) for the cyanine CN9, with N around ±1 of the actual electron number of the system. The

Figure 3. Plots of J2 of eq 2 versus the range-separation parameter γ, for CN9, with LC-PBE and LC-PBE0 RSE hybrid functionals. The optimal value for γ, minimizing J2, is given in figure key. Optimal γ for the other cyanines were obtained from similar plots (not shown).

Figure 2. Total energy of CN9 as a function of fractional electron number, ΔN, relative to the actual electron number of the system (ΔN = 0). The numerical values in the plot key correspond to the (ΔN)2 coefficient of quadratic fits of E(ΔN) in the electron-deficient and electron-rich segments, ΔN < 0 and ΔN > 0, respectively. For similar plots for the other cyanines see the Supporting Information.

system. The optimal values, minimizing J2, for the cyanine series are shown in Figure 4 (numerical data are in SI Figures S1 and S2). The corresponding functional parametrizations are indicated by an asterisk (e.g. LC-PBE*). In the delocalized cyanines, DFT exchange is seen to be important because the optimal γ parameters for the longer cyanines are significantly below values of 0.3 to 0.4, which are used in commonly applied global functional parametrizations. A γ of 0.3 atomic units translates to a switching distance of about 1.8 Å. As the chain length increases, the optimal γ decreases, which is in accord with

curvature values for the electron-deficient and electron-rich segments of the plots, extracted from quadratic fits of E(N), are listed separately. Table 3 provides similar curvature values for the full set of cyanines. The results are typical: the positive (and large) curvature for BP, PBE indicated too much delocalization, and the negative curvature of HF indicates quantifies that there is too little delocalization. The various hybrid functionals perform

Table 3. Delocalization Error and E(N) Curvature: Coefficients of (ΔN)2 of Quadratic Fits of E(ΔN)/eV for ΔN < 0 and ΔN > 0, Respectively, Extracted from Data Such as Shown in Figure 2 CN5

CN7

CN9

CN11

ΔN













PBE PBE0 LC-PBE LC-PBE0 LC-PBE* LC-PBE0* HF M06-HF M0 M06-2X M08-HX M08-SO

3.06 2.01 0.10 −0.27 0.01 0.00 −1.68 −0.46 1.97 1.05 1.15 0.96

2.99 1.99 0.06 −0.26 −0.03 0.00 −1.47 −0.35 1.96 1.12 1.19 0.98

2.59 1.69 −0.11 −0.42 0.00 −0.01 −1.62 −0.48 1.67 0.88 0.95 0.80

2.54 1.68 −0.12 −0.40 −0.01 −0.01 −1.59 −0.39 1.64 0.92 0.99 0.80

2.26 1.47 −0.25 −0.53 0.00 −0.01 −1.66 −0.53 1.44 0.73 0.80 0.66

2.23 1.45 −0.25 −0.52 0.00 0.00 −1.72 −0.44 1.41 0.78 0.83 0.67

2.02 1.30 −0.36 −0.63 −0.01 −0.02 −1.69 −0.59 1.28 0.61 0.68 0.54

2.00 1.28 −0.35 −0.62 −0.01 −0.02 −1.81 −0.48 1.24 0.66 0.71 0.56

E

dx.doi.org/10.1021/ct400649r | J. Chem. Theory Comput. XXXX, XXX, XXX−XXX

Journal of Chemical Theory and Computation

Article

previous optimal-tuning studies on conjugated π systems.24,36,47−50 Since γ is the inverse of a distance where the range-separation takes action, smaller γ means that DFT exchange acts over larger interelectronic distances. As shown previously,36 this goes along with a stronger (physically meaningful) delocalization of the π-system. Regarding the DE furnished by the optimally tuned functionals, the cyanine series is also no apparent exception: The DEs are the smallest among all functionals, and reasonably close to zero such that further explicit curvature minimization,24,35 by simultaneously adjusting γ and α = 1 − β of eq 1, appears unnecessary. The stock parametrization of LC-PBE also gives rather small DEs for the shorter cyanines, for the reason that γ obtained from the optimal tuning is close to 0.3. 3.3. Longest-Wavelength Excitations. Table 4 collects the results from calculations of the lowest-energy (longest-wavelength) singlet electronic excitations of the CNn series, along with a selection of data from the literature. For consistency with

Figure 4. Optimally tuned γ* as a function of cyanine dye length for MP2 and BP geometries.

Table 4. Lowest-Energy (Longest-Wavelength) Singlet Excitation Energies ΔEn (eV) for the CNn Cyanines of Figure 1, Calculated with Different Methodsa excitation energies and HOMO−LUMO gaps ΔE5

b

HF(RPA) HF(TDA) CIS PBE PBE(TDA) BP PBE0 M06 M06-HF M06-2X LC-PBE LC-PBE* LC-PBE0 LC-PBE0* CIS(D) cPBE+HFh

5.79 6.08 6.08 5.30 6.00 5.29 5.40 5.31 5.15 5.27 5.26 5.26 5.38 5.37 4.85 5.81

PBEe PBE0e M05e M05-2Xe M06-Le M06-HFe M06e M06-2Xe M08-HXe M08-SOe B2PLYPe,f CC2g CC3f,g CASPT2f,g DMCf,g

5.23 5.34 5.30 5.33 5.40 5.13 5.23 5.23 5.23 5.16 5.05 4.97 4.84 4.69 5.03

Δε5c

ΔE7

11.73

4.51 9.64 4.77 4.77 4.15 2.54 4.92 2.54 4.15 2.55 4.22 4.14 4.17 4.24 3.99 8.52 4.11 5.81 4.09 7.63 4.09 7.45 4.18 8.09 4.19 7.41 3.65 9.64 4.52 9.64 Previously Calculated Datad 4.12 4.19 4.15 4.17 4.25 3.98 4.10 4.09 4.09 4.04 3.92 3.79 3.65 3.53 3.83

3.56 3.56 3.56 5.48 5.57 10.48 7.40 9.28 9.44 9.86 9.39 11.73 11.75

Δε7

ΔE9

Δε9

ΔE11

Δε11

3.75 3.98 3.98 3.47 4.27 3.47 3.53 3.46 3.29 3.42 3.39 3.40 3.47 3.49 2.95 3.74

8.39

3.23 3.45 3.45 3.01 3.56 3.01 3.05 3.00 2.82 2.95 2.92 2.94 2.98 3.02 2.48 3.22

7.55

3.45 3.50 3.47 3.47 3.56 3.30 3.43 3.41 3.41 3.37 3.25 3.10 2.96 2.81 3.09

1.98 1.98 1.98 3.37 3.47 7.33 4.86 6.61 6.20 7.00 6.16 8.39 8.37

1.62 1.62 1.63 2.86 2.96 6.54 4.23 5.91 5.38 6.27 5.33 7.55 7.52

3.00 3.04 3.00 3.00 3.11 2.83 2.98 2.95 2.95 2.91 2.80 2.64 2.53 2.46 2.62

a TZVP basis, MP2/cc-pVQZ structures. HF(RPA) = time-dependent HF linear response calculation. TDA = Tamm−Dancoff approximation. All TDDFT calculations use the full (RPA-like) formalism except when TDA is indicated. CC2 and CC3 are approximate coupled-cluster models. DMC − = diffusion Monte Carlo. B2PLYP is a double hybrid functional. HOMO−LUMO gaps Δεn (eV) bFirst excitation energies in eV. cΔεn = εLUMO n (eV) dReference 8. eGround state PBE0/cc-pVQZ structures from corresponding citation. fReference 7. gMP2/cc-pVQZ structures from εHOMO n corresponding citation. Basis sets for first excitations in ref 8 are ANO-L-VQZP for n = 1 − 3 and ANO-L-VTZP for n = 4. In Citation 7 the ANO-LVTZP basis set is used to yield first excitations. hPBE correlation functional with 100% exact exchange, discussed in Section 3.4.

F

dx.doi.org/10.1021/ct400649r | J. Chem. Theory Comput. XXXX, XXX, XXX−XXX

Journal of Chemical Theory and Computation

Article

estimation of the CT excitation energy, and negative curvature (as typical with HF theory) goes along with an overestimation of the CT energy. An underestimation of the energies of ‘CT-like’ π → π* transitions with standard functionals that typically afford a large DE has also been noted for conjugated π systems.13 For the cyanines, this behavior is evidently not produced. The functionals affording the largest DE for the cyanine series (BP, PBE) produce singlet excitation energies that are very close to the optimally tuned RSE functionals, which afford the smallest DEs. Among the latter, LC-PBE* performs somewhat better than LCPBE0*. Both give similar excitation energies as the nontuned versions. CIS, and HF(RPA), with large negative E(N) curvature in the HF ground state, produce the highest excitation energies. However, M06-HF gives the lowest energies despite this functional also producing sizable negative curvature (albeit not as large as HF). As discussed in Section 3.1, the CN∞ limit of the excitation energy is weakly dependent on the extent to which the DE affects the optimized cyanine structures. The excitation energies also vary with the extent of the DE, but there is no obvious simple correlation. In Section 3.4, we provide expressions for the excitation energies from two-level models, where the orbital energy gap is the leading term. The optimally tuned RSE functionals produce, by construction, HOMO−LUMO energy gaps Δε that are very close to the fundamental gap, IP-EA (‘gap’ being used here as a molecular property not its band-structure analog for infiniteperiodic systems). The excitation energy is the ‘optical gap’, ΔE. Conceptually, one may consider the optical gap as being related to creating a noninteracting electron−hole (e-h) pair, with an associated energy of IP-EA, plus a negative binding energy of the e-h pair in the molecule. The tuned RSE functionals accordingly have Δε significantly larger than ΔE, by roughly a factor of 2. The linear response hybrid XC kernel in the excitation energy calculation needs to act in order to overcome this energy difference, but evidently it falls somewhat short. With HF, the orbital gap is even larger in comparison, for CN11 more than twice the calculated ΔE and almost three times the DMC reference value. The difference between Δε and ΔE is seen to be larger for HF than it is for the optimally tuned RSE functionals. For the pure functionals BP and PBE, the situation is reverse: Δε < ΔE, which is not untypical. The XC response in a GGA functional such as PBE and BP therefore must result in a positive adjustment of the leading Δε term toward the excitation energy and, apparently, overcorrects. The exact form of the GGA functional matters little; the numerical values for the orbital energies and the excitation energies are within 0.01 eV for PBE and BP. Due to the self-consistent nature of the ground state and the response calculations, it is not straightforward to make simple connections between HF, the RSE functionals, and the underlying GGA functions (PBE in this case). In most cases, the untuned versions of the RSE functionals have larger, in some cases significantly larger, Δε than the tuned versions, which goes along with more of a HF component overall due to the (in most cases) larger γ in the untuned versions. Upon optimal tuning, the decrease (in most cases) of γ means the functionals have DFT exchange active over longer ranges, which causes the orbital energy gaps to drop, but the excitation energies change very little. One may attribute this to the aforementioned positive ‘pure DFT’ XC response to the excitation energy, which becomes more important as γ decreases. The orbital energy gaps for M06-HF are in between the tuned RSE functionals and HF, due to this being a hybrid functional

the optimal tuning, with the fractional electron number calculations, and with the analysis for CN9 in the following subsection, the excitation energies were calculated with the TZVP basis. Additional excitation energies calculated with the larger aug-cc-pVTZ basis set can be found in SI Table S1. The aug-cc-pVTZ excitation energies are between 0.15 (CN5) and 0.08 eV (CN11) eV lower in energy. Comparing with Table 2, these changes are of similar magnitude as the variations of the LC-PBE excitation energies for the different optimized structures. We take the DMC and CC3 results to be close to exact complete-basis limit vertical excitation energies for the MP2-cc-pVQZ structures. Linear response (LR) calculations using HF theory and the Tamm-Dancoff approximation51 (TDA) are formally and numerically equivalent to configuration-interaction singles (CIS). HF(RPA) is equivalent to TDDFT-based LR calculations of excitation energies if in TDDFT the exact exchange functional for the KS single reference determinant is used and no correlation functional is present. The CIS excitation energies are much too high. HF(RPA) improves somewhat over CIS; that is, the TDA introduces unwelcomed errors that, for CN5, amount to as much as 0.3 eV. For the PBE-based TDDFT calculation the TDA produces even larger blue-shifts. The poor performance of TDA for the cyanine excitations has been noted before;6,7 the CNn series is one of the rare cases where deexcitation parameters for the linear response of the density matrix appear to be highly important. Other situations where the TDA may give significantly higher, but more correct, excitation energies is the case of triplet excitations where near triplet instabilities are present.52 Indeed, there are increasingly pronounced differences between RPA and TDA for the calculated HF triplet excitation energies with increasing chain length (Table S4), and an actual triplet instability is present for CN11 with HF(RPA) (Table S5). It is unclear whether nearinstabilities have any bearing on the smaller RPA-TDA difference in the PBE singlet excitations for CN11 as compared to the shorter chains. The corresponding RPA-TDA triplet energy differences are quite small with PBE for all systems (SI, Table S4). CIS(D) includes electron correlation from doubles substitutions perturbatively on top of a CIS calculation. The CIS(D) results vastly improve over CIS, HF(RPA), and most TDDFT calculations. Evidently, electron correlation is of paramount importance for the cyanine excitations. When considering the effects from differences in the basis set, that is, TZVP giving slightly too high excitation energies, CIS(D) is seen to overcorrect CIS. The improvement with CIS(D) for cyanine excitations has been noted previously by Grimme and Neese.6 A ‘double hybrid’ functional that included some of the same wave function-based correlation correcctions was accordingly found to improve the results, but it did not fully resolve all of the problems.6 The excitation energies calculated with all types of density functionals (i.e., not including HF or the double hybrid) differ by, at most, ∼0.3 eV for the set of cyanines. As the chain length gets longer the relative error (compared to DMC/CC3) in the excitation energy increases. In our data set, the highest excitation energies are produced by the PBE0 functional and the lowest energies are furnished by M06-HF. For systems with low-energy charge-transfer (CT) excitations, the performance of the various KS response methods is closely linked to the DE (for an example, see ref 36): A large DE, with positive curvature of E(N) goes along with a strong underG

dx.doi.org/10.1021/ct400649r | J. Chem. Theory Comput. XXXX, XXX, XXX−XXX

Journal of Chemical Theory and Computation

Article

with 100% global exact exchange. The trend in the orbital gaps correlates with the extent of negative curvature calculated for E(N) if one compares the various functionals. As can be seen from Table 4, M06 itself gives Δε ∼ ΔE (with Δε slightly higher for CN5 and CN7). The 100% exact exchange in the response is evidently providing a desired overcompensation for the increase in the orbital gap, relative to M06. The decreased error of overlocalization/negative curvature relative to HF, while keeping full exchange in the functional, works in favor of low excitation energies and better agreement with the DMC reference values. 3.4. Detailed Analysis for CN9. A detailed analysis has been performed for the CN9 molecule as a representative of the cyanine set, by following a similar protocol as used in ref 11. In particular, we were interested in how the TDDFT calculations compare with DFT-based ΔSCF calculations, where determinants for the ground state and excited electron configurations are optimized separately and the excitation energy is explicitly calculated from a difference of the total energies.53 Some authors have noted previously11,54−56 that DFT-based ΔSCF with approximate functionals sometimes yield rather accurate excitation energies in situations where linear-response (LR) TDDFT with XC response kernels derived from the same functionals produce large errors. Ziegler and co-workers have been motivated by this finding to justify a variational alternative for LR TDDFT for calculations of excitation energies,57−59 which has shown promising results, for instance for cases plagued by the CT problem.60 The term ΔD is used in the present work to indicate a ΔSCF-like procedure where the orbitals for the excited determinant have not been optimized separately but are the same as for the ground state. In other words, ΔD corresponds to ΔSCF without relaxation of the excited KS determinant. As prescribed by Ziegler et al.53 and others, the singlet excitation energy is obtained from 2E[S*] − E[T*] − E[GS], where the energies are calculated for excited single-determinant configurations representing the singlet ground state GS, and excited determinants without [S*(H → L)] and with spin-flip [T*(H → L)], in order to compensate for the fact that the open-shell singlet excited state cannot be represented by a single determinant. Following the analyses of refs 56 and 11, expressions for the LR TDDFT, LR TDHF, and for ΔD singlet (S) and triplet (T) excitation energies are taken from a two-level model for a transition between the HOMO φH of a closed-shell ground state and the LUMO φL in the Tamm−Dancoff approximation as

ESΔD = (εL − εH) +

αα αβ − [LL|f XC − f XC |HH]

E TΔD = (εL − εH) +

αα αβ E TTDDFT = (εL − εH) + [LH|f XC − f XC |LH]

(6b)

In the above expressions, the leading-order term in each case is the difference between the orbital energies, εL − εH. It is important to keep in mind that orbital gap in the HF calculations has a different meaning than in DFT and is much larger. Further, [pq|r−1 12 |rs] is a two-electron repulsion integral (ERI) involving the MOs. There are corresponding integrals in DFT involving the XC linear response kernel f στ XC with σ,τ referring to orbital spin labels α or β. We assume a local adiabatic XC response kernel, in which case the integrals are over the coordinates of one electron and independent of the excitation energy. The quantity Δρ = |φ2L| − |φ2H| is the ΔD density change upon excitation. The differences between singlet and triplet excitations is for TDHF given by 2[LH|r−1 12 |LH], a positive exchange integral, leading to ES > ET. This term is also present in the TDDFT expressions. The αα f στ XC integrals tend to be negative, with contributions from f XC > f αβ . Therefore, the linearized ΔD expressions should also XC produce ES > ET. Casida et al.56 pointed out that the second term on the righthand side of eq 6, which was referred to as a ‘charge transfer’ term (but explicitly cautioning that actual CT may not be involved), is an artifact of approximations in DFT. With the exact XC functional and the exact response kernel, the expressions for the linearized TDDFT and linearized ΔD two-level excitation energy αβ should be the same. For a semilocal functional, [LL|f αα XC − f XC| αα αβ HH] = [LH|f XC − f XC|LH]. With Δρ ≈ 0, representing an evident absence of CT, the linearized ΔD and TDDFT triplet excitation energies would then be the same. One may consider the magnitude of the ‘CT term’, ΔCT =

1 −1 αα [Δρ|r12 + f XC |Δρ] 2

(7)

relative to the excitation energy11 as a measure of the approximate character of the functional. Further, we can compare the magnitudes of the Coulomb and exchange integrals in the TDDFT versus the HF expressions to gauge the potential impact of correlation effects. Data for the model calculations for CN9 are collected in Table 5. For the XC response kernel, the model calculations utilize the VWN functional. Table 5 therefore also compares the singlet excitation energies obtained with NWChem and the full PBE response kernel, and ADF using PBE for the ground state and the VWN kernel, with basis sets of comparable flexibility. The difference is only 0.1 eV. For HF-based calculations, ADF employs an auxiliary basis set, which leads to slightly larger deviations between the excitation energies and orbital gaps obtained with the two codes, but the analysis data are sufficiently accurate to draw meaningful conclusions. We note that the ERIs listed in Table 5 are calculated numerically for a given orbital pair and do not use an auxiliary basis approximation. Figure 5 displays isosurfaces of the relevant orbitals (PBE), the HOMO−LUMO product, the ΔD density change Δρ, and the density difference between the excited singlet KS determinant and the ground state. Contributions from orbitals other than the HL pair in the excitation vectors, that is, expanding from the 2-level models to

(4a) (4b)

Adopting the TDA is equivalent to a linearization of the two-level model expressions. In time-dependent HF, the linearized eigenvalues are equal to the corresponding ΔD expressions and given by −1 −1 ESTDHF = (εL − εH) − [LL|r12 |HH] + 2[LH|r12 |LH]

(5a) −1 E TTDHF = (εL − εH) − [LL|r12 |HH]

(6a)

1 −1 αα [Δρ|r12 + f XC |Δρ] 2

αα αβ + [LL|f XC − f XC |HH]

αα αβ ESTDDFT = (εL − εH) + [LH|f XC + f XC |LH] −1 + 2[LH|r12 |LH]

1 −1 αα [Δρ|r12 + f XC |Δρ] 2

(5b)

DFT-based ΔD singlet and triplet energy expressions are obtained as in eq 6 after linearizing XC terms that depend on the density change: H

dx.doi.org/10.1021/ct400649r | J. Chem. Theory Comput. XXXX, XXX, XXX−XXX

Journal of Chemical Theory and Computation

Article

Table 5. Numerical Analysis Data for the Lowest-Energy Excitation of CN9, from ADF Calculations (PBE Orbitals, Selected Results from HF Orbitals)a GS orbitals LR (NW) LR (ADF) LR/TDA (NW) Δε (NW) Δε (ADF) EΔSCF S EΔD S Lin. EΔD S Lin. ETDDFT S Lin. ETDHF S αα (1/2)[Δρ|r−1 12 + f XC|Δρ] (1/2)[Δρ|f αα |Δρ] XC αβ [LH|f αα XC + f XC|LH] αα [LH|f XC − f αβ XC|LH] [LL|r−1 12 |HH] [LH|r−1 12 |LH] [LL|r−1 12 |LL] [HH|r−1 12 |HH] Lin. ΔD S/T Lin. TDDFT S/T Lin. TDHF S/T ⟨|φH|||φL|⟩

PBE 3.47 3.48 4.27 1.98 1.98 2.00 2.44 2.45 5.17b

HF 3.75 3.72 3.98 8.39 8.30

4.73 4.73c

0.30 −0.23 −0.24 −0.18 5.68 1.72 6.27 6.14 0.35 3.37 0.74

5.73 1.08 6.37 6.42

Figure 5. CN9 HOMO and LUMO, and the density change Δρ for the longest-wavelength excitation (HOMO→LUMO) from ΔSCF (PBE calculations). The ΔD density change (LUMO2 − HOMO2) and the HOMO−LUMO product are also shown. Isosurface values (±) as indicated.

2.17 0.68

a

All energies in eV. The DMC reference value for CN9 is 3.09 eV. Full linear response (LR) excitation energies and orbital gaps included to compare NWChem (NW) and ADF results. MP2 structure. In the ADF calculations, f XC corresponding to the VWN functional was used. ΔD = singlet excitation energy obtained from KS determinant energy differences with ground state orbitals as 2E[S*] − E[T*] − E[GS] (see text). The straight energy difference between excited ‘singlet’ and ground state determinants, E(S*) − E(GS), is 2.22 eV. ΔSCF = same as ΔD but with relaxation of the excited KS determinants. E(S*) − E(GS) = 1.86 eV. Δε = (εL − εH), Δρ = φ2L − φ2H. S/T = singlet− triplet excited state energy difference. b5.15 eV from a PBE(TDA) calculation with NWChem with all orbitals but HOMO and LUMO frozen, forcing the linearized 2-level model. c4.78 eV from a HF(TDA) is equal to Lin. calculation 2-level calculation with NW. Lin. ETDHF S EΔD S , eq 5.

99% for BP and PBE, to 96 to 97% for the tuned RSE functionals, to 94 to 93% for HF(RPA) and HF(TDA). A similar range of HL percentages in the transition density is found for the corresponding triplet excitations. Therefore, even with HF the excitation is described as a practically pure HL transition. Compared to analyses of excitations that we have performed in the past, for the cyanines the HL contribution in the HF calculations is rather large. However, energetically, the mixing with other orbital pairs is important. As noted, going beyond the TDA further lowers the excitation energy significantly, in particular for the DFT-based calculation. The difference between (linearized = TDA-based) TDDFT S vs T excitations for the two-level model is predominantly the term 2[LH|r−1 12 |LH]. The actual difference amounts to 3.37 eV, including a 3.44 eV contribution from twice the exchange ERI. In TDHF, two times the exchange ERI is exactly the S/T excitation energy difference, more than an eV smaller than with TDDFT because of the smaller exchange integral. It is worthwhile pointing out that the linearized TDDFT S/T excitation energy difference is larger than the actual singlet excitation energy from DMC, and comparable to the PBE TDDFT excitation energy. The full TDDFT excitation energy is also affected by this large positive contribution. The actual S/T excitation energy difference for CN9 calculated with NWChem is 1.79 eV (PBE), highlighting again the effects from orbital mixing in the transition vectors despite the fact that the transition is 99% HL. The DFT-based ΔSCF, ΔD, and linearized ΔD singlet excitation energies are all within a range of 0.5 eV from each other, within 2.00 and 2.45 eV, close to the HOMO−LUMO gap Δε (1.99 eV) and much below the DMC target. The full and linearized forms of ΔD give very similar results around 2.45 eV. Relaxation of the excited determinant in the ΔSCF procedure brings the excitation energy back close to the value for Δε. The qualitative appearance of the density change going from the GS to S* is shown in Figure 5 and seen to be quite similar to the ΔD

the full space of orbitals to describe the transition density, accounts for a significant lowering of the excitation energies. The linearized (TDA) TDDFT singlet energy from the 2-level model is 5.17 eV. A corresponding linearized HF value is 6.15 if one uses the HF HL gap, but otherwise the electron repulsion integrals calculated with PBE orbitals. With HF orbitals, the linearized HF excitation energy drops to 4.73 eV, in large parts because of the significantly smaller [LH|r−1 12 |LH] exchange integral. As explained in footnotes b and d of Table 5, TDA calculations with NWChem were performed with all orbitals but HOMO and LUMO ‘frozen’ in order to enforce the same linearized two-level models (using the full PBE response kernel for the DFT calculation), giving excellent agreement with the data assembled from the integrals in Table 5. This is another numerical test showing that the integrals used in the analysis extracted from ADF are sufficiently accurate. The linearized singlet excitations energies for PBE and HF are 0.90 and 0.75 eV above the full TDDFT and HF results (TDA), respectively. Relaxation of the excitation vector to include MO pairs other than HOMO−LUMO has therefore a significant effect both for HF and TDDFT. The contribution of the HOMO−LUMO orbital pair to the transition vector ranges from I

dx.doi.org/10.1021/ct400649r | J. Chem. Theory Comput. XXXX, XXX, XXX−XXX

Journal of Chemical Theory and Computation

Article

density change |φ2L|−|φ2H|. However, the effects from relaxation are readily apparent. As with the TDDFT and TDHF procedures discussed, a small amount of relaxation can produce significant changes in the energy. However, the ΔSCF vs ΔD vs linearized ΔD energy change is not as large as those seen in TDDFT and TDHF when going from the two-level models to the full calculations. ΔSCF and ΔD exhibit very poor performance for the CN9 cyanine. Even with ΔD, where the excitation energy is somewhat larger, the singlet excitation energy is 0.65 eV below the DMC target, and with ΔSCF it is more than 1 eV too low. Since the linearized ΔD excitation is very close to the ΔD energy difference, an analysis based on eq 6 is appropriate. The ‘CT term’ is rather small, only 0.30 eV, with small negative contributions from f αα XC being overpowered by the aggregate positive (but also not particularly large) (1/2)[Δρ|r−1 12 |Δρ]. The remaining term in the ΔD excitation energy, apart from the αβ αα αβ orbital energy gap, is −[LL|f αα XC − f XC|HH] = −[LH|f XC − f XC| LH], which is also small, 0.18 eV. This and the CT term are much too small in magnitude to push the excitation energy much above the HOMO−LUMO energy gap of 1.98 eV, resulting in the observed severe underestimation of the ΔD excitation. We assert that the ΔSCF procedure suffers from similar problems. The energy lowering of the excited determinants from the SCF procedure aggravates the problem. It is interesting in this context to consider the energy difference between singlet and triplet excitations (S/T). As pointed out, for the two-level models of TDDFT and TDHF, this is twice as large as the exchange integral between HOMO and LUMO. The term originates in the Coulomb potential of the system, not the XC potential, which is the reason why it occurs both in the TDDFT and TDHF expressions. For the various full calculations, the S/T difference is 1.8 to 2.4 eV indicating again the importance of even small mixing of different MO pairs in the transition vector for the resulting energy (see SI Table S3; the HF(RPA) triplet excitation is deemed unreliable because the energy is very low, even though the ground state is stable). We calculated triplet excitation energies of 1.73 eV for ΔSCF, 2.01 eV for ΔD and 2.10 for the linearized ΔD version. The corresponding S/T differences are very much smaller than those obtained from the LR calculations. In the linearized ΔD version, the S/T difference is −2[LL|f αα XC − αα αβ f αβ XC|HH] = −2[LH|f XC − f XC|LH] = 0.35 eV from Table 5. Evidently, these integrals are too small to be able to produce the magnitudes of the S/T gaps obtained in the LR calculations. Based on the available data, the finding is that for the HOMO− LUMO transition in the cyanine, the integrals resulting from a στ local XC response kernel of the type [LL|f στ XC|HH] = [LH|f XC| LH] are small in magnitude, and most likely much too small. Further numerical evidence is provided by a set of calculations with the PBE correlation functional and 100% HF exchange (Table 4). The HOMO−LUMO gaps and the excitation energies are virtually identical to HF(RPA); that is, the presence of the correlation functional component matters very little. Moreover, a calculation for CN9 with the PBE functional, but with the exchange component scaled to 1/2 gives results that are within about 0.2 eV from the full PBE calculations (SI Table S3). What, then, are the consequences of too small f XC contributions for the performance of TDDFT and ΔSCF with a local functional such as BP or PBE? The TDDFT singlet excitation energy includes a large exchange ERI, 2[LH|r−1 12 |LH], αβ while negative contributions from [LH|f αα XC ± f XC|LH] are likely severely underestimated. Relaxation from mixing contributions from other orbital pairs in the TDDFT excitation vector lowers

the excitation energy relative to the 2-level model, but it does apparently not qualitatively change the balance of XC versus ERI terms. A severe underestimation of the magnitude of negative στ [LL|f στ XC|HH] = [LH|f XC|LH] terms (even with small HOMO− LUMO energy gaps to begin with) then rationalizes why TDDFT with local functionals overestimates the cyanine singlet excitations. Here, the XC contribution to the excitation energy is negative. The ΔSCF and ΔD approaches fail by producing much too low excitation energies. The leading term, the orbital energy gap, is small, and the CT term is also small (Table 5). In the two-level model, the S/T difference is given by an overall positive but small αβ −2[LH|f αα XC − f XC|LH]. As in the TDDFT case, this term is likely much too small in magnitude. Incidentally, the TDDFT and ΔSCF triplet excitation energies calculated with PBE are 1.79 and 1.73 eV, respectively, which agree rather well much unlike the singlet excitations. One of the conclusions of this section is therefore that the shortcomings of TDDFT with standard local functionals in the DFT part of the XC functional is a too small magnitude of [LH| στ f XC |LH] contributions. These shortcomings place TDDFT singlet excitation energies of the cyanines with pure functionals on the same side as HF theory, both above the correct excitation energies. In the well-known CT problem, the HOMO and LUMO are spatially distant and correspondingly, the exchange integral is small. In this case, too small contributions from f XC in conjunction with small HOMO−LUMO energy gaps place the calculated CT energy below the correct value. However, in such a case, ΔSCF may perform well in comparison because the CT term tends to become important. An optimal tuning procedure, trading approximate DFT for HF contributions in the exchange functional, does evidently not improve the cyanine excitation energies, even though the process can eliminate large delocalization errors from the calculations. As already noted, going from CIS or time-dependent HF linear response to CIS(D) gives a large improvement of the excitation energies by lowering them toward DMC and CC3 (though overcorrecting). Therefore, another way to look at the changes introduced when going from a pure to an optimally tuned RSE is that introducing the long-range HF components in the RSE functional, to reduce the DE, may simply eliminate too much electron correlation simultaneously in order to give an improvement. As evident from Table 4, global hybrids such as PBE0 fare even worse. An analysis of what exactly happens in the Minnesota functionals is beyond the scope of this paper. One may note that with M06-HF the 100% HF exchange component is going to behave similar to HF(RPA), but the response calculation starts out with significantly smaller orbital energy gaps than HF (Table 4). This simple picture predicts the desired lowered excitation energies. However, the self-consistent nature of the ground state and TDDFT calculations complicates the analysis. One may speculate that for the symmetric cyanines a lack of CT character in the longest-wavelength transitions along with the large HL exchange integral exposes approximations made in TDDFT with common functionals. With a partial CT character, the known tendency of pure and global hybrid functionals to underestimate CT energies would likely provide some error compensation, by lowering the excitation energies toward accurate reference data. The ‘CT term’ in the ΔSCF and ΔD expressions would then presumably become larger, driven by −1 larger differences between [LL|r−1 12 |LL] and [HH|r12 |HH] versus −2[LL|r−1 |HH], and give improved results. 12 J

dx.doi.org/10.1021/ct400649r | J. Chem. Theory Comput. XXXX, XXX, XXX−XXX

Journal of Chemical Theory and Computation

Article

The spatial overlap, given by ⟨|φH|||φL|⟩ for the PBE orbitals, is 0.74 (Table 5) out of a range from 0 to 1. (|···| indicates the modulus of an orbital.) The self-repulsion integrals [HH|r−1 12 | HH] and [LL|r−1 12 |LL] are nearly identical (6.14 and 6.27 eV) and altogether not that different from the Coulomb repulsion [HH| r−1 12 |LL] between the two orbitals (5.68 eV). The similarity may be taken to support a lack of CT character,61 which is also evident from the MOs and the density change shown in Figure 5. There is a large transition dipole moment in the direction of the long axis of the cyanine, but no change of the dipole moment in the same direction upon excitation. Kuritz et al. invoked an additional numerical criterion in order to detect ‘CT-like’ character for excitations of π-chromophores:13 Linear combinations of the HOMO and LUMO were formed as ϕ1/2 = (φH ± φL)/21/2, and a very small value of ⟨|ϕ1|||ϕ2|⟩ was argued to indicate a CT-like problem. For the PBE orbitals of CN9, this spatial overlap is 0.50 and not particularly small (although, as pointed out previously by Masunov,18 there is a large amount of charge shifted upon the excitation). An interesting aspect related to the spatial overlap, and noted already, is the magnitude of the [LH|r−1 12 |LH] exchange integral, which is much larger when calculated with the PBE orbitals than with the HF orbitals. The spatial overlap of HOMO and LUMO listed in Table 5 is also noticeably larger for the PBE orbitals than for HF. Given that two times the exchange integral is a leading term entering the S/T gap, there is a possibility that this contribution to the excitation energy is simply very strongly overestimated in the PBE calculations, pushing the singlet excitation energy above the correct value. Supposedly, with sufficient balance from larger f XC related terms there would not be such problem. There is some support for this viewpoint in the data collection of the SI. For the triplet excitation energy (eV), HF(TDA) = CIS = 2.03, CIS(D) = 2.17; that is, the correlation effect is quite small. PBE gives 1.68, now underestimating with respect to CIS(D), with a PBE(TDA) value of 1.72 being close. Interestingly, the ‘two-level’ calculation with NWChem gives almost the same excitation energy. The CIS(D) estimate for the S/T gap is 0.78 eV, less than half of the 1.79 eV from PBE (but still ∼0.4 eV larger than those obtained from ΔSCF and ΔD). If the magnitude of the f XC contributions were the only factor influencing the excitation energy, the PBE triplet calculation should in fact benefit from an underestimation of these terms: Δε is 1.98, not too far to the CIS(D) triplet excitation energy. Small αβ negative terms of the type [LH|f αα XC − f XC|LH] keep the PBE result close, but still too low. Therefore, another problem is simply a too small HOMO−LUMO gap itself in the PBE calculation, negatively affecting all ΔSCF results and the TDDFT triplet excitation. On the other hand, the singlet excitation energy in the PBE calculation benefit from a too small Δε because here small-magnitude f XC contributions do not counterbalance the large positive exchange ERI. The trends regarding the triplet energies and the S/T difference noted for CN9 are found as well for the other cyanines. CC2 and CIS(D) singlet and triplet excitation energies are collected in the SI (Table S6). In comparison with CIS(D), one finds that TDDFT with PBE underestimates the triplet excitation energies by 0.4 to 0.6 eV and overestimates the S/T gap by 1.5 to 1.8 eV. HF(TDA) produces triplet excitation energies that are close to the CIS(D) results, at most 0.15 eV higher, but overestimates the S/T difference by 1.1 to 1.3 eV. Lastly, we discuss the orbital energies for CN9. The frontier orbital energies in HF theory are, per Koopmans’ theorem, approximations for -IP and -EA (the latter one typically being a

poor approximation). Since we have the orbital energies from the optimally tuned RSE functionals as well as ionization and electron detachment data from the tuning procedure available, a direct comparison with theoretical values for IP and EA is possible. The LUMO lies in the valence space and it is rather compact, as indicated by the [LL|r−1 12 |LL] integrals in Table 5 that are just as large as those for the HOMO. Therefore, basis set convergence related to an electron attachment should not be a serious concern. Selected relevant energetic data are collected in Table 6. The first point to note is that the HF Koopmans’ IP Table 6. CN9: Selected HOMO and LUMO Energies (eV) Times −1, and the HOMO−LUMO Gap, Calculated with Different Methods PBE PBE0 LC-PBE* LC-PBE0* HF

−εH

−εL

εL − εH

8.61 9.36 10.94 10.85 10.97

6.63 5.99 4.74 4.69 2.59

1.98 3.37 6.20 6.16 8.39

agrees within 0.1 eV with those from the tuned RSE functionals. However, the LUMO energy, as an estimate of the EA, is over 2 eV too small in magnitude with HF. (The EA is large because CN9 is treated as an isolated cation in gas phase.) Per the tuned RSE functionals, the EA is in the range of 4.7 eV, and the fundamental gap is about 6.2 eV. In other words, even when accounting for the fact that −εL ≈ EA is usually not a good approximation in HF theory, a difference of over 2 eV seems large. An intuitive explanation is that the electron-attached state, via occupation of the LUMO, affords too much electron repulsion. In calculations that include electron correlation to some degree, such as the DFT calculations with the RSE functionals, the anionic state is significantly stabilized. On the other hand, the effect from ionization is well described with HF, compared to the DFT calculations. While the ionized and electron-attached states are of course different from the excited states, we can nonetheless infer information about the interaction of the HOMO and the LUMO with other orbitals in the molecule from these calculations. In other words, for the cyanines, electron repulsion at the uncorrelated level is large if the LUMO is occupied, and this finding may well be transferable to the excited states. The argument relies on the fact that in each set of calculations the analysis of the linear response transition vectors indicates a relatively pure HOMO to LUMO transition.

4. CONCLUSIONS Electronic excitations of the cyanine dyes are certainly very interesting. In the present work, the performance of TDDFT has been investigated in relation to the DFT delocalization error (DE) and other factors that influence the calculated singlet and triplet excitations. A detailed analysis was performed for system CN9 of Figure 1. For an excitation with significant, or at least some, CT character, ‘pure’ DFT and HF often bracket the excitation energy, with HF being too high and DFT being too low. Tinkering with exact exchange in a hybrid functional is then likely to produce the desired result eventually. For the cyanines, HF and pure DFT produce too high singlet excitation energies, and the global hybrid PBE0 performs even worse. As it has been shown in recent years for a sizable number of cases (see citations in the Introduction and Computational Details), often a rational K

dx.doi.org/10.1021/ct400649r | J. Chem. Theory Comput. XXXX, XXX, XXX−XXX

Journal of Chemical Theory and Computation

Article

and systematic way of finding a well-performing functional for TDDFT is ‘optimal tuning’ of RSE hybrids, where the parameters are determined in order to satisfy a selection of basic properties of (generalized) KS systems. Tuned RSE functionals also tend to produce functionals with very small DEs (albeit moleculedependent). The class of problems that long-range corrected RSE functionals have been designed to fix are related to the asymptotic behavior of the XC potential and associated features. The present analysis suggests that such problems that are not primary issues for the cyanines. Moreover, minimizing the DE with RSE functionals turns out to be uncorrelated with better performance in TDDFT calculations, despite the fact that electron delocalization is evidently very important for the cyanines. Flexible Minnesota functionals that were recently shown to perform reasonably well for cyanine excitations afford significant DEs, although the DEs are smaller than for functionals such as BP, PBE, or PBE0. The better performance must be attributed to a better treatment of correlation while at the same time there is not much of a penalty for the presence of a DE. Grimme and Neese previously suggested6 that the differential correlation between the ground and excited state, and the lack of its proper treatment by many functionals, lies at the heart of the problem. Indeed, electron correlation, and the crucial balance of short-range exchange and correlation in the energy functional, the XC potential, and in the XC response kernel, are crucially important while the DE takes a secondary role. The longestwavelength excitations of symmetric cyanines have a low CT character, which exposes problems with TDDFT calculations in a somewhat different manner than what has usually been found for CT and ‘CT-like’ excitations. A numerical analysis for CN9 has shown that small (evidently too small) contributions from the DFT XC response kernel f XC cannot counter-balance a very large HOMO−LUMO exchange electron repulsion integral, giving too high singlet excitation energies with TDDFT. For an excitation with significant CT character and a large weight from a single occupied-unoccupied MO pair, the exchange integral would be much smaller. The lack of large f XC integrals at the same time causes a ΔSCF calculation to underestimate the singlet excitation energy severely. A comparison of singlet and triplet excitations further indicates that the HOMO−LUMO energy gap Δε is also to blame, since the TDDFT triplet excitation energy with the PBE GGA functional are too low. If the CIS(D) triplet excitation energy is reliable, the doubles correction furnishes small correlation effects for the triplet energyfor CN9 as well as for the other cyanines in the set. Interestingly, a too small Δε apparently provides some error compensation for the PBE based TDDFT singlet excitation while the global hybrid functionals expose the correlation problem even more. Given the observed trends, for molecules where electron correlation plays a similar role as in the cyanine series but where the excitations afford some (not too much) CT character, the results of a TDDFT calculation with a standard functional may end up close to the correct excitation energy. Future work will determine whether TDDFT problems with rhodamine and rosamine dyes mentioned in the Introduction have a similar origin as those for the cyanines.



behavior for MP2 and BP geometries. Singlet and triplet excitation energies for various functionals with percent contributions from the HOMO−LUMO orbital pair. Singlet and triplet energies calculated with CC2 and CIS(D). Eigenvalues from SCF instability calculations. This material is available free of charge via the Internet at http://pubs.acs.org.



AUTHOR INFORMATION

Corresponding Author

*E-mail: jochena@buffalo.edu. Notes

The authors declare no competing financial interest.



ACKNOWLEDGMENTS This research has been supported by the National Science Foundation, Grant Nos. CHE-0952253 and CHE-1265833. We thank the Center for Computational Research at the University at Buffalo for providing computational resources. J.A. thanks Profs. Boris Le Guennic, Denis Jacquemin, Leeor Kronik, and Tom Ziegler for comments on the manuscript.



REFERENCES

(1) Jacquemin, D.; Wathelet, V.; Perpète, E. A.; Adamo, C. J. Chem. Theory Comput. 2009, 5, 2420−2435. (2) Jacquemin, D.; Perpète, E. A.; Ciofini, I.; Adamo, C. Acc. Chem. Res. 2009, 42, 326−34. (3) Jacquemin, D.; Perpète, E. A.; Scuseria, G. E.; Ciofini, I.; Adamo, C. J. Chem. Theory Comput. 2008, 4, 123−135. (4) Schreiber, M.; Buß, V.; Fülscher, M. P. Phys. Chem. Chem. Phys. 2001, 3, 3906−3912. (5) Champagne, B.; Guillaume, M.; Zutterman, F. Chem. Phys. Lett. 2006, 425, 105−109. (6) Grimme, S.; Neese, F. J. Chem. Phys. 2007, 127, 154116. (7) Send, R.; Valsson, O.; Filippi, C. J. Chem. Theory Comput. 2011, 7, 444−455. (8) Jacquemin, D.; Zhao, Y.; Valero, R.; Adamo, C.; Ciofini, I.; Truhlar, D. G. J. Chem. Theory Comput. 2012, 8, 1255−1259. (9) Fabian, J. Dyes Pigm. 2010, 84, 36−53. (10) Calitree, B.; Donnelly, D. J.; Holt, J. J.; K.Gannon, M.; Nygren, C. L.; Autschbach, J.; Detty, M. R. Organometallics 2007, 26, 6248−6257. (11) Rudolph, M.; Ziegler, T.; Autschbach, J. Chem. Phys. 2011, 391, 92−100. (12) Autschbach, J.; Jorge, F. E.; Ziegler, T. Inorg. Chem. 2003, 42, 2867−2877. (13) Kuritz, N.; Stein, T.; Baer, R.; Kronik, L. J. Chem. Theory Comput. 2011, 7, 2408−2415. (14) Kuhn, H. J. Chem. Phys. 1949, 17, 1198−1212. (15) Autschbach, J. J. Chem. Educ. 2007, 84, 1840−1845. (16) Griffiths, J. Colour and Constitution of Organic Molecules; Academic Press: London, 1976, pp 246−247. (17) Jacquemin, D.; Perpète, E. A.; Ciofini, I.; Adamo, C.; Valero, R.; Zhao, Y.; Truhlar, D. G. J. Chem. Theory Comput. 2010, 6, 2071−2085. (18) Masunov, A. E. Int. J. Quantum Chem. 2010, 110, 3095−3100. (19) Meguellati, K.; Ladame, S.; Spichty, M. Dyes Pigm. 2011, 90, 114− 118. (20) Zhao, Y.; Truhlar, D. G. Acc. Chem. Res. 2008, 41, 157−167. (21) Cohen, A. J.; Mori-Sánchez, P.; Yang, W. Science 2008, 321, 792− 794. (22) Hiberty, P. C.; Shaik, S. Phys. Chem. Chem. Phys. 2004, 6, 224− 231. (23) Baer, R.; Livshits, E.; Salzner, U. Annu. Rev. Phys. Chem. 2010, 61, 85−109. (24) Refaely-Abramson, S.; Sharifzadeh, S.; Govind, N.; Autschbach, J.; Neaton, J. B.; Baer, R.; Kronik, L. Phys. Rev. Lett. 2012, 109, 226405. (25) Stein, T.; Autschbach, J.; Govind, N.; Kronik, L.; Baer, R. J. Phys. Chem. Lett. 2012, 3, 3740−3744.

ASSOCIATED CONTENT

S Supporting Information *

BP and MP2 geometries optimized with cc-pVQZ basis sets. Basis set convergence for the CNn series with LC-PBE functional. Tuning parameters and fractional occupation L

dx.doi.org/10.1021/ct400649r | J. Chem. Theory Comput. XXXX, XXX, XXX−XXX

Journal of Chemical Theory and Computation

Article

(26) Körzdörfer, T.; Parrish, R. M.; Sears, J. S.; Sherrill, C. D.; Brédas, J.-L. J. Chem. Phys. 2012, 137, 124305. (27) Frisch, M. J.; Trucks, G. W.; Schlegel, H. B.; Scuseria, G. E.; Robb, M. A.; Cheeseman, J. R.; Scalmani, G.; Barone, V.; Mennucci, B.; Petersson, G. A.; Nakatsuji, H.; Caricato, M.; Li, X.; Hratchian, H. P.; Izmaylov, A. F.; Bloino, J.; Zheng, G.; Sonnenberg, J. L.; Hada, M.; Ehara, M.; Toyota, K.; Fukuda, R.; Hasegawa, J.; Ishida, M.; Nakajima, T.; Honda, Y.; Kitao, O.; Nakai, H.; Vreven, T.; Montgomery, Jr., J. A.; Peralta, J. E.; Ogliaro, F.; Bearpark, M.; Heyd, J. J.; Brothers, E.; Kudin, K. N.; Staroverov, V. N.; Kobayashi, R.; Normand, J.; Raghavachari, K.; Rendell, A.; Burant, J. C.; Iyengar, S. S.; Tomasi, J.; Cossi, M.; Rega, N.; Millam, J. M.; Klene, M.; Knox, J. E.; Cross, J. B.; Bakken, V.; Adamo, C.; Jaramillo, J.; Gomperts, R.; Stratmann, R. E.; Yazyev, O.; Austin, A. J.; Cammi, R.; Pomelli, C.; Ochterski, J. W.; Martin, R. L.; Morokuma, K.; Zakrzewski, V. G.; Voth, G. A.; Salvador, P.; Dannenberg, J. J.; Dapprich, S.; Daniels, A. D.; Farkas, O.; Foresman, J. B.; Ortiz, J. V.; Cioslowski, J.; Fox, D. J. Gaussian 09, Revision A.02; Gaussian Inc.: Wallingford CT, 2009. (28) Dunning, T. H., Jr. J. Chem. Phys. 1989, 90, 1007. (29) Bylaska, E. J.; de Jong, W. A.; Govind, N.; Kowalski, K.; Straatsma, T. P.; Valiev, M.; van Dam, J. J.; Wang, D.; Apra, E.; Windus, T. L.; Hammond, J.; Autschbach, J.; Aquino, F.; Nichols, P.; Hirata, S.; Hackler, M. T.; Zhao, Y.; Fan, P.-D.; Harrison, R. J.; Dupuis, M.; Smith, D. M. A.; Glaesemann, K.; Nieplocha, J.; Tipparaju, V.; Krishnan, M.; Vazquez-Mayagoitia, A.; Jensen, L.; Swart, M.; Wu, Q.; Van Voorhis, T.; Auer, A. A.; Nooijen, M.; Crosby, L. D.; Brown, E.; Cisneros, G.; Fann, G. I.; Fruchtl, H.; Garza, J.; Hirao, K.; Kendall, R.; Nichols, J. A.; Tsemekhman, K.; Wolinski, K.; Anchell, J.; Bernholdt, D.; Borowski, P.; Clark, T.; Clerc, D.; Dachsel, H.; Deegan, M.; Dyall, K.; Elwood, D.; Glendening, E.; Gutowski, M.; Hess, A.; Jaffe, J.; Johnson, B.; Ju, J.; Kobayashi, R.; Kutteh, R.; Lin, Z.; Littlefield, R.; Long, X.; Meng, B.; Nakajima, T.; Niu, S.; Pollack, L.; Rosing, M.; Sandrone, G.; Stave, M.; Taylor, H.; Thomas, G.; van Lenthe, J.; Wong, A.; Zhang, Z. NWChem, A Computational Chemistry Package for Parallel Computers, Version 6.1 (2012 developer’s version); Pacific Northwest National Laboratory: Richland, WA, 2012. (30) Valiev, M.; Bylaska, E. J.; Govind, N.; Kowalski, K.; Straatsma, T. P.; Van Dam, H. J. J.; Wang, D.; Nieplocha, J.; Apra, E.; Windus, T. L.; de Jong, W. A. Comput. Phys. Commun. 2010, 181, 1477−1489. (31) Srebro, M.; Autschbach, J. J. Chem. Theory Comput. 2012, 8, 245− 256. (32) Yanai, T.; Tew, D. P.; Handy, N. C. Chem. Phys. Lett. 2004, 393, 51−57. (33) Perdew, J. P.; Burke, K.; Ernzerhof, M. Phys. Rev. Lett. 1998, 80, 891. (34) Kronik, L.; Stein, T.; Refaely-Abramson, S.; Baer, R. J. Chem. Theory Comput. 2012, 8, 1515−1531. (35) Srebro, M.; Autschbach, J. J. Phys. Chem. Lett. 2012, 3, 576−581. (36) Sun, H.; Autschbach, J. ChemPhysChem 2013, 14, 2450−2461. (37) Autschbach, J. ChemPhysChem 2009, 10, 1757−1760. (38) Levy, M.; Perdew, J. P.; Sahni, V. Phys. Rev. A 1984, 30, 2745− 2748. (39) Weigend, F.; Ahlrichs, R. Phys. Chem. Chem. Phys. 2005, 7, 3295− 3305. (40) Livshits, E.; Baer, R. Phys. Chem. Chem. Phys. 2007, 9, 2932−2941. (41) Karolewski, A.; Kronik, L.; Kümmel, S. J. Chem. Phys. 2013, 138, 204115. (42) Shao, Y.; Molnar, L. F.; Jung, Y.; Kussmann, J.; Ochsenfeld, C.; Brown, S. T.; Gilbert, A. T.; Slipchenko, L. V.; Levchenko, S. V.; O’Neill, D. P.; DiStasio, R. A., Jr.; Lochan, R. C.; Wang, T.; Beran, G. J.; Besley, N. A.; Herbert, J. M.; Yeh Lin, C.; Van Voorhis, T.; Hung Chien, S.; Sodt, A.; Steele, R. P.; Rassolov, V. A.; Maslen, P. E.; Korambath, P. P.; Adamson, R. D.; Austin, B.; Baker, J.; Byrd, E. F. C.; Dachsel, H.; Doerksen, R. J.; Dreuw, A.; Dunietz, B. D.; Dutoi, A. D.; Furlani, T. R.; Gwaltney, S. R.; Heyden, A.; Hirata, S.; Hsu, C.-P.; Kedziora, G.; Khalliulin, R. Z.; Klunzinger, P.; Lee, A. M.; Lee, M. S.; Liang, W.; Lotan, I.; Nair, N.; Peters, B.; Proynov, E. I.; Pieniazek, P. A.; Min Rhee, Y.; Ritchie, J.; Rosta, E.; David Sherrill, C.; Simmonett, A. C.; Subotnik, J. E.; Lee Woodcock, H., III; Zhang, W.; Bell, A. T.; Chakraborty, A. K.;

Chipman, D. M.; Keil, F. J.; Warshel, A.; Hehre, W. J.; Schaefer, H. F., III; Kong, J.; Krylov, A. I.; Gill, P. M. W.; Head-Gordon, M. Phys. Chem. Chem. Phys. 2006, 8, 3172−3191. (43) Head-Gordon, M.; Rico, R. J.; Oumi, M.; Lee, T. J. Chem. Phys. Lett. 1994, 219, 21−29. (44) Baerends, E. J.; Ziegler, T.; Autschbach, J.; Bashford, D.; Bérces, A.; Bickelhaupt, F. M.; Bo, C.; Boerrigter, P. M.; Cavallo, L.; Chong, D. P.; Deng, L.; Dickson, R. M.; Ellis, D. E.; van Faassen, M.; Fan, L.; Fischer, T. H.; Fonseca Guerra, C.; Ghysels, A.; Giammona, A.; van Gisbergen, S. J. A.; Götz, A. W.; Groeneveld, J. A.; Gritsenko, O. V.; Grüning, M.; Gusarov, S.; Harris, F. E.; van den Hoek, P.; Jacob, C. R.; Jacobsen, H.; Jensen, L.; Kaminski, J. W.; van Kessel, G.; Kootstra, F.; Kovalenko, A.; Krykunov, M. V.; van Lenthe, E.; McCormack, D. A.; Michalak, A.; Mitoraj, M.; Neugebauer, J.; Nicu, V. P.; Noodleman, L.; Osinga, V. P.; Patchkovskii, S.; Philipsen, P. H. T.; Post, D.; Pye, C. C.; Ravenek, W.; Rodríguez, J. I.; Ros, P.; Schipper, P. R. T.; Schreckenbach, G.; Seldenthuis, J. S.; Seth, M.; Snijders, J. G.; Solà, M.; Swart, M.; Swerhone, D.; te Velde, G.; Vernooijs, P.; Versluis, L.; Visscher, L.; Visser, O.; Wang, F.; Wesolowski, T. A.; van Wezenbeek, E. M.; Wiesenekker, G.; Wolff, S. K.; Woo, T. K.; Yakovlev, A. L. Amsterdam Density Functional, SCM, Theoretical Chemistry; Vrije Universiteit: Amsterdam, The Netherlands; URL http://www.scm.com (accessed 09/2013). (45) Vosko, S. H.; Wilk, L.; Nusair, M. Can. J. Phys. 1980, 58, 1200− 1211. (46) Jacquemin, D. J. Phys. Chem. A 2011, 115, 2442−2445. (47) Körzdörfer, T.; Sears, J. S.; Sutton, C.; Brédas, J.-L. J. Chem. Phys. 2011, 135, 204107. (48) Stein, T.; Eisenberg, H.; Kronik, L.; Baer, R. Phys. Rev. Lett. 2010, 105, 266802. (49) Refaely-Abramson, S.; Baer, R.; Kronik, L. Phys. Rev. B 2011, 84, 075144. (50) Salzner, U.; Aydin, A. J. Chem. Theory Comput. 2011, 7, 2568− 2583. (51) Hirata, S.; Head-Gordon, M. Chem. Phys. Lett. 1999, 314, 291− 299. (52) Peach, M. J. G.; Williamson, M. J.; Tozer, D. J. J. Chem. Theory Comput. 2011, 7, 3578−3585. (53) Ziegler, T.; Rauk, A.; Baerends, E. J. Theor. Chim. Acta 1977, 43, 261−271. (54) Rosa, A.; Ricciardi, G.; Gritsenko, O.; Baerends, E. J. Excitation energies of metal complexes with time−dependent density functional theory. In Principles and Applications of Density Functional Theory in Inorganic Chemistry I, Vol. 112; Kaltsoyannis, N.; McGrady, J. E., Eds.; Springer: Heidelberg, 2004; pp 49−115. (55) Ziegler, T. J. Chem. Soc., Dalton Trans. 2002, 642−652. (56) Casida, M. E.; Gutierrez, F.; Guan, J.; Gadea, F.-X.; Salahub, D. R.; Daudey, J.-P. J. Chem. Phys. 2000, 113, 7062−7071. (57) Cullen, J.; Krykunov, M.; Ziegler, T. Chem. Phys. 2011, 391, 11− 18. (58) Ziegler, T.; Seth, M.; Krykunov, M.; Autschbach, J. J. Chem. Phys. 2008, 129, 184114. (59) Ziegler, T.; Seth, M.; Krykunov, M.; Autschbach, J.; Wang, F. J. Chem. Phys. 2009, 130, 154102. (60) Ziegler, T.; Krykunov, M. J. Chem. Phys. 2010, 133, 074104. (61) Peach, M. J. G.; Benfield, P.; Helgaker, T.; Tozer, D. J. J. Chem. Phys. 2008, 128, 044118.

M

dx.doi.org/10.1021/ct400649r | J. Chem. Theory Comput. XXXX, XXX, XXX−XXX