Benchmark of Simplified Time‐Dependent Density Functional Theory for UV–Vis Spectral Properties of Porphyrinoids

Time‐dependent density functional theory is thoroughly benchmarked for the predictive calculation of UV–vis spectra of porphyrin derivatives. With the aim to provide an approach that is computationally feasible for large‐scale applications such as biological systems or molecular framework materials, albeit performing with high accuracy for the Q‐bands, the results given by various computational protocols, including basis sets, density‐functionals (including gradient corrected local functionals, hybrids, double hybrids and range‐separated functionals), and various variants of time‐dependent density functional theory, including the simplified Tamm–Dancoff approximation, are compared. An excellent choice for these calculations is the range‐separated functional CAM‐B3LYP in combination with the simplified Tamm–Dancoff approximation and a basis set of double‐ζ quality def2‐SVP (mean absolute error [MAE] of ≈0.05 eV). This is not surpassed by more expensive approaches, not even by double hybrid functionals, and solely systematic excitation energy scaling slightly improves the results (MAE ≈0.04 eV).

www.advancedsciencenews.com www.advtheorysimul.com computationally quite expensive for tackling a conjugated molecular system beyond the basic PP. Other prominent approaches, such as symmetry adapted cluster-configuration interaction [17] (SAC-CI) and similarity transformed equation-of-motion coupled-cluster [18] (STEOM-CC) are more accurate than CASPT2, but limited to molecular systems up to 50 atoms due to their high computational cost. To overcome the limits of CASPT2, secondorder N-electron valence state perturbation theory [19,20] (NEVPT2) has been introduced, which is more efficient than CASPT2, size consistent, and intruder-state-free, but like all multi-reference approaches, the computational cost of NEVPT2 is still high for larger systems. Overall, previously mentioned approaches yield reliable absorption energies for PPs, but suffer from a high computational cost for increased molecular systems, making them practically unsuitable for large molecules and materials.
To find an alternative for the prediction of excited state properties of large systems at a moderate cost, time-dependent density functional theory (TD-DFT) appears to be a promising candidate. TD-DFT is an extension of Kohn-Sham (KS) DFT and based on the almost 35-year-old Runge-Gross theorem, [21] which has been thoroughly reviewed in the literature. [21][22][23][24][25][26][27][28] Almost 25 years ago, Casida developed a constructive linear-response formalism for TD-DFT, known as Casida equations [23] but which we will refer to as the random-phase approximations (RPA), allowing to efficiently determine the solution of the TD-DFT equations, which are formulated in matrix equation involving the excitation and de-excitation matrices.
A popular approximation to the Casida equations is the Tamm-Dancoff approximation [29] (TDA), which simplifies the algebra and associated algorithms to obtain the electronic excitations, yet it typically yields electronic excitations close to those obtained by TD-DFT. [30,31] Unfortunately, TD-DFT on the basis of popular exchange correlation functionals suffers accuracy limitations which are most evident in the failure to correctly describe Rydberg and charge transfer (CT) states. [32,33] These drawbacks can usually be overcome by range-separated hybrid (RSH) functionals, which employ a large amount of Hartree-Fock (HF) exchange at large electron-electron distances and, therefore, reflect the correct asymptotic exchange potential.
TD-DFT in numerous variants has been applied to PPs, and in the following, we are summarizing the state of the art as found in the current literature.
In 1996, Bauernschmitt et al. [34] employed TD-DFT to compute the first four electronic excitations of PP to validate exchange correlation functionals, including the local density approximation (LDA: S-VWN), the generalized gradient approximation (GGA: BP86), and hybrid functional (B3LYP). Their results showed that TD-DFT excitation using the BP86 functional are in better accordance with experiment than configuration interaction singles (CIS) and TD-HF. Also, CASPT2 possesses an error of more than 0.3 eV for the Q-bands compared to experimental results. However, the employed basis set was overall small for post-HF approaches.
In 2010, Tian et al. [35] examined the performance of global hybrids (PBE0, B3LYP, M06, M06-2X, M06HF) and long-range corrected (LC) hybrid functionals ( B97X-D, B97X, B97, LC-PBE, and CAM-B3LYP) in TD-DFT calculations to predict the spectral properties of PP analogues. Among the many functionals tested, the LC functional B97X-D results in an error of 0.05 eV for Q y band. Moreover, they concluded that the results are robust with respect to subtle geometry changes resulting from the functional choice for geometry optimization and showed that diffuse functions have only a minor effect on calculated absorption spectra. However, the quite general study included only two porphyrinoids.
Eriksson et al. [36] (2011) investigated the ability of LC hybrid functionals B97, B97X and B97X-D within the TD-DFT framework. They found that B97X reproduces the experiment best with an error of up to 0.09 eV. Additionally, it was confirmed that the applied functional for geometry optimization has only a small influence on the calculated spectra.
Lee et al. [37] (2012) benchmarked five DFT functionals (B3LYP, LC-PBE, LC-BLYP, CAM-B3LYP, and B97X-D) using TD-DFT for PP derivative. It was found that B97X-D yields the best agreement to the reference for the LC functionals (Q ave bands: 0.055 eV). Overall, better results were obtained for B3LYP for the Soret and Q-bands. However, it was not recommended due to the susceptibility for CT excitations.
Fang et al. [39] (2014) compiled a subset of 96 excitations of 79 different organic and inorganic molecules, including basic PP. [40] They have assessed diverse DFT functionals (BP86, B3LYP, PBE0, M06-2X, M06-HF, CAM-B3LYP, and B97XD) and two wave-function based approaches (CIS and CC2). Overall, the lowest error was produced by CC2 with MAE of 0.19 eV. However, it was found that CC2 approach did not perform well for inorganic systems (MAE: 0.31 eV), while the MAE of B3LYP is only 0.22 eV.
Despite the many successful applications of TD-DFT on a wide range of molecular systems, it is often challenging in TD-DFT to calculate a sufficient number of excited states for a complete spectrum or spectra of extended biological systems. To overcome this challenge, the Grimme group presented two highly efficient approaches, the simplified Tamm-Dancoff approximation [42] (sTDA) and simplified time-dependent density functional theory approach [43] (sTD-DFT). In both approaches, the computational resources needed to tackle a targeted system is solely determined by the ground state DFT calculation. This is achieved by approximating Coulomb and exchange interactions of the electrons by monopole interactions. Additionally, the CI space is truncated with a screening based on second-order perturbation theory. Note, the central concepts of the sTDA approach to increase the computational efficiency can be also employed in tight-binding approaches and, thus, allows fast access to excited state properties of systems with an amazing size. [44] Computational studies validating sTDA or sTD-DFT for PPs are missing in the literature so far.
With the goal to identify a computational feasible approach to investigate the absorption properties of PP-containing materials and extended biological systems, we compare the semi-empirical sTDA, sTD-DFT, and canonical TD-DFT (RPA and TDA) for UVvis spectra calculations of porphyrinoids. After a short summary of computational details, we assess diverse DFT functionals regarding their performance for the calculations of absorption energies of the Q-bands. This includes the (for sTDA unsuccessful and unnecessary) attempt to improve the results by an empirical scaling of excitation energies that can improve results significantly.

