Minimal auxiliary basis set for time-dependent density functional theory and comparison with tight-binding approximations: Application to silver nanoparticles Open Access

The Time-Dependent Density Functional Theory (TDDFT) is the most used first-principles tool for the calculation of excitation energies in electronic systems.^1–5 TDDFT is nowadays implemented in many free or commercial codes and can, in principle, predict the exact excitation energies. In practice, its accuracy depends on the exchange–correlation (XC) functional. For molecular systems, TDDFT has an accuracy of about 0.2 eV–0.5 eV, depending on the accuracy of the XC functional.⁶ More recently, TDDFT has been applied to plasmonic systems,^7–14 e.g., silver and gold nanoparticles (NPs), which play a key role for applications in nano-optics, photovoltaics, surface-enhanced Raman scattering, and imaging.^15–19 The absorption spectrum of metal NPs is mainly characterized by a localized surface plasmon resonance (LSPR), which can be accurately described using classical electromagnetism,^20,21 but only for large systems. TDDFT, instead, includes all quantum effects, which are not considered in classical models or only approximated in other theoretical approaches for quantum plasmonics.^12,22–27

Thanks to the computational efficiency of TDDFT, calculations for systems with more than 100 atoms have been reported for silver^28–39 and gold^{29,31,34,35,39–45} nanosystems. In one of the largest application reported so far,⁴² the TDDFT absorption spectrum of Au₁₄₁₄ has been computed, but using more than 60 000 cores. Despite its efficiency, thus, TDDFT can be hardly applied to systems with thousands of atoms. To this end, different methods and algorithms have been developed^46–50 in order to reduce the computational scaling. TDDFT implementations can be distinguished into two main classes: (i) real-time TDDFT and (ii) linear-response TDDFT. In the real-time approach,⁵¹ the evolution of the time-dependent Kohn–Sham (KS) states is computed following a perturbation. Absorption spectra are, then, obtained as the Fourier transform of the time-dependent dipole moment of the system. In the linear-response approach, a Dyson equation is solved in the frequency space,^3,5 and excitation energies are obtained as poles of the interacting density–density response, or alternatively, the excitation energies can be obtained directly solving an eigenvalue problem as first pointed out by Casida.² The spectra of both methods are equivalent.^51,52

In this work, we will focus on the Casida formalism and, in particular, on the construction of the response matrix, which requires the calculation of Coulomb and XC kernel integrals. A standard approach to reduce the TDDFT computational cost is the resolution-of-the-identity (RI) method, in which an auxiliary basis set is used to expand the electronic density for a faster computation of the Coulomb integrals.^53–61 More recently, following ideas of Time-Dependent Density Functional Tight Binding (TD-DFTB) methods,^62–66 approximated TDDFT approaches have been proposed: in the simplified Tamm-Dancoff approximation (sTDA),⁶⁷ sTDDFT,⁶⁸ and TDDFT + TB⁶⁹ methods, the computation of kernel integrals is avoided completely using semiempirical formulas. Although the applications of those methods for organic molecules are well established, investigations for more complex systems such as metal NPs have appeared only recently.^70–72 

In this work, we introduce a modification of the first-principles RI-TDDFT approach, in which the auxiliary basis set is composed of a single s-type function per atom. In the conventional RI approach, the size of the auxiliary basis set is very large as it is usually optimized to reproduce the total Coulomb energy with chemical accuracy.^60,61,73

Which is the impact of such minimal auxiliary basis set on the absorption spectra?

We show that our method, named TDDFT-as, for auxiliary s-type, can predict excitation energies for different silver nanoparticles with an average accuracy of only 12 meV, which is really negligible considering the much larger errors related to the orbital basis set and the choice of the XC functional.

The TDDFT-as method resembles the TDDFT + TB approach. In this work, we analyze and compare these two approaches in details. We found that the accuracy of TDDFT-as for silver nanoparticles is distinctively superior, and the origin on the TDDFT + TB lower accuracy is clarified.

