Accurate excitation energies of molecules and oligomers from a semilocal density functional

Excitation energy is the energy absorbed to transit one electron from the ground state to an excited state. It plays a predominant role in energy conversion, biological processes, chemiluminescence, and controllable preparation of molecular devices.^1,2 Contrary to their ground state counterparts, excited states are short-lived, highly reactive and often coupled together, leading to the quest to master their manipulation being an extremely involved demand. However, accurate description of molecular excitation remains a major challenge for chemists and physicists.^1,2 Time-dependent density functional theory (TDDFT)^3–9 based on the linear-response theory is a powerful and most widely used tool for efficiently evaluating dynamical properties of molecules.^4–9 

In TDDFT, a critical issue is to approximate the time-dependent dynamical effective potential,^3–5 which consists of the kinetic part, the classical Hartree potential, the exchange-correlation potential, as well as the external part. Among these ingredients, the exchange-correlation part is essential, because other contributions can be expressed explicitly in terms of the instantaneous Kohn-Sham orbitals. There are two general ways to approximate the time-dependent dynamical exchange-correlation potential. One way and also the simplest way is to approximate it from the ground-state potential but with the ground-state electron density or orbitals replaced by the instantaneous electron density or orbitals. This is the adiabatic approximation. The other way is to go beyond the adiabatic approximation by developing a frequency-dependent or memory part. The former is simple and easy to be implemented, while the latter is more complicated but more useful.

Many benchmark studies^10–31 have been focused on the calculation of low-lying vertical excitation energies within the adiabatic approximation. Among these studies, most of them have been devoted to singlet valence excited states,^10–13 with some considerations on triplet valence^14–16 and singlet and triplet Rydberg excited states.^14–31 In many previous studies^{10,11,14–16} the geometry optimization was performed using a single level of theory (e.g., PBE0 functional³² or Møller Plesset second-order perturbation theory). Based on the optimized geometry, one calculated excitation energies with various density functionals. Alternatively, the same density functional can be used for both the ground-state geometry optimization and excitation energy calculation. The latter approach appears more consistent. At energy minimum, density functionals should be able to describe ground-state properties, such as geometries. Few benchmark studies have been focused on the performance of meta-generalized gradient approximation (meta-GGA) density functionals for the calculation of excitation energies within TDDFT.^18–24 It is therefore very interesting to assess the performance of meta-GGA density functionals on various excitation energies such as singlet, triplet, valence, and Rydberg excited states within TDDFT formalism.^25–32 Such benchmark tests directly reflect the reliability of exchange-correlation (xc) density functionals, because low-lying excitation is closely related to the performance of the ground-state DFT.

Recently, Tao and Mo (TM)³³ have developed a meta-GGA functional. The basic idea is based on an observation that the conventional exchange hole is fully nonlocal, but it can be made localized under a general coordinate transformation.^34,35 With this property, we can perform a density matrix expansion of the transformed exchange hole. This essentially leads to the TM exchange functional. Many exact or nearly exact properties on the exact exchange energy as well as on the exact exchange hole have been automatically built into this functional. Unlike many other semilocal functionals, for which the underlying exchange-correlation holes were unknown when they were originally proposed, TM functional is directly obtained from its associated exchange-correlation hole. The TM correlation part is based on a modification of the original TPSS correlation functional. In this work, we make an assessment of the TM functional on low-lying molecular excitation energies. To have a better understanding of the exchange part, we also include the combination TMTPSS (TMx + TPSSc)^36–39 in this evaluation. Our previous studies show that both functionals can achieve high accuracy for proton affinities, harmonic vibrational frequencies, hydrogen-bonded dissociation energies, bond lengths, low-lying atomic excitation energies, and lattice constants.^36–39 Here we show that TM and TMTPSS can achieve good performance on the excitation energies of small molecules and conjugated oligomers. They are comparable to or even more accurate than computationally more expensive hybrids for excitation energies of small molecules, but less accurate than hybrids for conjugated oligomers, due to the nonlocality problem.