Computational Methods
All geometries have been fully optimized using the Turbomolesuite, [45] employing the BLYP functional with Grimme's D3 correction for London dispersion (BLYP-D3) [46][47][48] in combination with the resolution of identity (RI) approximation, [49][50][51] and the TZVP [52] split-valence basis set of triple-quality with polarization functions. The convergence criterion for the selfconsistent field approach was increased to 10 −8 Hartree. This approach is very fast for the molecules studied here, but also easily affordable in periodic calculations that suffer severe performance loss for hybrid functionals. As hybrid functionals are known to produce more accurate structures, we assessed the influence of the geometry on the excited state properties. For that purpose, we reoptimized all geometries using the B3LYP-D3 hybrid functional, [46][47][48]53] again employing the TZVP basis set and the RIJK.
To speed up the calculations, we employed the RI approximation throughout, including its variant for double hybrid functionals, [72,73] and the RIJCOSX approximation [74] was employed for the global hybrid and RSH functionals. For comparison, we have also included the post-HF methods CIS [75] and CIS(D). [76,77] The performance of each approach was assessed by calculation of the mean error (ME), the mean absolute error (MAE), and the absolute maximum error (MAXE) to the experimental reference values.

Benchmark Set
We have included diverse variants of porphyrinoids starting from basic PP to the extension of conjugated -system of the central core followed by ring functionalization and modification of metal atoms. The molecules included in our benchmark set are given in Figure 1, while Table 1 lists the experimental references.

