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

Yvon Maday
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. 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, H₂ 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 A–E that follow Sec. 4.

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}=|{𝒓}_{𝒊}-{𝒓}_{𝒋}|$ is the distance between electron $i$ (at position ${𝒓}_{𝒊}$) and electron $j$ (at position ${𝒓}_{𝒋}$). Finally, the external potential $V$ is written like

$$V=\sum _{i=1}^{N}v({𝒓}_{𝒊}).$$ (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 =\mathrm{\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 =\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 )=E-E(\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 )=⟨\mathrm{\Psi }\left|\left(W-W(\mu )\right)\right|\mathrm{\Psi }⟩,\text{for }\mu \to \mathrm{\infty }.$$ (2.9)

By changing the integration variables ${𝒓}_{𝒊}$ to $\mu {𝒓}_{𝒊}$ 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,

$${\tilde{\chi }}_j(\mu )=1-j\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 ${\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.

Fig. 1: Basis functions ${\tilde{\chi }}_j$ of Eq. (2.11), continuous curves with the color corresponding to $j$; and an (unknown) function to be approximated by linear combination on this basis (dot-dashed, gray). The unknown function in this figure is proportional to $\overline{E}(\mu )$ of harmonium.

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 $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 $\mathrm{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 $\overline{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\ne H$ obtained with the model wave functions, $\mathrm{\Psi }(\mu )$,

$$A(\mu )=⟨\mathrm{\Psi }(\mu )|A|\mathrm{\Psi }(\mu )⟩.$$ (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

$$⟨\mathrm{\Psi }(\mu )|A|\mathrm{\Psi }(\mu )⟩={\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 $⟨\mathrm{\Psi }(\mu )|A|\mathrm{\Psi }(\mu )⟩$ 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 =\mathrm{\infty })$. The plots show the error in the estimate of $E$ made when considering approximations that use information only for ${\tilde{\mu }}_m\le \mu$. From the plots, we read off how small $\mu$ can be and still have “reasonable” accuracy. In thermochemistry, $\mathrm{kcal}{\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\mathrm{eV}\approx 2.3\mathrm{kcal}{\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 .

The (maximal) parameter $\mu$ is considered between $0$ and $3{\mathrm{bohr}}^{-1}$. The model without correction (blue curve) reaches chemical accuracy for $\mu \approx 3{\mathrm{bohr}}^{-1}$ for H^– and harmonium, but only at $\mu \approx 5{\mathrm{bohr}}^{-1}$ for H₂ in its ground state.

3.2 General behavior of errors

Fig. 2: Errors for harmonium (top), H₂ (middle), and H^– (bottom) using FLEIM with one to four points (1: blue curve, 2: brown curve, 3: green curve, 4: red curve). The abscissa represents the biggest $\mu$ allowed for use in the FLEIM algorithm. The error of the model without correction (blue curve) does not show up in the figure for the H₂ molecule because it is larger than the domain covered by the plot.

The plots in Fig. 2 for harmonium, H₂, 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_{K-1}|$ 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: $⟨r_{12}⟩$ and $⟨r_{12}^2⟩$ for harmonium

Fig. 3: Errors made for the expectation value of the distance between electrons (top) and the distance squared (bottom) for harmonium, by using a model wave function, $\mathrm{\Psi }(\mu )$, and after correcting with FLEIM (the different curves correspond to the number of points used). The insets zooms in.

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\mathrm{pm}$ $\approx$ $0.02\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, $⟨r_{12}^2⟩$, 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, $\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 $⟨r_{12}^2⟩$. 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 $⟨r_{12}^2⟩$ 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 $⟨r_{12}^2⟩$. In the case of harmonium, this can be shown analytically. Changing ${𝒓}_{\mathit{1}}$ and ${𝒓}_{\mathit{2}}$ to center-of-mass, $𝑹$, and inter-particle distance, ${𝒓}_{\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 ${𝒓}_{\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 $⟨r_{12}^2⟩$. 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 =\mathrm{\infty }$.

Fig. 4: Absolute errors for the harmonium molecule at equilibrium distance (top), for the H₂ molecule at equilibrium distance (middle), and in the first excited state of the same symmetry as the ground state (bottom): a $\mu$-dependent LDA (black dashed curve), combined with the a $\mu$-dependent Perdew–Burke–Ernzerhof approximation (PBE, gray dashed curve), FLEIM (3 points) (green curve), combined with a $\mu$-dependent local density approximation. The abscissa represents the biggest $\mu$ allowed for use in the FLEIM algorithm. The insets zoom into the regions of small errors, the dotted line corresponding to the value of “chemical accuracy.” Note the different ranges for $\mathrm{\Delta }E$.

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

$$F[\rho ]\approx {\int }_{{ℝ}^3}f(\rho (𝒓),|\mathbf{\nabla }\rho (𝒓)|,\mathrm{\dots })d𝒓.$$ (A.1)

The limitation of such an approach can be seen for a simple density functional, the Hartree term of the energy,

$$E_H[\rho ]={\int }_{{ℝ}^6}\frac{\rho ({𝒓}_{\mathit{1}})\rho ({𝒓}_{\mathit{2}})}{|{𝒓}_{\mathit{1}}-{𝒓}_{\mathit{2}}|}$$ (A.2)

when $\rho (𝒓)={\rho }_{𝑨}(𝒓)+{\rho }_{𝑩}(𝒓)$ and ${\rho }_A(𝒓){\rho }_{𝑩}(𝒓)\approx \mathit{0}$ (i.e., ${\rho }_A$ and ${\rho }_B$ are spatially separated; their overlap decreases much faster than $1/|{𝒓}_{\mathit{1}}-{𝒓}_{\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 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 ({𝒓}_{\mathit{1}}-{𝒓}_{\mathit{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,

$$⟨\mathrm{\Psi }(\mu )\left|\left(W-W(\mu )\right)\right|\mathrm{\Psi }(\mu )⟩=\frac{1}{2}{\int }_{{ℝ}^6}\rho ({𝒓}_{\mathit{1}})\rho ({𝒓}_{\mathit{2}})\left(\frac{\mathit{1}}{{𝒓}_{\mathit{12}}}-𝒘({𝒓}_{\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 one-particle basis functions, so $\rho$ is a linear combination of Gaussian functions. We consider a generic term,

$${\int }_{{ℝ}^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$.

Fig. 5: Errors for harmonium using FLEIM with up to four points (1: blue curve, 2: brown curve, 3: green curve, 4: red curve), using the basis sets given in (B.3) (top) and (B.4) (bottom). The abscissa represents the biggest $\mu$ allowed for use in the FLEIM algorithm.

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,\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

$${\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 $𝒞$, e.g., a norm $𝒞[\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 $𝒞[\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$,

$${\tilde{\chi }}_0≔{\chi }_0.$$ (C.1)

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

$${\tilde{\mu }}_0≔\underset{j\in J}{\max }{\mu }_j.$$ (C.2)

We can then define the first normalized interpolation function as

$$q_0≔\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

$${ℐ}_0[f]≔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,\mathrm{\dots },K-2$. We define the ${ℐ}_{K-2}$ Lagrangian interpolation function as

$${ℐ}_{K-2}[f]≔\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 }}_{\mathrm{\ell }})=f({\tilde{\mu }}_{\mathrm{\ell }})\mathit{\hspace{1em}}\text{for }\mathrm{\ell }=0,\mathrm{\dots },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}≔\underset{{\chi }_i,i\in \tilde{I}}{\arg \max }\left\{𝒞\left[{\chi }_i-{ℐ}_{K-2}[{\chi }_i]\right]\right\},$$ (C.7)

and the next interpolation point as

$${\tilde{\mu }}_{K-1}≔\underset{{\mu }_j,j\in \tilde{J}}{\arg \max }\left\{\left|{\tilde{\chi }}_{K-1}({\mu }_j)-{ℐ}_{K-2}[{\tilde{\chi }}_{K-1}]({\mu }_j)\right|\right\}.$$ (C.8)

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

$$q_{K-1}≔\frac{{\tilde{\chi }}_{K-1}-{ℐ}_{K-2}[{\tilde{\chi }}_{K-1}]}{{\tilde{\chi }}_{K-1}({\tilde{\mu }}_{K-1})-{ℐ}_{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})≔\underset{{\chi }_i,i\in \tilde{I}}{\arg \max }\underset{{\mu }_j,j\in \tilde{J}}{\arg \min }\left\{𝒞\left[{\chi }_i({\mu }_j)-{ℐ}_{K-2}[{\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 )$

Fig. 6: Test of EIM using model $E(\mu )=1+{\chi }_j(\mu )$.

Fig. 7: Test of FLEIM using model $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(\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{\mathrm{bohr}}^{-1}$, the FLEIM results are only slightly changed, as shown in Fig. 8.

Fig. 8: Results using a denser grid (lighter colors), for harmonium, H₂ at equilibrium distance, and H^–; the curves obtained with the denser grid are shown in lighter colors. The error of the model without correction (blue curve) does not show up in the figure for the H₂ molecule because it is larger than the domain covered by the plot.

D.3 Systems

Two-electron systems studied are

1. harmonium, having

   $$v(𝒓)=\frac{\mathit{1}}{\mathit{2}}{\omega }^{\mathit{2}}{𝒓}^{\mathit{2}}$$ (D.1)

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

2. H^– anion,

   $$v(𝒓)=-\frac{\mathit{1}}{𝒓};$$ (D.2)

3. H₂ molecule,

   $$v(𝒓)=-\frac{\mathit{1}}{|𝒓-{𝑹}_{𝑨}|}-\frac{\mathit{1}}{|𝒓-{𝑹}_{𝑩}|}$$ (D.3)

   where the nuclei are in the equilibrium position, ${𝑹}_{𝑨}=-{𝑹}_{𝑩}$ with $|{𝑹}_{𝑨}|=\mathit{0.7}\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 H₂. 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\mathrm{bohr}$.

For H₂ 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 H₂, as it has an important contribution of ionic states (H⁺$\mathrm{\cdots }$H^–). For H₂ the equilibrium distance of $1.4\mathrm{bohr}$ was chosen for the ground state. For the excited state, the distance of $4.2\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 ${𝒓}_{\mathit{1}},{𝒓}_{\mathit{2}}$ can be changed to those corresponding to the center of mass and the distance between particles,

$𝑹$ $=$ ${\displaystyle \frac{1}{2}}\left({𝒓}_{\mathit{1}}+{𝒓}_{\mathit{2}}\right),$

${𝒓}_{\mathit{12}}$ $=$ $\left({𝒓}_{\mathit{1}}-{𝒓}_{\mathit{2}}\right)$ (E.1)

yielding for the potential energy

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

Note that $⟨\mathrm{\Psi }(\mu )|r_{12}^2|\mathrm{\Psi }(\mu )⟩$ 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,

${\displaystyle \frac{\partial }{\partial {\omega }^2}}⟨\mathrm{\Psi }(\mu )|H(\mu )|\mathrm{\Psi }(\mu )⟩$ $=$ $⟨\mathrm{\Psi }(\mu )|{\displaystyle \frac{1}{2}}(r_1^2+r_2^2)|\mathrm{\Psi }(\mu )⟩$ (E.4)

$=$ $⟨\mathrm{\Psi }(\mu )|(R^2+{\displaystyle \frac{1}{4}}{r}_{12}^2)|\mathrm{\Psi }(\mu )⟩$

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

$$\mathrm{\Psi }(\mu )=\mathrm{\Psi }(𝑹,{𝒓}_{\mathit{12}},\mu )=𝜱(𝑹)\varphi ({𝒓}_{\mathit{12}},\mu )$$ (E.5)

$\mathrm{\Phi }$ and $\varphi$ both normalized to one, and

$$⟨\mathrm{\Psi }(\mu )|\frac{1}{2}(r_1^2+r_2^2)|\mathrm{\Psi }(\mu )⟩={\int }_{{ℝ}^3}\rho (𝒓,\mu ){𝒓}^{\mathit{2}}\mathrm{d}𝒓$$ (E.6)

where $\rho (𝒓,\mu )$ is the density of the model system, we have

$$\frac{1}{4}⟨\mathrm{\Psi }(\mu )|r_{12}^2|\mathrm{\Psi }(\mu )⟩={\int }_{{ℝ}^3}\rho (𝒓,\mu ){𝒓}^{\mathit{2}}\mathrm{d}𝒓-{\int }_{{ℝ}^{\mathit{3}}}|𝜱(𝑹){|}^{\mathit{2}}{𝑹}^{\mathit{2}}\mathrm{d}𝑹$$ (E.7)

Note that the change of $⟨r_{12}^2⟩$ with $\mu$ is only due to the change of the density with $\mu$. Thus, the difference between $\frac{1}{4}⟨\mathrm{\Psi }(\mu )|r_{12}^2|\mathrm{\Psi }(\mu )⟩$ and $\frac{1}{4}⟨\mathrm{\Psi }|r_{12}^2|\mathrm{\Psi }⟩$ 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).