All-electron real-time and imaginary-time time-dependent density functional theory within a numeric atom-centered basis function framework

Many physical, chemical, and biological processes of interest in modern science and technology ultimately originate from non-equilibrium dynamics of electrons under external stimuli. Electron dynamics are at the heart of phenomena including photo-excitation, hot carrier relaxation, charge recombination, and electronic stopping. Starting with the time-dependent Hartree–Fock method in the late 1970s and early 1980s,¹ various theoretical methodologies including correlated wavefunction methods² have been developed for studying non-equilibrium response properties with the real-time propagation approach.³ As a computationally simpler alternative to wavefunction-based methods, the real-time propagation approach to time-dependent density functional theory (TDDFT) has gained great popularity in the last few decades. Since TDDFT can be rigorously justified due to the Runge–Gross theorem,⁴ it offers a powerful extension to density functional theory (DFT) for studying dynamical phenomena from first principles.^5–9 Through its formulation for the linear response (LR) regime (i.e., LR-TDDFT)^10–12 in the frequency domain, the TDDFT approach has become extremely popular for studying optical excitation of molecules and materials.¹³ At the same time, the real-time (RT) propagation approach to TDDFT (i.e., RT-TDDFT)¹⁴ has become a computationally powerful simulation method for studying non-equilibrium dynamics of electrons, in general, not limited to the realm of linear response theory. Since some of the first uses of the real-time propagation approach in the 1990s by Yabana and Bertsch,^14,15 RT-TDDFT has gained great popularity in various areas of chemistry and condensed-matter physics.^16–22 RT-TDDFT can be used to describe both linear and nonlinear responses of matter to perturbations of various types and strengths. For large systems and certain regimes of electronic excitation, RT-TDDFT can be also more computationally efficient than the commonly used linear response TDDFT method even for calculating optical absorption spectra.²³ In recent years, there has been a surge in applications of RT-TDDFT methodologies to a wide range of excited state phenomena, such as electronic stopping,^24–28 optical absorption,^{15,16,29–31} core electron excitations,^32–36 electronic circular dichroism spectra,³⁷ exciton dynamics in nanostructures,^38,39 atom-cluster collisions,^40,41 high harmonic generation,^42,43 laser-induced water splitting,⁴⁴ and ionization processes.^45–47 The promulgation of RT-TDDFT as a means for simulating dynamic phenomena of electrons has led to its implementation in a variety of electronic structure codes, including NWChem,^16,48 SIESTA,^49,50 CP2K,^51,52 SALMON,⁵³ Octopus,^54–56 Q-Chem,^57–59 GAUSSIAN,^20,60 MOLGW,^61,62 Qbox/Qb@ll,^63–66 Exciting,⁶⁷ and ELK.⁶⁸ Numerical details of the implementation vary widely among these codes, especially in the underlying set of basis functions used. In the present work, we describe an implementation of RT-TDDFT using numeric atom-centered orbital (NAO) basis functions in the electronic structure code FHI-aims,⁶⁹ a framework that (for ground-state DFT) combines high numerical precision^70,71 with computational efficiency and scalability to large system sizes (e.g., Fig. 11 in Ref. 72). Importantly, both non-periodic and periodic boundary conditions are supported, i.e., both molecular and condensed-phase simulations can be efficiently pursued. For comparison with RT-TDDFT, FHI-aims already includes an implementation of LR-TDDFT, as well as of the Bethe–Salpeter equation for neutral excitations of the electronic system, for non-periodic systems.⁷³ 

Although both RT-TDDFT and it-TDDFT are well-defined problems mathematically, numerical integration of coupled non-linear differential equations requires attention to detail; the question of which numerical method is most appropriate depends on the particular implementation. As mentioned above, we here provide these details for the NAO-based framework employed in the FHI-aims code.⁶⁹ Importantly, FHI-aims treats core and valence electrons on equal footing within the full potentials of the electrons and of the nuclei. In this sense, our RT-TDDFT/it-TDDFT implementation has the potential to treat both molecular and condensed-phase physics and can simulate both valence and core level physics. The use of an all-electron NAO basis set offers the potential of efficient RT-TDDFT simulations in a very general treatment.

This paper is structured as follows: we first present implementation details of RT-TDDFT/it-TDDFT in FHI-aims. We then compare the RT-TDDFT implementation to the linear-response theory formulation of TDDFT (LR-TDDFT) for a molecular benchmark set. Subsequently, several applications of RT-TDDFT are presented for molecular and condensed-matter systems and valence and core level excitations. The it-TDDFT is then shown to work among systems of different boundary conditions as well. We additionally provide a short overview of the numerical scaling.

We included support for local-density approximation (LDA), generalized gradient approximation (GGA), and hybrid functionals (the latter currently only for finite systems). For hybrid functionals, the computation of the exchange matrix (which is naturally complex in the time-dependent case¹⁶) uses the existing efficient resolution-of-identity methodology.⁸⁷ 

In terms of the computational resources required, it should also be noted that real-space operations (density update, integration of the Hamiltonian matrix, etc.) still dominate the overall computational cost for typical everyday tasks, i.e., semilocal DFT involving several tens or perhaps somewhat above a hundred atoms. For small systems (several tens of atoms), the real-space operations can be several orders of magnitude more costly in comparison to matrix-related operations like diagonalization. The matrix-related operations will nevertheless become comparatively costly for very large systems composed of hundreds of atoms and basis functions.

In all benchmarks and examples provided in this paper, the nuclei remain clamped, i.e., we focus on assessing electron dynamics only. However, the overall RT-TDDFT implementation described here also supports Ehrenfest dynamics for both non-periodic and periodic systems, a topic that will be covered elsewhere.

