Meta-generalized gradient approximations in time dependent generalized Kohn–Sham theory: Importance of the current density correction Open Access

Time-dependent density functional theory (TDDFT) is frequently used for calculating electronic excitations. The accuracy of a TDDFT calculation hinges on the approximation that is employed for the exchange–correlation (xc) functional. The local density approximation (LDA) and Generalized Gradient Approximations (GGAs) come at a low computational cost but suffer from systematic shortcomings that limit their usefulness. Functionals such as range-separated hybrids can be much more accurate, rivaling many-body perturbation theory and wavefunction methods in accuracy,^1,2 yet they also come at a much higher numerical cost and conceptual complexity than semilocal functionals. In the context of ground-state density functional theory (DFT), meta-GGAs that depend on the non-interacting kinetic energy density have emerged as an attractive class of xc approximations with an accuracy better than GGAs,^3–18 while still being semilocal in terms of the computational expense. Recently, it has been shown that a non-empirical construction strategy can give meta-GGAs that show (ultra-)nonlocality and can describe static long-range charge-transfer in great similarity to Fock exchange,¹⁷ while also yielding reasonable results for many ground-state properties.¹⁹ This raises the hope that using meta-GGAs in TDDFT may allow us to overcome the systematic, qualitative failures of LDA and GGAs, without having to face the substantial increase in computational cost that is associated with fully non-local functionals.

Using a meta-GGA in TDDFT, however, is less straightforward than using explicit density functionals because the dependence on the kinetic energy density introduces a gauge dependence.^20,21 In this paper, we revisit this issue and address it for two modern meta-GGAs, the r²SCAN¹⁸ and TASK¹⁷ functionals. In doing so, we rely on the real-time propagation formalism, which has been developed in different technical realizations in the past years.^22–36 The real-time approach is finding increasing interest because, on the one hand, it allows us to go beyond the linear response regime and, on the other hand, it scales well with the system size³⁷ and can reach outstanding computational efficiency.

In Sec. II, we summarize the basic facts about meta-GGAs from the ground-state perspective, motivating why we focus on the r²SCAN and TASK functionals and discussing the technically relevant issue of using pseudopotentials in meta-GGA calculations. In Sec. III, we review the use of meta-GGAs in TDDFT and discuss the gauge variance problem and the current density correction that has been suggested to address it. Sections IV and V are devoted to the discussion of TDDFT meta-GGA results obtained with TASK and r²SCAN, respectively. We find that the current density correction can sometimes be small but can also have a significant impact and, therefore, should generally be taken into account when meta-GGAs are used in TDDFT. We close this paper with a summary and conclusions.

In our work, we focus on two relatively recently developed meta-GGAs. One is the r²SCAN functional,¹⁸ which was derived from first-principles as a broadly applicable xc approximation that describes many different types of electronic bonding well, while defusing the numerical intricacies of the original SCAN construction¹⁴ to some extent. SCAN and r²SCAN are considered to be among the most accurate, non-empirical general-purpose meta-GGAs. However, here we can only report TDDFT results for r²SCAN and not SCAN because we could not reach a numerically long-term stable propagation with SCAN. This latter finding is no surprise (also see the comments in the supplementary material) in view of earlier reports about the numerical challenges associated with SCAN.^18,38–41

The second functional that we study is the TASK meta-GGA,¹⁷ which is also derived non-empirically from first-principles, fulfilling the same constraints on the exchange energy as SCAN, but with a focus on incorporating a realistic derivative discontinuity into the energy functional. This leads to a reliable prediction of bandgaps^17,41,42 and, with a suitable correlation functional, also a reasonable accuracy for atomization energies and barrier heights.¹⁹ 

As for any orbital-dependent density functional,⁴³ one has the choice of using mGGAs either in the Kohn–Sham scheme, i.e., with the potential calculated as the functional derivative with respect to the density in the optimized effective potential formalism, or in the generalized Kohn–Sham (GKS) scheme, i.e., calculating functional derivatives with respect to the orbitals. In TDDFT, meta-GGAs have been explored within the Kohn–Sham scheme, e.g., in the context of linear-response calculations that used the Sternheimer approach.⁴⁴ Here, however, we follow Refs. 21, 45, and 46 and focus on the GKS approach to avoid the serious difficulties^26,47 that are associated with the solution of the time-dependent optimized effective potential equation, which could so far only be solved for one-dimensional model systems.^48,49

