Post-processing of the planewave approximation of Schrödinger equations. Part I: linear operators

Abstract

Abstract In this article we prove a priori error estimates for the perturbation-based post-processing of the plane-wave approximation of Schrödinger equations introduced and tested numerically in previous works (Cancès, Dusson, Maday, Stamm and Vohralík, (2014), A perturbation-method-based a posteriori estimator for the planewave discretization of nonlinear Schrödinger equations. C. R. Math., 352, 941--946; Cancès, Dusson, Maday, Stamm and Vohralík, (2016), A perturbation-method-based postprocessing for the planewave discretization of Kohn–Sham models. J. Comput. Phys., 307, 446--459.) We consider here a Schrödinger operator ${{\mathscr{H}} \,}= -\frac{1}{2}\varDelta +{\mathscr{V}}$ on $L^2(\varOmega )$, where $\varOmega $ is a cubic box with periodic boundary conditions and where ${\mathscr{V}}$ is a multiplicative operator by a regular-enough function ${\mathscr{V}}$. The quantities of interest are, on the one hand, the ground-state energy defined as the sum of the lowest $N$ eigenvalues of ${{\mathscr{H}} \,}$, and, on the other hand, the ground-state density matrix that is the spectral projector on the vector space spanned by the associated eigenvectors. Such a problem is central in first-principle molecular simulation, since it corresponds to the so-called linear subproblem in Kohn–Sham density functional theory. Interpreting the exact eigenpairs of ${{\mathscr{H}} \,}$ as perturbations of the numerical eigenpairs obtained by a variational approximation in a plane-wave (i.e., Fourier) basis we compute first-order corrections for the eigenfunctions, which are turned into corrections on the ground-state density matrix. This allows us to increase the accuracy by a factor proportional to the inverse of the kinetic energy cutoff ${E_{\textrm{c}}}^{-1}$ of both the ground-state energy and the ground-state density matrix in Hilbert–Schmidt norm at a low computational extra cost. Indeed, the computation of the corrections only requires the computation of the residual of the solution in a larger plane-wave basis and two fast Fourier transforms per eigenvalue.