What are the core numerical routines used in density functional theory?

It looks to me to that many of the things one would do within the framework of density functional theory (DFT) ultimately boil down to solving a generalized eigenvalue problem $$Av = \lambda Bv$$ where $$A$$ and $$B$$ are matrices, $$\lambda$$ is the eigenvalue and $$v$$ is the corresponding generalized eigenvector.

What are the core numerical routines used in DFT, besides solving generalized eigenvalue problems? What about optimization routines?

• For Kohn-Sham DFT there is also the optimization of orbitals! – Nike Dattani May 6 '20 at 14:01
• @NikeDattani Can you write an answer? :D – edwinksl May 6 '20 at 14:05
• I know almost nothing about DFT, so I will have to do some reading. I know KS-DFT has an SCF procedure similar to Hartree-Fock, which means that a KS-DFT code could have a numerical routine akin to DIIS. But I don't know all the other numerical routines in a DFT calculation so my answer would not be thorough. Perhaps if the question specified that it was seeking 1 numerical routine per answer, I could give it a try, similar to this question which seeks 1 topic per answer: materials.stackexchange.com/questions/17/… – Nike Dattani May 6 '20 at 15:23

DFT is a nonlinear eigenvalue problem. In the community of chemistry and material science, it is almost always solved via the self-consistent-field method (SCF).

