# Fast computation of the Kohn-Sham susceptibility of large systems

Dietrich Foerster

► To cite this version:

Dietrich Foerster. Fast computation of the Kohn-Sham susceptibility of large systems. Physical Review B: Condensed Matter and Materials Physics (1998-2015), 2005, 72 (7), pp.073106 (1-4). 10.1103/Phys-RevB.72.073106. hal-00004344

HAL Id: hal-00004344

https://hal.science/hal-00004344v1

Submitted on 1 Mar 2005

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons CC BY-SA 4.0 - Attribution - ShareAlike - International
License
Fast computation of the Kohn-Sham susceptibility of large systems

D. Foerster

CPMOH, Université Bordeaux 1

351, cours de la Libération

33405 TALENCE Cedex, France

(March 1, 2005)

For hybrid systems, such as molecules grafted onto solid surfaces, the calculation of linear response in time dependent density functional theory is slowed down by the need to calculate, in ~ N⁴ operations, the susceptibility of N non interacting Kohn-Sham reference electrons. We observe that this susceptibility can be calculated N times faster, with a precision of ~ (Δω/ε)², where Δω is the width of the frequency intervals and ε a displacement off the real frequency axis. By itself or in combination with previously known procedures of accelerating such calculations, our procedure should significantly facilitate the calculation of TDDFT response and optical spectra of hybrid systems.

*Motivation.* Time dependent density functional theory (TDDFT) has proved to be very successful in determining optical properties of molecules and clusters, for a review see [1]. Therefore it would be interesting, to study within TDDFT linear response, the simplest version of this theory, also hybrid materials, such as dyes bound to semiconductor or metal surfaces. Such materials exhibit interesting phenomena, for example charge transfer upon photo excitation in the Grätzel cell [2] or increased fluorescence for organic molecules chemisorbed onto noble metal clusters [3].

Density functional algorithms based on localized orbitals, such as the siesta package [4] have been used to study the ground state of systems as large as DNA chains [5], but the computation of optical spectra and properties of excited states of hybrid systems in such schemes appears still to be very difficult.

The present note is motivated by recent progress in the calculation of optical spectra in TDDFT linear response [6] using localized orbitals and the observation that even this recent advance is not sufficient to deal with hybrid systems. To model such systems one must deal with a slab of a few layers of solid onto which a molecule may be adsorbed, without the benefit of any symmetries, and with N, the number of orbitals, at least of the order of a few thousand. In the local orbital approach of [6] the needed CPU time grows rather steeply, as ~ N⁴, with N and paradoxically, most of the CPU time is actually spent on the charge susceptibility of the non interacting Kohn-Sham reference fermions.

Here we focus on the Kohn-Sham charge susceptibility and describe a simple way to calculate it N times faster than by straightforward computation, in ~ N³ rather than ~ N⁴ operations, with a relative precision of ~ (Δω/ε)², where Δω is the spacing of frequencies,

and we also study the performance of this technique by applying it to a random hamiltonian.

The equations of TDDFT linear response [8] are based on the definition of the time dependent Kohn-Sham potential $V_{KS}(r, t)$ of non interacting reference electrons:

