Time-dependent density functional theory using atomic orbitals and the self-consistent Sternheimer equation

We present the implementation of linear-response time-dependent density functional theory based on the self-consistent Sternheimer equation and employing a basis set of numerical pseudo-atomic orbitals. We demonstrate this method by presenting test calculations on systems of increasing size ranging from benzene to chlorophyll $a$, and by comparing our results with those obtained within Casida's formalism and with previous calculations. We provide a detailed assessment of the accuracy of this method, both in relation to the use of local orbitals for describing electronic excitations and to the handling of the frequency response using Padé approximants. We establish a simple criterion for estimating a priori the accuracy of the basis set in the calculation of optical spectra. We show that the computational cost of this method scales quadratically with the system size.

Figure 1

(a) Isotropic average $\overline{\alpha }$ of the polarizability tensor of benzene (in units of bohr${}^3$) as a function of photon energy. Our calculation (blue solid curve) is compared to experiment (dashed black line) [61] and to previous calculations using plane waves (dashed-dotted red curve) [13] or real-space methods (dashed-double dotted purple curve) [11]. All curves are scaled so as to obey the $f$-sum rule, except the calculation of Ref. [11], which was scaled in order to match the $\pi$-${\pi }^{*}$ peak. Our calculations include a Lorentzian broadening of 0.3 eV. (b) Joint density of states of benzene calculated using a DZP basis (blue solid curve), a 5Z5P basis (purple dotted curve), and a plane-wave basis set (dashed red curve). The inset shows the corresponding density of states for the DZP and plane-wave basis. The vertical dashed line indicates the HOMO.

Figure 2

Isotropic average $\overline{\alpha }$ of the polarizability tensor of benzene (in units of bohr${}^3$) as a function of photon energy. We compare our calculations performed using a DZP basis and the self-consistent Sternheimer approach of Sec. 1a2 (black solid line) with the Casida approach of Sec. 1a1 (blue dashed line).

Figure 3

(a) Isotropic average $\overline{\alpha }$ of the polarizability tensor of benzene from experiment (dotted black line) [61] as shown in Fig. 1 and calculated using various basis sets: DZP (blue solid line), TZP (purple dashed line), TZDP (dashed-dotted red line), and 5Z5P (dashed-double-dotted green line). (b) Joint density of states of benzene calculated using the same basis sets as in (a). The arrow indicates that the onset of the JDOS is aligned to that of the absorption spectrum in (a). The key for the curves is the same as in (a). A Lorentzian broadening of 0.3 eV has been used in both cases.

Figure 4

(a) Calculation of the component ${\alpha }_{xx}\left(\omega \right)$ of the polarizability of benzene using the self-consistent Sternheimer equation: spectrum evaluated explicitly near the real axis (black solid line), and approximate analytic continuation using Padé functions (blue dashed line). Both spectra correspond to complex frequencies with a small imaginary component $\hslash \kappa =0.3$ eV. Inset: the blue dots represent the seed frequencies for the Padé functions in the complex plane. (b) Result obtained with different seed frequencies for the Padé functions. The different values for the imaginary part of the seed frequencies $\kappa$ are indicated in eV. The bottom-most (blue) curve is the same as the Padé curve in (a). The middle (purple) curve uses the same imaginary part, but the real part of the seed frequencies is shifted by 0.1 eV.

Figure 5

Number of iterations required to achieve convergence in the self-consistent Sternheimer approach for benzene. The blue disks correspond to calculations at complex frequencies $\tilde{\omega }=\omega +i\kappa$ with $\hslash \kappa =3$ eV. The black diamonds correspond to $\hslash \kappa =0.3$ eV. The resulting spectra are those shown in Fig. 4.

Figure 6

Component ${\alpha }_{xx}\left(\omega \right)$ of the polarizability of benzene obtained by solving the Sternheimer equation for different numbers of seed frequencies equally spaced in the interval $\hslash \omega =0$–40 eV. The spectra were calculated for $\tilde{\omega }=\omega +i\kappa$ with $\hslash \kappa =3$ eV and continued near the real axis by means of Padé approximants. The number of seed frequencies used in the approximants is indicated next to each curve.

Figure 7

The complex zeros ${\tilde{\omega }}_0$ of the denominator of Padé approximants corresponding to the spectra shown in Fig. 3. The black diamonds were obtained by using Padé functions containing only even powers, the blue circles were obtained by imposing only Eq. (28) on the seed frequencies. Inset: zoom of the figure around $\tilde{\omega }=0$.

Figure 8

(a) Calculated isotropic average $\overline{\alpha }$ of the polarizability tensor of C${}_{60}$ fullerene. The solid blue line is our calculation using the self-consistent Sternheimer equation and Padé approximants. The purple dashed-dotted line is the result of Ref. [30] using propagation in real time and the same DZP basis adopted by us. The dashed red line is a plane-waves calculation from Ref. [14]. (b) Comparison between the JDOS of C${}_{60}$ fullerene as calculated using local orbitals and a DZP basis set (solid blue line) and a plane-waves code (dashed red line) [62]. Both joint densities of states were shifted in order to match the onset of the absorption spectrum in (a), as indicated by the arrow.

Figure 9

Isotropic average $\overline{\alpha }$ of the polarizability tensor of chlorophyll $a$. (a) Comparison between the spectra calculated using the Sternheimer equation and Padé approximants (black dashed line) and the Casida method (blue solid line). For the Padé approximants, we used 200 seed frequencies $\tilde{\omega }=\omega +i\kappa$ with $\hslash \omega$ uniformly distributed between 0 and 10 eV, and $\hslash \kappa =2$ eV. Both spectra correspond to an artificial broadening $\hslash \kappa =0.3$ eV. Inset: schematic representation of the structure of chlorophyll $a$. (b) Comparison of our calculation based on the Casida method (black solid line) with previous calculations and with experiment: a calculation using plane waves (dash-dotted red line) [13], one using onetep (dashed-double dotted purple line) [15], and experiment (dashed blue line) [67].

Figure 10

(a)–(d) Isosurfaces of the frontier orbitals of chlorophyll $a$: (a) $\mathrm{HOMO}-1$, (b) HOMO, (c) LUMO, (d) $\mathrm{LUMO}+1$. The phytyl side chain C${}_{20}$H${}_{39}$ has been truncated for clarity. (e)–(h) Variation of the electronic charge density $\Delta n(\mathbf{r},\omega )$ corresponding to the optical transitions indicated in Fig. 9. (e) Q${}_x$ transition, (f) Q${}_y$ transition, (g), (h) the first Soret transition obtained with an external field that is polarized in the $x$ direction (g) and in the $y$ direction (h). [Note: Figures were rendered using the Visual Molecular Dynamics (VMD) program [68].]

Figure 11

Timing of the inversion of the linear system in Eq. (14) for the model polythiophene chains, as a function of the number of atoms (purple disks). The left arrow indicates that values are shown on the left $y$ axis. Timing of the potential update in Eq. (25) is shown as blue diamonds, for the same models. The right arrow indicates that values are shown on the right $y$ axis. The dashed curves are parabolic fits to the data points, showing the $\mathcal{O}\left(N^2\right)$ scaling. The timing refers to calculations on a single core. Schematic representations of the chains are given for clarity.