Results and Discussion
In this section, the performance of diverse functionals is presented. We will start the validation with the computationally least costly DFT approach, the GGA functionals, and will finish with the most expensive one, the double hybrid functionals. We would like to add here that the present TD-DFT studies involve only transitions in the frozen ground state of the molecules. Therefore, only the 0 → 0 transitions of the Q-bands can be obtained from the calculations. We will focus on the Q-bands in the following because they absorb in the visible light range. Furthermore, these transitions can be clearly distinguished from other transitions in the excited state calculations.

GGA and meta-GGA Functionals
GGA and meta-GGA functionals do not require four-center integrals as the Coulomb interaction can be calculated directly via the electron density, which is particularly beneficial for periodic calculations and for codes employing different basis functions than Gaussian-type orbitals. On the other hand, pure KS DFT is very prone to errors originating from the self-interaction error (SIE). Thus, an overall poor general performance for excited state calculations can be expected due to weakly bound electrons.
As can be seen in Table 2 and Figure 2a, the performance of TD-DFT for the GGA functionals PBE, BP86, and BLYP is very similar. Compound 8 shows a large error resulting in an outlier of about 0.5 eV. This can be attributed to a significant contribution of CT excitations to the Q-bands (see Table S2 and Figure  S2a,b, Supporting Information, for a detailed description). Employing the meta-GGA functionals TPSS and M06-L reduces significantly the MAXE. However, the MAE (denoted by "+" sign in Figure 2a) is not improved and still exceeds 0.12 eV. Furthermore, meta-GGAs tend stronger to over-estimate the absorption energies compared to GGAs. Employing the RPA-approach does not result in significant improvements in comparison to the TDA-approach (see Table 2 and Figure 2b). For instance, the calculated MAE from both approaches in combination with the GGA-functionals is nearly similar. Only for M06-L, the MAE is reduced by about 0.03 eV. Additionally, the RPA-approach tends to lower absorption energies than the TDA-approach.
The computationally cheaper approaches, sTDA and sTD-DFT, show a comparable MAXE for GGAs (see Table 2 and Figure  3a,b). However, the MAE of the GGA-functionals is strongly affected by the selected approach and basis set. For example, the ME and MAE of sTDA with the functional BLYP increases by about 0.1 eV when the larger def2-TZVP basis set is employed (see Table 2 and Figure 4a). Changing from sTDA to sTD-DFT results also in large MAE values. The ME and MAE of a given approach basis-set combination is overall close to each other in all cases and, thus, highlighting systematic deviations. Please note that the sTDA approach includes empirical parameters optimized for functionals with a HF exchange contribution between 20% and 60% and not for pure DFT functionals. [42] TPSS shows an improvement to the GGAs but is still significantly affected by the selected basis set. The MAE of M06-L for sTDA and sTD-DFT is only slightly affected by the choice of basis set which is in contrast to the GGAs. To sum up, the MAE of GGAs and meta-GGAs exceeds 0.08 eV, while the MAXE is reduced to 0.30 eV only for the M06-L functional. Finally, a systematic underestimation of absorption energies, especially for the sTDA and sTD-DFT in combination with a large basis set, is observed suggesting a global scaling of the obtained energies to match better the experimental reference. This approach is rather semi empirical but allows to access excited state calculations of extended system sizes due to the low cost of pure KS DFT.
After scaling of energies, significant improvements of the MAE and MAXE are only obtained for the sTDA and sTD-DFT approaches in combination with the large basis set (see Figure 4b for BLYP). All functionals with scaled error value are listed in Table S3, Supporting Information, and a graphical illustration can be seen in Figure S5, Supporting Information. Nonetheless, the MAE is still above 0.05 eV and high MAXE is obtained as well. Thus, GGAs and meta-GGAs cannot be recommended for calculations of UV-Vis spectra of porphyrinoids.