In plane wave and real space calculations, the external potential in Eq. (12) is typically a pseudopotential that represents the interaction between the atomic cores and the valence electrons. In our calculations, we rely on the real-space code BTDFT,^31,50 which uses norm-conserving pseudopotentials of Troullier–Martins type.⁵¹ For orbital functionals, constructing consistent pseudopotentials is considerably more demanding than for semilocal functionals, such as LDA or GGAs, as demonstrated for exact exchange.^43,52 First steps toward constructing pseudopotentials for meta-GGAs have been taken,⁵³ but including the differential operator part of the GKS potential, i.e., Eq. (11), into the pseudopotential construction is challenging, as the usual inversion from the pseudo-orbitals to the screened pseudopotential maps to a local potential. Thus, consistent pseudopotentials for r²SCAN and TASK are, unfortunately, not available yet. Using a non-consistent pseudopotential is always unsatisfactory from a formal point of view. In the context of exact exchange calculations, it has been demonstrated that the practical consequences of using a non-consistent pseudopotential need not necessarily be large.⁵⁴ However, in the context of meta-GGAs, special attention should, nevertheless, be devoted to the pseudopotential question because α of Eq. (3) (or the ratio τ^W/τ) is typically used to detect and adjust the meta-GGA to iso-orbital regions. A simple example demonstrates that the pseudopotential can have a significant influence on the iso-orbital character: In a sodium atom described by a usual LDA or GGA pseudopotential, which only treats the 3s electron as a valence electron, all of the valence density is iso-orbital and the meta-GGA treats xc effects in the Na atom as in the hydrogen atom. Meanwhile, in the all-electron calculation, the core electrons are present and the Na atom density is mostly not of iso-orbital character; thus, the meta-GGA does not switch to its iso-orbital limit. Consequently, a meta-GGA calculation for sodium based on a GGA pseudopotential may lead to very different results than a meta-GGA calculation based on all electrons. Thus, depending on the system that one studies, meta-GGAs can have a larger sensitivity to pseudopotential inconsistency than, for example, GGAs.

For the calculations presented in Secs. IV and V, we can, therefore, be sure that the results for molecules that only consist of light atoms, such as N, C, O, and Li, e.g., CO and Li₂, are not at all affected by the pseudopotential iso-orbital problem. For molecules also containing atoms with more than just 1s core electrons, pseudopotentials with the core correction of Eq. (15) are used for the atoms with only a 1s core, whereas standard pseudopotentials are used for the heavier atoms. The comparison with all-electron calculations has shown that this pragmatic strategy leads to overall reliable results for the systems that we study in this paper. For general future work on meta-GGAs, consistent pseudopotentials are, of course, very desirable.

Due to the inherently local Kohn–Sham potential, TDDFT calculations using meta-GGAs that were performed in the Kohn–Sham instead of the GKS approach⁴⁴ did not need to worry about the challenges that arise from the operator-potential Eq. (11). They were, however, limited to the Krieger–Li–Iafrate approximation. TDDFT calculations for meta-GGAs in the GKS scheme with the operator potential and including the current density correction have been performed by Bates et al.^21,60 for the Tao–Perdew–Staroverov–Scuseria⁹ functional. They found the effect of the current density correction to be relatively small. Other TD-GKS meta-GGA calculations did not take the correction into account.⁶¹ As it is clear from Eqs. (28) and (31) that a more pronounced τ-dependence will lead to more pronounced effects, the development of meta-GGAs, such as TASK, which have a strong τ-dependence, and to some extent also r²SCAN, calls for a re-evaluation of the importance of the current density correction. In the following, we, therefore, present a systematic comparison of excitation spectra calculated with and without the use of the current density correction. We here focus on just a few model systems for which accurate calculations can be performed and for which accurate reference data are available. The chosen systems, though few, nevertheless, span a range of paradigmatic electronic structures, from few electron diatomic molecules to conjugated aromatic molecules and organic semiconductors.

