The goal of quantum-chemical calculations is the simulation of materials and molecules. In density-functional theory (DFT) the first step along this line is obtaining the electron density minimising an energy functional. However, since energies and the density are usually not very tractable quantities in an experimental setup, comparison to experiment and scientific intuition also requires the computation of properties. Important properties include the forces (i.e. the energetic change due to a displacement of the structure) polarisabilities (change in dipole moment due to an external electric field) or phonon spectra (which can be measured using infrared spectroscopy). Therefore an efficient and reliable property computation is crucial to make quantum-chemical simulations interpretable and to close the loop back to experimentalists.

In DFT property calculations are done using density-functional perturbation theory (DFPT), which essentially computes the linear response of the electronic structure to the aforementioned changes in external conditions (external field, nuclear displacements etc.). Solving the equations underlying DFPT can become numerically challenging as (especially for metallic systems) the equations are ill-conditioned.

In a collaboration with my former PostDoc advisor Benjamin Stamm and my old group at the CERMICS at École des Ponts, including Eric Cancès, Antoine Levitt, Gaspard Kemlin, we just published an article, where we provide a more mathematical take on DFPT. In our work we provide an extensive review of various practical setups employed in main-stream codes such as ABINIT and QuantumEspresso from a numerical analysis point of view, highlighting the differences and similarities of these approaches. Moreover we develop a novel approach approach to solve the so-called Sternheimer equations (a key component of DFPT), which allows to make better use of the byproducts available in standard SCF schemes (the algorithm used to obtain the DFT ground state). With our approach we show savings up to 40% in the number of matrix-vector products required to solve the response equations. Since these are the most expensive step in DFPT this implies a similar saving in computational cost overall. Naturally our algorithm has been implemented as the default response solver in our DFTK code, starting from version 0.5.9.

Most of this work was done during a two-month visit of Gaspard Kemlin with Benjamin and myself here in Aachen. I think I speak for the both of us when I say that it has been a great pleasure to have Gaspard around, both on a professional as well as a personal level.

The full abstract of the paper reads

Response calculations in density functional theory aim at computing the change in ground-state density induced by an external perturbation. At finite temperature these are usually performed by computing variations of orbitals, which involve the iterative solution of potentially badly-conditioned linear systems, the Sternheimer equations. Since many sets of variations of orbitals yield the same variation of density matrix this involves a choice of gauge. Taking a numerical analysis point of view we present the various gauge choices proposed in the literature in a common framework and study their stability. Beyond existing methods we propose a new approach, based on a Schur complement using extra orbitals from the self-consistent-field calculations, to improve the stability and efficiency of the iterative solution of Sternheimer equations. We show the success of this strategy on nontrivial examples of practical interest, such as Heusler transition metal alloy compounds, where savings of around 40% in the number of required cost-determining Hamiltonian applications have been achieved.

Article
Eric Cancès, Michael F. Herbst, Gaspard Kemlin, Antoine Levitt and Benjamin Stamm. Numerical stability and efficiency of response property calculations in density functional theory. Letters in Mathematical Physics, 113, 21 (2023). DOI 10.1007/s11005-023-01645-3 [code]