Abstract

We introduce a new implementation of time-dependent density-functional theory which allows the entirespectrum of a molecule or extended system to be computed with a numerical effort comparable to that of a single standard ground-state calculation. This method is particularly well suited for large systems and/or large basis sets, such as plane waves or real-space grids. By using a superoperator formulation of linearized time-dependent density-functional theory, we first represent the dynamical polarizability of an interacting-electron system as an off-diagonal matrix element of the resolvent of the Liouvillian superoperator. One-electron operators and density matrices are treated using a representation borrowed from time-independent density-functional perturbation theory, which permits us to avoid the calculation of unoccupied Kohn–Sham orbitals. The resolvent of the Liouvillian is evaluated through a newly developed algorithm based on the nonsymmetric Lanczos method. Each step of the Lanczos recursion essentially requires twice as many operations as a single step of the iterative diagonalization of the unperturbed Kohn–Sham Hamiltonian. Suitable extrapolation of the Lanczos coefficients allows for a dramatic reduction of the number of Lanczos steps necessary to obtain well converged spectra, bringing such number down to hundreds (or a few thousands, at worst) in typical plane-wave pseudopotential applications. The resulting numerical workload is only a few times larger than that needed by a ground-state Kohn–Sham calculation for a same system. Our method is demonstrated with the calculation of the spectra of benzene, fullerene, and of chlorophyll a.

Received 09 January 2008Accepted 27 February 2008Published online 16 April 2008

Acknowledgments:

S.B. wishes to thank Renata Wentzcovitch for hospitality at the University of Minnesota Department of Chemical Engineering and Materials Science and vLab, where some of the work reported in this paper was started in summer 2005 and the Leverhulme trust for a visiting professorship at the Department of Earth Sciences of the University College London, where this paper was written in 2007. R.G. and S.B. thank Brent Walker and Marco Saitta for participating in the early developments of the ideas on which this paper is founded. D.R. gratefully acknowledges useful conversations with Giuseppe Pastori Parravicini and correspondence with Liana Martinelli, the former also bearing the responsibility for making S.B. aware of the power and beauty of Lanczos methods. Y.S. acknowledges support from NSF under Grant No. NSF/ITR-0428774, the U.S. Department of Energy under Grant No. DE-FG-03ER25585, and from the Minnesota Supercomputer Institute.

30.The batch representation of response functions has been rediscovered several times, and given several different names, in the quantum chemistry community, since it was introduced in the context of DFPT (Ref. 14). See, e.g., Refs. 17, 29, and 18.

44.Ultrasoft pseudopotentials are taken from the public QUANTUM ESPRESSO pseudopotential library. Pseudopotential files can be downloaded from URL http://www.pwscf.org/pseudo/1.3/UPF/xxx, where xxx = C.pbe-rrkjus.UPF, H.pbe-rrkjus.UPF, C.pw91-van_ak.UPF, H.pw91-van_ak.UPF, Mg.pw91-np-van.UPF.N.pw91-van_ak.UPF.o.pw91-van_ak.UPF.

Abstract

We introduce a new implementation of time-dependent density-functional theory which allows the entirespectrum of a molecule or extended system to be computed with a numerical effort comparable to that of a single standard ground-state calculation. This method is particularly well suited for large systems and/or large basis sets, such as plane waves or real-space grids. By using a superoperator formulation of linearized time-dependent density-functional theory, we first represent the dynamical polarizability of an interacting-electron system as an off-diagonal matrix element of the resolvent of the Liouvillian superoperator. One-electron operators and density matrices are treated using a representation borrowed from time-independent density-functional perturbation theory, which permits us to avoid the calculation of unoccupied Kohn–Sham orbitals. The resolvent of the Liouvillian is evaluated through a newly developed algorithm based on the nonsymmetric Lanczos method. Each step of the Lanczos recursion essentially requires twice as many operations as a single step of the iterative diagonalization of the unperturbed Kohn–Sham Hamiltonian. Suitable extrapolation of the Lanczos coefficients allows for a dramatic reduction of the number of Lanczos steps necessary to obtain well converged spectra, bringing such number down to hundreds (or a few thousands, at worst) in typical plane-wave pseudopotential applications. The resulting numerical workload is only a few times larger than that needed by a ground-state Kohn–Sham calculation for a same system. Our method is demonstrated with the calculation of the spectra of benzene, fullerene, and of chlorophyll a.