While the linear-response theory formulation of TDDFT (LR-TDDFT) is widely used to calculate optical absorption spectra, RT-TDDFT, when restricted to small electric fields, is a formally equivalent alternative to obtain the same information. Here, we briefly discuss the pros and cons from a computational perspective. While LR-TDDFT is computationally much cheaper than RT-TDDFT for small systems, there exists a limit for which RT-TDDFT becomes more efficient, as LR-TDDFT (without restrictions to the occupied and unoccupied orbital spaces considered) formally scales as $O(Nat6)$⁠, with N_at being the number of atoms in the system. In addition, the number of excitations included in a LR-TDDFT simulation becomes prohibitive⁹⁹ at some point, typically around (10³–10⁴).²³ In contrast, RT-TDDFT simulations yield the whole excitation spectrum when a broadband excitation, i.e., a delta-pulse, is used to excite the system. Another advantage of RT-TDDFT is that in contrast to LR-TDDFT, no exchange–correlation kernel (the second variational derivate of the XC functional with respect to the charge density) is required—a quantity that needs to be computed additionally and that may introduce numerical problems. Similar conclusions can be drawn from recent studies comparing LR- and RT-TDDFT approaches in the linear-response regime.^16,23,100

Despite their formal equivalence, the technical differences between the LR- and RT-TDDFT approaches deserve further study, especially the appropriate choice of basis sets used to represent Kohn–Sham wave functions and the exchange–correlation functional (and the exchange–correlation kernel in LR-TDDFT). In addition, the techniques used to extract and compare the same information (i.e., oscillator strengths and excitation energies) from both approaches need to be discussed.^101,102

To assess these points and also to validate our implementation for this regime, we chose to perform benchmark calculations for a popular electronic excitations ab initio simulation test set, known as Thiel’s set by Schreiber et al.,¹⁰³ consisting of 28 small to medium-sized organic molecules, built from the chemical elements H, C, N, and O.

We emphasize that we do not attempt to evaluate the accuracy of our method with regard to excitation energies, whose “best estimates” were defined by Schreiber et al. but rather take a route similar to Liu et al.⁷³ where the numerical precision of both a G₀W₀-based Bethe–Salpeter equation (BSE@G₀W₀) and LR-TDDFT implementation in the FHI-aims code was analyzed. This approach offers a direct comparison between different levels of theory within a single consistent numerical framework, especially with respect to basis set convergence. An important result of the study by Liu et al. was that combining standard NAO basis sets⁷³ with additional, Gaussian-type augmentation basis functions as pioneered by Kendall et al.¹⁰⁴ provides a significant accuracy increase at reasonable cost when performing absorption spectrum calculations with LR-TDDFT and BSE approaches using NAO basis sets. Since these small-to-medium-sized calculations impose a much higher computational demand when done via RT-TDDFT as compared to LR-TDDFT, we chose to perform calculations only for two basis sets, called “tight” and “tight + aug2” in Table I—which are probably of most practical relevance. Specifically, these basis sets combine the so-called “tier2” NAO basis sets⁶⁹ with “tight” numerical settings for integration grids, basis function cutoff radii, and Hartree potential expansion order in FHI-aims. In the case of “tight + aug2,” called “tier2 + aug2” in Ref. 73, s and p Dunning-type augmentation functions¹⁰⁴ are additionally included. We extracted singlet and triplet excitation energies, the latter based on an approach where the spin symmetry is broken by applying a spin-dependent external field, as described in Ref. 105.

The excitation energies {ω_n} selected for analysis are refined by a least-squares fit that simultaneously yields the oscillator strengths {f_n}. The resolution of the ω-axis and, thus, the accuracy of the fit can be increased by adding zeros to the dipole time series (zero-padding) before performing the discrete Fourier transforms in order to obtain the DSF. Schelter et al. further discussed another accuracy-increasing approach by Bruner et al.¹⁰²—based on a Padé approximant for the electronic dipole’s Fourier transform, allowing to arbitrarily set the resolution for the Fourier transform of the absorption strength function. While Schelter et al. found their method preferable, they suggested that more data could be extracted by combining both methods. The key idea is to first use the Padé approximant to simulate highly resolved, very sharp spectra to visually locate peak positions, which are then passed to Schelter’s fitting method for refinement. In the Padé approach, the Fourier transform of the signal is represented as a ratio of power series whose coefficients are fitted to the data; as these are not depending on the frequency, the spectrum can be calculated for arbitrary frequencies, i.e., it is interpolated in frequency space—this corresponds to an extrapolation of the signal in time. We refer the interested reader to Ref. 102 for details of the Padé method.

We decided to implement and compare both approaches because it became quite clear that they complement each other: the fitting method often fails in practice for close lying and/or very weak excitations and is increasingly difficult for spectra with dense excitations; the Padé method can help one to detect also weak and close lying excitations but may introduce artifacts for very short simulation times, as discussed in Ref. 102. See Fig. 2 for a depiction of the methods. It should be noted that the approach by Bruner et al. is further based on applying the Padé method to specific occupied–unoccupied transition dipole moments in order to accelerate spectra convergence by utilizing the sparse spectral density of those individual contributions to the total spectrum.

We here only verified this technique but did not apply it to all calculations since the effort required for this analysis is significant. Nevertheless, this approach should be very practical when evaluating “difficult” spectra.

We used custom python scripts to perform post-processing, making use of the Numpy¹⁰⁷ and SciPy¹⁰⁸ libraries for numerical treatment, especially fast Fourier transform (FFT) (scipy.fft) and non-linear least-squares fitting for the DSF fit (curve_fit from scipy.optimize).

