Qualitatively Incorrect Features in the TDDFT Spectrum of Thiophene-Based Compounds.

Ab initio molecular electronic structure computations of thiophene-based compounds constitute an active field of research prompted by the growing interest in low-cost materials for organic electronic devices. In particular, the modeling of electronically excited states and other time-dependent phenomena has moved toward the description of more realistic albeit challenging systems. We demonstrate that due to its underlying approximations, time-dependent density functional theory predicts results that are qualitatively incorrect for thiophene and thienoacenes, although not for oligothiophene chains. The failure includes spurious state inversion and excitation characters, wrong distribution of oscillator strengths and erroneous potential energy surfaces. We briefly analyze possible origins of this behavior and identify alternative methods that alleviate these problems.

Thiophene-based materials play a central role in the field of organic electronics.[1−3] Derivatives of thiophene, oligothiophenes and oligothienoacenes are extensively used in organic photovoltaics,[4−10] light emitting diodes,[11−15] field effect transistors,[16−21] and so forth. Typically high extinction coefficients make them particularly suitable for solar cell materials. As a matter of fact, organic dyes, which do not contain thiophene motifs, seem to be in minority. This is part of the reason why excited states of simple thiophene compounds have drawn interest from numerous theoretical perspectives, including spectroscopy,[22−36] excited state geometries,[24,30,34] and nonadiabatic molecular dynamics.[37,38] In this Letter, we discuss the computational conundrum associated with the low-lying bright ππ* singlet excited states of thiophene, short oligothiophenes, and oligothienoacenes (planar fused oligomers). Figure 1 compares the absorption spectra of thiophene computed with linear response time-dependent density functional theory (TDDFT) and a post-Hartree–Fock (post-HF) formalism. The overall band structure appears rather similar, yet decomposition into individual states (vide infra) reveals an inversion that could dramatically impact the computational prediction of photochemical and photophysical processes of thiophene-based compounds. We believe that this dichotomy and especially its consequences have been overlooked in the literature. The following analysis discusses the underlying origin for this contrasting behavior and identifies methods that alleviate the problem.

Figure 1

Photoabsorption spectra computed from the Wigner distribution (red line) of thiophene. The spectra were decomposed into contributions from two lowest excited states, S1 and S2. The red peak is mainly due to the excitations from highest occupied to the lowest unoccupied molecular orbital (HOMO → LUMO, B2), whereas the blue peak is dominated by HOMO-1 → LUMO, A1, character. The cc-pVTZ basis set was used.

TDDFT,[39] or more precisely linear response TDDFT within the adiabatic approximation,[40] has emerged as the most widely used method for molecular excited states.[41] Its exalted status certainly arises from the combination of low computational cost and, in many cases, amazing predictive power for both excitation energies[42] and excited state properties[43] of fairly large systems. This contrasts sharply with highly accurate multireference methods, such as CASPT2 (complete active space second order perturbation theory)[44] and MCQDPT2 (multiconfigurational quasi-degenerate second order perturbation theory),[45] which can be applied only to small systems that appear particularly challenging (of relevance to the present study, we point out the controversy with CASPT2 results for bithiophene[32]). Although TDDFT is a formally exact theory, standardly employed approximations result in several limitations. These include the underestimation of charge-transfer and Rydberg excitations,[46] underestimation of triplet excitation energies,[47] and a general inapplicability for high-lying excited states.[48] The ubiquitous adiabatic approximation is responsible for the inaccurate description of conical intersections as well as the absence of doubly and multiply excited states[49,50] (though it also underlies the charge-transfer problem[51]). Valence ππ* states with mainly single excitation character (which are studied here) are usually considered less problematic, although exceptions do exist. Grimme and Parac studied the two lowest singlet ππ* states of oligoacenes,[52] short-axis polarized L (HOMO → LUMO) and long-axis polarized L (HOMO-1 → LUMO, HOMO → LUMO+1). They found L to be strongly underestimated by hybrid and GGA (generalized gradient approximation) functionals, with spurious inversion of states in case of naphthalene. L was later described as an “ionic”,[52] “charge-transfer in disguise”[53] and “charge-transfer-like”[54] excitation. Substantial improvements were found with range separated hybrid functionals.[55,56] Most recently, by employing a range separated hybrid functional with optimally tuned parameter,[57] Baer, Kronik et al.[54] reported excellent agreement with reference CC2 (approximate coupled cluster singles and doubles)[58] results, thereby restoring the predictive power of TDDFT for the oligoacene series.