This article is organized as follows: In Sec. II, we review the RI-TDDFT and the TDDFTB + TB theories and discuss the conditions under which these approaches can be considered equivalent; in Sec. III, we describe the systems under investigation, the computational details of the TDDFT-as approach, and the procedure we follow to compare absorption spectra computed by two different methods; in Sec. IV, we discuss the results of the exponent optimization and the comparison of the TDDFT-as results with reference TDDFT and the approximated TDDFT + TB ones; finally, in Sec. V, conclusions and future perspectives are drawn.

The TDDFT Casida eigenvalue equation^2,4,5 is (considering for simplicity, closed-shell systems and singlet excitations)

where ω_I are excitation energies and

In Eqs. (2) and (3), a, b (s, t) denote occupied (unoccupied) KS orbitals with eigenvalues ϵ, the single-particle gaps are ω_as = ϵ_s − ϵ_a, and $fxc=δ2Excδρ↑δρ↑+δ2Excδρ↑δρ↓$ is the XC kernel. The integrals of the Coulomb kernel, $(as|bt)=∫∫ϕarϕsrϕbr′ϕtr′/r−r′d3rd3r′$⁠, appearing in Eq. (3) can be simplified using the RI approach, where an orbital basis function product can be approximated as

where μ, ν, … indicate orbital basis functions, P, Q indicate the auxiliary basis functions, and the coefficients $ΔμνP$ minimizing the error in the Coulomb norm are given by^56,59

Here and hereafter, we use the notation $(a|b)=∫∫χarχbr′/|r−r′|d3rd3r′$⁠. Orbital (and auxiliary) basis functions are atom-centered functions of the type $χμr=χAlmr=RAl(|r−RA|)Ylm((r−RA)/|r−RA|)$⁠, where Y_lm are real spherical harmonics and R_A is the position of atom A.

In the RI approach, the Coulomb integrals appearing in Eq. (3) become

where $Cμa$ are the coefficients of the KS orbital ϕ_a and

is the contribution of the KS orbital product ϕ_aϕ_s to the auxiliary basis function χ_P.

The eigenvectors of Eq. (1) are related to the transition density, i.e., the first-order change of the electronic density: for excitation I, the transition density is

It is interesting to show that Eq. (7) can be recast to the approximation, which is used in TB schemes.^62–69 In fact, if in Eq. (5) we do a Mulliken approximation [i.e., $χμr→⟨χμ|χν⟩χνr=Sμνχνr$⁠, and symmetrizing], we obtain

where in the second line, we assume that each $|χμr|2=|χAlmr|2$ corresponds to the auxiliary basis function $χAμ$ [i.e., by keeping only the monopolar term of $|χμr|2$ so that every orbital basis function on the atom A corresponds to the same auxiliary basis function $χAμ$]. Using Eq. (11) in Eq. (8), we obtain (using the atom index A to name the auxiliary basis function as there is just one auxiliary basis function per atom)

which equals the TD-DFTB Mulliken atomic transition charges.⁶² 

It is important to observe that atomic transition charges in Eq. (12) are “population analysis,” i.e., they depend only on the overlap between atomic orbitals and do not depend at all on the auxiliary basis function. This is very different from Eq. (8), where $ΔμνP$ clearly depends on the auxiliary function P, i.e., on the exponent of the Gaussian function, see Eq. (5). Recently, another population analysis (Löwdin) has been used in STDA and TDDFT + TB^67,69 instead of the Mulliken population.

Finally, combining Eq. (12) with Eq. (7), we obtain

which is the TD-DFTB approximation for the Coulomb kernel integrals.^62–69 

Actually, Eq. (13) corresponds to the full TDDFT kernel, i.e., with Coulomb and XC effects, if the latter are treated with semilocal XC functionals. In fact, in this case, the last term in Eq. (3) can be approximated as