$$V_{KS}(r,t) = V_{ext}(r,t) + \int dr' \frac{n(r',t)}{|r-r'|} + V_{xc}(n(r,t)) \quad (1)$$

where $V_{ext}$ is the external potential and $V_{xc}$ the exchange correlation potential which reduces to an ordinary function of the density $n(r,t)$ in the local and adiabatic approximation. Differentiating both sides of this equation with respect to the charge density $n(r,t)$ one finds a Dyson like equation for the interacting susceptibility $\chi_{int} = \frac{\delta n}{\delta V_{xc}}$:

$$\chi_{int}(rt, r't') = \frac{1}{\chi^{-1}(rt, r't') - \frac{\delta(t-t')}{|r-r'|} - f_{xc}(rt, r't')} \quad (2)$$

with a local exchange kernel $f_{xc} = \frac{dV_{xc}}{dn}$. The poles of $\chi_{int}$ as a function of frequency provide information on the excited states. The preceding equation requires as input the susceptibility $\chi$ of independent Kohn-Sham reference fermions

$$\chi(\omega, r, r') = \sum_{i,j} \frac{n_i - n_j}{\omega + i\epsilon - (E_i - E_j)} \phi_i^*(r) \phi_i(r') \phi_j(r) \phi_j^*(r') \quad (3)$$

where {$\phi_i(r), E_i$}, $i = 1..N$ is a set of N Kohn-Sham eigenfunctions and energies. From eq (3) we see that $\chi$ is a sum of ~ N² terms that must be computed for a set of $N_\omega$ frequencies and with $r, r'$ each ranging over $V/\Delta V$ cells in space, with V the volume of the system and $\Delta V$ an elementary grid volume. Because the number of orbitals N grows linearly with the volume V, the calculation of Kohn-Sham susceptibility $\chi$ then requires ~ N⁴$N_\omega$ operations, more than the remaining operations needed to calculate the interacting susceptibility $\chi_{int}$, the number of which grows with N only as ~ N³$N_\omega$.

*Accelerated Computation.* In LCAO schemes such as, for example, the siesta approach [4], one develops the fermion operators in terms of a set of orbitals {$\psi_a(r)$} of finite range. This implies a corresponding expansion for the Kohn-Sham density response function $\chi$ and one must now find

$$\chi_{abcd}(t) = \sum_{E,F} \psi_a^* \psi_b^E \psi_c^{*F} \psi_d^F \frac{n(E) - n(F)}{\omega + i\epsilon - (E - F)}. \quad (4)$$
Here the $\psi_a^E$ are eigenvectors of the Kohn-Sham hamiltonian in the orbital basis and the indices of $\chi_{abcd}$ each range over $N$ orbitals (which we take to be real, so we drop the complex conjugation). There appear to be far too many independent quantities $\chi_{abcd}$, but for orbitals of finite range locality limits the distance between, respectively, the orbitals $(a, c)$ and $(b, d)$ to a few atomic distances which is also quite natural as these pairs of orbitals correspond to coincident points in the continuum limit. Therefore, in a basis of localized orbitals, only a subset of $\sim N^2$ susceptibilities $\chi_{abcd}$ contributes in the linear response equations, see [6] for a more detailed discussion.

The Kohn-Sham susceptibility of free fermions is just a product (or convolution in frequency) of propagators,

$$ \chi_{abcd}(t) = -iG_{ab}(t)G_{cd}(-t) \quad (5) $$

so its calculation is essentially finished once we know the propagator of the free Kohn-Sham fermions. The determination of $G_{ab}(\omega)$ requires of the order of $\sim N^3$ operations per frequency (even if the eigenfunctions are already available) or $\sim NN_\omega$ operations for a given set of indices $\{a, b, c, d\}$ at all frequencies, and these operations totally dominate the low $\sim N_\omega \log N_\omega$ calculational cost of a fast convolution.

For this reason, and also because convolutions of Green functions require great care to avoid finite size effects, we give here a way of calculating $\chi_{abcd}(\omega)$ that is comparable in overall speed, but more accurate than a convolution. We decompose the susceptibility into positive and negative frequency components, $\chi_{abcd}(\omega) = \chi_{abcd}^{-}(\omega) + \chi_{cdab}^{-}(\omega)$ with

$$ \chi_{abcd}^{-}(\omega) = \sum_{E<0<F} \frac{\rho_{ab}^{E} \rho_{cd}^{F}}{\omega - i\epsilon - (E - F)} \quad (6) $$

where $\rho_{ab}^E = \psi_a^E \psi_b^E$ and rewrite $\chi_{abcd}^{-}(\omega)$ as

$$ \chi_{abcd}^{-}(\omega) = -\sum_{E<0} \rho_{ab}^{E} G_{cd}^{+}(-\omega + E) \quad (7) $$

$$ G_{cd}^{+}(\omega) = \sum_{F>0} \frac{\rho_{cd}^{F}}{\omega + i\epsilon - F} \quad (8) $$

where $G_{cd}^{\pm}(\omega)$ denotes the positive/negative frequency part of the propagator. We also discretize the frequency axis to a finite mesh of $N_\omega$ points which we choose to be uniform, for simplicity, and we interpolate $G_{cd}^{+}(\omega)$ linearly between the mesh points. Because $G_{cd}^{+}(\omega)$ varies on a frequency scale of $\epsilon$ this interpolation introduces a relative error in $\chi_{abcd}$ of the order of $(\Delta\omega/\epsilon)^2$ where $\Delta\omega$ is the frequency spacing. The sum in eq (7) can be represented as a convolution, $\chi_{abcd}^{-} = \frac{1}{\pi} \text{Im} G_{ab}^{-} \otimes G_{cd}^{+}(-\omega)$ that can be calculated much faster, but the result is less accurate and only useful in conjunction with an accelerated calculation of $G_{ab}^{\pm}$ itself.

After reconstituting $\chi_{abcd}$ from $\chi_{abcd}^{-}$ we have spent $\sim NN_\omega$ operations on $\chi_{abcd}$ for a given set of indices

$\{a, b, c, d\}$. As mentioned before, only $\sim N^2$ combinations of such indices enter into the equation of the interacting susceptibility, and so we have indeed reduced the computational cost of the Kohn-Sham susceptibility from $\sim N^4 N_\omega$ to $\sim N^3 N_\omega$ operations.

For a molecule bound to a slab of surface layers and similar systems that involve thousands of orbitals, this speed up by a factor $N$ could turn out to be essential for the feasibility of the computation.

**Test of the procedure.** To illustrate the quality of our approximate procedure as a function of the number of orbitals $N$, the number of frequency points $N_\omega$ and the smoothing parameter $\epsilon$, we choose a symmetric hamiltonian $h$ at random with a probability:

$$ \mathrm{prob}(h) \sim \exp(-N \sum_{i,k=1..N} h_{ik}^2) \quad (9) $$

from the so called "orthogonal ensemble" of random matrices. For $N \to \infty$, the density of eigenvalues of this matrix ensemble tends to a semicircular distribution in the interval $(-1, 1)$, see [9]. We therefore choose units in which the band ranges between $(-1, 1)$, set $E_f = 0$ and calculate $\chi_{abcd}$ (i) exactly and (ii) approximately in, respectively $\sim N^2 N_\omega$ and $\sim NN_\omega$ operations and compare our results.

With a reasonable choice of the values of $\epsilon$ and the number of mesh points $N_\omega$, the exact and approximate values of $\chi$ coincide to within the linewidth of the plot, as seen in figure (1). Figure (2) confirms that the error scales like $(\Delta\omega/\epsilon)^2$ with the mesh spacing. For $N = 2000$, the approximate calculation is about a thousand times faster than the exact one.

**Conclusion.** In summary, we have described an accelerated procedure for computing the charge response of Kohn-Sham reference fermions in order to facilitate the calculation of TDDFT linear response for systems that are intrinsically of large size, such as molecules bound to solid surfaces. Our speedup from $N^4 N_\omega$ to $N^3 N_\omega$ operations for calculating the subset of elements of $\chi_{abcd}$ that contribute in linear response comes at the price of a finite precision of the order of $\sim (\Delta\omega/\epsilon)^2$ where $\Delta\omega$ is the frequency spacing.

A reduction from $N^4 N_\omega$ to $N^3 N_\omega$ operations in TDDFT linear response was achieved before by limiting the density response to a finite manifold of fit functions [7]. Because our procedure is very simple it might be combined with such or similar methods to further reduce the growth of operations with $N$.

We expect our algorithm to facilitate the study, by TDDFT linear response, of optical properties of hybrid systems such as molecules docked on surfaces, systems that are difficult to treat by existing computational methods and which exhibit many interesting features that deserve further study.

**Acknowledgements.** I wish to thank Andrei Postnikov (Osnabrück) for useful discussions and continued correspondence and Xavier Blase (Lyon) for sending his
work on TDDFT linear response prior to publication. I
am indebted to Daniel Sanchez-Portal (San Sebastian)
for discussions on the siesta code and to Mark Casida
(Grenoble) for interesting comments on his approach to
linear response. Computer time was provided by IDRIS
(Orsay) and Regatta (Bordeaux) and financial support
by GDR-DFT. I thank Alois Würger for discussions on
the manuscript.

[1] See the reviews of E.K.U. Gross, M. Casida and other authors at the meeting on "Time-Dependent Density-Functional Theory and the Dynamics of Complex Systems", Santa Fe, 2004, http://campus.umr.edu/tddft/tddft.

[2] O. Regan and M. Grätzel, Nature **353** (1991) 737, M. Grätzel, Nature **414**, (2002) 338.

[3] G.K.Thomas, P.V. Kamat, J. Am. Chem. Soc. **122**, 2655 (2000).

[4] J.M. Soler, E. Artacho, J.D. Gale, A. Garcia, J. Junquera, P. Ordejon, and D. Sanchez-Portal, J. Phys. C14 (2004) 2745.

[5] D. Sanchez-Portal, P. Ordejon, E. Artacho, K. Soler, Int. J. Phys. Chem **65** (1997) 453.

[6] X. Blase and P. Ordejon, Phys. Rev. B69 (2004) 85111.

[7] M.E. Casida in "Recent Developments and Applications of Modern Density Functional Theory", edited by J.M. Seminario, Amsterdam 1996; S.J.A. van Gisbergen, J. G. Snijders and E. J. Baerends, Comp. Phys. Comm. **118** (1999) 119.

[8] M. Petersilka, U. J. Gossmann, and E. K. U. Gross, Phys. Rev. Lett. **76** (1996) 1212. For a review, see also M. A. L. Marques and E. K. U. Gross, Ann. Rev. Phys. Chem. **55** (2004) 427.

[9] E. Wigner, Ann. of Math. **67** (1958) 325.

**Figure Captions:**

Figure (1) : Reχ<sub>abcd</sub>, Imχ<sub>abcd</sub> as a function of fre-
quency for N = 2000 orbitals, ε = 0.03, N<sub>ω</sub> = 256
frequencies and for a random choice of {a,b,c,d}. The
difference between the exact and the approximate val-
ues of χ<sub>abcd</sub> is smaller than the linewidth of the figure,
therefore one can only see two distinct curves.

Figure (2) The relative error in χ_{abcd} for N = 2000 orbitals and ε = 0.03 as a function of the number of frequency points N_ω. The relative error is defined as the ratio of the root mean squares of error and signal and, as expected, it scales as ~N_ω^2 or ~(Δω)^2.