Simulations were carried out for the two basis sets denoted as “tight” and “tight + aug2” in Table I. We performed geometry relaxations of all molecules and for the two basis sets using the Perdew–Burke–Ernzerhof (PBE) exchange–correlation functional.¹⁰⁹ These geometries were then used for any further TDDFT calculations.

Linear-response TDDFT calculations were started from PBE self-consistent field solutions and carried out using the Perdew–Wang local-density approximation (PW-LDA)¹¹⁰ exchange–correlation kernel as provided by the LibXC library.¹¹¹ Analogous to Ref. 73, we use the notation “LR-TDDFT-LDA@PBE” for this approach.

All RT-TDDFT simulations were carried out for a simulation time of 2500 atomic units (a.u.) = 60.48 fs, applying a delta-kick with a full width at half-maximum (FWHM) of 2.8 a.u. = 0.07 fs after 10 a.u. = 0.24 fs and using the EM propagator with a time step of 0.4 a.u. = 0.0096 fs. The chosen FHWM excites modes until well above 50 eV, and we verified the approach by using a weak instant δ-like excitation only applied at the first time step, which yielded the same spectra in this region relevant for our study. The numerical stability was judged by conservation of energy, where we set a maximum of |ΔE| = 0.001 eV. The field amplitude was set to E₀ = 0.01 a.u. = 5.1 × 10⁹ V/m.

We start the validation of the results by inspecting the sensitivity of the lowest-lying singlet and triplet excitation to the basis set size. Figure 3 shows the singlet values for ΔE₁ = E₁(tight) − E₁(tight + aug2) both for RT- and LR-TDDFT. As noted in Refs. 73, the results obtained with the “tight + aug2” basis set are generally well converged when compared with much larger benchmark sets. One can directly observe that both the RT and LR method display the same sensitivity of the molecules to the choice of basis set. The average basis-set induced energy difference for all molecules, ΔE_1,RT − ΔE_1,LR, is around −0.001 eV, with minimum and maximum deviations of −0.022 and 0.016 eV, respectively. The overall characteristic shows significant differences, and we group the results as follows: small changes below 0.02 eV are observed for ten molecules, intermediate-level changes between 0.02 and 0.15 eV are seen in eight molecules, and high changes above 0.15 eV are visible in ten molecules with a maximum of around 0.6 eV for pyrrole. We suspect that molecules with their lowest excitation coming already close to ionization are most prone to basis-set errors. This is in line with previous studies^73,104 that identified the importance of augmentation basis functions to correctly describe electron affinities, i.e., to converge the lowest unoccupied molecular orbital (LUMO) level. To visualize possible correlations, we also included the LUMO energies (formally equivalent to negative electron affinities) for all molecules in Fig. 3 on a separate scale; these values were obtained from our PBE calculations with the “tight + aug2” basis set and should not be taken as accurate representations of electron affinities. The visible correlation between the basis-set sensitivity of the excitations and of the LUMO levels confirms the hypothesis that excited-state wave functions close to the ionization threshold, reaching out far into the vacuum, are the common cause for the requirement of the “tight + aug2” basis set to converge both the excitation energies and the electron affinities. We can conclude at this point that both the RT and LR methods show similar convergence properties with respect to basis set change from the “tight” to the “tight + aug2” set.

We present analogous data for triplet excitations in Fig. 4. It is found that the energy range is approximately halved compared to the singlet case; thus, the basis effect is smaller. We note that this observation only holds for the lowest triplet excitation, while higher-lying triplets, as can be seen in Fig. 5 showing the (unrefined) DSF spectra for ethene, are more strongly affected by the basis set. Some comparably larger differences between LR- and RT-TDDFT are found for formamide (0.03 eV) and acetamide (0.07 eV). In all these cases, the basis set sensitivity is slightly higher for LR-TDDFT compared to RT-TDDFT. The differences for the other 22 molecules are lower than 0.03 eV, and the average sensitivity of the LR-TDDFT peaks is only 0.004 eV larger than in RT-TDDFT. In summary, both methods show excellent agreement in both the lowest-lying singlet and triplet excitations. Again, we emphasize that the “tight + aug2” set was shown to correspond to good overall basis set convergence, compared to much larger benchmark basis sets, in Ref. 73.

Next, we show that the RT-TDDFT approach, in conjunction with the methods to extract peak positions and oscillator strengths described before, is capable of reproducing a whole range of excitations of a molecule from a single run. For this purpose, we selected the first 8–18 lowest-energy clearly identifiable excitations from the DSF spectra, each verified by both the fitting and Padé approximant techniques, if possible. The differences between both evaluation methods are generally in the meV range, and we do not provide a direct comparison. Nevertheless, we note again that combining both methods can be very helpful when dealing with dense and/or weak excitations in the spectra. For very close-lying excitations, assessment of transition dipole moments (e.g., via the method in Ref. 102) could further help identify excitations that correspond to one another in LR-TDDFT spectra.

Where peaks could be clearly resolved, we computed the oscillator strengths via the fitting technique. The differences between LR- and RT-TDDFT are found to be small on average (2.2% for the “tight” set and 9.15% for the “tight + aug2” set), but caution is warranted when excitations are nearly degenerate, e.g., occurring in benzene, as these can hardly be resolved via RT-TDDFT. Additionally, the fitting uncertainty is larger for smaller peaks, as can be seen in Fig. 6, where the relative differences between RT- and LR-TDDFT oscillator strengths (OS) are shown in dependence of absolute oscillator strengths. While most recorded excitations have low OS, i.e., between 0.005 and 0.2, the relative deviations are large in this range, whereas RT-TDDFT faithfully reproduces the oscillator strengths computed via LR-TDDFT for strong peaks in the DSF spectra. For smaller values, RT-TDDFT OS are on average higher than LR-TDDFT OS and more so for the “tight + aug2” basis set.