where W_A is an onsite parameter, which can be, thus, incorporated into Eq. (13). In fact, in the conventional TD-DFTB approaches and related methods, the two-center integrals (A|B) are usually not computed from the basis set but with an analytic function of the internuclear distance (R_AB) and Hubbard parameters (U) of the atom types. For example, the sTDA/sTDDFT methods^67,68 (for the same atom type, as in the present work) employ

where a is a parameter. Thus, we have that U = γ(0) can be computed from isolated atoms and includes both the Coulomb and XC effects. Note also that for the metallic nanoparticles considered in this work, the role of the XC kernel is very small.⁷⁴ 

In summary, we have that the TB integrals in Eq. (13) employ three approximations:

• a single and spherical function for each atom is used to expand the transition densities;

• an atomic Hubbard parameter U is used to define the single and spherical function;

• a population analysis (Mulliken or Löwdin) approximation is used to compute the atomic transition charges $Γ̃asA$⁠.

On the other hand, no approximations (beside the size of the auxiliary basis set) are used for the general expression in Eqs. (5), (7), and (8).

Despite the quite good accuracy of the conventional TD-DFTB approaches and related methods,^62–69 no detailed study of the role of approximations (i)–(ii)–(iii) is present in the literature.

In this work, we propose to keep only one Gaussian s-type orbital function per atom in the auxiliary basis set, i.e., approximation (i), but the Gaussian exponent is fitted to reproduce the absorption spectra of test systems, i.e., we do not use an atomic Hubbard parameter as in approximation (ii). Moreover, the atomic transition charges are computed exactly from Eqs. (5) and (8), without using the population approximation (iii). In this way, we can test directly the accuracy of using only approximation (i). We named such a method with one s-type function per atom in the auxiliary basis set as “TDDFT-as.”

In TDDFT-as, the two-center integral between two normalized s-type functions, i.e., $gsr=(α/π)3/2⁡exp(−αr2)$ with $∫gsrd3r=1$⁠, centered on atoms A and B is^75,76

where R_AB is the distance between atoms A and B. When R_AB = 0, we have

which relates the exponent of the Gaussian s-type function to the Hubbard parameter. Note, however, that, in general, the auxiliary basis set is Coulomb normalized, i.e., $gCsr=Uα−1/2gsr=(23/4α5/4)/(2π5/4)exp(−αr2)$⁠. If α = 0, all the kernel contributions vanish and a single-particle absorption spectrum is obtained. In Sec. IV A, we will verify if an optimal α exists for silver NPs.

In this work, we considered eight silver NPs of different sizes and shapes: Ag₂₀, $Ag35+$⁠, and Ag₁₂₀ (tetrahedron, T_d symmetry); Ag₃₂ (cubic rod, D_4d symmetry); $Ag55−1$ and $Ag555+$ (cuboctahedron, O_h symmetry); and $Ag553−$ and $Ag555+$ (icosahedron, I_h symmetry).

All the simulations including geometry relaxations and ground state and excited state calculations have been performed using the TURBOMOLE program package.⁷⁷ For all these tasks, we used the Perdew–Burke–Ernzerhof (PBE)⁷⁸ functional for the XC energy, the default effective core potential (ECP) with 28 core electrons,⁷⁹ and the split valence plus polarization (def-SVP) basis set.⁸⁰ This choice represents a good compromise between accuracy and computational efficiency.

For each system, we first performed a reference TDDFT calculation with the conventional auxiliary⁸¹ basis set, i.e., 7s4p3d3f2g, which corresponds to 97 Cartesian basis functions per atom. Then, we ran TDDFT-as calculations with just one s-type auxiliary basis function and neglecting the XC kernel.

All TDDFT calculations have also been done freezing the 4s and 4p occupied orbitals (i.e., with eigenvalues below 2 a.u.) and a number of virtual orbitals corresponding to the number of atoms (i.e., with eigenvalues above 50 a.u.). The number of considered occupied and virtual orbitals is reported in Table I. Freezing these orbitals has no effect on the spectra in the considered range (i.e., up to about 5 eV).