We first calculate the absorption spectra of the diatomic molecules Li₂ and CO, because for these small systems, all numerical parameters can be chosen to yield high accuracy and accurate reference results are available. Our calculations focus on TASKx, where “x” denotes that we are using exchange only without a correlation functional, because in this way, we can focus on the effect of the pronounced τ-dependence that is part of the construction strategy of the TASK meta-GGA for exchange. However, we also do an exemplary check of the effect of correlation. For ground-state energetics, reasonable accuracy is achieved by combining TASKx with the cc-functional of Ref. 19, i.e., a self-interaction corrected LDA correlation. However, for the closed shell systems that we study here, the cc energy correction vanishes, reducing to LDAc, and therefore, we use TASKx+LDAc in our TDDFT calculation.

Following Ref. 62, we choose a bond length of 2.1421a₀ for carbon monoxide (CO). In a first set of calculations, we check the relative accuracy of our real-time real-space calculations using a pseudopotential with a core correction as described in Sec. II B in comparison with explicitly linearized all-electron TDDFT with different Gaussian basis sets as implemented in TURBOMOLE.⁶³ The real-time real-space calculations use the program package BTDFT and the propagation evaluation as described in Ref. 31. Concerning TURBOMOLE, we note that discrepancies that we observed with earlier versions of TURBOMOLE disappeared after the recently made update.⁶⁰ As a guide to the eye, in the following, we plot schematic excitation spectra with a Gaussian broadening $exp(−ℏ(ω−ω0)/η)2$⁠, where η = 0.025 and ℏω₀ is the excitation energy. Tables listing all the excitation energies that are discussed here and in the following are given in the supplementary material.

The upper panel of Fig. 1 shows the schematic spectra obtained for CO with TASKx calculated with BTDFT and TURBOMOLE, in straightforward GKS TD calculations, i.e., without the current density correction. In BTDFT, we used a grid spacing of 0.2a₀ and a spherical grid of radius 15a₀. Using a bigger grid of 22a₀ changes the excitation energy by less than 0.01 eV. In TURBOMOLE, the results are noticeably influenced by the basis set, and with increasing size of the basis set, the excitation energy systematically shifts toward the value that we obtained with BTDFT. The high sensitivity toward the basis set and requirement for a very large basis is in line with earlier findings reported in the context of Kohn–Sham linear-response calculations for CO.⁶⁴ We, therefore, strongly assume that the remaining small discrepancy between the TURBOMOLE and BTDFT spectra is due to the different numerical representations and would vanish in calculations that could reach the numerical limits. The lower panel of Fig. 1, which shows the excitation energy for TASKx combined with LDA correlation, confirms these findings. The LDA correlation shifts the excitation energy up by ∼0.1 eV, and the basis-set dependence in TURBOMOLE is a bit less pronounced, but qualitatively, the results are the same. Therefore, we do not discuss the combination with LDA correlation for other systems.

As we have thus established an estimate of the numerical accuracy for the excitation energies, we investigate in our next step the influence of the current density correction. Figure 2 shows the comparison of the results from our real-time real-space calculations without (dashed line) and with (solid line) the current density correction for TASKx (upper panel) and TASKx+LDAc (lower panel). The difference between GKS and CGKS is considerable, as the current density correction lowers the excitation energy by about 0.46 eV in both cases, shifting the excitation to 8.54 eV (TASKx) and 8.64 eV (TASKx+LDAc).⁶⁵ This is an interesting observation in view of the experimental value of 8.51 eV for the first excitation energy in CO⁶⁶ and the third-order response coupled cluster (CC3) results of Ref. 62 that find the prominent excitation line at 8.47 eV (calculated with Dunning’s basis set aug-cc-pVQZ⁶⁷): The current density correction obviously brings the TD-meta-GAA results in much better agreement with the reference data.