Global Hybrid Functionals
A hybrid functional is defined as an approximate KS density functional where a part or all the semi-local DFT exchange expression E DFT X is replaced by exact HF exchange E HF X . The amount of HF exchange for typical hybrid functionals lies in the range of 10-25% but can be as high as 50-55% like in the BHLYP and the M06-2X functional. Increasing the amount of HF exchange also increases the likelihood of encountering triplet instabilities (i.e., imaginary triplet excitation energies which indicate that the ground state is unstable with respect to symmetry breaking). Moreover, incorporation of exact HF exchange reduces the SIE, which is a significant troublemaker in TD-DFT. However, it also introduces a four-index integral into the Hamiltonian which leads to higher computational cost in comparison to the GGA and meta-GGA functionals.
As shown in Tables 2 and 3, similar trends are observed for global hybrids as for pure DFT functionals; for example, RPA tends to lower absorption energies than TDA. In contrast, the large error which stems from the outliers are significantly reduced for the global hybrids with respect to GGAs (see Table 3 and Figure 2a,b). TDA in combination with global hybrids tends strongly to overestimate absorption energies of Q-bands which can be reduced by employing the RPA approach, especially in combination with the BHLYP and M06-2X functional. Nonetheless, MAE is still above 0.10 eV and thus, this approach cannot be recommended.
The approximate sTDA and sTD-DFT approaches fail significantly for the global hybrids with a large amount of HF exchange like BHLYP and M06-2X. However, global hybrids with HF-exchange contributions in the range of 20-25% show significant improvement compared to the pure DFT functionals (see Table 3 and Figure 3a,b). For example, sTDA in combination with the B3LYP functional and def2-SVP basis set possesses an MAE of 0.06 eV, while the MAXE is 0.18 eV. Increasing the basis set increases the deviation, similarly as for the GGA functionals (see Figure 4a for B3LYP and Figure S5, Supporting Information, for all the tested functionals).
The MAE and MAXE can be reduced by energy scaling due to systematic deviations. A global hybrid with large HF-exchange contribution works best: sTDA/sTD-DFT in combination with the functionals BHLYP or M06-2X and the def2-SVP basis set possess MAEs of only 0.06 eV (see Table S4 and Figure S5, Supporting Information). Moreover, MAXE is reduced to 0.12 eV   for M06-2X functional in combination with the RPA-approach and def2-SVP basis set. Therefore, employing global hybrids with large HF-exchange contribution and energy scaling might be a suitable, albeit somewhat empirical approach to estimate the absorption energies of porphyrinoids.

