After about a year of work on the density-functional toolkit (DFTK) we finally started using the code for some new science. Given the mathematics background of the CERMICS the first two DFTK-related articles were not so much about large-scale applications, but rather deal some fundamental questions of Kohn-Sham density-functional theory (DFT). The first submission provides a detailed analysis of self-consistent iterations and direct-minimisation approaches for solving equations such as DFT. My main focus in the past weeks, however, was rather the second project, an invited submission for the Faraday discussions New horizons in density functional theory, which are to take place in September this year.

In this article we deal with numerical error in DFT calculations. More precisely we want to address the question: What is the remaining error in the result obtained by a particular DFT simulation, which has already been performed?

This is naturally a rather broad statement, which cannot realistically be addressed in the scope of a single article. As we detail we only address the question of bounding the error of a particular quantity of interest, namely band energies near the Fermi level. Other quantities, such as the forces or the response to an external perturbation, might benefit from our ideas, but will need extensions on top. Also one needs to keep in mind that there are plenty of sources for error in a DFT calculation, including:

1. The model error due to replacing the (almost) exact model of the many-body Schrödinger equation by a reduced, but more feasible model like DFT.
2. The discretisation error due to employing only finitely many basis functions instead of solving analytically in an infinite-dimensional Hilbert space.
3. The algorithmic error resulting from using non-zero convergence tolerances in the eigensolvers as well as the SCF procedure.
4. The arithmetic error obtained by doing computation only in finite-precision floating-point arithmetic.

Of course in practice people often have a good ballpark idea of the error of DFT as a method or the errors obtained from a particular basis cutoff. Rightfully one might ask, why one should go through the effort of deriving provable bounds to the error in a DFT result? While there are surely many takes on this question, I only want to highlight two aspects in this summary:

• Educated guesses for convergence parameters taken from experience can fail. Typically they fail exactly when interesting things happen in a system and thus the usual heuristic breaks down. In other words converging a simulation to investigate what's going on becomes difficult when it's most needed. Detailed error analysis splitting up the error during an SCF iteration into contributions like the errors 1 to 4 (or even finer) can help to hint at the parameters worth tweaking or can provide insight into which error term behaves unusual in an SCF.
• Thinking such aspects one step further, a good bound to the individual terms even allows to equilibrate sources of error during a running calculation. The outlook of this idea would be a fully black-box scheme for DFT calculations where typical convergence parameters are automatically while the calculation progresses in order to yield the cheapest path to a desired target accuracy.

With respect to the errors 1 to 4 mentioned above one would of course like to be able to provide an upper bound to each of them. Unfortunately especially obtaining an estimate for the model error is a rather difficult task. In our work we have therefore concentrated on the errors 2 to 4 and moreover we only focused on Kohn-Sham models without self-consistency, i.e. where none of the terms in the Hamiltonian depend on the density / orbitals. It goes without saying that especially this last restriction needs to be lifted to make our results useful in practice. This angle we have left for future work.

Already at the stage of the current reduced model, however, there are a few aspects to consider when finding an upper bound:

• We wanted our bound to be fully-guaranteed, which means that we wanted to design a bound where we are able to prove that the exact answer must lie inside the bounds we give. This means when we provide error bars for band energies it is mathematically guaranteed to find the exact answer of the full-dimensional (complete-basis set limit) calculation at infinite accuracy and at zero convergence tolerances inside our bound.
• To be useful our bound should be (cheaply) computable, because otherwise plainly checking the calculation at vastly increased precision might end up being the better option. Notice that our bound does require to use a finer discretisation for some terms, but this is only a one-shot a posteriori step and not an iterative one.
• Ideally the bound should be sharp meaning that the upper bound to the error we report should not be too far off the true error. Even better would be an optimal bound, where we are as close to the true error as possible (given proposed structure in the error estimate). Such considerations are very important to not end up with completely correct but also useless statements like: "The error in the band energy is smaller than 10 million Hartree".

Finding a balance between these three aspects is not always easy and in our work we often take the pragmatic route to obtain a simpler, albeit less sharp error. Still, our bound is fully computable and allowed us, for the first time, to report band structure diagrams of silicon annotated with fully-guaranteed error bars of combined discretisation, algorithm and arithmetic errors. Details of our methodologies are given in the paper. Its full abstract reads

We address the problem of bounding rigorously the errors in the numerical solution of the Kohn-Sham equations due to (i) the finiteness of the basis set, (ii) the convergence thresholds in iterative procedures, (iii) the propagation of rounding errors in floating-point arithmetic. In this contribution, we compute fully-guaranteed bounds on the solution of the non-self-consistent equations in the pseudopotential approximation in a plane-wave basis set. We demonstrate our methodology by providing band structure diagrams of silicon annotated with error bars indicating the combined error.

Article
Michael F. Herbst, Antoine Levitt and Eric Cancès. A posteriori error estimation for the non-self-consistent Kohn-Sham equations. Faraday Discussions, 224, 227 (2020). DOI 10.1039/D0FD00048E [code]