Next, we turn to the energies of the singlet and triplet excitations obtained from both methods. Figure 7 shows the distribution of the deviations for singlet values {E_i(RT) − E_i(LR)} of each molecule, i.e., the sets of energy differences between corresponding excitations, obtained via RT- and LR-TDDFT. As noted before, 8–18 excitations per molecule were considered for the evaluation. In the statistical analysis, the red bars indicate the mean value and the box-plot captures the range of the observed values (see figure caption for details). While for the total average over all molecules and excitations, RT-TDDFT peaks are smaller only about 0.007 eV than LR-TDDFT peaks (for both basis sets), a few larger individual differences are visible in the plot. Only ethene, butadiene, pyridine, pyrazine, s-tetrazine, and formaldehyde exhibit deviations above 0.05 eV. For the “tight” basis set, one transition in formaldehyde shows up as a noticeable outlier at −0.08 eV. When comparing results for both basis sets, we note that the overall deviations are comparable in size. Assessing the deviations for individual molecules, we observe that both the average deviations as well as their spread among the different electronic transitions appear correlated in the “tight” and “tight + aug” basis sets. This points to methodological differences unrelated to the basis set as the likely cause of the deviations.

The analogous evaluation for triplet excitations can be found in Fig. 8. The overall deviations are approximately three times larger compared to the singlet case. Averaged over all molecules, the RT-TDDFT results are 0.022 eV (0.020 eV) higher for the “tight” (“tight + aug2”) basis set. The energy range of the differences in the “tight” basis set is noticeably larger compared to the singlet case, 15 of 28 molecules exhibit differences larger than 0.05 eV. The largest average differences are seen for ethene and formaldehyde, both around 0.8 eV with individual differences up to 0.17 eV. These molecules show smaller deviations with the “tight + aug2” basis set, hinting toward a basis set effect. Overall, the characteristic is very similar, although larger than for the singlet excitations, for both basis sets, again indicating a methodological cause.

To validate the implementation of the velocity-gauge formalism, we choose the standard model of bulk silicon. We used the two-atom unit cell of bulk crystalline silicon with a lattice parameter of a = 5.402 Å. The “tight” basis set was employed to expand the electronic wave functions. In the case of a periodic system, the significance of diffuse augmentation functions is much diminished since there are no basis function tails that extend into a vacuum; we therefore do not need to consider augmentation functions here. The external potential was switched on at time t = 0 with a strength of 0.001 a.u. = 5.1 × 10⁸ V/m. The strength falls within the linear response regime. Such step function-like external potential will excite all the modes in the system. For the time propagation, we choose the exponential midpoint integrator with a time step of 0.1 a.u. = 0.0024 fs. In total, we propagated the system for 3000 a.u. = 72.6 fs. The electronic current density was saved during the simulation, and the dielectric function was obtained via the procedure described in Sec. IV A. The imaginary part of the frequency-dependent dielectric function corresponds to the optical absorption spectrum in a condensed-matter system.¹¹⁴ 

We also observed an artificial peak feature around 0 eV as has been discussed in the literature, which is believed to be a consequence of the velocity-gauge.^{49,67,115–117} In order to remove this artifact from the simulations, the data analysis was done in the following way: first we performed RT-TDDFT simulations with a step function like “kick” in the external vector potential as shown in Fig. 9 (top left), and the resulting current is shown in Fig. 9 (middle left). Then, Fourier transform of the current according to Eqs. (25) and (26) yielded the imaginary part of the dielectric function as shown in Fig. 9 (top right). In order to remove the artificial peak, we took the Fourier transform of the background current with the step function that, at time t = 0 fs, switches to the value of the RT-TDDFT current at infinity as shown in Fig. 9 (bottom left). The imaginary part of the dielectric function from this background current is plotted in Fig. 9 (middle right). The subtraction of the background current from that of the RT-TDDFT simulation removes the artificial peak around 0 eV in the imaginary part of the dielectric function as shown in Fig. 9 (bottom right).

We now discuss the frequency-dependent imaginary part of the dielectric function for bulk silicon, which was generated with the aforementioned techniques. The lattice parameter is 5.402 Å. The calculation is non-relativistic, although for production calculations, a scalar-relativistic formalism could be used without problems and should be used. The dielectric function calculated using the Lindhard formula at the RPA level^118,119 with a k-point grid of 48 × 48 × 48 is also considered for comparison. For the RT-TDDFT calculations, we used different k-grids of 16 × 16 × 16, 32 × 32 × 32, and 48 × 48 × 48 resolution to check the k-grid integration convergence.

As shown in Fig. 10 (left), with the Γ-centered k-point mesh, we observe two major artificial peaks around 4 and 6 eV for the 16 × 16 × 16 k-grid. For the dense 48 × 48 × 48 k-grid, those peaks disappear and one major peak around 3.8 eV is seen with a shoulder like peak at around 2.9 eV. As suggested by Ref. 67, we also used a shifted k-point mesh for the dielectric function calculation with an offset vector of (0.01, 0.45, 0.37) in units of the reciprocal cell. In this manner, the symmetry of the k-point mesh is broken. The calculated imaginary part of the dielectric function is shown in Fig. 10 (right) for k-grids of 16 × 16 × 16 and 32 × 32 × 32. With the shifted k-grid, the calculated dielectric function with 32 × 32 × 32 k-grid has similar quality as the non-shifted k-grid of 48 × 48 × 48. The shifting of the k-grids thus allows us to achieve a faster k-grid convergence by breaking the symmetry in RT-TDDFT calculations for the dielectric function.