The second small molecule that we study as a test case is Li₂. For our calculations, we adopt the experimental bond length of 5.051a₀ from Ref. 68. With a grid-spacing of 0.2a₀ and a grid radius of 15a₀, our real-time real-space calculations show two prominent excitation lines at 2.63 and 3.33 eV for TASKx-GKS, i.e., without the current density correction, as depicted in Fig. 3. In the TURBOMOLE calculations, we again observe a noticeable basis set dependence, in particular, for the second excitation. Whereas the TZVPP result deviates non-negligibly from the real-space real-time result for the second excitation, the QZVPP⁶⁹ and QZVPPD results are close to each other and to the real-time results. The remaining small discrepancy of 0.09 eV is in the range of what is to be expected on numerical grounds for such different techniques. With the aug-cc-pVQZ basis set, the second excitation shifts further in energy, but the calculation appears sensitive to numerical details and a previously unseen third excitation appears. We, therefore, regard the result with the QZVPPD basis set as the most reliable TURBOMOLE spectrum.

With the numerical trust range thus established, we again check the impact that the current density correction has on the spectra. Figure 4 shows the comparison of the spectra obtained with BTDFT for TASKx without (dashed line) and with (solid line) the current density correction. Again, we observe a substantial difference. The two prominent excitations are at 2.63 and 3.33 eV for TASKx-GKS and at 1.90 and 2.46 eV for TASKx-CGKS. Thus, the excitations shift by 0.73 and 0.87 eV, respectively.

Once more it is instructive to compare these values to reference data. We used the experimental and coupled cluster (CC3) results from Ref. 68 to calculate the respective values for the vertical excitation energies with the help of the harmonic approximation. This results in 1.86 and 2.57 eV for the experimental values and 1.83 and 2.57 eV for the CC3 values. Therefore, taking the current density correction into account reduces the relative error significantly from 40.9% to 2.2% for the first excitation and from 29.6% to 4.3% for the second excitation, compared to the experimental data. Thus, the results for Li₂ confirm the previous findings for CO: The current density correction can have a large impact on the excitation energies, and the TASKx-CGKS results are in good agreement with the experimental and theoretical reference values.

We continue our study with a few paradigm organic molecules that have a conjugated π-electron system and thus an electronic structure that shows other features than in diatomic molecules. Our first test case is benzene (C₆H₆), with Fig. 5 showing the schematic spectrum that we obtain with the real-space real-time approach without and with the current density correction. Again, the excitation shifts to a lower energy when going from GKS to CGKS. The magnitude of the shift, however, is with 0.15 eV smaller than what we observed for CO and Li₂. A possible explanation for this observation is seen in Eq. (28): It shows that the derivative $∂excmGGA/∂τ$ is decisive for how strongly the continuity equation is violated, and this is also an indication for the strength of the current density correction that one expects. The importance of this derivative has also been pointed out in the context of meta-GGA calculations of the optical response of semiconductors.⁷⁰ It is one of the construction principles of TASK to generally have this derivative larger than previous meta-GGAs to increase the (ultra-)nonlocality of the meta-GGA exchange. However, the magnitude of this derivative for a given system depends on which values of n, s, and α are realized by the system’s electronic structure and which parts of the enhancement factor are thus probed by the given system.⁴¹ Organic molecules, such as benzene, that have delocalized electrons are more “electron gas like” than a diatomic molecule. Consequently, one can expect that such molecules probe the enhancement factor closer to the electron-gas limit. Since both TASK and r²SCAN are designed to respect the electron gas limit, one expects a reduced derivative and thus a smaller impact of the current density correction, in agreement with the numerical observation.

We further note that in contrast to the energetic position of the excitations, for which we always observe a downshift, there does not seem to be a similarly systematic trend in the changes in the oscillator strengths of the main excitations: Whereas for CO and Li₂, the oscillator strength of the dominant excitation decreased when going from GKS to CGKS, it slightly increases for benzene.

As a test case of a certain practical relevance, we next look at a model organic semiconductor molecule (NDI-1). In NDI-1, a naphthalene diimide core serves as an electron acceptor and two thiophene molecules are attached to it, one to the left and the other to the right, serving as electron donors. This system was designed with the aim of being used in organic solar cells, and it has been used in several TDDFT studies as a test case,^44,71 because its lowest excitation is of charge-transfer character, with electron density being transferred from the thiophenes to the NDI core. This is the type of charge-transfer excitation that typically is of relevance in organic electronics, with a spatial separation of electron and hole densities, yet carrying some oscillator strength.