Range Separated Hybrid Functionals
RSH functionals posses a different contribution of HF exchange in short and long interelectronic distances. Short range corrected functionals, like HSE06, possess a medium amount of HF exchange in the short range while it drops commonly to zero at long interelectronic distances. This allows a faster calculation of solid-state properties compared to global hybrid functionals, but an improvement for excited states cannot be expected for this type of functional. In contrast, LC-RSH exhibits a large amount of HF exchange at long interelectronic distances. This significantly reduces errors originating from Rydberg states and CT excitations in TD-DFT calculations. [66,[90][91][92] Therefore, we will focus solely on LC-RSH. The performance of TDA for the investigated RSH functionals is overall comparable to that of global hybrids, as visible in Tables 3 and 4 and Figure 2a. Additionally, TDA tends stronger to overestimate the absorption energies. The RPA approach for CAM-B3LYP results in overestimated absorption energies, while other LC-RSH functionals tend to underestimate these energies. Overall, RPA has a strong dependency on the type of functional and produces lower absorption energies than TDA (see Table 4 and Figure 2b).  On the other hand, as can be seen in Table 4 and Figure 3a,b, the performance of sTDA and sTD-DFT is improved over global hybrids only for CAM-B3LYP, which amounts up to 46% HFexchange. In contrast to CAM-B3LYP, LC-RSH functionals incorporating up to 85-100% HF-exchange result in remarkably large errors. This is in agreement with the original work where sTDA was developed for LC hybrids where best results of excitation energies were obtained with CAM-B3LYP for CT-free systems. [93] Employing large basis sets does not improve the re-sults based on simplified approaches (see Figure 4a for CAM-B3LYP case and Figure S5, Supporting Information, for all the tested functionals). Thus, results obtained with the def2-SVP in combination with simplified approaches are most reliable and reasonable.
To sum up and as shown in Figure 5a, GGA functionals produce large errors in the form of outliers, for example, BLYP in combination with def2-SVP basis set. Global hybrid functionals, for example, B3LYP, produces MAE of 0.06 eV in combination  with sTDA and the def2-SVP basis set and appears to produce overall reliable results. However, best results, and indeed excellent ones, are obtained with the RSH, CAM-B3LYP in combination with sTDA and an overall small def2-SVP basis set which yields an MAE of about 0.05 eV. This can be barely improved by energy scaling (see Figure 5b and scaled error values are listed in Table S5, Supporting Information).

Double Hybrid Functionals and Post-Hartree-Fock Approaches
In addition to the exact HF exchange, double hybrid functionals include a second-order perturbation theory correction term (MP2) for the correlation part of the functional. This improves mainly the consideration of dispersion forces. However, the com-putational time is comparable to MP2. Therefore, we have also included some traditional post-HF approaches with comparable computational cost, CIS, and CIS(D) in our study. As can be seen in Table 5 as well as in Figure 6a, double hybrid functionals produce large errors comparable to CIS, including perturbative double corrections results in even larger errors of the CI approach. However, we would like to highlight that the employed basis set is only of double-quality due to the system size. The strong systematic overestimation of absorption energies for the double hybrid functionals and post-HF methods suggests a scaling of the obtained absorption energies. Indeed, the results are significantly improved and an accuracy comparable to CAM-B3LYP can be reached (see Figure 6b and Table S6, Supporting Information). Thus, CIS(D) with scaled absorption energies might be a suitable approach to verify results obtained with sTDA, def2-SVP, and CAM-B3LYP.

Influence of Diffuse Basis Set Functions and Ground State Structure
Commonly, diffuse basis sets are recommended for weakly bound electrons found in anions or in excited states. Therefore, we have selected three functionals, BLYP, B3LYP, and CAM-B3LYP, and extended the employed Ahlrichs basis set by diffuse functions. Independent of the employed functional type, including diffuse basis set functions provides poorer results compared to the def2-SVP double-basis set, (see Figure S7a and Table  S7, Supporting Information). The worse performance cannot be explained by the ϵ (HOMO) criterion [94] (see Table S8, Supporting Information). Also scaling of energy does not improve results, since including diffuse basis sets increases the scattering of the calculated absorption energies in most cases (see Figure  S8a, Supporting Information). Thus, unintuitively, the smallest basis set provides the most accurate results.
Finally, we investigated the influence of the electronic structure method during structure optimization. Instead of the GGA BLYP, the more expensive hybrid functional B3LYP was selected for structure optimization. The influence is overall negligible for absorption energies obtained by BLYP and B3LYP (compare Tables S7 and S9, Supporting Information). In the case of the RSH CAM-B3LYP, the errors without energy scaling are even increased pointing to some error compensation for the most reliable approach (see Figure S7b, Supporting Information). Nonetheless, global scaling of energy provides nearly identical results. Thus, as long as the correct combination of scaling factor, structure optimization setup, and absorption energy calculation approach are selected, results can be barely improved.

Conclusions
We have presented a detailed validation of the simplified timedependent DFT method developed by Grimme et al. for the calculation of UV-vis spectra of porphyrinoids, including free base and metal containing PPs. The original RPA approach tends to smaller absorption energies than TDA, which is also visible for the simplified versions. Local GGA functionals produce large errors and therefore cannot be recommended. In contrast to local DFT functionals, global hybrids yield significantly improved results only after energy scaling, especially BHLYP and M06-2X. We can recommend as global hybrid B3LYP, which produces MAE of 0.06 eV in combination with sTDA, and the def2-SVP basis set, which can be barely improved by energy scaling. The best results without energy scaling are obtained with the RSH CAM-B3LYP in combination with sTDA and the def2-SVP basis set yielding an MAE of about 0.05 eV. Significantly more expensive perturbative corrected double hybrid functionals tend to yield results comparable to CAM-B3LYP solely when energies are scaled. Apart from that, employing a hybrid instead of a GGA functional for geometry optimization has less significant effect on the calculated absorption bands, whereas increasing the basis set does not improve the calculated absorption bands. Most notable, including diffuse basis functions even leads to worse results. Thus, employing a cheap GGA like BLYP for structure optimization, selecting an overall small basis set of double-quality in combination with a CAM-B3LYP, and the sTDA approach provides a cost-efficient approach to estimate the absorption spectra of porphyrinoids which can be barely improved by more expensive approaches. Unfortunately, none of the local functionals has sufficient predictive power, which is an obstacle in particular for periodic calculations.

Supporting Information
Supporting Information is available from the Wiley Online Library or from the author.