All calculations were performed using Gaussian 09,⁴⁰ in which the TM and TMTPSS functionals have been implemented^36–39 via a local modification. To make our assessment reliable, a large basis set 6-311++G(3df, 3pd) was used. The ground-state geometries of molecules and 3-hexyl-thiophene oligomers were optimized with respective density functionals LSDA,⁴¹ PBE,⁴² TPSS,⁴³ TM,^33,36,37 TMTPSS,^38,39 B3LYP,⁴⁴ and PBE0.³² Then we evaluated the vertical excitation energies based on the optimized ground-state geometries. Since our calculations involve the treatment of both ground-state (via optimization of molecular geometries) and excited states (within the adiabatic approximation in TDDFT), for consistency, we chose the same basis in all the calculations. This large basis set was previously used to perform a comprehensive assessment⁴⁵ of the TPSS functional for molecules. The ultrafine grid (Grid = UltraFine) in numerical integration and the tight self-consistent field convergence criterion (SCF = Tight) were used. The TDDFT calculations of 3-hexyl-thiophene oligomers in chloroform solvent were performed using the PCM (polarizable continuum model) method. To make our comparison consistent, we performed our own calculations with all reference functionals, rather than attempting to extract data from the literature for 3-hexyl-thiophene oligomers. Throughout the paper, we calculated the mean error (ME) using the sign convention: error = theory-experiment. A mean error shows us whether excitation energies are underestimated or overestimated on the average for a particular density functional, while the mean absolute error (MAE) shows us how far a density functional theoretical estimate is from experiment.

Excitation energies of 17 small molecules consisting of 134 excited states have been widely used to benchmark various DFT methods,^10–31 due to their experimental excitation energies in the literature. Table I shows the molecules considered and their excited states. Here the entire benchmark test set is examined with the LSDA,⁴¹ GGA (PBE),⁴² meta-GGA (TPSS,⁴³ TM,³³ and TMTPSS^38,39), and hybrid (B3LYP⁴⁴ and PBE0³²) functionals. The analysis of density functional performance on excited states will be focused on the general trend. Individual errors obtained for each molecule and each DFT method used in the benchmark are available in the supplementary material.

The benchmark set of 134 excited states can be classified into 76 singlet and 58 triplet states.^20,31 Table II lists the ME and MAE for singlet and triplet excited states of 17 molecules. To show the trend, we make a graphical representation of the MAE of these molecules, as shown by Fig. 1. The data corresponding to Table II and Fig. 1 are summarized in Table S1 of the supplementary material. From Fig. 1, we find that the MAEs of the seven adiabatic density functionals in predicting the singlet and triplet excitation energies of molecules can be ordered as PBE < LSDA < TPSS < TM < TMTPSS < B3LYP < PBE0 and PBE < LSDA < TM = B3LYP < PBE0 < TMTPSS < TPSS, respectively. For singlet excited states, the PBE0 gives the best MAE (0.25 eV), followed by B3LYP (0.36 eV), TMTPSS (0.37 eV) and TM (0.40 eV). The PBE functional gives the poorest performance among the DFT methods considered, suggesting the inadequacy in the ground-state exchange-correlation potential. Among semilocal functionals, TMTPSS gives the best singlet excitation energy, and TM also offers the good performance for this excited state. Negative MEs indicate that the calculated excitation energies tend to be red-shifted with respective to the experimental values. Tozer⁴⁶ also shows that the PBE GGA functional underestimates the excitation energies of molecules, while mixing some amount of exact exchange into PBE can increase the excitation energies.

From Table II, we find that for triplet excited states, the PBE0 (MAE = 0.38 eV), TPSS (MAE = 0.36 eV), and TMTPSS (MAE = 0.37 eV) functionals offer the best performances among all the DFT methods, while the PBE (MAE = 0.53 eV) functional is the worst. These results for the triplet excited states are slightly different from the work of Jacquemin et al.,¹⁵ in which they found that the MAEs of PBE0 and B3LYP are slightly larger than 0.40 eV. The difference between the current benchmark study and that of Jacquemin et al.¹⁵ is due to the fact that the performance of the functionals in the present work is compared with experimentally determined excitation energies, while the previous analysis is based on the comparison of density functional performances with wave function methods.