The NDI-1 molecule contains sulfur (S), which requires a pseudopotential with more than just 1s electrons in the core. Consequently, the core-correction described in Sec. II B cannot be used for S and we fall back to a regular Troullier–Martins pseudopotential for this atom (see the supplementary material for details), while using the core correction for all first-row atoms. To check whether this leads to inaccuracies, we calculated the GKS spectrum of NDI-1 also with all electrons using TURBOMOLE (spectra shown in the supplementary material). The two calculations agree closely. We are thus confident that for NDI-1, the pseudopotential inconsistency for the S atoms does not deteriorate the quality of the results.

Figure 6 shows the comparison of the result from the real-time real-space GKS and CKGS calculations. The first striking observation is that, now, the spectra with and without the current density correction are very similar. For TASKx-GKS, the first line appears at 2.05 eV, and for TASKx-CGKS, the first line appears at 2.03 eV, i.e., the difference is minute. Furthermore, it is instructive to compare the energy of the first peak, which corresponds to the above-mentioned charge-transfer excitation, with the energies that one obtains for this excitation with other xc approximations: With TD-LDA (and similarly with typical GGAs), this excitation is found at 1.69 eV, whereas an optimally tuned range-separated hybrid places it at 2.52 eV, close to the experimental result.⁷¹ Thus, this excitation paradigmatically shows the well-known problem that usual (semi-)local approximations severely underestimate charge-transfer excitation energies.⁷² In light of this comparison, the findings for TD-TASK appear as a mixed result: On the one hand, there is a clear improvement compared to TD-LDA, as about half of the gap between TD-LDA and the optimally tuned result is closed. This is plausible in view of the construction strategy of the TASK functional that stresses ultra-nonlocality and leads to field-counteracting terms.¹⁷ On the other hand, we see that the TD-TASK result is still noticeably lower than the optimally tuned reference, which shows that the TDDFT charge-transfer problem is not fully solved by TD-TASK.

These conclusions are confirmed by our last test, in which we calculate the excitations of two aggregated bacteriochlorophyll (BChl) chromophores, BChl-302 and BChl-303, from the B850 ring of the light harvesting complex 2 (LH2) of purple bacteria Rhodoblastus (Rbl.) acidophilus. This is a challenging test because one BChl molecule by itself shows two excitations, the so-called Q_y and Q_x excitations, but here, the two molecules are coupled. Consequently, one expects that there should be four excitations in the combined BChl-302 and BChl-303 system: The Q_y excitations typically show the characteristic J-aggregate type coupling, i.e., a symmetric and an antisymmetric coupling³¹ that leads to one excitation of high and one of low oscillator strength, and the Q_x excitations couple just weakly and mostly retain their individual molecular character. A calculation using the optimally tuned range-separated hybrid ωPBE, which we adapted from Ref. 73 and depict in the upper panel of Fig. 7, shows exactly these characteristics: There is one line of Q_y character with a high oscillator strength and one with low oscillator strength slightly below 1.8 eV, and the two Q_x excitations appear between 2.2 and 2.3 eV with a small energetic separation. The only further excitation in the shown energy range is a very weak line that appears above the coupled Q_x lines at 2.33 eV. A TD-LDA calculation, on the contrary, also shown in the upper panel, shows the characteristic problems that one has with spurious charge-transfer excitations in this system: Instead of two well-defined coupled Q_y excitations, there are now four excitations around 1.8 eV, because in addition to the two coupled Q_y lines, there are two low-lying spurious CT excitations. We marked the low-intensity excitations that lie close to the high-intensity ones by little arrows below the energy axis for clarity of the visualization, and the numerical values of all excitation energies can be found in the supplementary material.

The lower panel of Fig. 7 shows the TASKx results. As for NDI-1, the influence of the current density correction is negligible. In addition, as for NDI-1, the results with respect to the charge-transfer excitations are mixed: On the one hand, the coupling of the excitations is more correct than with the TD-LDA, as we observe the qualitatively correct picture of one high-strength and one low-strength excitation for the coupled Q_y and Q_x excitations. On the other hand, there are additional excitations slightly below 1.8 eV, at 2.04 eV, and at about 2.2 eV that are spuriously low, showing that the TDDFT charge-transfer problem is not fully fixed by TASKx.

