Fast construction of the Kohn--Sham response function for molecules

The use of the LCAO (Linear Combination of Atomic Orbitals) method for excited states involves products of orbitals that are known to be linearly dependent. We identify a basis in the space of orbital products that is local for orbitals of finite support and with a residual error that vanishes exponentially with its dimension. As an application of our previously reported technique we compute the Kohn–Sham density response function $\chi_{0}$ for a molecule consisting of $N$ atoms in $N^{2}N_{\omega}$ operations, with $N_{\omega}$ the number of frequency points. We test our construction of $\chi_{0}$ by computing molecular spectra directly from the equations of Petersilka–Gossmann–Gross in $N^{2}N_{\omega}$ operations rather than from Casida’s equations which takes $N^{3}$ operations. We consider the good agreement with previously calculated molecular spectra as a validation of our construction of $\chi_{0}$. Ongoing work indicates that our method is well suited for the computation of the GW self-energy $\Sigma=\mathrm{i}GW$ and we expect it to be useful in the analysis of exitonic effects in molecules.

💡 Research Summary

The paper presents a novel algorithm for constructing the Kohn–Sham non‑interacting response function χ₀ of molecular systems with a computational cost that scales as N² N_ω, where N is the number of atoms and N_ω the number of frequency points. The authors start from the well‑known problem that, within a linear‑combination‑of‑atomic‑orbitals (LCAO) framework, the set of orbital products φ_i(r) φ_j(r) is highly linearly dependent. This redundancy makes a direct representation of χ₀ prohibitively expensive, leading to the traditional O(N³) scaling of Casida’s TDDFT formulation.

To overcome this, the authors construct a compact, localized product basis. Each atomic orbital is assumed to have finite spatial support; the products of two such orbitals are projected onto a set of auxiliary functions that are orthogonalized and truncated by singular‑value decomposition. The truncation criterion is chosen so that the residual error decays exponentially with the size of the auxiliary basis, guaranteeing rapid convergence while preserving locality. Because the auxiliary functions are confined to the overlap region of the original orbitals, the resulting basis size grows only linearly with N, and matrix elements can be assembled with O(N²) effort.

With this basis, the non‑interacting response function is expressed as a matrix χ₀(μ,ν; ω) in the auxiliary space. The frequency dependence resides solely in the familiar denominator ω – (ε_a – ε_i) + iη, while the spatial part is captured by the pre‑computed product‑basis integrals. Consequently, evaluating χ₀ at all required frequencies costs O(M² N_ω), where M ∝ N is the dimension of the product basis, yielding the advertised O(N² N_ω) scaling.

The authors then bypass Casida’s eigenvalue problem entirely. Instead, they insert the newly built χ₀ into the Petersilka‑Gossmann‑Gross (PGG) linear‑response equation,