Using crystalline silicon as a representative model system, we also examined convergence of some key numerical settings in the FHI-aims code. For silicon atoms, in the 2020 default settings, the tier1 basis set (including the minimal basis of free-atom functions, as well as one s, p, d, and f basis function each) is used for both the “light” and “intermediate” settings; the “intermediate” settings have a denser integration grid with 17 918 grid points per atom compared to the “light” settings of 5604 grid points per atom. The “tight” settings in FHI-aims for silicon atoms have the same integration grid as the “intermediate” settings, while the basis set includes an additional d and g basis function each. We calculated the imaginary part of the dielectric function for silicon with these three levels of precision and show them in Fig. 11. No significant dependence on the numerical settings is observed. This indicates that even the basis set and numerical integration used in the default (2020) “light” settings for Si in FHI-aims are well converged for the RT-TDDFT based dielectric function of crystalline silicon.

In the FHI-aims code, the core and valence electrons are treated on the same footing without employing approximations such as pseudopotentials or different types of basis functions in the core region. Similar to those RT-TDDFT methods that are based on Gaussian-type orbitals (GTOs),⁶⁰ the RT-TDDFT methodology in the FHI-aims code allows us to calculate core electron excitations without further approximations. While codes for GTO-based RT-TDDFT simulation often do not implement periodic boundary conditions for studying extended systems, our approach here can be applied also to solid-state systems like crystalline silicon carbide as well as to isolated molecules like water on an equal footing, as demonstrated below for calculations of x-ray absorption spectra.

The same computational procedure as for the (“valence”) optical absorption spectrum calculations can be used to obtain x-ray absorption spectra; the system is excited by a weak delta-like external field that induces broadband excitations. At the same time, the integration step size needs to be sufficiently small so that the Fourier transform from time to frequency space is able to capture high-energy core–valence excitations. Particularly for the core-electron excitations, one faces the problem of a large sensitivity of the excitation energies on the exchange–correlation (XC) approximation, and this often makes it necessary to artificially shift the resulting excitation spectra to compare to experimental ones. For interpreting experimental x-ray absorption near edge structure (XANES), the oscillator strength and relative excitation energies with respect to the lowest core-electron excitation peak are usually of greater interest, justifying the significant value of these calculations even with the ad-hoc energy shift necessary due to the XC approximation.

For the K-edge spectrum (i.e., excitation from the oxygen 1s orbital), we provide comparison to previous theoretical results³³ and to experimental data.¹²⁰ The theoretical results come from a real-time TDDFT simulation using the GTO basis set code NWChem.⁴⁸ Here, Lopata et al. used the B3LYP hybrid XC functional in combination with the Sapporo-QZP basis set,¹²¹ which includes additional diffuse basis functions. As FHI-aims is also able to employ hybrid XC functionals^87,122,123 (based on the adiabatic approximation), B3LYP is equivalently used in our work for comparison. We provide results not only for the FHI-aims “tight + aug2” basis set but also for the Sapporo-QZP set (adapted from Ref. 121 via the Basis Set Exchange,¹²⁴ using dense integration grids without cutoff), as used in the theory reference. The theoretical simulation results by Lopata et al. include ad hoc relativistic corrections to shift the overall spectrum, whereas we probed relativistic corrections via the specific atomic zero-order regular approximation (atomic ZORA) approach in FHI-aims as described in Ref. 69.

The RT-TDDFT simulations were performed for ∼72 fs in total, with a time step of 0.1 a.u. = 0.0024 fs and using an instant electric field kick at t = 0 of 0.001 a.u. = 5.14 × 10⁸ V/m. The spectrum was obtained by averaging the three RT-TDDFT simulations with the perturbing field for each of the three Cartesian directions.

Comparing to a non-relativistic RT-TDDFT calculation, we found that the relativistic shift of the excitation energy amounts to ∼0.8 eV, equally for both employed basis sets and in line with findings in Ref. 125.

Furthermore, RT-TDDFT simulations with the PBE XC approximation were found to result in a much higher shift of 28.9 eV but left the spectral structure nearly unchanged. Each theoretical spectrum has been shifted such that the lowest peak matches its experimental analog. The shift for our “tight + aug2” spectra is +15.4 eV, which is a bit smaller than the shift of +15.55 eV used in Ref. 33 (in the figure, we added another +0.15 eV such that the peaks align slightly better, i.e., 15.4 eV were originally provided). The peak positions and shift from our calculations with the Sapporo-QZP basis set shown in Fig. 12 are in excellent agreement with the reference, confirming the validity of our implementation. The “tight + aug2” basis set yields very comparable results to the Sapporo-QZP basis with some deviations for the third and fourth peak. In particular, the fourth peak appears at 0.3 eV lower energy with the “tight + aug2” basis set. Overall, a very good agreement is found, and the “tight + aug2” basis set yields excellent results in this regime.

As an example of an x-ray absorption spectroscopy (XAS) calculation for extended systems, we study crystalline silicon carbide (Si) in the 2H–SiC polytype.

The 2H polytype of SiC forms a hexagonal unit cell with four atoms. Following the reference RT-TDDFT calculation in Ref. 49, we use the experimentally determined structure with lattice parameters of a = 3.076 Å and c = 5.048 Å. A k-point mesh of 20 × 20 × 12 was used for the Brillouin zone sampling. Within FHI-aims, we use the intermediate setting with its default basis set; for silicon, this is a tier1 basis set (including up to f functions); for carbon, this is a tier1 basis set (including up to d functions) plus an extra f-orbital basis function from tier2. As mentioned previously, in FHI-aims, the free-atom derived “minimal” basis functions to describe core electrons are always included.

