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:
- 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.
- The discretisation error due to employing only finitely
many basis functions instead of solving analytically
in an infinite-dimensional Hilbert space.
- The algorithmic error resulting from using non-zero convergence
tolerances in the eigensolvers as well as the SCF procedure.
- 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.