In the following, we investigate the importance of the current density correction for the r²SCANx functional by looking at the same paradigm test systems in Sec. IV. We again focus on only the exchange part in order to allow for a transparent comparison with the TASKx results, and also because the r²SCAN meta-GGA correlation is expected to have a τ-dependence of opposite sign to meta-GGA exchange,⁴¹ obscuring the non-locality effects in r²SCAN for x and c.

In the spectrum that we obtain for CO in GKS and CGKS with the real-time real-space propagation, we again observe a considerable difference between GKS and CGKS, similar to what we found for TASKx in Sec. IV B. The excitation is shifted by 0.36 eV from 8.82 eV (r²SCANx-GKS) to 8.46 eV (r²SCANx-CGKS), i.e., the shift is smaller than for TASKx but still very noticeable. Once more taking the current density correction into account improves the agreement with the experimental reference.

The situation is similar for Li₂. As for TASKx, we see a red shift of the excitation spectrum when switching from GKS to CGKS. The first excitation is shifted from 2.15 to 1.93 eV, and the second excitation is shifted from 2.76 to 2.49 eV. So again, the shift is noticeable but smaller than for TASK, and again, the current density correction improves the agreement with the experimental values of 1.86 and 2.57 eV. The oscillator strength shows a drop from GKS to CGKS for r²SCANx as well as for TASKx.

Turning toward the organic molecules, we again start with benzene. When going from r²SCANx-GKS to r²SCANx-CGKS, the dominant excitation line is once more red shifted but only by 0.07 eV from 7.24 to 7.17 eV, i.e., to a much lesser extent than for CO or Li₂. The trend thus corresponds to the trend that is seen with the TASK functional, but the absolute scale is smaller. As discussed in Sec. IV B, these findings are explained by the fact that the magnitude of the current density correction is related to the electronic structure of the system, and the different meta-GGAs pick these differences up with different strengths, related to the different magnitude of ∂e_x/∂τ in the different functionals.

Figure 8 shows the r²SCANx results for NDI-1, and again, we see the trend confirmed that we already saw for the TASK functional: There are only tiny differences between GKS and CGKS, i.e., the current density correction does not have any significant effect. With respect to the energy of the charge-transfer excitation, we observe some differences between r²SCANx-CGKS and TASKx-CGKS: With the former, the charge-transfer excitation shows up at 1.92 eV, whereas it is at 2.03 eV with TASKx. Thus, TASKx does somewhat better than r²SCAN in working against the spurious downshift of charge-transfer excitations that is common for semilocal functionals. This is in line with the construction strategy of TASK, which put emphasis on the ultra-nonlocality.

Finally, we look at the description of the coupling of the BChl-molecules of the combined BChl-302-BChl-303 system. The corresponding spectrum is shown in Fig. 9. One immediately sees that the effect of the current density correction is very small. In addition, the coupling of the Q_y excitations is similarly described with r²SCANx as with TASKx: There is one spurious low-strength excitation near the coupled Q_y-lines. The situation in the energy range of the Q_x excitations is more difficult with r²SCANx: Whereas for TASKx, one could identify the two coupled Q_x lines, as the spurious CT excitations were energetically separated from them, with r²SCANx, there are four excitations that lie so close to each other that one cannot count on the correct description of the coupling. So, overall, the results for r²SCANx are qualitatively similar to the ones obtained with TASKx, but the TASKx functional is doing somewhat better on the charge-transfer excitations. Both meta-GGAs, however, do not reach the accuracy of the tuned range-separated hybrid.

We have calculated the absorption spectra of paradigm molecular model systems with the TASK and r²SCAN meta-GGA exchange functionals in the GKS scheme with and without a current density correction that restores gauge invariance and the continuity equation. The linear response was obtained by real-time propagation on a real-space grid, and we discussed how the reliability of pseudopotential results can be increased via a core-correction when meta-GGAs are used with usual norm-conserving pseudopotentials that have been constructed based on LDAs or GGAs.