Here we focus upon the ππ* state inversion in related categories of π-conjugated thiophenes and oligothienoacenes that are, nevertheless, distinct in nature due to the presence of third row heteroatoms (i.e., low-lying d-orbitals close in energy) and lower symmetry. The problem illustrated for thiophene in Figure 1 is further examined by computing the vertical excitation energies of the two lowest ππ* states at various TDDFT and post-Hartree–Fock levels (Figure 2a). TDDFT generally predicts that two states, one dominated by the HOMO → LUMO (B2) and another by the HOMO-1 → LUMO (A1) transition, appear in a reversed order (or nearly degenerate) when compared to more accurate post-HF methods. The dependence on the specific exchange-correlation kernels and especially on the fraction of exact exchange is assessed via consideration of different categories of functionals: PBE (GGA), B3LYP and PBE0 (global GGA hybrid), M06 (meta-GGA hybrid), M06-2X (meta-GGA hybrid with 54% Hartree–Fock exchange), long-range corrected hybrid, ωB97x-D as well as optimally tuned long-range corrected hybrid LC-PBE*. Despite their conceptual differences, each of these functionals gives the same qualitative results: two transitions of similar intensities with B2 being slightly lower than A1. Including an implicit polar solvent does not alter the ordering of the transitions and only slightly affects the excitation energies (see Supporting Information). This state ordering is amplified with the CIS (configuration interaction singles) approach, and its response analogue TD-HF (time dependent Hartree–Fock).[59] However, the energy of the HOMO-1 → LUMO (A1) state drops by more than 1 eV, and the sequence reverses when correlation is added to CIS by means of doubles and approximate triples using the CIS(D) method.[60] This indicates that A1 is highly sensitive to the correlation effects introduced by CIS(D). CC2 and the second order algebraic diagrammatic construction, ADC(2),[61] predict the same ordering as CIS(D). This post-HF picture is further solidified by the computationally more involved EOM-CCSD (equation of motion coupled cluster singles doubles),[62] SAC–CI (symmetry adapted cluster–configuration interaction)[63] and multi state MS-CASPT2 of ref (34). Additionally, the agreement achieved by the double hybrid TDA-B2PLYP with wavefunction-based methods is excellent. Similar in spirit to CIS(D), the double hybrid approach of Grimme and Neese[64] adds perturbative second-order corrections to the CIS-like representation of TDDFT (i.e., TDDFT in Tamm-Dancoff approximation; TDA). Although less dramatic, an inversion similar to CIS/CIS(D) is found for the two excited states of thiophene with the double hybrid after introducing the perturbative correction (for details, see Table 1 in Supporting Information).

Figure 2

Excitation energies of the two lowest ππ* states (S1 and S2) of (a) thiophene and (b) dithieno[2,3-b:2′,3′-d]thiophene. Radii of the circles are proportional to the oscillator strengths, and the colors indicate the dominant character. Relevant M06 orbitals are shown. Note that there is no qualitative difference between Kohn–Sham and Hartree–Fock orbitals. The cc-pVQZ basis set was employed. *EOM-CCSD indicates that the energies of dithienothiophene were computed with the cc-pVTZ basis set. The **CASPT2 results for thiophene are taken from ref (34) and correspond to experimental geometry of Bak et al.[65]