For the exponent optimization (see Sec. IV A), we considered all excited states up to the energy ω_max reported in Table I, with the corresponding number of (optically allowed) states considered. In this way, we considered not only the main plasmon peak but also structures in the high-energy region.

We also performed TDDFT + TB calculations using a locally modified version of the STDA code,⁸² in which the kernel integrals have been replaced by Eq. (17). The STDA code, in contrast to TURBOMOLE, does not use any symmetry: in order to reduce the computational cost, we have removed from the transition space all single-particle transitions, which are optically dark by symmetry.

To estimate the accuracy of TDDFT-as results with respect to the reference TDDFT ones, we considered the following root mean square (rms) error indicators for excitation energies and oscillator strengths (f_i):

where N is N_states in Table I and the function j(i) associates the excited state i of the reference TDDFT calculation with the excited state j in the TDDFT-as calculation. Assuming j(i) = i means that the order of excited states is the same in TDDFT-as and TDDFT, which is often not the case. The association is done considering the quantity

and, for a given i, selecting the j-state with the largest O_ij. Note that 0 < O_ij < 1. In Eq. (21), the first part is the overlap between the eigenvectors of Eq. (1), whereas the second exponential factor is included to select states with closer oscillator strengths. In fact, in some cases, the overlap between different eigenvectors can be similar, and the second exponential factor helps to select the right state.

The indicators E_ene[α] and E_osc[α] are very accurate, but also quite complicated and require that the two sets of excitation energies and oscillator strengths are close each other otherwise the association will fail (i.e., O_ij will be always too small). A simpler indicator with a more general applicability can be defined as

where σ^ref(ω) and σ(ω)[α] are the photoabsorption cross sections for TDDFT and TDDFT-as, respectively, and ω_max is defined in Table I.

The photoabsorption cross section is obtained from the computed excitation energies and oscillator strengths using a Lorentzian broadening, namely,

The function E_spe[α] does not include any direct information about the error on excitation energies; moreover, it depends on the broadening g. However, it is very simple to compute and, as we will show in this article, can give quite accurate information.

To optimize the exponent of the auxiliary basis set, we considered seven silver NPs: all the ones in Table I but Ag₁₂₀, which has been considered only for timings, see Sec. IV D.

In Fig. 1(a), we report the error E_ene[α] for all systems considered. The first interesting result is that all the curves have a similar shape with a single, quite flat minimum at around α = 0.03, more precisely in the range from α = 0.026 for $Ag55−3$(I_h) to α = 0.040 for Ag₃₂(D_4h). Although the values of optimal α change a little among different systems, most of the curves are quite flat, meaning that accurate results can be obtained within a significant range of α values.

The second important result is that the value of E_ene[α] at the minima is very small, in the range 7 meV–17 meV, meaning that TDDFT-as predicts all excitation energies with an average rms error less than 20 meV, which means a very high accuracy (less than the orbital basis set error).

Then we considered, see Fig. 1(b), the quantities ⟨E_ene⟩ and ⟨E_osc⟩, which are the average among all systems. The ⟨E_ene⟩ minimum is for α = 0.036 with an average rms error among all systems of only 12 meV.

For the oscillator strength, the ⟨E_osc⟩ curve does not show any well-defined minimum: it is very flat, with an average rms error of 0.25. Thus, we can consider the optimized value for ⟨E_ene⟩,

as a fixed exponent for all TDDFT-as calculations of silver NPs. An error distribution analysis (see Fig. S1 of the supplementary material) shows that all energies are underestimated by TDDFT-as. The results for the E_spe indicator (see Fig. S2 of the supplementary material) give similar conclusion but with a smaller exponent (α = 0.034). Although the absolute value of E_spe depends on the broadening parameter, it turned out that the optimized α is not.