Our results show that the current density correction can have a substantial influence on the excitations, changing the excitation energy by 40% and more than 0.8 eV in the most pronounced case of Li₂. The current density correction shifts the excitations to lower energies and, in the cases studied here, improves the agreement with experiment. For larger conjugated organic molecules, however, we found that the correction is negligible. The differing importance of the current density correction is a consequence of differences in the electronic structure of the different molecules to which the meta-GGA enhancement factor is sensitive. That is, the magnitude of the current density correction is system dependent because different molecules probe different values of n, ∇n, and α, and these are mapped to differing degrees of non-locality by the enhancement factor F of Eq. (1). The conjugated molecules have more delocalized densities than small molecules, such as CO and Li₂, i.e., they are somewhat more “electron gas like.” TASK and r²SCAN respect the homogeneous electron gas limit and thus reduce their non-locality for systems that become closer to the electron-gas. As the magnitude of the current density correction is related to the degree of non-locality of the meta-GGA, one expects a reduced importance of the current density correction for systems with more delocalized, electron gas-like densities. It is, however, hard to a priori quantitatively predict the extent to which the non-locality of a meta-GGA is relevant in a given system. Furthermore, attempts to do so would require a very detailed analysis of the electronic structure, and while examples of such an analysis are known, e.g., in the context of understanding bandgaps,⁴¹ analyzing the interplay of the three variables n, ∇n, and α is not anything that one wants and can do routinely. Therefore, it is advisable not only for fundamental reasons but also for practical reasons to use the current density correction in all TDDFT calculations that employ meta-GGAs.

A comparison of TASKx and r²SCANx showed that the overall trends are very similar for the two functionals, but the effects of non-locality are more pronounced with TASKx than with r²SCANx. This is in line with the different non-empirical construction strategies of the two functionals, with r²SCAN focusing on ground-state energetics, while a special emphasis was given to the ultra-nonlocality in TASK. This reasoning also explains the somewhat better performance of TASKx for charge-transfer excitations. However, neither TASK nor r²SCAN reaches the quantitative accuracy for charge-transfer excitations that one sees in optimally tuned range-separated hybrid functionals. This finding is plausible and, to some extent, expected: The derivative discontinuity physics that TASK captures is necessary for a proper description of the field-counteracting term^74,75 and charge-transfer, as repeatedly discussed in the past,^17,72,76 yet other properties of a functional are also very important for the proper description of charge transfer excitations, especially being not or only weakly affected by self-interaction. Neither TASK nor r²SCAN is free from self-interaction, as their semilocal potentials cannot cancel the non-local self-Hartree potential. Designing a non-local self-interaction correction for meta-GGAs is thus a relevant direction of future research.

See the supplementary material for a detailed discussion of the continuity equation for meta-GGAs, the real-space representations of the GKS and CGKS operators, pseudopotential parameters, further numerical parameters and computational details, tables and figures listing further excitation details, comments on propagating SCANx, and the author contribution statement.

We acknowledge the financial support by the Deutsche Forschungsgemeinschaft (DFG) under Grant No. KU 1410/4-1, the Elite Study Program “Biological Physics” of the Elite Network of Bavaria, and the Collaborative Research Network “Solar Technologies go Hybrid” from the Bavarian State Ministry of Science, Research, and the Arts.

The authors have no conflicts to disclose.

Rian Richter: Conceptualization (supporting); Data curation (equal); Formal analysis (equal); Investigation (equal); Methodology (equal); Software (supporting); Validation (equal); Visualization (equal); Writing – original draft (equal); Writing – review & editing (supporting). Thilo Aschebrock: Conceptualization (supporting); Investigation (equal); Methodology (equal); Software (lead); Writing – original draft (equal); Writing – review & editing (supporting). Ingo Schelter: Software (supporting); Supervision (supporting); Validation (supporting); Writing – review & editing (supporting). Stephan Kümmel: Conceptualization (lead); Formal analysis (equal); Funding acquisition (lead); Investigation (equal); Methodology (equal); Project administration (lead); Resources (lead); Supervision (lead); Writing – original draft (equal); Writing – review & editing (lead).

The data that support the findings of this study are available within the article and its supplementary material. Raw numerical data are available from the corresponding author upon reasonable request.