It was shown⁴⁶ that, for triplet states, PBE underestimates the excitation energies, but improvement via the increase of the exact exchange is much less pronounced than for the singlet states. The reason for this difference is clear from an analysis of individual excitations. In many cases the triplet excitation energy does increase via the increase of the exact exchange, while in other cases it drops, leading to a degradation in accuracy. The relatively poor performance of the hybrid functionals in Fig. 1(b) can be understood in terms of the ground state triplet instability problem, consistent with Refs. 25 and 47.

The benchmark set can be also divided into 62 valence and 72 Rydberg states.^20,31 The MEs and MAEs for valence and Rydberg excited states of small molecules are listed in Table III, while Fig. 2 gives the intuitive graphical representation of the statistical errors on the same excited states. The corresponding raw data are available in Tables S1 of the supplementary material. As shown in Table III and Fig. 2, the accuracy of the density functionals in predicting the valence and Rydberg excited states can be ordered as LSDA < PBE < TPSS < TM < TMTPSS < PBE0 < B3LYP and PBE < LSDA < TPSS < TM < TMTPSS < B3LYP < PBE0, respectively. From Table III, we can observe that for valence excited states, the TMTPSS gives an MAE of 0.30 eV, which is almost the same as those of hybrid functionals B3LYP (MAE = 0.27 eV) and PBE0 (MAE = 0.29 eV). These two hybrids give the best overall performances among the DFT methods, but with higher computation cost. The LSDA (MAE = 0.46 eV) performs the worst among the DFT methods studied. The TM functional (MAE = 0.34 eV) offers better performance than those of TPSS and PBE for excited valence states. For Rydberg states, the nonlocal PBE0 (MAE = 0.30 eV) functional gives the best performance among all DFT methods tested here, while the PBE (MAE = 0.69 eV) is the worst functional. TMTPSS (MAE = 0.42 eV) gives the best overall performances among semilocal DFT methods, while the TM functional (MAE = 0.48 eV) provides a better performance than those of TPSS and PBE for Rydberg states.

As shown by Li and Truhlar,⁴⁸ the hybrid functionals PBE0 and B3LYP are quite good for valence excitations but they underestimate Rydberg excitations. The reason is that Rydberg excited states are usually higher relative to the ground state than valence excited states. Therefore, a density tail potential correction is more important for Rydberg excited states, but usual hybrid functionals only add the exact exchange globally. A way to solve this problem is through the development of range-separated density functionals. There are two ways to construct range-separated functionals. One way is to replace the semilocal exchange part with the Hartree-Fock exchange.^49,50 This kind of nonlocal functionals are developed primarily for solid-state calculations. The other is to replace the semilocal long-range part with the Hartree-Fock exchange. These functionals with the long-range correction (LC) scheme are usually for molecular calculations because for molecules, the LC scheme is more important, as demonstrated by many range-separated functionals such as CAM-B3LYP,⁵¹ LC-BOP,^52–54 LC-BLYP,^52,53 and LC-PBE.⁵⁵ From the comparison of BLYP with LC-BLYP and B3LYP with CAM-B3LYP, one can also see that the long-range or density-tail potential correction is important for Rydberg excitations.

Figure 3 summarizes the calculated MAE for each density functional with respect to the total number of excited states calculated using the TDDFT approach. The detailed results for each molecule and overall statistical results can be found in Table S1 of the supplementary material. From Fig. 3, the order of accuracy of the density functionals in predicting the overall excited energies can be placed as PBE < LSDA < TPSS < TM < TMTPSS < B3LYP < PBE0. The PBE0 functional gives the best performance among all DFT methods. From Table S1, we can see that PBE (MAE = 0.69 eV) is the worst performing functional for Rydberg states, as pointed out above. TMTPSS (MAE = 0.37 eV) gives the best overall performance among semilocal functionals. The TM functional (MAE = 0.42 eV) offers a better performance than those of LSDA, PBE, and TPSS. All density functionals possess a negative ME value. The underestimation of excitation energies has been ascribed to the inadequacy in the long-range behavior of the exchange-correlation potential.³¹ A reduction in ME is observed as one goes from LSDA and GGA functionals to meta-GGA and hybrid functionals, suggesting the effect of nonlocality.

