## Chapter 1 FLEIM: A stable, accurate and robust extrapolation method at infinity for computing the ground state of electronic Hamiltonians

Étienne Polack
Université Bourgogne Franche-Comté and CNRS,
Laboratoire de Mathématiques de Besançon, F-25030 Besançon, France

Sorbonne Université, CNRS, Université de Paris,
Laboratoire Jacques-Louis Lions (LJLL), F-75005 Paris,
and Institut Universitaire de France

Andreas Savin
CNRS and Sorbonne Université,
Laboratoire de Chimie Théorique (LCT), F-75005 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 exchange-correlation 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 short-range behavior of the wave function, as present in the Kato cusp condition [14, 7, 26, 6], and in higher-order terms [23, 17]. We further note that approximating numerically the short-range 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 short-range 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 long-range 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 many-body, 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 two-electron systems that are numerically (and sometimes even analytically) easily accessible: the harmonium, the hydrogen anion, H, and the hydrogen molecule, H2 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.

Some conclusions and perspectives are presented in Sec. 4. Finally, in order to facilitate reading the manuscript, various details are given in Appendices AE that follow Sec. 4.

### 2 Approach

#### 2.1 The model Schrödinger equation

We study a family of Schrödinger equations,

 $H(\mu)\Psi(\mu)=E(\mu)\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,

 $W(\mu)=\sum_{i (2.3)

choosing

 $w(r_{ij},\mu)=\frac{\operatorname{erf}(\mu r_{ij})}{r_{ij}}$ (2.4)

where $r_{ij}=\lvert\vecfont{r}_{i}-\vecfont{r}_{j}\rvert$ is the distance between electron $i$ (at position $\vecfont{r}_{i}$) and electron $j$ (at position $\vecfont{r}_{j}$). Finally, the external potential $V$ is written like

 $V=\sum_{i=1}^{N}v(\vecfont{r}_{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 one-particle operators are denoted by lower case letters.

Note also that for $\mu=0$ we have a trivial noninteracting system, while for $\mu=\infty$ we recover the Coulomb system. The operator $w$ is long-ranged: 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 long-range 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 one-particle 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=\infty$, e.g., $E=E(\mu=\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=\infty$. We thus need to estimate the difference in eigenvalues:

 $\bar{E}(\mu)=E-E(\mu).$ (2.6)

Since $\bar{E}(\mu)$ tends to zero at infinity, the idea is to first expand this difference $\bar{E}(\mu)$ in a basis (of functions that tend to zero at infinity), retaining $M$ terms,

 $\bar{E}(\mu)\approx\bar{E}_{M}(\mu)=\sum_{j=1}^{M}c_{j}\chi_{j}(\mu),$ (2.7)

 $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,\dotsc,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\leq j\leq M$, and the values $\mu_{m}$, for $m=0,\dotsc,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

 $\bar{E}(\mu)=\left\langle\Psi\middle|\mathopen{}\bigl{(}W-W(\mu)\bigr{)}% \middle|\mathopen{}\Psi\right\rangle,\quad\text{for \mu\rightarrow\infty}.$ (2.9)

By changing the integration variables $\vecfont{r}_{i}$ to $\mu\vecfont{r}_{i}$ we see that

 $\bar{E}(\mu)\propto\mu^{-2}\quad\text{as \mu\rightarrow\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,

 $\tilde{\chi}_{j}(\mu)=1-j\;\mu(1+j^{2}\;\mu^{2})^{-1/2},\quad j=1,\dotsc,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 $\tilde{\chi}_{2}$ — for small $\mu$ — and basis function $\tilde{\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,\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,\dots,M$ known as “magic points.”

In the $K$th iteration of EIM, one starts from a set of $K-1$ basis functions (for us, $\tilde{\chi}_{j}$) and $K-1$ points (for us, $\tilde{\mu}_{j}$ belonging to some (discretized) interval, say, close to zero, to benefit at most of the regularization of the $\operatorname{erf}$ function). One then chooses the $K$th function $\tilde{\chi}_{K-1}$ (among the remaining $M-K$ basis functions) as being the one that is most poorly approximated by the current interpolation (based on the $K-1$ basis functions and the $K-1$ 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 $\tilde{\mu}_{K-1}$ 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 $\bar{E}$). It is thus a cheap step compared with the calculation of $E(\tilde{\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\neq H$ obtained with the model wave functions, $\Psi(\mu)$,

 $A(\mu)=\langle\Psi(\mu)\rvert A\lvert\Psi(\mu)\rangle.$ (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\rightarrow 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

 $\langle\Psi(\mu)\rvert A\lvert\Psi(\mu)\rangle=\partial_{\lambda}E(\lambda,\mu% )\bigr{|}_{\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 $\langle\Psi(\mu)\rvert A\lvert\Psi(\mu)\rangle$ is not possible in DFT, without having a property-specific 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=\infty)$. The plots show the error in the estimate of $E$ made when considering approximations that use information only for $\tilde{\mu}_{m}\leq\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 $\pm0.1\text{\,}\mathrm{eV}\approx2.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 $\tilde{\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(\tilde{\mu}_{m})$, chosen by the algorithm, are used with $\tilde{\mu}_{m}<\mu$.

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\approx3\text{\,}{\mathrm{bohr}}^{-1}$ for H and harmonium, but only at $\mu\approx5\text{\,}{\mathrm{bohr}}^{-1}$ for H2 in its ground state.

#### 3.2 General behavior of errors

The plots in Fig. 2 for harmonium, H2, 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 ${\lvert E_{K}-E_{K-1}\rvert}$ as an asymptotically valid error estimate for $E_{K-1}$. One can notice in the above figures that, when the difference between, say, the 2- and the 3-point approximation error is larger than “chemical accuracy,” so is the error in the 2-point approximation.

#### 3.4 Expectation values with FLEIM: $\langle r_{12}\rangle$ and $\langle r_{12}^{2}\rangle$ for harmonium

We look at the average distance between the electrons in harmonium. Figure 3(top) shows the error made by using $\Psi(\mu)$ instead of $\Psi(\mu=\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, $\langle r_{12}^{2}\rangle$, in harmonium. While for computing the energy we explored correcting the missing short-range part of the interaction, we now ask whether it is possible to correct the error of using the model wave function, $\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=\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, $\Psi(\mu)$, in computing $\langle r_{12}^{2}\rangle$. Figure 3(bottom) shows the error made by using $\Psi(\mu)$ instead of $\Psi(\mu=\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 $\langle r_{12}^{2}\rangle$ 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 one-particle operator, while $r_{12}^{2}$ is a two-particle operator. However, changing the one-particle operator also modifies the wave function and this affects the value of $\langle r_{12}^{2}\rangle$. In the case of harmonium, this can be shown analytically. Changing $\vecfont{r}_{1}$ and $\vecfont{r}_{2}$ to center-of-mass, $\vecfont{R}$, and inter-particle distance, $\vecfont{r}_{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 $\vecfont{r}_{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 $\langle r_{12}^{2}\rangle$. 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 center-of-mass 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=\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 3-point FLEIM. This is confirmed when comparing the results with DFAs and for the H2; 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 ground-state 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 H2 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.

### 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(\vecfont{r})$, while possible (see, e.g., Ref. [4]), is exceedingly time-consuming. Computationally convenient DFAs exploit the knowledge of the density around a given point $\vecfont{r}$ in space. Typically,

 $F[\rho]\approx\int_{\mathbb{R}^{3}}f\left(\rho(\vecfont{r}),\lvert\boldsymbol{% \nabla}\rho(\vecfont{r})\rvert,\dotsc\right)\,\mathrm{d}\vecfont{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}\vecfont{r}_{1}\,\mathrm{d}\vecfont{r}_{2}$ are omitted when the context is unambiguous.

 $E_{H}[\rho]=\int_{\mathbb{R}^{6}}\frac{\rho(\vecfont{r}_{1})\rho(\vecfont{r}_{% 2})}{\lvert\vecfont{r}_{1}-\vecfont{r}_{2}\rvert}$ (A.2)

when ${\rho(\vecfont{r})=\rho_{A}(\vecfont{r})+\rho_{B}(\vecfont{r})}$ and ${\rho_{A}(\vecfont{r})\rho_{B}(\vecfont{r})\approx 0}$ (i.e., $\rho_{A}$ and $\rho_{B}$ are spatially separated; their overlap decreases much faster than $1/|\vecfont{r}_{1}-\vecfont{r}_{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 one-electron 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 (long-ranged), but short-ranged. For instance, if the interaction is Dirac’s $\delta(\vecfont{r}_{1}-\vecfont{r}_{2})$ function, $E_{H}$ becomes exactly of the form of Eq. (A.1). For other short-range 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 long-range 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 well-defined (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,

 $\left\langle\Psi(\mu)\middle|\mathopen{}\bigl{(}W-W(\mu)\bigr{)}\middle|% \mathopen{}\Psi(\mu)\right\rangle=\frac{1}{2}\int_{\mathbb{R}^{6}}\rho(% \vecfont{r}_{1})\rho(\vecfont{r}_{2})\left(\frac{1}{r_{12}}-w(r_{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 $\bar{E}(\mu)$. Most molecular codes use Gaussian one-particle 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 $\tilde{\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.

 $\tilde{\chi}_{1,j}(\mu)=\left(1+(j\,\mu)^{2}\right)^{-1}$ (B.3)

where we choose again $j=1,2,\dotsc,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

 $\tilde{\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 3-point approximations, as well as a change of sign of the error in the 4-point approximation. This direction of investigation deserves to be pursued.

### C The empirical interpolation method

The empirical interpolation method is a model-order 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\lVert\cdot\rVert$, or, if we are, as in this contribution, only interested in correctly approximating the value for $\mu\equiv\infty$, we choose the cost function as the absolute value at that extrapolation point $\mathcal{C}[\varphi]=|\varphi(\infty)|=\lim_{\mu\rightarrow\infty}|\varphi(\mu)|$.

First, select one of the basis functions. We can choose to add the constant function $\chi_{0}$,

 $\tilde{\chi}_{0}\coloneqq\chi_{0}.$ (C.1)

We then select the first interpolation point as the largest admissible $\mu$ available,

 $\tilde{\mu}_{0}\coloneqq\max_{j\in J}\mu_{j}.$ (C.2)

We can then define the first normalized interpolation function as

 $q_{0}\coloneqq\frac{\tilde{\chi}_{0}}{\tilde{\chi}_{0}(\tilde{\mu}_{0})}.$ (C.3)

We can then create the first approximation with an interpolation scheme, for instance Lagrangian interpolation as

 $\mathcal{I}_{0}[f]\coloneqq f(\tilde{\mu}_{0})q_{0}.$ (C.4)

We now assume to have chosen $K-1$ functions $\tilde{\chi}_{k}$, normalized functions $q_{k}$. Let us also assume that we have selected $K-1$ interpolation points $\tilde{\mu}_{k}$, $k=0,\dotsc,K-2$. We define the $\mathcal{I}_{K-2}$ Lagrangian interpolation function as

 $\mathcal{I}_{K-2}[f]\coloneqq\sum_{k=0}^{K-2}\beta_{k}q_{k},$ (C.5)

where the coefficients $\beta_{k}$ are determined by solving the system

 $\sum_{k=0}^{K-2}\beta_{k}q_{k}(\tilde{\mu}_{\ell})=f(\tilde{\mu}_{\ell})\quad% \text{for\ }\ell=0,\dotsc,K-2.$ (C.6)

If we denote $\tilde{I}$ as the set of indices of remaining basis functions and $\tilde{J}$ as the set of indices of remaining interpolation points, we choose the next function as

 $\tilde{\chi}_{K-1}\coloneqq\operatorname*{arg\,max}_{\chi_{i},\,i\in\tilde{I}}% \bigl{\{}\mathcal{C}\bigl{[}\chi_{i}-\mathcal{I}_{K-2}[\chi_{i}]\bigr{]}\bigr{% \}},$ (C.7)

and the next interpolation point as

 $\tilde{\mu}_{K-1}\coloneqq\operatorname*{arg\,max}_{\mu_{j},\,j\in\tilde{J}}% \bigl{\{}\left\lvert\tilde{\chi}_{K-1}(\mu_{j})-\mathcal{I}_{K-2}[\tilde{\chi}% _{K-1}](\mu_{j})\right\rvert\bigr{\}}.$ (C.8)

We can then define the $K$th normalised interpolation function as

 $q_{K-1}\coloneqq\frac{\tilde{\chi}_{K-1}-\mathcal{I}_{K-2}[\tilde{\chi}_{K-1}]% }{\tilde{\chi}_{K-1}(\tilde{\mu}_{K-1})-\mathcal{I}_{K-2}[\tilde{\chi}_{K-1}](% \tilde{\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

 $(\tilde{\chi}_{K-1},\tilde{\mu}_{K-1})\coloneqq\operatorname*{arg\,max}_{\chi_% {i},\,i\in\tilde{I}}\operatorname*{arg\,min}_{\mu_{j},\,j\in\tilde{J}}\bigl{\{% }\mathcal{C}\bigl{[}\chi_{i}(\mu_{j})-\mathcal{I}_{K-2}[\chi_{i}](\mu_{j})% \bigr{]}\bigr{\}}.$ (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(\leq 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},\dotsc$ 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

Two-electron systems studied are

1. 1.

harmonium, having

 $v(\vecfont{r})=\frac{1}{2}\omega^{2}r^{2}$ (D.1)

where for $\omega=1/2$ the exact energy is known ($E=2\text{\,}\mathrm{hartree}$);

2. 2.

H anion,

 $v(\vecfont{r})=-\frac{1}{r};$ (D.2)
3. 3.

H2 molecule,

 $v(\vecfont{r})=-\frac{1}{\lvert\vecfont{r}-\vecfont{R}_{A}\rvert}-\frac{1}{% \lvert\vecfont{r}-\vecfont{R}_{B}\rvert}$ (D.3)

where the nuclei are in the equilibrium position, ${\vecfont{R}_{A}=-\vecfont{R}_{B}}$ with $\lvert\vecfont{R}_{A}\rvert=0.7\text{\,}\mathrm{bohr}$.

#### D.4 Obtaining the model energy

In order to simplify the test of FLEIM, $E(\mu)$ was pre-calculated for a dense range of values $\mu$ and interpolated. The values for $E(\mu)$ were obtained with the program Molpro [1] for H and H2. 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 center-of-mass 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 H2 the V5Z basis set of [11] was used; for H the aug-V5Z 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 aug-V5Z basis set was also used for the excited state of H2, as it has an important contribution of ionic states (H+$\cdots$H). For H2 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 $\vecfont{r}_{1},\vecfont{r}_{2}$ can be changed to those corresponding to the center of mass and the distance between particles,

 $\displaystyle\vecfont{R}$ $\displaystyle=$ $\displaystyle\frac{1}{2}\left(\vecfont{r}_{1}+\vecfont{r}_{2}\right),$ $\displaystyle\vecfont{r}_{12}$ $\displaystyle=$ $\displaystyle\left(\vecfont{r}_{1}-\vecfont{r}_{2}\right)$ (E.1)

yielding for the potential energy

 $\frac{1}{2}\omega^{2}\left(\vecfont{r}_{1}^{1}+\vecfont{r}_{2}^{2}\right)=% \omega^{2}\left(\vecfont{R}^{2}+\frac{1}{4}\vecfont{r}_{12}^{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 $\vecfont{R}$) that is independent of the model, and one (in $\vecfont{r}_{12}$) that through 2.1 depends on the model ($\mu$).

Note that $\langle\Psi(\mu)\rvert r_{12}^{2}\lvert\Psi(\mu)\rangle$ 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{\operatorname{erf}(% \mu r_{12})}{r_{12}}$ (E.3)

Indeed, by using the Hellmann–Feynman theorem,

 $\displaystyle\frac{\partial}{\partial\omega^{2}}\langle\Psi(\mu)\rvert H(\mu)% \lvert\Psi(\mu)\rangle$ $\displaystyle=$ $\displaystyle\langle\Psi(\mu)\rvert\frac{1}{2}(r_{1}^{2}+r_{2}^{2})\lvert\Psi(% \mu)\rangle$ (E.4) $\displaystyle=$ $\displaystyle\langle\Psi(\mu)\rvert(R^{2}+\frac{1}{4}r_{12}^{2})\lvert\Psi(\mu)\rangle$

where $\vecfont{R}=(\vecfont{r}_{1}+\vecfont{r}_{2})/2$. As it is possible to separate variables $\vecfont{R}$ and $\vecfont{r}_{12}$ in the Schrödinger equation,

 $\Psi(\mu)=\Psi(\vecfont{R},\vecfont{r}_{12},\mu)=\Phi(\vecfont{R})\phi(% \vecfont{r}_{12},\mu)$ (E.5)

$\Phi$ and $\phi$ both normalized to one, and

 $\langle\Psi(\mu)\rvert\frac{1}{2}(r_{1}^{2}+r_{2}^{2})\lvert\Psi(\mu)\rangle=% \int_{\mathbb{R}^{3}}\rho(\vecfont{r},\mu)r^{2}\,\mathrm{d}\vecfont{r}$ (E.6)

where $\rho(\vecfont{r},\mu)$ is the density of the model system, we have

 $\frac{1}{4}\langle\Psi(\mu)\rvert r_{12}^{2}\lvert\Psi(\mu)\rangle=\int_{% \mathbb{R}^{3}}\rho(\vecfont{r},\mu)\,r^{2}\,\mathrm{d}\vecfont{r}-\int_{% \mathbb{R}^{3}}|\Phi(\vecfont{R})|^{2}R^{2}\,\mathrm{d}\vecfont{R}$ (E.7)

Note that the change of $\langle r_{12}^{2}\rangle$ with $\mu$ is only due to the change of the density with $\mu$. Thus, the difference between $\frac{1}{4}\langle\Psi(\mu)\rvert r_{12}^{2}\lvert\Psi(\mu)\rangle$ and $\frac{1}{4}\langle\Psi\rvert r_{12}^{2}\lvert\Psi\rangle$ 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 ISITE-BFC (contract ANR15-IDEX-0003) (É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] R. D. Amos, A. Bernhardsson, A. Berning, P. Celani, D. L. Cooper, M. J. O. Deegan, A. J. Dobbyn, F. Eckert, C. Hampel, G. Hetzer, P. J. Knowles, T. Korona, R. Lindh, A. W. Lloyd, S. J. McNicholas, F. R. Manby, W. Meyer, M. E. Mura, A. Nicklass, P. Palmieri, R. Pitzer, G. Rauhut, M. Schütz, U. Schumann, H. Stoll, A. J. Stone, R. Tarroni, T. Thorsteinsson, and H.-J. Werner (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] M. Barrault, Y. Maday, N. C. Nguyen, and A. T. Patera (2004-11) An ‘empirical interpolation’ method: application to efficient reduced-basis discretization of partial differential equations. Comptes Rendus Mathematique 339 (9), pp. 667–672 (en). External Links: ISSN 1631-073X, Document Cited by: §1.2, §C.
• [3] G. E. W. Bauer (1983) General operator ground-state expectation values in the hohenberg-kohn-sham density-functional formalism. Phys. Rev. B 27, pp. 5912. Cited by: §2.3.
• [4] F. Colonna and A. Savin (1999-01) Correlation energies for some two- and four-electron systems along the adiabatic connection in density functional theory. The Journal of Chemical Physics 110 (6), pp. 2828–2835. External Links: ISSN 0021-9606, Document Cited by: §A.
• [5] W. E. and D. H. (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] H. Flad, G. Flad-Harutyunyan, and B. Schulze (2020-11) Explicit Green operators for quantum mechanical Hamiltonians. II. Edge-type singularities of the helium atom. Asian-European Journal of Mathematics 13 (07), pp. 2050122 (en). External Links: ISSN 1793-5571, 1793-7183, Link, Document Cited by: §1.1.
• [7] S. Fournais, M. Hoffmann-Ostenof, Th. Hoffmann-Ostenhof, and Th. O. Sorensen (2009) Analytic Structure of Many-Body Coulombic Wave Functions. Commun. Math. Phys. 289, pp. 291. Cited by: §1.1.
• [8] E. Goll, T. Leininger, F. R. Manby, A. Mitrushchenkov, H. Werner, and H. Stoll (2008-06) Local and density fitting approximations within the short-range/long-range hybrid scheme: application to large non-bonded complexes. Physical Chemistry Chemical Physics 10 (23), pp. 3353–3357. External Links: ISSN 1463-9084, Document Cited by: §3.5.
• [9] E. Goll, H. Werner, H. Stoll, T. Leininger, P. Gori-Giorgi, and A. Savin (2006-10) A short-range gradient-corrected spin density functional in combination with long-range coupled-cluster methods: Application to alkali-metal rare-gas dimers. Chemical Physics 329 (1), pp. 276–282 (en). External Links: ISSN 0301-0104, Document Cited by: §2.2.2.
• [10] E. Goll, H. Werner, and H. Stoll (2005-11) A short-range gradient-corrected density functional in long-range coupled-cluster calculations for rare gas dimers. Physical Chemistry Chemical Physics 7 (23), pp. 3917–3923. External Links: ISSN 1463-9084, Document Cited by: §3.5.
• [11] D. H. (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] J. Harris and R. O. Jones (1974) The surface energy of a bounded electron gas-solid. J. Phys. F 4, pp. 1170–1186. Cited by: §1.1.
• [13] P. Hohenberg and W. Kohn (1964-11) Inhomogeneous Electron Gas. Physical Review 136 (3B), pp. B864–B871. External Links: Document Cited by: §A, §B.
• [14] T. Kato (1957) On the eigenfunctions of many-particle systems in quantum mechanics. Communications on Pure and Applied Mathematics 10 (2), pp. 151–177 (en). External Links: ISSN 1097-0312, Link, Document Cited by: §1.1.
• [15] H. F. King (1996-12) The electron correlation cusp. Theoretica chimica acta 94 (6), pp. 345–381 (en). External Links: ISSN 1432-2234, Link, Document Cited by: §3.4.
• [16] W. Kohn and L. J. Sham (1965-11) Self-Consistent Equations Including Exchange and Correlation Effects. Physical Review 140 (4A), pp. A1133–A1138. External Links: Document Cited by: §A.
• [17] Y. I. Kurokawa, H. Nakashima, and H. Nakatsuji (2016) General coalescence conditions for the exact wave functions: higher-order relations for coulombic and non-coulombic systems. Adv. Quantum Chem. 73, pp. 59. Cited by: §1.1.
• [18] M. Levy (1979-12) Universal variational functionals of electron densities, first-order density matrices, and natural spin-orbitals and solution of the v-representability problem. Proceedings of the National Academy of Sciences 76 (12), pp. 6062–6065 (en). External Links: ISSN 0027-8424, 1091-6490, Document Cited by: §A.
• [19] E. H. Lieb (1983) Density functionals for Coulomb systems. International Journal of Quantum Chemistry 24 (3), pp. 243–277 (en). External Links: ISSN 1097-461X, Document Cited by: §A.
• [20] Y. Maday, N. C. Nguyen, A. T. Patera, and S. H. Pau (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] H. Nakashima and H. Nakatsuji (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, $\Sigma$, $\Pi$, $\Delta$, and $\Phi$. J. Chem. Phys. 149, pp. 244116. Cited by: §D.4.
• [22] S. Paziani, S. Moroni, P. Gori-Giorgi, and G. B. Bachelet (2006-04) Local-spin-density functional for multideterminant density functional theory. Physical Review B 73 (15), pp. 155111. External Links: Document Cited by: §3.5.
• [23] V. A. Rassolov and D. M. Chipman (1996) Behavior of electronic wave functions near cusps. J. Chem. Phys. 104, pp. 9908. Cited by: §1.1.
• [24] A. Savin (1996-01) On degeneracy, near-degeneracy 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: Cited by: §A, §3.5.
• [25] A. Savin (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] H. Yserentant (2010) Regularity and approximability of electronic wave functions. Springer. Cited by: §1.1.