For the periodic system of 2H–SiC, we used the velocity-gauge formula of RT-TDDFT, calculating and recording the current during the simulation. We again use a step function-like external potential switched on at t = 0 to study the dielectric response of the system in the linear-response regime. The strength of the external potential is 0.001 a.u. = 5.1 × 10⁸ V/m. We used a time step of 0.01 a.u. = 0.0024 fs for a total simulation length of 400 a.u. = 9.7 fs. The system was only excited in the z-direction in order to capture the zz-component of the dielectric function. Similar to the case of crystalline silicon discussed above, we again remove the artifact background due to the Fourier transform by subtracting the dielectric function calculated from the step function current as described in Sec. IV B.

We compare the zz-component of the dielectric function calculated by Ref. 49 via the RT-TDDFT functionality of the ELK code,⁶⁸ as plotted in Fig. 13. Note that in Ref. 49 the results calculated by the SIESTA code are compared to the results from ELK. The two codes in the reference generate quite similar results. We choose the results from the ELK code for the comparison here since ELK is another all-electron code. Since FHI-aims offers the possibility to calculate the dielectric function from the RPA-level Lindhard formula, we included this result for comparison in Fig. 14. We plotted the zz-component of the dielectric function at three energy windows, each 20 eV wide: the energy window of valence electron excitation (0–20 eV), the energy window of Si L-edge electron excitation (87–107 eV), and the energy window of the C K-edge electron excitation (258–278 eV).

For the valence electron excitation region, as shown in Fig. 13 (left) and Fig. 14 (left), the imaginary part of the dielectric function has two peaks at the energy of 6.4 and 7.4 eV; the peak strengths and peak positions are consistent between the FHI-aims code and the reference data from the ELK code; in addition, the results from RT-TDDFT and from the Lindhard formula are consistent with each other within the FHI-aims code. For the Si L-edge electron excitation region, as shown in Fig. 13 (middle) and Fig. 14 (middle), the main feature in the imaginary part of the dielectric function starts at 92 eV and remains a plateau-like feature until 107 eV; comparing to ELK’s result, FHI-aims has some features below the edge at 92 eV, while the RT-TDDFT and Lindhard formula in FHI-aims are consistent with each other. The features below the edge may reflect transitions in a different energy range altogether [valence to high-lying empty conduction states (see below)], and such transitions might also account for other differences seen between ELK and FHI-aims for the Si L-edge. For the C K-edge electron excitation region, as shown in Fig. 13 (right) and Fig. 14 (right), the feature in the imaginary part of the dielectric function starts at 264 eV and reaches a peak at around 272 eV. The results are consistent among the FHI-aims code with RT-TDDFT, the Lindhard formula, and the reference result by the ELK code.

We further evaluate the numerical settings convergence for the 2H–SiC core level excitation calculations. Similar to our discussion of crystalline silicon, we tested the default “light,” “intermediate,” and “tight” settings in FHI-aims. As shown in Fig. 15, for the valence electron excitation region and the Si L-edge excitation region, little difference is observed for the imaginary part of the dielectric function when comparing the different numerical settings. For the C K-edge region, the imaginary parts of the dielectric functions calculated by the “light” and “intermediate” settings are consistent with each other. However, for the “tight” setting, the imaginary part of the dielectric function has features below 264 eV and extra peaks show up at around 268 eV. As noted by Pemmaraju et al.,³⁴ this is caused by a decaying background of high-energy valence excitations. To confirm this explanation, we performed Lindhard formula dielectric function calculations with “intermediate” and “tight” settings. An additional Lindhard formula dielectric function calculation was therefore performed with the “tight” settings, summing up only the orbital pairs excited from the C K-shell. From Fig. 16 it is apparent that different backgrounds exist for “intermediate” and “tight” settings, especially at the energy region above 110 eV. In the C K-shell-only Lindhard calculation, the imaginary part of the dielectric function is zero below 264 eV and its shape is identical with the one from the “intermediate” settings. This confirms that the extra feature seen in the “tight” settings originates from the decaying background of the high-energy valence excitations. Thus, some caution is warranted regarding straightforward RT-TDDFT core level absorption simulations with large basis sets that also introduce high-lying unoccupied states in energy ranges comparable with core level binding energies. Given that “light” and “intermediate” settings already yield a converged imaginary part of the dielectric function, one qualitative remedy is to use these smaller basis sets for the study of core level spectra to avoid the interference from the decaying background of the valence to high-energy unoccupied-state excitations.

In order to validate our implementation of it-TDDFT, we first examine two typical systems—crystalline silicon and the benzene molecule.

For these systems, a conventional SCF method can be used to obtain the ground state easily in order to compare with the SCF solution obtained using the it-TDDFT approach. We used the “light” setting with tier1 basis sets and the LDA functional for these tests. The crystalline silicon has a unit cell with two silicon atoms, and the k-point grid was set to 8 × 8 × 8. As shown in Fig. 17, for both systems, the absolute energy differences between the fully converged it-TDDFT method and the SCF method are within 10⁻⁸ eV. For crystalline silicon, the total energy converges in 3 fs of imaginary time. For the benzene molecule, the total energy converges in 1 fs of imaginary time. These two examples show that, within the all-electron NAO framework, it-TDDFT can be employed for both isolated molecules and periodic systems.