To study the influence of electron delocalization upon the performance of DFT, we calculate the oscillator strength and transition energies of the first low-lying excited state transitions for conjugated 3-hexyl-thiophene oligomers in chloroform within the adiabatic TDDFT formalism. The equilibrium structures of 3-hexyl-thiophene oligomers (T)_n (n = 1-8) from TM calculations at 6-311++G(3df, 3pd) basis set are shown in Fig. S1 of the supplementary material. Table IV and Fig. S2 of the supplementary material show the lowest vertical excitation energies calculated for 3-hexyl-thiophene oligomers (T)_n (n = 1-8), in comparison with other DFT methods. As shown in Fig. S2, the lowest excitation energy decreases with repeat units for all DFT methods. From Table IV and Fig. S2, we can see that PBE0 (MAE = 0.532 eV) yields the best performance, while PBE (MAE = 0.898 eV) is the worst performing functional, suggesting the significance of nonlocality via the exact exchange mixing. TMTPSS (MAE = 0.848 eV) gives the best overall performance among semilocal DFT methods. The TM functional (MAE = 0.852 eV) offers a better performance than those of TPSS and PBE for excited states of 3-hexyl-thiophene oligomers, in consistency with molecular excitations.

The performance of newly developed density functionals for the calculation of vertical excitation energies within the linear response TDDFT formalism has been assessed on a benchmark set consisting of 134 experimental excited state energies. The PBE0 (MAE = 0.30 eV), TMTPSS (MAE = 0.42 eV), B3LYP (MAE = 0.43 eV), and TM (MAE = 0.48 eV) functionals are recommended for Rydberg excited states. For valence states, the B3LYP (MAE = 0.27 eV), PBE0 (MAE = 0.29 eV), TMTPSS (MAE = 0.30 eV), TM, and TPSS functionals offer equivalent performance. In particular, the TMTPSS functional is an attractive option for the calculation of valence excited states since it is a semilocal functional, without requiring the exchange integrals needed for the admixture of HF exchange in hybrid functionals. For singlet states, the PBE0 functional offers the best performance, while B3LYP and TMTPSS functional offer a similar good performance. For triplet states, TM and TPSS yield equivalent performance. Overall, the best performing semilocal DFT is TMTPSS (MAE = 0.37 eV), followed by TM (MAE = 0.42 eV). The best overall performing hybrid functional is PBE0 (MAE = 0.29 eV). This finding is consistent with the work of Ref. 57. The best overall density functional based on this benchmark study is PBE0, which offers a balanced treatment of singlet (MAE = 0.25 eV), triplet (MAE = 0.32 eV), valence (MAE = 0.30 eV), and Rydberg (MAE = 0.25 eV) excited states. For conjugated oligomers, TPSS, TM, and TMTPSS are at the same level of accuracy, while hybrids B3LYP, and particularly PBE0 give better accuracy. In view of the good overall performance of the TMTPSS and TM functional for diverse systems and a wide class of properties, we conclude that TMTPSS and TM are attractive as efficient universal density functionals. An obvious advantage is that we are able to use the same method and the same basis set to simultaneously describe different class of problems in chemistry, physics, and materials science. A limitation of TMTPSS and TM is that, like other semilocal DFT methods, they are unable to describe properties, in which nonlocality is significantly important, such as charge transfer and reaction barrier height.

See supplementary material for all calculated properties of individual complex.

G.T. was supported by the National Natural Science Foundation of China under Grant No. 51264021, by Back-up Personnel Foundation of Academic and Technology Leaders of Yunnan Province under Grant No. 2011HR013, and China Scholarship Council. Y.M. acknowledges support from the NSF under Grant No. CHE−1640584. J.T. acknowledges support from the NSF under Grant No. DMR−1607868.