Similar patterns are seen for the A1 and B2 representations of fused dithienothiophene (Figure 2b), which belongs to the same symmetry point group, C2, as the thiophene example discussed above. Transiting from CIS to CIS(D) inverses the energy ordering of the two ππ* excitations. The TDDFT computations generally match the CIS trend, with the exception of pure PBE, which still gives too low excitation energies. On the other hand, wavefunction-based methods, along with TDA-B2PLYP, slightly improve the CIS(D) results through reduction of the energy gap between the two states. The possibility of double excitations playing a role can be excluded as diagnostic tools associated with ADC(2) and CC2 indicate that both states are dominated by single excitations from the single configurational ground state. The origin of the TDDFT problems is more subtle than that of CIS (i.e., correlation effects). In fact, from the energetic perspective the HOMO-1 → LUMO state is reasonably well positioned by both the reference wavefunction methods and TDDFT, especially when functionals with moderate amount of exact exchange are chosen. By carefully examining Figure 2a and b and drawing a parallel with the oligoacene example discussed above,[52−56] it appears that the HOMO → LUMO state is underestimated by TDDFT. Inclusion of exact exchange shifts the energy of this transition; however, it simultaneously shifts the other state, which is affected by correlation effects. The overall result being that no qualitative improvement is observed. Optimally tuned long-range corrected functional[54] performs no better than hybrid and long-range corrected hybrid functionals. The state inversion in the TDDFT context is certainly caused by interplay of different effects that influence each of the two states. Because none of the standard functionals (within the strict linear response TDDFT) that are tested here provides a suitable description and the “single state engineering” through increasing the amount of exact exchange is ineffective, the error is inherently connected with subtler effects in the description of correlation for the excited state. The trends obtained with the double hybrid are qualitatively correct, emphasizing the importance of perturbative doubles corrections applied to TDA. Upon inclusion of virtual orbitals, the CIS(D)-like correction captures nonlocal correlation effects missed by conventional approximate functionals in the adiabatic approximation, as advocated by Grimme et al.[64,66,67] We note, however, that the (D) correction term is only a pragmatic fix in the context of TDDFT. The mechanisms leading to an improvement between uncorrected TDDFT and TDDFT/CIS(D), in the framework of exact TDDFT, are surely more complex to grasp than those between CIS and CIS(D). In this context, we notice that even the Tamm–Dancoff approximation applied to certain functionals (for instance PBE0) may improve the results of the parent TDDFT approach.

C2 thienothiophene (Figure 3a) is a somewhat more complicated example. CIS(D) again inverts CIS, and so does the double hybrid with respect to its underlying TDA energies. Yet the TDA-B2PLYP excitation energies do not align well with those from the post-Hartee–Fock methods. The TDDFT picture contrasts with the wavefunction-based methods with the intense state lying slightly below the second weaker state. Becaues of symmetry reasons (two states of Bu representation close in energy), our analysis of relaxed density differences reveals a large degree of character mixing, rather than the clear-cut inversion implied by the simplified orbital representation. The character ambiguity arises from the different extents of density increase/depletion at sulfur atoms and the central CC bond (see Supporting Information). Nevertheless, the problematic behavior of TDDFT concerning both the excitation energies and the corresponding oscillator strengths of thienothiophene ππ* states is fully in line with the two former examples (Figure 2). In sharp contrast, nonplanar bithiophene (Figure 3b) and terthiophene (Supporting Information) are rather unproblematic, at least qualitatively. The small dependence on the exact exchange fraction and the energy lowering of the second bright state when moving from CIS and CIS(D) still holds, yet the energy gap between the two ππ* states is consistently large and no inversion occurs. Such a separation arises from the limited conjugation between flexible oligothiophene units. Unlike in thiophene and thienoacenes, where we always refer to S1 and S2, for bithiophene and terthiophene, the two bright transitions are typically S1 and S3/S7, respectively (see Supporting Information).

Figure 3