One advantage of using it-TDDFT over other conventional iterative SCF methods is its guaranteed convergence.⁸¹ Regarding the difficult-to-converge situations/systems, chemical bond formation/dissociation presents a particular challenge because the excited state and the ground state often come quite close in energy. Here, we apply it-TDDFT to two challenging cases: one is a fluoromethane molecule with an elongated C–F bond of 4 Å and the other is a Na–Cl dimer with the Na–Cl distance of 6 Å, as shown in Fig. 18. Such configurations may be encountered in chemical reactions, and they are also employed as model systems for the theoretical development of ground state methods. The default setting for SCF in FHI-aims is the Pulay density mixing scheme^126,127 with the mixing coefficient adjusted by an estimated presence or absence of a bandgap at the end of the second SCF iteration, after initialization by a superposition of atomic densities. These default SCF settings usually converge well for a wide range of materials. However, with the default setting of SCF in FHI-aims, the two test systems mentioned above do not converge to a stationary state. It is quite possible that other SCF settings (away from the default) would yield one or more stationary solutions, but the search for non-default SCF convergence parameters is typically time-consuming and not desirable. With the it-TDDFT approach, the ground state is obtained successfully, although longer imaginary times are required for convergence than in the simple test systems of crystalline silicon and the benzene molecule. For the fluoromethane molecule with an elongated C–F bond, as shown in Fig. 18(a), the absolute energy difference is converged after an imaginary time of 60 fs. For the Na–Cl dimer with the Na–Cl distance of 6 Å as shown in Fig. 18(b), the absolute energy difference is converged after the imaginary time of 4 fs. Nevertheless, convergence can be achieved, which (in some circumstances, especially for more complex systems that are difficult to treat outrightly by multireference approaches) may offer an alternative pathway to circumvent technical difficulties that are otherwise difficult to overcome.

To validate it-TDDFT for a hard-to-converge periodic system, we test the examples of a boron vacancy and a nitrogen vacancy, respectively, in hexagonal boron nitride single layer, discussed, e.g., in Ref. 128. With a zero occupation broadening, a straightforward SCF cycle converges for the boron vacancy but not for the nitrogen vacancy. For the it-TDDFT method, as shown in Fig. 19, both systems converge. For the boron vacancy, the absolute energy difference between the SCF and it-TDDFT method is also around 10⁻⁸ eV, similar to the difference observed for crystalline silicon and benzene. For the hard-to-converge boron vacancy case, the it-TDDFT takes 60 fs imaginary time to converge. Interestingly, we observe a plateau in the imaginary time range of 0–20 fs. This might be an indication that the electronic solution is intermittently trapped in a meta-stable electronic state, although we did not pursue this point.

For the three cases where the conventional SCF calculation converged, namely, benzene, crystalline silicon, and hexagonal boron nitride with nitrogen vacancy, we further compared the convergence of both methods by means of the individual Kohn–Sham eigenvalues. For it-TDDFT, the final imaginary-time propagated Hamiltonian was diagonalized, whereas the standard SCF procedure always yields the KS eigenvalues by nature of it. The individual final SCF/it-TDDFT KS eigenvalues only differ by no more than 2.0 × 10⁻⁹ eV in all cases, underlining the desired convergence behavior.

In conclusion, the it-TDDFT approach implemented within the present NAO framework offers an alternative pathway to achieve stationary electronic solutions for otherwise difficult molecular and periodic systems, in line with the literature.^81–83 

As stated earlier, LR-TDDFT is the method of choice for the calculation of optical spectra for small and medium sized systems, as the common implementations are very efficient in such cases. RT-TDDFT is thus not a frequently used alternative. However, due to the formal $O(N6)$¹²⁹ scaling of the LR-TDDFT approach (with N being the number of atoms) and its potentially huge memory demand, RT-TDDFT is a potentially more viable approach for larger systems (although several approximate LR-TDDFT methods exist to tackle larger systems, see, e.g., Refs. 129–133).

The computational scaling of our RT-TDDFT implementation mainly depends on real-space operations, such as the electron density update or matrix element integrations, and the computation of the propagator, i.e., matrix operations (such as multiplication and diagonalization). For small systems like those considered in the molecular benchmark section, real-space operations take about 97% of the total computation time as the involved matrices are of small size, e.g., 582 × 582 (basis functions) for the largest molecule in the test set, naphtalene. The $O(N3)$ matrix operations become much more dominant with system size, compared to real-space operations that asymptotically scale as $O(N)$⁠.⁶⁹ With regard to possible applications of our RT-TDDFT approach for the calculation of optical absorption spectra (or other useful applications) for large systems, we here provide an analysis of the numerical scaling.

To provide a systematic benchmark, we employ oligoacenes of different chain lengths, i.e., the number of benzene rings for our calculations (see Table II). We used the “tight” basis set in combination with the PBE exchange–correlation functional for the calculations, as these choices may reflect rather prototypical settings. The calculations were not spin polarized, but the cost for this modification can be estimated to roughly double. We measure the effort by the wall clock time for one real-time propagation step, which consists of two applications of the propagator and real-space operations (updating density, Hartree potential and Hamiltonian) using the semi-implicit predictor-corrector scheme.

The overall simulation time to calculate a sufficiently well resolved absorption spectrum depends on the setup and on the possible use of refining techniques as presented before, but we assume that a minimum propagation time of 40 fs = 1666 a.u. is typical. Further assuming a time step of 0.3 a.u. = 0.0072 fs, this estimate would require a total amount of around 5550 real-time propagation steps, which we will use as reference below. Figure 20 shows single-step timings for our benchmark set, computed with 32, 64, and 128 parallel tasks (Intel Xeon E5-2670, 2.6 GHz CPUs) on the OCuLUS cluster of the Paderborn Center for Parallel Computing.