The accuracy of TDDFT-as in reproducing the TDDFT reference spectra can also be directly recognized in Fig. 2, where three structures (the smallest one and the two extreme cases for the optimal α) have been reported. Figure 2 clearly shows a very high accuracy over the whole energy range. More in detail, we note that the intensity of the main plasmon peak is a bit reduced in TDDFT-as. We found that the intensity on the main peak in TDDFT-as can be increased reducing the α value; however, in this case, also the intensity of all other peaks at high-energy is increased, thus decreasing the overall accuracy.

In Fig. 3, we report the reference transition densities, i.e., computed using KS orbitals, see Eq. (9), and the ones using the auxiliary basis set, see Eq. (10), considering just one optimized s-type function.

The reference transition densities show a very rich structure due to atomistic details, and it is not very simple to distinguish the direction of the dipole. On the other hand, the TDDFT-as approximated transition densities are much smoother and without any atomistic structure, but both having similar (deviations about 10%, see Table S1 of the supplementary material) oscillator strength of the reference one. Thus, a much simpler interpretation of dipole moment and the charge redistribution in the excited state is obtained; moreover, the calculation of transition density plots is much cheaper.

However, we note that the transition dipole moments [$μI=∫rδρIrd3r$], which are used to compute the TDDFT-as absorption spectrum, are computed directly from orbitals, i.e., using Eq. (9), and not from Eq. (10). In fact, the auxiliary basis set is used in TDDFT-as only to approximate the Coulomb integrals in the TDDFT kernel matrix. An alternative approach is to use directly the expansion in Eq. (10) also for the computation of the transition dipole moments. This is, indeed, very efficient when only one s-type basis set is involved; in fact, in this case, the transition dipole moment is

where the sum in A is over all atoms and

Thus, the calculation of transition dipole moment between each pairs of atomic orbitals is not required anymore. A comparison of the oscillator strength for the main plasmon peaks for all systems is reported in the supplementary material (Table S1). We have also recomputed E_osc for all systems with this alternative definition of the transition dipole moment. Results show (see Fig. S3 of the supplementary material) that a well-defined minimum is present, in contrast to Fig. 1(b). In this case, in fact, there are stronger constraints (both Coulomb integrals and dipole moments) on the optimization of the α exponent. We found that the optimized value is α = 0.039, thus very close to the one obtained in Sec. IV A. However, the rms error on the oscillator strength increases by about 40% (⟨E_osc⟩ = 0.034), and the resulting spectra are not very accurate at high energies. Thus, while transition density plots at the plasmon peak can be safely used, Eq. (9) is more accurate for the calculation of the absorption spectrum, and it is used by default in TDDFT-as calculations.

The optimized α value for the TDDFT-as approach (α = 0.036) for the silver atom corresponds to an onsite Hubbard parameter U = 0.151 a.u., see Eq. (18). This value is much smaller than the one used in DFTB parameterization (U = 0.2191 a.u.) or in the sTDA and sTDDFT methods (U = 0.250 a.u.).

In Fig. 4, we report the absorption spectra computed with TDDFT + TB with different values of Hubbard parameters, for three selected systems Ag₂₀(T_d), Ag₃₂(D_4h), and $Ag55−3$(I_h).

Recall that while in TDDFT-as atomic transition charges are computed from Eqs. (5) and (8), in TDDFT + TB atomic transition charges are computed with the Löwdin population approximation.⁶⁹ For two-center integrals, we also use Eq. (17) in TDDFT + TB so that our implementation of TDDFT + TB is a little bit different from the original scheme.⁶⁹ A modification of the two-center formula has some effects on the spectrum, but much less than the modification of the Hubbard parameter: recall, in fact, that all two-center formulas have γ(0) = U and γ(R) = 1/R for R large; thus, modifications can occur only for first neighboring atoms. In both TDDFT-as and TDDFT + TB, transitions dipole moments are computed directly from KS orbitals, i.e., computing the dipole moment of Eq. (9). Thus, the only difference between TDDFT-as and our TDDFT + TB is in the atomic transition charge computational scheme, besides the actual value of the Hubbard parameter (or α exponent). The top panel of Fig. 4 clearly shows that if the Hubbard parameter from our TDDFT-as approach is combined with a population analysis approximation for the atomic transition charges, very bad absorption spectra are obtained for Ag₂₀ (T_d). The TDDFT + TB spectrum using U = 0.151 a.u. is also very bad for the other two systems, and it is not reported for graphical reasons.