If the DFT is solved on the level of GGA or LDA (referring to Jacob's Ladder), then, within each SCF iteration, the most expensive numerical routine is indeed (generalized) eigenvalue problem. In addition to the linear (generalized) eigenvalue problem, the density/orbital mixing also involves cubic scaling numerical routines. To my knowledge, they are still much cheaper than the eigensolver.

However, if you further move upon Jacob's ladder to hybrid-functionals or above, then you will encounter numerical routines of much higher order scaling than the eigensolver. In practice, some approximation and numerical tricks can be used to reduce the computational cost of hybrid-functional to the same level as eigensolver. Hence the discussion is more complicated in this case.

Jumping out of the chemistry and material science community, some mathematicians have tried to solve the nonlinear eigenvalue problem directly from an optimization point of view and not uses SCF+mixing. These methods have only been applied to LDA-level functionals. Hence the core numerical routine varies from method to method.

• Very nice answer! Since you mentioned hybrids, I'll also mention that there's double-hybrids like B2PLYP, which for example involves MP2. – Nike Dattani May 6 '20 at 21:03

The answer to your question depends a little on the context you are talking about, most importantly the type of basis function you use. As I will describe especially this aspect crucially influences the mathematical structure of the problem you try to solve.

Structure of Kohn-Sham

Being a bit sloppy on the mathematical details (for more see e.g. section 4.6 of my thesis Kohn-Sham DFT and related problems such as Hartree-Fock are concerned with solving the following minimisation problem

$$\min_{\varphi_i}\left( E[\{\varphi_i\}_i] \right) \quad \text{under the constraint} \quad \int \varphi_i(r) \varphi_j(r) \text{d}r = \delta_{ij}$$

where one seeks to find the minimising set of orthonormal orbitals for the Kohn-Sham functional $$E[\{\varphi_i\}_i]$$. Its precise form depends, for example, on the type of exchange-correlation functional employed.

Now finding the minimum of this problem means as usual that we need to take the functional derivative of that equation and set it to zero. The resulting Euler-Lagrange equations are the famous Kohn-Sham equations

$$\hat{F}[\{\varphi_i\}_i] \varphi_i = \epsilon_i \varphi_i$$

involving the orbital-dependent Kohn-Sham operator $$\hat{F}[\{\varphi_i\}_i]$$. Notice that we are still in the infinite-dimensional regime of functions and operators at this point and notice further that $$\hat{F}$$ is actually the derivative of the energy functional with respect to the orbitals.

Self-consistent field

As Yingzhou Li already mentioned there are two ways to tackle this problem. The most common approach is the self-consistent field (SCF) procedure, where our angle of attack is the second equation. Introducing a basis the spectral problem of the operator becomes a generalised matrix eigenvalue problem

$$F[\{v_i\}_i] v_i = \epsilon_i S v_i.$$

Since $$F$$ depends on the eigenvectors $$\{v_i\}_i$$, the procedure is to first guess some $$\{v_i\}_i$$, build an initial $$F$$, diagonalise it to get new $$\{v_i\}_i$$ and repeat until convergence. Because the diagonalisation is done many times (order of 20 to 50) for a single ground-state DFT calculation this is where most time is spent. I should mention that I am glossing over details here. In practice you never follow this simple procedure but additional measures like "mixing", "damping" or various forms of preconditioning are employed, but they are not the time-dependent steps.

The structure of $$F$$ and $$S$$ depend on the basis you employ. If your basis is orthogonal (like properly normalised plane-waves) than $$S$$ is just the identity and can be dropped and you get a standard eigenvalue problem. If you use for example Gaussians, than $$S$$ needs to be kept. Also for plane waves you typically need lots of them, say like 1 million or more, whereas for Gaussians typically large bases are around 1000 functions. On the other hand the matrices $$F$$ and $$S$$ for Gaussians are pretty densely populated, whereas for plane-wave-discretised problems $$F$$ is very sparse and contains many zeros. (I am focusing on Gaussians versus plane waves here, but the picture is similar for other types of bases).

It is easy to imagine that storing a $$10^6$$ times $$10^6$$ matrix in memory is going to be a problem, whereas storing a $$1000$$ times $$1000$$ matrix is completely fine. As a result one usually does not use the same diagonalisation procedures for both type of basis sets: For Gaussians one assembles the matrix and uses a dense diagonalisation procedure from LAPACK, whereas for plane-waves one uses iterative methods like Davidson or LOBPCG, where the key step is to be able to apply $$F$$. That is form matrix-vector products $$F x$$ for trial eigenvectors $$x$$ that come up during the iterative procedure.

To summarise: For an SCF-based approach for Gaussians building $$F$$ is the key step, which primarily involves the computation of the two-electron integrals as the most expensive constituent of $$F$$. For plane waves, however, the application of $$F$$ is most crucial. Here the most costly steps is the computation of fast-fourier transforms, which are needed to transform between Fourier space (the space in which the elements of $$v$$ and $$F$$ live) and real space (the space in which the potential is usually stored).

Direct minimisation

The second approach to tackle Kohn-Sham is not to use a self-consistent field procedure, but instead to directly tackle the minimisation problem I mentioned first. For this one typically uses Newton or Quasi-Newton minimisation procedures to directly minimise the energy by varying the input orbitals (hence the name). These algorithms all require the derivative of the energy wrt. the variations, hence again $$F$$. But similar to iterative diagonalisations $$F$$ only needs to be applied during the Newton minimsation. Therefore these procedures can be used both for Gaussian as well as plane-wave basis sets and the most costly step is, as before, building $$F$$ or the application of it.

Regarding what Yingzhou Li said, there is no limit to using direct minimisation LDA I am aware of. In fact in DFTK we can easily do direct minimisation for both LDA and GGA functionals and I do not immediately see a reason why an extension to hybrids would not be possible, but we never tried it.

It should be said, however, that some restrictions do apply to direct minimisation. In the typical implementations it is assumed that only the orbitals need to be optimised and not the occupation numbers. For systems with zero or small gaps (such as metals) direct minimisation is therefore not suitable.

• Very good answer! Thanks so much Michael! – Nike Dattani May 10 '20 at 14:10
• We use direct minimisation in CASTEP as a robust fall-back for cases where SCF mixing methods are fragile. It is typically much more expensive per iteration, but also much more robust. There are no particular issues when moving from LDA to GGA, meta-GGA, Hubbard U, Fock exchange etc. – Phil Hasnip May 15 '20 at 1:05
• @PhilHasnip can you expand a bit on that? Why do you say that it is more expensive per iteration? It's about the same cost (applying the hamiltonian to all bands), isn't it? In my (admittedly limited) experience it converges in about the same number of total hamiltonian applications as SCF, for systems with a bandgap. Also I couldn't find in the castep docs what scheme it uses for direct minimization with fractional occupancies, I assume it's the one from Marzari/Vanderbilt/Payne '97? – Antoine Levitt May 18 '20 at 12:59
• @AntoineLevitt Yes, we use the ensemble DFT method of Marzari, Vanderbilt & Payne, and the expense there comes from the additional Fourier transforms (for small systems and/or many processes) and the additional matrix diagonalisations (for large systems). For insulators the method is much quicker, although density mixing is still usually faster (assuming it converges!). It is complicated, because the direct minimisation methods need more Fourier transforms per wavefunction update, but typically require fewer wavefunction updates. – Phil Hasnip May 19 '20 at 13:27
• OK, thanks! So the problem is the bad convergence of the M/V/P scheme for metals? (it's something I plan to work on soon) – Antoine Levitt May 20 '20 at 14:41

Well, arguably the core numerical routine in density functional theory is the evaluation of the exchange-correlation functional itself, and the evaluation by quadrature of the integrals where the exchange-correlation functional sits.