Excitation energies of the two lowest bright ππ* states of (a) thieno[3,2-b]thiophene (S1 and S2) and (b) 2,2′-bithiophene (S1 and S3/S4). As in Figure 2, radii of the circles are proportional to the oscillator strengths, and colors denote the orbital excitations with largest coefficients. The cc-pVQZ basis set was employed. The *EOM-CCSD and SAC–CI energies of bithiophene were computed at the cc-pVTZ level. The **CASPT2 results of bithiophene are taken from 32 and correspond to the ground state structure optimized at CASPT2 level with C2 symmetry.

The consequence of the state inversion is best illustrated by the absorption spectra of thiophene highlighted earlier (Figure 1) and its derivatives (Figure 4), which are computed using the geometries sampled from semiclassical Wigner distribution (see Computational Methods) at three illustrative levels, PBE0, M06-2X, and ADC(2). Despite the use of an identical set of 500 initial structures generated from a Wigner distribution in the ground electronic state, the TDDFT spectra of thiophene and dithienothiophene contrast with the ADC(2) spectra. The decomposition of the spectrum into the relevant state contributions (color coding is the same as in Figures 1 and 2) reveals that the inversion is not just an artifact of the optimized geometries: the overall peak intensities largely follow the trends given by the oscillator strengths in the static computations. For thiophene, both functionals predict two overlapping peaks, which results in a single unstructured hill, whereas the ADC(2) spectrum shows both a different state ordering and a longer tail at higher energies, consistent with experimentally measured cross sections.[36] The situation is similar for dithienothiophene, where ADC(2) places two states close in energy resulting in an intense peak differing from TDDFT predictions. The seemingly reliable TDDFT spectra of thienothiophene actually emerge from a peculiar mixing of states. Although we intentionally do not attribute the dominant character from the orbital contributions, the contrast between the methods is clearly visible from the small TDDFT intensities of S2, which is consistent with the oscillator strengths computed at the optimized geometry (Figure 3a). On the other hand, the dichotomy does not apply to bithiophene, for which all three methods give similar pictures with only minor difference in the position and peak intensity. Note that the two functionals give similar results in all four cases.

Figure 4

Photoabsorption spectra computed from the Wigner distribution (red line) of dithienothiophene (first row), thienothiophene (second row), and bithiophene (third row). The dithienothiophene and thienothiophene spectra were decomposed into contributions from S1 and S2 and S1 and S3+S4 in the case of bithiophene. The color code reflects the main excitation character (in line with Figures 1 and 2). The same color is used for the bands of thienothiophene to reflect the character ambiguity. The cc-pVTZ basis set was used.

Finally, it is instructive to examine the impact of the preceding observations on the excited state geometries. The inversion of the states seen for ground state equilibrium geometries has little meaning as molecules that relax in the excited state can adopt conformations far from the Franck–Condon region. The PBE0, M06-2X, and ADC(2) optimized minima of the excited state with the initial A1 symmetry (Figure 2) of thiophene are displayed in Figure 5a. ADC(2) predicts the S-puckered minimum, whereas TDDFT converged to the C-puckered structures. The C-symmetric S-puckered geometry is a transition state at the TDDFT level, as pointed out by Marian et al.[30] On the other hand, both DFT-MRCI[30] and CASPT2[34] (δCCCS = 26.7°) optimizations confirm the S-puckered structure as being a minimum. Slightly different S-puckered minima are found when the B2 state is optimized with both PBE0 and M06-2X (Figure 5b), albeit not with ADC(2). Actually, reoptimization of the TDDFT structures with ADC(2) leads to an intersection with the ground state, consistent with CASPT2 computations of Stenrup,[34] suggesting that B2 may be unbound. Even more surprising are the S1 geometries of thienothiophene and dithienothiophene (Figure 5c and d), for which large discrepancies exist for both the geometries and the vertical excitation energies when the two functionals are compared to ADC(2). For both systems, ADC(2) leads to tilted structure, whereas the two functionals prefer nearly planar frameworks. These systems were recently investigated by Jacquemin et al., who looked at how solvation models affect excited state geometries.[68] Despite the different purpose of that work, we noticed that the quality of the results obtained with the usually reliable M06-2X functional should be taken with care. In contrast to previous examples, good agreement is achieved for the S1 minimum of bithiophene. All three methods predict the well-known planar quinoidal structure that shows the characteristic alternation of single and double bonds with only slight disagreement in the lengths of the terminal C–C bonds (see Supporting Information).