Using a larger Hubbard parameter, improved spectra are obtained for Ag₂₀(T_d). Yet, the default TDDFT + TB still overestimates the intensity of the main plasmon peak. A better accuracy is obtained using the optimized value (U = 0.337 a.u.), obtained minimizing E_spe, which improved the description of the plasmon peak. Nevertheless, the intensities of high-energy part (>4 eV) of the spectrum are overestimated.

For $Ag55−3$(I_h), bottom panel of Fig. 4, the situation is similar. In this case, the optimized Hubbard parameter is even larger (U = 0.4 a.u.), and the agreement with the reference TDDFT is better. Instead, the behavior of Ag₃₂(D_4h) is very different. With default TDDFT + TB, very bad absorption spectra are obtained, with a big peak at around 3.8 eV, not present in the reference TDDFT. In this case, increasing the Hubbard parameter does not help because in this way, the main plasmon peak at 2.9 eV disappears.

In Ref. 72, TDDFT and TDDFT + TB absorption spectra for silver and gold clusters have been reported. Silver clusters are similar to the ones of this work, but the TDDFT reference employs a different XC functional⁸³ and TDDFT + TB a different formula for two-center integrals. In any case, results of Ref. 72 indicate a similar accuracy and overestimation of the oscillator strengths of TDDFT + TB with respect to TDDFT at high energy, which is consistent with our findings.

Comparing Fig. 4 with Fig. 2, the much higher accuracy of TDDFT-as with respect to TDDFT + TB is, thus, evident. A more quantitative analysis is reported in Fig. 5 for all systems. The default TDDFT + TB yields an error for E_spe in the range 0.70–1.08. Optimizing the Hubbard parameter for each system, the E_spe can be reduced down to 0.25–0.69. Clearly, these values of E_spe are much larger than the ones that can be obtained using the TDDFT-as with α = 0.036, with an error of E_spe in the range 0.09–0.20.

This direct comparison between TDDFT + TB and TDDFT-as clearly indicates that the origin of the reduced accuracy of TDDFT + TB is related to the population analysis approximation (as all other quantities are the same in the two methods).

In Fig. 6, we compared the atomic transition charges computed with our TDDFT-as method, i.e., using Eq. (8), with the ones from the TDDFT + TB approach, for two different single-particle transitions of Ag₂₀. The plots clearly show that the TDDFT + TB atomic transition charges are much smaller than the TDDFT-as ones; thus, a larger Hubbard parameter is required in the former to obtain similar Coulomb integrals. This finding explains the results in Fig. 4 where an incorrect spectrum is obtained with U = 0.151 a.u.

Finally, we considered a larger nanoparticle, Ag₁₂₀ (T_d), not included in the fitting procedure. The computed photoabsorption cross sections from TDDFT and TDDFT-as are reported in Fig. 7. We used a log-scale for the y axis and a small broadening (g = 20 meV) so that the accuracy for all peaks can be well distinguished. Clearly, the agreement with the reference TDDFT calculation is excellent as the difference between the two curves can be hardly distinguished, for all the energy range considered.

We also considered the timings for the Ag₁₂₀(T_d) nanoparticle. We report in Table II the time (in seconds) for the Coulomb, XC, and RI part, i.e., the calculation in memory of the three-center integrals (Q|μν), see Eq. (5). The TDDFT-as computational cost for all these steps is less than 0.5% of the full TDDFT results, which is consistent with the reduction in the auxiliary basis size of two orders of magnitude and the reduced number of Davidson iteration steps required.

