Chapter 1 FLEIM:
A stable, accurate and robust extrapolation method at
infinity for computing the ground state of electronic Hamiltonians
Étienne Polack
Université
Bourgogne FrancheComté and CNRS,
Laboratoire de Mathématiques de Besançon, F25030 Besançon,
France
Yvon Maday
Sorbonne Université, CNRS, Université de Paris,
Laboratoire JacquesLouis Lions (LJLL), F75005 Paris,
and Institut Universitaire de France
Andreas Savin
CNRS and Sorbonne Université,
Laboratoire de Chimie Théorique (LCT),
F75005 Paris, France
Abstract
The Kohn–Sham method uses a single model system, and corrects it by a density functional the exact user friendly expression of which is not known and is replaced by an approximated, usable, model. We propose to use instead more than one model system, and use a greedy extrapolation method to correct the results of the model systems. Evidently, there is a higher price to pay for it. However, there are also gains: within the same paradigm, e.g., excited states and physical properties can be obtained.
1 Introduction
1.1 Motivation
Density functional theory (DFT) has a weak point: its approximations (DFAs). First, the Hohenberg–Kohn theorem tells us that there is a density functional for electronic systems, $F[\rho ]$, that is universal (that is, independent of the potential of the nuclei), but does not give us a hint on how systematic approximations can be constructed. In practice, models are produced to be fast in computations, typically by transferring properties from other systems, like the uniform electron gas. Second, the most successful approximations are using the Kohn–Sham method (introducing a fermionic wave function) that decomposes $F[\rho ]$ into the kinetic energy, the Hartree energy and an exchangecorrelation energy contribution although the question of how and what part of $F[\rho ]$ should be approximated is, in principle, open.
In the present contribution we totally change the paradigm in the following way still led by the issue of universality. Let us start with a physical consideration. When electrons are close, the Coulomb repulsion is so strong that some of its features dominate over the effect of the external potential. This is also reflected mathematically in the shortrange behavior of the wave function, as present in the Kato cusp condition [14, 7, 26, 6], and in higherorder terms [23, 17]. We further note that approximating numerically the shortrange part of the wave function needs special care, due to the singularity of the Coulomb interaction when the electrons are close.
The considerations above and the independence of the interaction between electrons from that between them and the external potential provides a basis for constructing approximations. Thus, we propose to solve accurately a Schrödinger equation with a Hamiltonian that is modified to eliminate the shortrange part of the interaction between the electrons which is one of the difficult parts in the numerical simulations. The way to do it is not unique, and we try to turn this to our advantage: we use several models, and from them we try to extrapolate to the physical system [25]. In other words, we follow an “adiabatic connection” (see [12]), without ever constructing a density functional. This new paradigm thus explores the possibility to replace the use of DFAs by mathematically controlled approximations: we make density functional theory “without density functionals.”
Our approach has introduced an additional difficulty nonexistent in the Kohn–Sham method: the longrange part of the interaction has to be treated accurately, and not only its electrostatic component. One may ask whether this additional effort is justified, and whether one gains anything with respect to a calculation where the physical (Coulomb) interaction is used. For a single calculation, the gain is due to the lack of singularity in the interaction expressed by a weak interaction potential allowing for simplified treatments, such as perturbation theory. However, as the extrapolation to the physical system needs more than one point, it is essential that the number of points stays very small, and the interaction weak.
1.2 Objective and structure of the paper
We first choose, in Sec.2.1, a family of model (parameter dependent) Hamiltonians that are more flexible than using only the Kohn–Sham (noninteracting) Hamiltonian.^{∗}^{∗}∗Note however that this is at the prize of working in ${\mathbb{R}}^{3N}$ instead of ${\mathbb{R}}^{3}$, and thus requiring accurate manybody, e.g., configuration interaction calculations. This is followed by a description of how universality is introduced, namely by analyzing how a nonsingular interaction approaches the Coulomb one, and not by transfer from other systems, as usually done in DFAs. The physical system of interest is one among the parameter dependent models corresponding to some precise value of the parameter; in Sec. 2.2 its solution is extrapolated from the solutions to the models for other values of the parameter, expected that these solutions are more simple to be approximated. This extrapolation is efficiently handled in the general framework of the model reduction methods and more precisely referring to a variation of the Empirical Interpolation Method [2].
We believe that such an approach can not only discuss what DFAs are really doing, but can evolve to being used in applications. Some argument supporting this statement is given. However, in this paper numerical examples (gathered in Sec. 3) are only presented for twoelectron systems that are numerically (and sometimes even analytically) easily accessible: the harmonium, the hydrogen anion, H^{–}, and the hydrogen molecule, H_{2} in the ground state, at the equilibrium distance.
As we do not use the Hohenberg–Kohn theorem, the technique can be applied without modification also to excited states. We provide in Sec. 3.5, as an example, the first excited state of the same symmetry as the ground state.
2 Approach
2.1 The model Schrödinger equation
We study a family of Schrödinger equations,
$$H(\mu )\mathrm{\Psi}(\mu )=E(\mu )\mathrm{\Psi}(\mu ),$$  (2.1) 
where $\mu $ is some nonnegative parameter. More precisely, in this paper, we use
$$H(\mu )=T+V+W(\mu ),$$  (2.2) 
where $T$ is the operator for the kinetic energy, $V$ is the external potential (in particular that of the interaction between nuclei and electrons) and $W(\mu )$ represents the interaction between electrons. Although not required by the general theory, in this paper we introduce the dependence on $\mu $ only by modifying the interaction between electrons,
$$  (2.3) 
choosing
$$w({r}_{ij},\mu )=\frac{\mathrm{erf}(\mu {r}_{ij})}{{r}_{ij}}$$  (2.4) 
where ${r}_{ij}={\bm{r}}_{\bm{i}}{\bm{r}}_{\bm{j}}$ is the distance between electron $i$ (at position ${\bm{r}}_{\bm{i}}$) and electron $j$ (at position ${\bm{r}}_{\bm{j}}$). Finally, the external potential $V$ is written like
$$V=\sum _{i=1}^{N}v({\bm{r}}_{\bm{i}}).$$  (2.5) 
where $v$ is the local one particle operator. Note that the $N$particle operators are denoted by upper case letters, while the oneparticle operators are denoted by lower case letters.
Note also that for $\mu =0$ we have a trivial noninteracting system, while for $\mu =\mathrm{\infty}$ we recover the Coulomb system. The operator $w$ is longranged: as $\mu $ increases, the Coulomb interaction $1/{r}_{12}$ starts being recovered from large distances. The first reason for this choice is that, as mentioned above, we expect a universal character for short range (this is related to the difficulty of common DFAs to correctly describe longrange contributions, cf. Appendix A). The second reason is that the solution of Eq. (2.1) is converging more rapidly with (conventional) basis set size when the interaction has no singularity at ${r}_{12}=0$.
In principle, introducing a dependence of the oneparticle operators ($T$ and $V$) on $\mu $ makes the formulas a bit more clumsy, but does not introduce important difficulties in its application. Using such a dependence might improve the results, but it is not discussed in this contribution. In the following, in order to simplify notation, we drop the argument $\mu $, when $\mu =\mathrm{\infty}$, e.g., $E=E(\mu =\mathrm{\infty})$.
2.2 The correction to the model
2.2.1 Using a basis set
Of course, solving the Schrödinger equation for the model, Eq. (2.1) with finite $\mu $s, does not provide the desired solution, i.e., the one that is obtained for $\mu =\mathrm{\infty}$. We thus need to estimate the difference in eigenvalues:
$$\overline{E}(\mu )=EE(\mu ).$$  (2.6) 
Since $\overline{E}(\mu )$ tends to zero at infinity, the idea is to first expand this difference $\overline{E}(\mu )$ in a basis (of functions that tend to zero at infinity), retaining $M$ terms,
$$\overline{E}(\mu )\approx {\overline{E}}_{M}(\mu )=\sum _{j=1}^{M}{c}_{j}{\chi}_{j}(\mu ),$$  (2.7) 
leading to
$$E(\mu )\approx E\sum _{j=1}^{M}{c}_{j}{\chi}_{j}(\mu ),$$ 
or, more precisely since $E$ is not known, we replace it by an approximation denoted as ${E}_{M}$,
$$E(\mu )\approx {E}_{M}\sum _{j=1}^{M}{c}_{j}{\chi}_{j}(\mu ).$$  (2.8) 
The idea then proceeds by determining the unknown ${E}_{M}$ values and the coefficients ${c}_{i}$ from $M+1$ values of $E({\mu}_{m})$, for $m=0,\mathrm{\dots},M$ for an appropriate choice of the parameter values ${\mu}_{m}$. Finally, taking into account that the functions ${\chi}_{j}$ tend to zero at infinity, the proposed approximation for $E$ is ${E}_{M}$. Of course, this extrapolation approach often fails if care is not enough taken in the choice of the functions ${\chi}_{j},1\le j\le M$, and the values ${\mu}_{m}$, for $m=0,\mathrm{\dots},M$.
First, one has to decide about their form. Second, one has to find a way to keep $M$ as small as possible to reduce computational cost while preserving a good accuracy.
2.2.2 Approaching the Coulomb interaction
As recalled above, we derive from the leading term of the Coulomb interaction between the electrons that, to leading order, the solutions of differential equations are determined at short range by the singularities. The interaction $w$ in Eq. (2.1) has no singularity at ${r}_{12}=0$, for any finite $\mu $. However, as the parameter $\mu $ increases, $w(\cdot ,\mu )$ approaches the singular Coulomb potential.
In order to see how this limit is approached, let us perturb the exact solution. To first order, the perturbation correction to the energy is given by
$$\overline{E}(\mu )=\u27e8\mathrm{\Psi}\left\left(WW(\mu )\right)\right\mathrm{\Psi}\u27e9,\text{for}\mu \to \mathrm{\infty}.$$  (2.9) 
By changing the integration variables ${\bm{r}}_{\bm{i}}$ to $\mu {\bm{r}}_{\bm{i}}$ we see that
$$\overline{E}(\mu )\propto {\mu}^{2}\mathit{\hspace{1em}}\text{as}\mu \to \mathrm{\infty},$$  (2.10) 
providing a leading behavior that we want the basis functions ${\chi}_{i}$ to reflect. It is possible to continue this analysis for higher order terms. In fact, the next term (in ${\mu}^{3}$) has a coefficient proportional to that of ${\mu}^{2}$, the proportionality coefficient being determined by the nature of the Coulomb singularity [9].
2.2.3 Choice of the basis functions
In the main part of this contribution we use a simple ansatz,
$${\stackrel{~}{\chi}}_{j}(\mu )=1j\mu {(1+{j}^{2}{\mu}^{2})}^{1/2},j=1,\mathrm{\dots},M,$$  (2.11) 
that respects indeed the condition of Eq. (2.10). The motivation for this specific choice, that is arbitrary to a certain degree, as well as some results obtained with other choices of basis functions, is given in Appendix B.
The first functions of this basis set are presented in Fig. 1, together with an example of a function it has to approximate. It illustrates that the function we want to describe is between basis function ${\stackrel{~}{\chi}}_{2}$ — for small $\mu $ — and basis function ${\stackrel{~}{\chi}}_{3}$ —for large $\mu $. However, a simple linear combination between these (only) two surrounding basis functions from the family in Eq. (2.11) does not improve much the accuracy, but of course, more (and more appropriate) functions in the family can (and should) be called.
2.2.4 Reducing the basis set
Using a large set of ${\chi}_{j}$ (a large $M$) can rapidly become computationally prohibitive (because it requires a large number of evaluations of $E({\mu}_{m})$, for $m=0,\mathrm{\dots},M$) and numerically unstable (because it is classically much more difficult to stabilize extrapolation than interpolation). In order to reduce their number and increase the stability of the extrapolation, we use a greedy (iterative) method, as in the Empirical Interpolation Method (EIM) leading to proper choices of ${\mu}_{m}$, for $m=0,\mathrm{\dots},M$ known as “magic points.”
In the $K$th iteration of EIM, one starts from a set of $K1$ basis functions (for us, ${\stackrel{~}{\chi}}_{j}$) and $K1$ points (for us, ${\stackrel{~}{\mu}}_{j}$ belonging to some (discretized) interval, say, close to zero, to benefit at most of the regularization of the $\mathrm{erf}$ function). One then chooses the $K$th function ${\stackrel{~}{\chi}}_{K1}$ (among the remaining $MK$ basis functions) as being the one that is most poorly approximated by the current interpolation (based on the $K1$ basis functions and the $K1$ points) in a sense dedicated to the final goal we want to achieve (that can be uniform error, error on some part of the domain, or even at some value) and the $K$th point ${\stackrel{~}{\mu}}_{K1}$ that, in the admissible set, brings the more information. In this contribution, as we are only interested at extrapolating the value of $\mu $ at infinity, so we chose the error as the absolute value of difference between the $K$th basis function and its interpolant at infinity as the final goal we want to achieve.
Note that the procedure selecting the next point and function does not make any use of the function to be approximated (here $\overline{E}$). It is thus a cheap step compared with the calculation of $E({\stackrel{~}{\mu}}_{m})$ on the system of interest.
To improve the results for the extrapolation, we have modified the EIM algorithm into what we call the Forward Looking EIM (FLEIM). While EIM tries to get the maximal improvement through a sequential choice of, first the new basis function, then the new point of interpolation, FLEIM tries to get the best pair for improvement in the selected goal. The method is explained in more detail in Appendix C. In what follows, we present the results of FLEIM as they are better and more stable than those of EIM, as is illustrated in App D.1.
2.3 Computing other physical properties
FLEIM can be used to approximate other physical properties, i.e., correct expectation values of operators $A\ne H$ obtained with the model wave functions, $\mathrm{\Psi}(\mu )$,
$$A(\mu )=\u27e8\mathrm{\Psi}(\mu )A\mathrm{\Psi}(\mu )\u27e9.$$  (2.12) 
This can be seen immediately by noting that the derivation in Sec. 2.2.1 is not specific for correcting $E(\mu )$, but can also be applied to $A(\mu )$.
For the choice of the basis functions, we point out that properties are obtained by perturbing the Hamiltonian with the appropriate operator, say, $A$,
$$H\to H(\lambda )=H+\lambda A.$$  (2.13) 
The expectation value of $A$ can be obtained as the derivative of $E(\lambda )$ w.r.t. $\lambda $, at $\lambda =0$. Of course, this procedure can be applied to model Hamiltonians, yielding $E(\lambda ,\mu )$ and
$$\u27e8\mathrm{\Psi}(\mu )A\mathrm{\Psi}(\mu )\u27e9={\partial}_{\lambda}E(\lambda ,\mu ){}_{\lambda =0}$$  (2.14) 
Thus, in this contribution, we use the same type of basis functions for $A(\mu )$ as for $E(\mu )$; see the results in Sec. 3.4. Note that computing $\u27e8\mathrm{\Psi}(\mu )A\mathrm{\Psi}(\mu )\u27e9$ is not possible in DFT, without having a propertyspecific density functional [3].
3 Numerical results
3.1 Guidelines
The quality of the corrections using Eqs. (2.2.1) and (2.11) is explored numerically. Technical details on the calculations are given in Appendix D.
The plots show the errors done by the approximations in the estimate of the energy: we choose a model, $\mu \mapsto E(\mu )$, and let the empirical interpolation method choose which easier models (with weaker interactions) to extrapolate and get an approximation for $E=E(\mu =\mathrm{\infty})$. The plots show the error in the estimate of $E$ made when considering approximations that use information only for ${\stackrel{~}{\mu}}_{m}\le \mu $. From the plots, we read off how small $\mu $ can be and still have “reasonable” accuracy. In thermochemistry, $\mathrm{kcal}\text{}{\mathrm{mol}}^{1}$ is a commonly considered unit, and is often considered as “chemical accuracy.” For electronic excitations, one often uses $\mathrm{eV}$ units, and one often indicates it with one decimal place. “Chemical accuracy” is marked in the plots by horizontal dotted lines. The plots show the errors in the range of $\pm 0.1\text{}\mathrm{eV}\approx 2.3\text{}\mathrm{kcal}\text{}{\mathrm{mol}}^{1}$.
We consider approximations using up to four points (thus chosen in $[0,\mu ]$). The first point ${\mu}_{0}$ is always the value chosen ${\mu}_{0}=\mu $ shown on the $x$axis of the plots, and the basis function associated to it is ${\stackrel{~}{\chi}}_{0}$, the constant function; note that using only this pair $({\chi}_{0},{\mu}_{0})$ corresponds to choosing $E\simeq E({\mu}_{0})$= the value provided by the model, i.e., no correction is applied. When the number of points is increased, further values of $E({\stackrel{~}{\mu}}_{m})$, chosen by the algorithm, are used with $$.
The (maximal) parameter $\mu $ is considered between $0$ and $3\text{}{\mathrm{bohr}}^{1}$. The model without correction (blue curve) reaches chemical accuracy for $\mu \approx 3\text{}{\mathrm{bohr}}^{1}$ for H^{–} and harmonium, but only at $\mu \approx 5\text{}{\mathrm{bohr}}^{1}$ for H_{2} in its ground state.
3.2 General behavior of errors
The plots in Fig. 2 for harmonium, H_{2}, and H^{–} have similar features and are discussed together. As the number of points used increases, the smallest value of $\mu $ for which the good accuracy is reached decreases. Note that FLEIM produces very small errors for values of $\mu $ larger than 2. However, with the chosen basis set, the algorithm presented in this contribution has difficulties correcting the errors for $\mu $ smaller than 1.
3.3 Possibility of error estimates
Some tests can be done to estimate the quality of the approximation. For example, we can compare how the approximations change when increasing the number of basis functions, $K$, in our approximation and consider ${E}_{K}{E}_{K1}$ as an asymptotically valid error estimate for ${E}_{K1}$. One can notice in the above figures that, when the difference between, say, the 2 and the 3point approximation error is larger than “chemical accuracy,” so is the error in the 2point approximation.
3.4 Expectation values with FLEIM: $\u27e8{r}_{12}\u27e9$ and $\u27e8{r}_{12}^{2}\u27e9$ for harmonium
We look at the average distance between the electrons in harmonium. Figure 3(top) shows the error made by using $\mathrm{\Psi}(\mu )$ instead of $\mathrm{\Psi}(\mu =\mathrm{\infty})$ in computing the expectation value of ${r}_{12}$, as well as the correction that can be achieved with FLEIM, using the same basis set as above (2.11). The inset in Figure 3(top) concentrates on the errors made in the region that could be considered chemically relevant ($1\text{}\mathrm{pm}$ $\approx $ $0.02\text{}\mathrm{bohr}$). We note the similarity with the behavior of in correcting $E(\mu )$.
Let us now examine the average square distance between the electrons, $\u27e8{r}_{12}^{2}\u27e9$, in harmonium. While for computing the energy we explored correcting the missing shortrange part of the interaction, we now ask whether it is possible to correct the error of using the model wave function, $\mathrm{\Psi}(\mu )$ for the expectation value of an operator that is important at long range.
For $\omega =1/2$, we know the exact values of the expectation value of ${r}_{12}^{2}$ at $\mu =0$ and $\mu =\mathrm{\infty}$; they are 6 and $(42\sqrt{\pi}+64)/(5\sqrt{\pi}+8)\approx 8.21$, respectively (see, e.g., Ref. [15]). Note the large effect of the model wave function, $\mathrm{\Psi}(\mu )$, in computing $\u27e8{r}_{12}^{2}\u27e9$. Figure 3(bottom) shows the error made by using $\mathrm{\Psi}(\mu )$ instead of $\mathrm{\Psi}(\mu =\mathrm{\infty})$ in computing the expectation value of ${r}_{12}^{2}$, as well as the effect of the correction that can be achieved with FLEIM, using the same basis set (2.11) as above. We note again the similarity with the behavior of in correcting $E(\mu )$ or the expectation value of the distance between electrons.
The expectation value $\u27e8{r}_{12}^{2}\u27e9$ also illustrates another aspect: the effect of a change of the external potential on the energy. At first sight this may seem surprising, as the external potential is a oneparticle operator, while ${r}_{12}^{2}$ is a twoparticle operator. However, changing the oneparticle operator also modifies the wave function and this affects the value of $\u27e8{r}_{12}^{2}\u27e9$. In the case of harmonium, this can be shown analytically. Changing ${\bm{r}}_{\mathit{1}}$ and ${\bm{r}}_{\mathit{2}}$ to centerofmass, $\bm{R}$, and interparticle distance, ${\bm{r}}_{\mathit{12}}$, cf. Appendix E, allows us to see explicitly that a modification of ${\omega}^{2}$, the parameter that specifies the external potential, affects the Schödinger equation in ${\bm{r}}_{\mathit{12}}$. It introduces a term proportional to ${\omega}^{2}{r}_{12}^{2}$. The first order change in the energy when we change the external potential (${\omega}^{2}$) is thus proportional to $\u27e8{r}_{12}^{2}\u27e9$. Our results in Fig. 3(bottom) show that our conclusions on model corrections are not modified by small changes in the external potential. Note that the centerofmass Schrödinger equation also depends on ${\omega}^{2}$, but it is independent of $\mu $ and thus does not affect our discussion on model correction.
3.5 Comparison with DFAs
Instead of using extrapolation with FLEIM, one can use DFAs. While up to now the external potential did not change with $\mu $, in DFA calculations a one particle potential that depends on $\mu $ is added in order to correct the density.
We consider here two DFAs, the local density approximation, LDA [24, 22], and one that reproduces that of Perdew, Burke and Ernzerhof (PBE) [10, 8] at $\mu =0$. Both approximations are modified to be $\mu $dependent. In particular, they vanish at $\mu =\mathrm{\infty}$.
As shown in Fig. 4(top) for harmonium, DFAs are clearly much better at small $\mu $. However, they are not good enough. The figure suggests the range of $\mu $ for which DFAs are within chemical accuracy is similar to that obtained with the 3point FLEIM. This is confirmed when comparing the results with DFAs and for the H_{2}; see Fig. 4(middle). Note that with FLEIM the errors at large $\mu $ are smaller.
Note also that the curves obtained with extrapolation are significantly flatter at large $\mu $ than those obtained with DFAs. This should not be surprising: DFAs transfer the large $\mu $ behavior, while extrapolation extracts it from information available for the system under study.
Furthermore, using groundstate DFAs for excited states does not only pose a problem of principle (questions its validity, as the Hohenberg–Kohn theorem is proven for the ground state), but can also show a deterioration of quality. However, there is no question of principle from the perspective of this contribution (of using a model and correcting it by extrapolation). Also, the error in the excited state seems comparable to that in the ground state, as seen in the example of the H_{2} molecule, in the first excited state of the same symmetry as the ground state; see Fig. 4(bottom).
4 Conclusion and perspectives
In this contribution we have illustrated with a few models how to simplify the Hamiltonians by smoothly getting rid of the singularities in the system and thus have more numerically tractable problems. This simplification is obtained by introducing a parameter that, when it is equal to infinity, it corresponds to the original, plain Hamiltonian. After numerically solving a few simplified problems, the solution of interest is obtained by extrapolation.
We present a new (in the field) method for extrapolating the quantities of interest from few finite values (hence easy to solve) of the parameter by a technique borrowed from reduced basis paradigm: the empirical interpolation method. In contrast to DFAs, no parameters are fitted, no transfer from different system are made: only extrapolation is used. Note that in contrast to DFAs, improvement can be envisaged by either adding further points or using more appropriate basis functions and error estimates are asymptotically accessible.
Appendices
A On density functional approximations
In DFT, the existence of a universal functional of the density, $F[\rho ]$, i.e., the Hohenberg–Kohn theorem [13], is rigorously proven [18, 19]. However, obtaining accurately the value of $F$ for a given density $\rho (\bm{r})$, while possible (see, e.g., Ref. [4]), is exceedingly timeconsuming. Computationally convenient DFAs exploit the knowledge of the density around a given point $\bm{r}$ in space. Typically,
$$F[\rho ]\approx {\int}_{{\mathbb{R}}^{3}}f(\rho (\bm{r}),\mathbf{\nabla}\rho (\bm{r}),\mathrm{\dots})d\bm{r}.$$  (A.1) 
The limitation of such an approach can be seen for a simple density functional, the Hartree term of the energy,^{∗}^{∗}∗The volume elements $\mathrm{d}{\bm{r}}_{\mathit{1}}\mathrm{d}{\bm{r}}_{\mathit{2}}$ are omitted when the context is unambiguous.
$${E}_{H}[\rho ]={\int}_{{\mathbb{R}}^{6}}\frac{\rho ({\bm{r}}_{\mathit{1}})\rho ({\bm{r}}_{\mathit{2}})}{{\bm{r}}_{\mathit{1}}{\bm{r}}_{\mathit{2}}}$$  (A.2) 
when $\rho (\bm{r})={\rho}_{\bm{A}}(\bm{r})+{\rho}_{\bm{B}}(\bm{r})$ and ${\rho}_{A}(\bm{r}){\rho}_{\bm{B}}(\bm{r})\approx \mathit{0}$ (i.e., ${\rho}_{A}$ and ${\rho}_{B}$ are spatially separated; their overlap decreases much faster than $1/{\bm{r}}_{\mathit{1}}{\bm{r}}_{\mathit{2}}$). Hohenberg and Kohn recognized the difficulty of approximating ${E}_{H}$ by expressions of the type given in Eq. (A.1) and suggested separating it from $F[\rho ]$. However, this does not fundamentally solve the problem, as one can immediately see in oneelectron systems, where ${E}_{H}$ has to be canceled by another term commonly expressed in DFAs by an ansatz of the forms given in Eq. (A.1). Note, that the problem would not exist for interactions that are not Coulomb (longranged), but shortranged. For instance, if the interaction is Dirac’s $\delta ({\bm{r}}_{\mathit{1}}{\bm{r}}_{\mathit{2}})$ function, ${E}_{H}$ becomes exactly of the form of Eq. (A.1). For other shortrange interactions one can approach such a form by using Taylor expansions. In recent years it has become popular to compensate for the limitation of the ansatz in Eq. (A.1) by adding “empirical” energy corrections to describe longrange effects.
Another problem is that the antisymmetry of the electronic wave function is hidden in $F[\rho ]$. As the most important effect of the antisymmetry is the Pauli repulsion, Kohn and Sham [16] proposed to consider the variational principle for a model in which particles do not interact. However, DFAs following the pattern of Eq. (A.1) are still in use. In fact, further separating terms from $F[\rho ]$ may even lead (for degenerate cases) to the question whether the limit of a noninteracting system is welldefined (see, e.g., Sec. 3.5 in [24]).
B Basis functions
In order to get an idea how the leading term of the correction behaves, we consider the missing part of the Hartree term,
$$\u27e8\mathrm{\Psi}(\mu )\left\left(WW(\mu )\right)\right\mathrm{\Psi}(\mu )\u27e9=\frac{1}{2}{\int}_{{\mathbb{R}}^{6}}\rho ({\bm{r}}_{\mathit{1}})\rho ({\bm{r}}_{\mathit{2}})\left(\frac{\mathit{1}}{{\bm{r}}_{\mathit{12}}}\bm{w}({\bm{r}}_{\mathit{12}},\mu )\right).$$  (B.1) 
Separating the Hartree part from $F[\rho ]$ was already proposed by Hohenberg and Kohn [13], and it is also the dominant part in $\overline{E}(\mu )$. Most molecular codes use Gaussian oneparticle basis functions, so $\rho $ is a linear combination of Gaussian functions. We consider a generic term,
$${\int}_{{\mathbb{R}}^{6}}{e}^{\alpha {r}_{1}^{2}}{e}^{\alpha {r}_{2}^{2}}\left(\frac{1}{{r}_{12}}w({r}_{12},\mu )\right).$$  (B.2) 
This integral is easily computed, e.g., by using Fourier transforms, and one obtains an expression that is proportional to the form of the basis function ${\stackrel{~}{\chi}}_{j}$ in Eq. (2.11), after arbitrarily relating $1/{\alpha}^{2}$ to the basis index $j$.
Let us now consider another basis set, constructed from the requirement that the functions decay as ${\mu}^{2}$ and are finite at the origin.
$${\stackrel{~}{\chi}}_{1,j}(\mu )={\left(1+{(j\mu )}^{2}\right)}^{1}$$  (B.3) 
where we choose again $j=1,2,\mathrm{\dots},10$. The results are only slightly worse than those obtained when using the basis set of Eq. (2.11), compare Fig. 2(top) and Fig. 5(top).
Let us consider
$${\stackrel{~}{\chi}}_{2,j}(\mu )={a}_{j}\mu {\left(1+{({a}_{j}\mu )}^{3}\right)}^{1}$$  (B.4) 
where we choose ${a}_{j}={2}^{j/2}$. The result for harmonium is shown in Fig. 5(bottom). We note a slight improvement of the results for the 2 and 3point approximations, as well as a change of sign of the error in the 4point approximation. This direction of investigation deserves to be pursued.
C The empirical interpolation method
The empirical interpolation method is a modelorder reduction method introduced in [2] as a way to efficiently find a reduced basis and approximate one particular function within a manifold of parameter dependent functions. The points at which to do the interpolation are referred as magic points [20].
For a family of basis functions ${\chi}_{i}$, $i\in I$ with discrete points ${\mu}_{j}$, $j\in J$ chosen on a regular grid close to zero (see Appendix D.2 for an analysis of the influence of the choice of the grid), we want to find a family of $K$ functions and interpolation points with which to interpolate a test function $f\in \mathrm{Span}\{{\chi}_{i},i\in I\}$.
C.1 Algorithm for EIM
We assume that we have chosen some cost function $\mathcal{C}$, e.g., a norm $\mathcal{C}[\cdot ]\equiv \parallel \cdot \parallel $, or, if we are, as in this contribution, only interested in correctly approximating the value for $\mu \equiv \mathrm{\infty}$, we choose the cost function as the absolute value at that extrapolation point $\mathcal{C}[\phi ]=\phi (\mathrm{\infty})={lim}_{\mu \to \mathrm{\infty}}\phi (\mu )$.
First, select one of the basis functions. We can choose to add the constant function ${\chi}_{0}$,
$${\stackrel{~}{\chi}}_{0}\u2254{\chi}_{0}.$$  (C.1) 
We then select the first interpolation point as the largest admissible $\mu $ available,
$${\stackrel{~}{\mu}}_{0}\u2254\underset{j\in J}{\mathrm{max}}{\mu}_{j}.$$  (C.2) 
We can then define the first normalized interpolation function as
$${q}_{0}\u2254\frac{{\stackrel{~}{\chi}}_{0}}{{\stackrel{~}{\chi}}_{0}({\stackrel{~}{\mu}}_{0})}.$$  (C.3) 
We can then create the first approximation with an interpolation scheme, for instance Lagrangian interpolation as
$${\mathcal{I}}_{0}[f]\u2254f({\stackrel{~}{\mu}}_{0}){q}_{0}.$$  (C.4) 
We now assume to have chosen $K1$ functions ${\stackrel{~}{\chi}}_{k}$, normalized functions ${q}_{k}$. Let us also assume that we have selected $K1$ interpolation points ${\stackrel{~}{\mu}}_{k}$, $k=0,\mathrm{\dots},K2$. We define the ${\mathcal{I}}_{K2}$ Lagrangian interpolation function as
$${\mathcal{I}}_{K2}[f]\u2254\sum _{k=0}^{K2}{\beta}_{k}{q}_{k},$$  (C.5) 
where the coefficients ${\beta}_{k}$ are determined by solving the system
$$\sum _{k=0}^{K2}{\beta}_{k}{q}_{k}({\stackrel{~}{\mu}}_{\mathrm{\ell}})=f({\stackrel{~}{\mu}}_{\mathrm{\ell}})\mathit{\hspace{1em}}\text{for}\mathrm{\ell}=0,\mathrm{\dots},K2.$$  (C.6) 
If we denote $\stackrel{~}{I}$ as the set of indices of remaining basis functions and $\stackrel{~}{J}$ as the set of indices of remaining interpolation points, we choose the next function as
$${\stackrel{~}{\chi}}_{K1}\u2254\underset{{\chi}_{i},i\in \stackrel{~}{I}}{\mathrm{arg}\mathrm{max}}\left\{\mathcal{C}\left[{\chi}_{i}{\mathcal{I}}_{K2}[{\chi}_{i}]\right]\right\},$$  (C.7) 
and the next interpolation point as
$${\stackrel{~}{\mu}}_{K1}\u2254\underset{{\mu}_{j},j\in \stackrel{~}{J}}{\mathrm{arg}\mathrm{max}}\left\{\left{\stackrel{~}{\chi}}_{K1}({\mu}_{j}){\mathcal{I}}_{K2}[{\stackrel{~}{\chi}}_{K1}]({\mu}_{j})\right\right\}.$$  (C.8) 
We can then define the $K$th normalised interpolation function as
$${q}_{K1}\u2254\frac{{\stackrel{~}{\chi}}_{K1}{\mathcal{I}}_{K2}[{\stackrel{~}{\chi}}_{K1}]}{{\stackrel{~}{\chi}}_{K1}({\stackrel{~}{\mu}}_{K1}){\mathcal{I}}_{K2}[{\stackrel{~}{\chi}}_{K1}]({\stackrel{~}{\mu}}_{K})}.$$  (C.9) 
This system is represented by a lower triangular matrix with ones on the diagonal, and hence has a unique solution. The algorithm ends when the desired target accuracy is reached.
C.2 The forward looking empirical interpolation method (FLEIM)
To better adapt the method for extrapolation, we propose a double loop alternative: Instead of selecting sequentially first for a new basis function and then a new interpolation point — Eqs. (C.7) and (C.8) — we select the best pair
$$({\stackrel{~}{\chi}}_{K1},{\stackrel{~}{\mu}}_{K1})\u2254\underset{{\chi}_{i},i\in \stackrel{~}{I}}{\mathrm{arg}\mathrm{max}}\underset{{\mu}_{j},j\in \stackrel{~}{J}}{\mathrm{arg}\mathrm{min}}\left\{\mathcal{C}\left[{\chi}_{i}({\mu}_{j}){\mathcal{I}}_{K2}[{\chi}_{i}]({\mu}_{j})\right]\right\}.$$  (C.10) 
D Numerical details of the calculations
D.1 Testing EIM and FLEIM with $E(\mu )=1+{\chi}_{j}(\mu )$
In this subsection, we compare on a simple function, the behaviour of EIM and FLEIM on a analytic test function $E(\mu )$ behaves like $1+{\chi}_{j}(\mu )$. The results for EIM are given in Fig. 6, and for FLEIM, in Fig. 7.
First, we note that FLEIM provides better approximation together with more stable results when $\mu $ varies. Second, we note that the “wall” for $j>1$ at small $\mu $ is also imporved. At large $\mu $, all ${\chi}_{j}$ have the same decay at large $\mu $: this makes both methods fit to work in this regime.
D.2 Discretization for FLEIM
The interval between $0$ and the value of $\mu $ under study was divided in $10$ equal intervals (producing $11$ points). FLEIM was used to select $K\phantom{\rule{veryverythickmathspace}{0ex}}(\le 4)$ points and basis functions on which $E(\mu )$ was calculated.
We now investigate the effect of changing the grid of $\mu $ and use, for comparison a fixed finer grid of values of ${\mu}_{1},{\mu}_{2},\mathrm{\dots}$ ranging from 0 to the value $\mu $ indicated in the plots, with a step of $0.01\text{}{\mathrm{bohr}}^{1}$, the FLEIM results are only slightly changed, as shown in Fig. 8.
D.3 Systems
Twoelectron systems studied are