Figure 5

Minima at the S1 adiabatic potential energy surface of: (a, b) thiophene, (c) thienothiophene, and (d) dithienothiophene. We report the vertical excitation energies (in parentheses), the dominant orbital excitations at given geometry (denoted by red and blue circles as in Figures 2 and 3), and the relevant dihedral angles. δCCCC of dithienothiophene denotes the dihedral angle between the middle and the side ring. The cc-pVTZ basis set was used.

To summarize, we demonstrated the problematic performance of TDDFT for the two lowest ππ* states of thiophene and short thienoacenes. The failures include incorrect state ordering, poor distribution of oscillator strengths, and erroneous descriptions of critical points on the excited state potential energy surfaces. Similar trends are also identified for furan and its fused derivative (see Supporting Information), but without the spurious inversion. What makes the thiophene example unique with respect to the furo- and oligoacene compounds is the large oscillator strengths of the two states involved and the experimental relevance. The incorrect state ordering and oscillator strengths thus notably impact the computed absorption spectra. In addition, serious discrepancies were found between TDDFT and ADC(2) optimized excited state geometries. These qualitative failures directly affect the possible prediction of adiabatic excitation energies and photoemission properties. Achieving accurate potential energy surfaces is also crucial for excited state molecular dynamics simulations, for which inexpensive TDDFT is often the method of choice. More generally, the present study emphasizes the importance of systematically carrying out careful comparisons with more accurate wavefunction-based methods.

Computational Methods

Ground state optimized structures and vibrational frequencies were obtained at the M06/DGDZVP[69] level with the Gaussian09[70] program package. Excited states were analyzed in their respective point groups: C2 for thiophene, furan, and dithienothiophene, C2 for thienothiophene and furofuran, C2 for bithiophene, and C for terthiophene. Excitation energies were computed with Gaussian09 (M06, M06-2X, ωB97x-D, LC-PBE, LC-PBE*, TD-HF, CIS/CIS(D), EOM-CCSD, SAC–CI), Turbomole 6.5[71] (PBE, PBE0, TDA(PBE0), B3LYP, ADC(2), CC2), and Orca 3.0.2[72] (TDA-B2PLYP), and numerical values are given in Supporting Information (Tables 1 and 2). In ADC(2) and CC2 computations, we apply the resolution of identity and the frozen core approximations. For the B2PLYP computations, the resolution of identity was employed. M06 and M06-2X computations use an ultrafine integration grid. SAC–CI is performed with the default parameters and convergence criteria of the Gaussian09 program. The CIS(D) oscillator strengths in Figures 2 and 3 are taken from CIS. LC-PBE* computations were performed by tuning the range separation parameter γ to match the HOMO energy and the difference between the total energies of the cation and neutral molecule. All the excitation energies were converged at the cc-pVQZ level, unless otherwise specified. Comparisons with smaller and augmented basis sets show a small and systematic deviation for two valence excitations. To reduce the computational burden, the photoabsorption spectra and excited state geometry optimizations were performed with the cc-pVTZ basis set. Excited state geometries were optimized with Turbomole 6.5 (PBE0, ADC(2)) and Gaussian09 (M06-2X) with tight convergence criteria and no symmetry constraints. The nuclear configurations used for spectral simulations were sampled by an uncorrelated Wigner distribution in the ground sate,[73,74] as implemented in Newton-X package.[75] A total of 500 structures were taken for each compound and vertical excitation energies and oscillator strengths were computed. Each transition was broadened by a Lorentzian using a phenomenological width of 0.05 eV. The spectra were decomposed into contributions from different states and the main character was assigned according to the dominant orbital excitations. Molecules, orbitals, and difference densities were visualized with VMD 1.9.1 program package.[76]