We further decompose the total single-step time into contributions due to the real-space operations, the linear solver, and the remaining matrix operations. Timings for output (e.g., total energy) are not explicitly mentioned as these are negligible. The presented data show several principal characteristics: while the real-space operations are here always dominant over the propagator-related operations (with a factor of >10 to ∼2 from smaller to larger molecules), the latter show a non-linear increase, while the former show a clear kink between n = 40 and n = 80. This can be explained by the code switching from an orbital-based to a density-matrix-based updating method of the electron density, the latter scaling as $O(N)$ for large systems.⁶⁹ A linear fit for the real-space integration scaling (n ≥ 40) and a cubic fit for the propagation scaling yield a crossover at ∼n = 136, corresponding to an oligoacene molecule with around 800 atoms. Nevertheless, the scaling with respect to parallel tasks is very good, i.e., doubling the number of parallel tasks approximately cuts the computation time in half in our benchmark. We note that for particularly costly real-space integrations (for dense and extended radial grids, e.g., when using augmented basis sets), the effort for time propagation will be comparatively much smaller.

Using a sufficient number of parallel tasks ≥128, one would be able to compute absorption spectra also for the largest molecule with nearly 500 atoms within a few days or less using the criteria defined above. It should be noted that memory demand is not a problem here due to the use of distributed storage.

We here describe and assess an all-electron, full-potential, real- and imaginary-time TDDFT implementation based on numeric atom-centered basis functions, building on the framework of the FHI-aims electronic structure code. The applicability of the approach to a range of applications is demonstrated, covering both molecular (non-periodic system geometries) and solid-state (periodic system geometries) examples and addressing the time-dependence of both core and valence excitations. We benchmark the method’s numerical precision, its efficiency, and basis-set related characteristics.

As a benchmark of valence electron excitations in non-periodic systems, we compare to results for Thiel’s test set¹⁰³ of small molecules, showing that the present RT-TDDFT implementation achieves numerical convergence with the same FHI-aims NAO basis sets that were shown to be appropriate for linear-response TDDFT calculations.⁷³ The analysis covers both singlet and triplet excitations, as well as the extraction of optical absorption spectra from real-time TDDFT simulations.

Periodic systems are covered by real-time velocity gauge simulations, which are shown to accurately reproduce reference data from the literature, specifically, the imaginary dielectric function of silicon.

The all-electron framework furthermore gives access to near edge x-ray absorption spectra via RT-TDDFT. In this context, we also show the capability of our implementation to incorporate hybrid exchange–correlation functionals. The oxygen K-edge spectrum of water is accurately reproduced compared to GTO-based RT-TDDFT reference data. We show that the native FHI-aims basis sets are well-suited for this task. We also provide core level absorption spectra for the 2H polytype of silicon carbide by performing real-time simulations with periodic boundary conditions and show that the results agree very well with earlier all-electron RT-TDDFT reference results.

We also implemented imaginary-time TDDFT functionality, applicable to both molecular and periodic systems, since the it-TDDFT approach offers an alternative way to obtain converged ground-state solutions of the time-independent Kohn–Sham equation. As shown in earlier work by Flamant et al.,⁸¹ the method is very easy to use and always converges, given an appropriate time step.

Finally, we delivered an analysis of our implementation’s computational scalability, showing that real-time TDDFT simulations can be carried out efficiently also for large systems with hundreds of atoms. While the established framework is already optimized, it is very likely that even better scaling can be achieved by incorporating customized numerical linear algebra methods.

We summarize the features and advantages of our implementations in the following condensed list.

• The all-electron NAO methodology facilitates the simulation of core and valence electron dynamics on equal footing.

• Both finite and periodic systems can be simulated using the highly accurate and efficient FHI-aims framework.

• As a further advantage, the all-electron description offers added flexibility in the choice of the gauge. Having a strictly local potential in the Hamiltonian is particularly convenient when performing velocity gauge simulations of periodic systems in radiation fields at optical wavelengths.

Even in view of this progress, one should still be aware of the remaining shortcomings of present-day implementations of the RT-TDDFT methodology. Some are rooted in inherent approximations of the functional itself and require basic research to overcome them, while others may be lifted by further improvement of numerical implementations. As an example for the first kind, it is well known that, for periodic systems with external electromagnetic fields, gauge invariance calls for a TD-current-DFT (TDCDFT) formulation rather than a functional using the instantaneous charge density only. This could also offer a route to lift the restrictions of the adiabatic exchange–correlation functionals that are implemented in present-day RT-TDDFT codes. As examples for the second type of issues that could and should be addressed in future code development, we list the following:

• Fully relativistic treatment of core electrons: While this is very important for heavy elements, the basic methods for its implementation are known and have been worked out for molecular systems.^134,135 A quasi-four-component relativistic implementation exists in the FHI-aims code,¹³⁶ although not yet for RT-TDDFT.

• Added flexibility of the basis set: Depending on the type of electronic excitations to be studied, localized basis sets come with inherent limitations, e.g., when electrons are promoted into continuum states. While the NAO basis functions in FHI-aims offer a very good description of the ground-state density and wavefunction, a higher spatial resolution (at least in some regions) may be required for some excited-state properties. A key point is the accurate description of the density difference between the ground state and the excited state. This could be achieved by introducing additional NAOs, but most generally by some unbiased basis functions, and ultimately by an additional real-space grid if required. This issue leaves ample room for future code improvement.

This work was partially funded by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) via the CRC1242 “Non-Equilibrium Dynamics of Condensed Matter in the Time Domain” (Project No. 278162697). Several calculations were performed on the OCuLUS cluster of the Paderborn Center for Parallel Computing. This work used the Extreme Science and Engineering Discovery Environment (XSEDE),¹³⁷ which is supported by the National Science Foundation (Grant No. ACI-1548562). This work used the Extreme Science and Engineering Discovery Environment (XSEDE) Expanse at the San Diego Supercomputer Center through Allocation No. TG-DMR200077.

V. B. is a member of the executive board of MS1P e.V., the non-profit organization that licenses the FHI-aims electronic structure code used in this work. V. B. does not receive any financial gains from this position. All other authors have no conflicts to declare.