However, we point out that we did not optimize the code for efficiency. In fact, the main focus of the present work was to verify the accuracy of a single s-type function in the auxiliary basis set.

In this article, we have shown that it is possible to reproduce TDDFT results for silver nanoparticles with an accuracy of 12 meV using only one s-type function per atom in the auxiliary basis set. In addition, the exponent α of this s-type function is optimized to model XC effects too. In this approach, named TDDFT-as, the computational cost for the calculation of the TDDFT kernel matrix elements is reduced by two orders of magnitude, with respect to TDDFT with a default auxiliary basis set. The TDDFT-as method is based on two observations: (i) to reproduce the absorption spectrum of silver NPs, there is no need to reproduce exactly all the atomistic details of the transition density: they can be averaged out for the calculation of Coulomb integrals and also for the dipole moment and (ii) the effects of a (semi-)local XC kernel are much smaller than the Coulomb one;⁷⁴ thus, the former can be included as an onsite correction of the latter.

In this work, we have optimized the exponent among seven silver NPs with different sizes and symmetries, showing the transferability of the α exponent. The TDDFT-as approach accurately reproduces the main plasmon peak as well as the high-energy part of the absorption spectrum.

The TDDFT-as approach resembles the TDDFT + TB scheme: both are based on the KS eigenvalues and eigenvectors and use a monopolar approximation to compute the Coulomb-XC kernel matrix elements. TDDFT + TB employs two other additional approximations, i.e., the population analysis to compute the atomic transition charges and a Hubbard atomic parameter. In this work, we have shown that TDDFT-as is by far more accurate than any TDDFT + TB approach with a population analysis approximation and optimized Hubbard parameter. For some systems, e.g., a cubic silver rod, TDDFT + TB even predicts a spurious peak in the high-energy part of the spectrum.

The TDDFT-as approach is very simple, and it is easily available in all RI-TDDFT codes (changing the auxiliary basis set and switching-off XC contributions).

There are many directions to be followed in the future:

In this work, we have, in fact, only considered weakly charged systems and the PBE XC functional. It will be interesting to study if the optimal α will change with the KS potential, which can be largely modified using a different XC functional (i.e., the LB94 functional⁸³) or by a large system charge. In this case, an alternative approach could be used to extract a system-dependent optimal α, using information from the ground-state KS orbitals.

Another path to follow is to further optimize the auxiliary basis set in order to verify if it is possible to obtain even better accuracy, without increasing the computational cost. It will be interesting, for example, to use a contraction of more s-type functions or different s-type functions for different atoms (e.g., on the surface or inside the nanoparticle).

Finally, TDDFT-as can be extended to other systems, e.g., other metals, metal alloys, or even molecular systems. We expect that TDDFT-as can be quite accurate for all these systems as it is closer to TDDFT than TD-DFTB and TDDFT + TB. Preliminary calculations show that TDDFT-as is quite accurate in reproducing the absorption spectra of gold NPs. However, the overall accuracy is lower for silver NPs; in fact, gold NPs have a more complex absorption spectrum with stronger contributions from the d-band,⁸⁴ which cannot be fully described with just one s-type function. When systems with different atom types are considered, it must be kept in mind that the optimized α for a given atom might also depend on the local environment as it aims to describe transition densities, which can be quite delocalized. Extensions of the TDDFT-as method (e.g., including two basis functions per atom) are also possible.

See the supplementary material for TDDFT-as error distribution for excitation energies and oscillator strengths, excitation energies and oscillator strengths at the plasmon peak for all systems, error E_spe for all systems and averaged value, and error E_osc for all systems and averaged value using transition dipole moments from the auxiliary basis set.

The data that support the findings of this study are available from the corresponding author upon reasonable request.