1.
harmonium, having
$$v(\bm{r})=\frac{\mathit{1}}{\mathit{2}}{\omega}^{\mathit{2}}{\bm{r}}^{\mathit{2}}$$ (D.1) where for $\omega =1/2$ the exact energy is known ($E=2\text{}\mathrm{hartree}$);

2.
H^{–} anion,
$$v(\bm{r})=\frac{\mathit{1}}{\bm{r}};$$ (D.2) 
3.
H_{2} molecule,
$$v(\bm{r})=\frac{\mathit{1}}{\bm{r}{\bm{R}}_{\bm{A}}}\frac{\mathit{1}}{\bm{r}{\bm{R}}_{\bm{B}}}$$ (D.3) where the nuclei are in the equilibrium position, ${\bm{R}}_{\bm{A}}={\bm{R}}_{\bm{B}}$ with ${\bm{R}}_{\bm{A}}=\mathit{0.7}\text{}\mathrm{bohr}$.
D.4 Obtaining the model energy
In order to simplify the test of FLEIM, $E(\mu )$ was precalculated for a dense range of values $\mu $ and interpolated. The values for $E(\mu )$ were obtained with the program Molpro [1] for H^{–} and H_{2}. This program was also used for the density functional calculations.
For harmonium, it is possible to separate the variables in the Schrödinger equation. The centerofmass equation can be solved exactly. The equation in ${r}_{12}$ was solved by discretization on a grid of ${10}^{4}$ points between $0$ and $10\text{}\mathrm{bohr}$.
For H_{2} the V5Z basis set of [11] was used; for H^{–} the augV5Z basis set of [5]. (The error in the energy of H^{–} is too large if we do not augment the V5Z basis set with a diffuse basis function.) The augV5Z basis set was also used for the excited state of H_{2}, as it has an important contribution of ionic states (H^{+}$\mathrm{\cdots}$H^{–}). For H_{2} the equilibrium distance of $1.4\text{}\mathrm{bohr}$ was chosen for the ground state. For the excited state, the distance of $4.2\text{}\mathrm{bohr}$ was chosen. It is close to a minimum of the potential energy curve. Furthermore, this value can be compared with the accurate calculation of Ref. [21].
E Change of coordinates in harmonium
Harmonium is characterized by the external potential given by Eq. (D.1). For two particles the variables ${\bm{r}}_{\mathit{1}},{\bm{r}}_{\mathit{2}}$ can be changed to those corresponding to the center of mass and the distance between particles,
$\bm{R}$  $=$  $\frac{1}{2}}\left({\bm{r}}_{\mathit{1}}+{\bm{r}}_{\mathit{2}}\right),$  
${\bm{r}}_{\mathit{12}}$  $=$  $\left({\bm{r}}_{\mathit{1}}{\bm{r}}_{\mathit{2}}\right)$  (E.1) 
yielding for the potential energy
$$\frac{1}{2}{\omega}^{2}\left({\bm{r}}_{\mathit{1}}^{\mathit{1}}+{\bm{r}}_{\mathit{2}}^{\mathit{2}}\right)={\omega}^{2}\left({\bm{R}}^{\mathit{2}}+\frac{\mathit{1}}{\mathit{4}}{\bm{r}}_{\mathit{12}}^{\mathit{2}}\right)$$  (E.2) 
The transformation of variables can be done also for the kinetic energy, and makes the model Schrödinger equation separable into a part (in $\bm{R}$) that is independent of the model, and one (in ${\bm{r}}_{\mathit{12}}$) that through 2.1 depends on the model ($\mu $).
Note that $\u27e8\mathrm{\Psi}(\mu ){r}_{12}^{2}\mathrm{\Psi}(\mu )\u27e9$ is also measuring the error due to the change of density with $\mu $, for harmonium
$$H(\mu )=T+\frac{1}{2}{\omega}^{2}({r}_{1}^{2}+{r}_{2}^{2})+\frac{\mathrm{erf}(\mu {r}_{12})}{{r}_{12}}$$  (E.3) 
Indeed, by using the Hellmann–Feynman theorem,
$\frac{\partial}{\partial {\omega}^{2}}}\u27e8\mathrm{\Psi}(\mu )H(\mu )\mathrm{\Psi}(\mu )\u27e9$  $=$  $\u27e8\mathrm{\Psi}(\mu ){\displaystyle \frac{1}{2}}({r}_{1}^{2}+{r}_{2}^{2})\mathrm{\Psi}(\mu )\u27e9$  (E.4)  
$=$  $\u27e8\mathrm{\Psi}(\mu )({R}^{2}+{\displaystyle \frac{1}{4}}{r}_{12}^{2})\mathrm{\Psi}(\mu )\u27e9$ 
where $\bm{R}=({\bm{r}}_{\mathit{1}}+{\bm{r}}_{\mathit{2}})/\mathit{2}$. As it is possible to separate variables $\bm{R}$ and ${\bm{r}}_{\mathit{12}}$ in the Schrödinger equation,
$$\mathrm{\Psi}(\mu )=\mathrm{\Psi}(\bm{R},{\bm{r}}_{\mathit{12}},\mu )=\bm{\Phi}(\bm{R})\varphi ({\bm{r}}_{\mathit{12}},\mu )$$  (E.5) 
$\mathrm{\Phi}$ and $\varphi $ both normalized to one, and
$$\u27e8\mathrm{\Psi}(\mu )\frac{1}{2}({r}_{1}^{2}+{r}_{2}^{2})\mathrm{\Psi}(\mu )\u27e9={\int}_{{\mathbb{R}}^{3}}\rho (\bm{r},\mu ){\bm{r}}^{\mathit{2}}\mathrm{d}\bm{r}$$  (E.6) 
where $\rho (\bm{r},\mu )$ is the density of the model system, we have
$$\frac{1}{4}\u27e8\mathrm{\Psi}(\mu ){r}_{12}^{2}\mathrm{\Psi}(\mu )\u27e9={\int}_{{\mathbb{R}}^{3}}\rho (\bm{r},\mu ){\bm{r}}^{\mathit{2}}\mathrm{d}\bm{r}{\int}_{{\mathbb{R}}^{\mathit{3}}}\bm{\Phi}(\bm{R}){}^{\mathit{2}}{\bm{R}}^{\mathit{2}}\mathrm{d}\bm{R}$$  (E.7) 
Note that the change of $\u27e8{r}_{12}^{2}\u27e9$ with $\mu $ is only due to the change of the density with $\mu $. Thus, the difference between $\frac{1}{4}\u27e8\mathrm{\Psi}(\mu ){r}_{12}^{2}\mathrm{\Psi}(\mu )\u27e9$ and $\frac{1}{4}\u27e8\mathrm{\Psi}{r}_{12}^{2}\mathrm{\Psi}\u27e9$ also indicates how much the density is affected by the model.
Funding
Part of this work was supported by the French “Investissements d’Avenir” program, project ISITEBFC (contract ANR15IDEX0003) (ÉP). This work has also received funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation program (grant agreement No 810367–project EMC2) (YM & AS).
References
 [1] (2002) MOLPRO, a package of ab initio programs designed by H.J. Werner and P. J. Knowles, version 2002.2. Birmingham, UK. Cited by: §D.4.
 [2] (200411) An ‘empirical interpolation’ method: application to efficient reducedbasis discretization of partial differential equations. Comptes Rendus Mathematique 339 (9), pp. 667–672 (en). External Links: ISSN 1631073X, Document Cited by: §1.2, §C.
 [3] (1983) General operator groundstate expectation values in the hohenbergkohnsham densityfunctional formalism. Phys. Rev. B 27, pp. 5912. Cited by: §2.3.
 [4] (199901) Correlation energies for some two and fourelectron systems along the adiabatic connection in density functional theory. The Journal of Chemical Physics 110 (6), pp. 2828–2835. External Links: ISSN 00219606, Document Cited by: §A.
 [5] (1994) Gaussian basis sets for use in correlated molecular calculations. iv. calculation of static electrical response properties. The Journal of Chemical Physics 100 (4), pp. 2975–2988. External Links: Document, Link, https://doi.org/10.1063/1.466439 Cited by: §D.4.
 [6] (202011) Explicit Green operators for quantum mechanical Hamiltonians. II. Edgetype singularities of the helium atom. AsianEuropean Journal of Mathematics 13 (07), pp. 2050122 (en). External Links: ISSN 17935571, 17937183, Link, Document Cited by: §1.1.
 [7] (2009) Analytic Structure of ManyBody Coulombic Wave Functions. Commun. Math. Phys. 289, pp. 291. Cited by: §1.1.
 [8] (200806) Local and density fitting approximations within the shortrange/longrange hybrid scheme: application to large nonbonded complexes. Physical Chemistry Chemical Physics 10 (23), pp. 3353–3357. External Links: ISSN 14639084, Document Cited by: §3.5.
 [9] (200610) A shortrange gradientcorrected spin density functional in combination with longrange coupledcluster methods: Application to alkalimetal raregas dimers. Chemical Physics 329 (1), pp. 276–282 (en). External Links: ISSN 03010104, Document Cited by: §2.2.2.
 [10] (200511) A shortrange gradientcorrected density functional in longrange coupledcluster calculations for rare gas dimers. Physical Chemistry Chemical Physics 7 (23), pp. 3917–3923. External Links: ISSN 14639084, Document Cited by: §3.5.
 [11] (1989) Gaussian basis sets for use in correlated molecular calculations. i. the atoms boron through neon and hydrogen. The Journal of Chemical Physics 90 (2), pp. 1007–1023. External Links: Document, Link, https://doi.org/10.1063/1.456153 Cited by: §D.4.
 [12] (1974) The surface energy of a bounded electron gassolid. J. Phys. F 4, pp. 1170–1186. Cited by: §1.1.
 [13] (196411) Inhomogeneous Electron Gas. Physical Review 136 (3B), pp. B864–B871. External Links: Document Cited by: §A, §B.
 [14] (1957) On the eigenfunctions of manyparticle systems in quantum mechanics. Communications on Pure and Applied Mathematics 10 (2), pp. 151–177 (en). External Links: ISSN 10970312, Link, Document Cited by: §1.1.
 [15] (199612) The electron correlation cusp. Theoretica chimica acta 94 (6), pp. 345–381 (en). External Links: ISSN 14322234, Link, Document Cited by: §3.4.
 [16] (196511) SelfConsistent Equations Including Exchange and Correlation Effects. Physical Review 140 (4A), pp. A1133–A1138. External Links: Document Cited by: §A.
 [17] (2016) General coalescence conditions for the exact wave functions: higherorder relations for coulombic and noncoulombic systems. Adv. Quantum Chem. 73, pp. 59. Cited by: §1.1.
 [18] (197912) Universal variational functionals of electron densities, firstorder density matrices, and natural spinorbitals and solution of the vrepresentability problem. Proceedings of the National Academy of Sciences 76 (12), pp. 6062–6065 (en). External Links: ISSN 00278424, 10916490, Document Cited by: §A.
 [19] (1983) Density functionals for Coulomb systems. International Journal of Quantum Chemistry 24 (3), pp. 243–277 (en). External Links: ISSN 1097461X, Document Cited by: §A.
 [20] (2009) A general multipurpose interpolation procedure: the magic points. Communications on Pure & Applied Analysis 8 (1), pp. 383 (en). External Links: Document Cited by: §C.
 [21] (2018) Solving the Schrödinger equation of hydrogen molecule with the free complement–local Schrödinger equation method: potential energy curves of the ground and singly excited singlet and triplet states, $\mathrm{\Sigma}$, $\mathrm{\Pi}$, $\mathrm{\Delta}$, and $\mathrm{\Phi}$. J. Chem. Phys. 149, pp. 244116. Cited by: §D.4.
 [22] (200604) Localspindensity functional for multideterminant density functional theory. Physical Review B 73 (15), pp. 155111. External Links: Document Cited by: §3.5.
 [23] (1996) Behavior of electronic wave functions near cusps. J. Chem. Phys. 104, pp. 9908. Cited by: §1.1.
 [24] (199601) On degeneracy, neardegeneracy and density functional theory. In Theoretical and Computational Chemistry, J. M. Seminario (Ed.), Recent Developments and Applications of Modern Density Functional Theory, Vol. 4, pp. 327–357. External Links: Link, Document Cited by: §A, §3.5.
 [25] (2011) Correcting model energies by numerically integrating along an adiabatic connection and a link to density functional approximations. J. Chem. Phys. 134, pp. 214108. Cited by: §1.1.
 [26] (2010) Regularity and approximability of electronic wave functions. Springer. Cited by: §1.1.