Re-using model results to determine materials properties: connector theory approach

TL;DR

This paper introduces a general connector theory that uses model results to accurately approximate real material properties, enabling new functional design and simplifying complex calculations.

Contribution

It presents a novel, exact theoretical framework for reusing model results in materials property predictions through the concept of connectors.

Findings

Connector theory provides a rigorous basis for approximations.

The approach can be used to design new density functional functionals.

Illustrative examples demonstrate practical applications.

Abstract

Computational materials design often profits from the fact that some complicated contributions are not calculated for the real material, but replaced by results of models. We turn this approximation into a very general and in principle exact theory by introducing the concept of a connector, which is a prescription of how to use the results of a model system in order to simulate a real system. We set the conditions that must be fulfilled for the existence of an exact connector. We demonstrate that, and why, this approach is a very convenient starting point for approximations. We also show that the connector theory can be used to design new functionals, for example for density functional theory. We illustrate our purposes with simple but realistic examples.

Figures

Equations

$$O(x;Q^R)=O(x;Q_x^c).$$

$$O(x;Q^R)={\mathcal{O}}_x({\mathcal{Q}}_x^c).$$

$${\mathcal{Q}}_x^c={\mathcal{O}}_x^{-1}(O(x;Q^R)).$$

$$O(x;Q^R)={\mathcal{O}}_x({\mathcal{Q}}_x^c).$$

$${\mathcal{Q}}_x^{c,\mathrm{approx}}=({\mathcal{O}}_x^{\mathrm{approx}}{)}^{-1}(O^{\mathrm{approx}}(x;Q^R)).$$

$$O(x;Q^R)\approx O_x^c\equiv \mathcal{O}({\mathcal{Q}}_x^{c,\mathrm{approx}}).$$

$$V(z)=\{\begin{array}{ll}\frac{1}{2}{\omega }_0^2{\left(z+\frac{L}{2}\right)}^2& z<-\frac{L}{2}(\mathrm{I})\\ 0& -\frac{L}{2}\le z\le \frac{L}{2}(\mathrm{II})\\ \frac{1}{2}{\omega }_0^2{\left(z-\frac{L}{2}\right)}^2& z>\frac{L}{2}(\mathrm{III}),\end{array}$$

$$\bigl{[}-\frac{1}{2}\frac{d^{2}}{dy^{2}}+v(y)\bigr{]}\phi(y)=\varepsilon\phi(y)\,,$$

$$E^{\mathrm{approx}}(j;{\omega }_0,L)={\mathcal{E}}_j(L)\left[1-\frac{2\Gamma (\frac{1}{4})}{\Gamma (\frac{3}{4})}\sqrt{\alpha ({\omega }_0,L)}\right],$$

$$L_{j{\omega }_{0}L}^c=\frac{\pi j}{\sqrt{2E(j;{\omega }_0,L)}}.$$

$${\mathcal{E}}_j^{\mathrm{approx}}(\tilde{L})$$

$$={\mathcal{E}}_j(L)\left[1-2\frac{\tilde{L}-L}{L}\right].$$

$$L_{{\omega }_{0}L}^{c,\mathrm{approx}}=L\left[1+\frac{\Gamma (\frac{1}{4})}{\Gamma (\frac{3}{4})}\sqrt{\alpha ({\omega }_0,L)}\right].$$

$$E^{c}(j;\omega_{0},L)=\mathcal{E}_{j}\bigl{(}L^{c,{\rm approx}}_{\omega_{0}L}\bigr{)}=\frac{\mathcal{E}_{j}(L)}{\left[1+\frac{\Gamma(\frac{1}{4})}{\Gamma(\frac{3}{4})}\sqrt{\alpha(\omega_{0},L)}\right]^{2}}.$$

$$G_{\mathbf{K}\mathbf{K}}(\mathbf{k},\omega ;[V])={\mathcal{G}}_{\mathbf{k}\omega \mathbf{K}}({\mathcal{V}}_{\mathbf{k}\omega \mathbf{K}}^c),$$

$${\mathcal{G}}_{\mathbf{k}\omega \mathbf{K}}({\mathcal{V}}^c)=\frac{1}{\omega -|\mathbf{k}+\mathbf{K}{|}^2/2-{\mathcal{V}}^c},$$

$${\mathcal{V}}_{\mathbf{k}\omega \mathbf{K}}^c=\omega -\frac{|\mathbf{k}+\mathbf{K}{|}^2}{2}-\frac{1}{G_{\mathbf{K}\mathbf{K}}(\mathbf{k},\omega ;[V])}.$$

$${\mathcal{V}}_{\mathbf{k}\omega \mathbf{K}}^{capprox}={\mathcal{V}}_{{\mathbf{k}}_0\omega \mathbf{K}}^c,$$

$$G_{\mathbf{K}\mathbf{K}}^{capprox}(\mathbf{k},\omega ;[V])=\frac{1}{-\frac{|\mathbf{k}+\mathbf{K}{|}^2}{2}+\frac{|{\mathbf{k}}_0+\mathbf{K}{|}^2}{2}+\frac{1}{G_{\mathbf{K}\mathbf{K}}({\mathbf{k}}_0,\omega ;[V])}}.$$

$$v_{\mathrm{xc}}(\mathbf{r},[n])=v_{\mathrm{xc}}^h(n_{\mathbf{r}}^c),$$

$$n_{\mathbf{r}}^{c,\mathrm{approx}}=(v_{\mathrm{xc}}^{h,\mathrm{approx}}{)}^{-1}(v_{\mathrm{xc}}^{\mathrm{approx}}(\mathbf{r};[n])).$$

$$E_{\rm xc}^{\rm WDA}[n]=\int d\mathbf{r}d\mathbf{r^{\prime}}\frac{n(\mathbf{r})n(\mathbf{r^{\prime}})}{2\mathbf{|r-r^{\prime}|}}C(\tilde{n}(\mathbf{r,r^{\prime}}))\Big{(}1-e^{\frac{-\lambda(\tilde{n}(\mathbf{r,r^{\prime}}))}{\mathbf{|r-r^{\prime}|^{5}}}}\Big{)},$$

$$v_{\mathrm{xc}}^c(\mathbf{r},[n])=v_{\mathrm{xc}}^h\left(n^h=n_{\mathbf{r}}^{c,\mathrm{approx}}([n])\right),$$

$$v_{\mathrm{xc}}^h(n_0)+\int d{\mathbf{r}}^{\prime }(n({\mathbf{r}}^{\prime })-n_0)f_{\mathrm{xc}}(|\mathbf{r}-{\mathbf{r}}^{\prime }|;n_0)=v_{\mathrm{xc}}^h(n_0)+(n_{\mathbf{r}}^{c,\mathrm{approx}}-n_0)f_{\mathrm{xc}}^h(n_0),$$

$$n_{\mathbf{r}}^{c,\mathrm{approx}}([n])=\frac{1}{f_{\mathrm{xc}}^h(n_0)}\int d{\mathbf{r}}^{\prime }n({\mathbf{r}}^{\prime })f_{\mathrm{xc}}(|\mathbf{r}-{\mathbf{r}}^{\prime }|;n_0),$$

$$v_{\mathrm{xc}}(\mathbf{r})$$

$$v_{\mathrm{x}}(\mathbf{r})$$

$$v_{\mathrm{c}}(\mathbf{r})$$

$$E_{xc}[n]\equiv \int d\mathbf{r}n(\mathbf{r}){ϵ}_{xc}(n_{\mathbf{r},e}^c).$$

$$v_{\mathrm{xc}}(\mathbf{r})$$

$$v_{xc}^c(\mathbf{r})={ϵ}_{xc}(n_{\mathbf{r},e}^c)+n_{\mathbf{r},e}^c\frac{d{ϵ}_{xc}(n_{\mathbf{r},e}^c)}{d{n}_{\mathbf{r},e}^c},$$

$$n_{\mathbf{r}}^c$$

Full text

I Introduction

Computational materials design Curtarolo et al. (2013); Saal et al. (2013); Marzari (2016) is complicated by the complexity of materials and by interaction effects. This hampers both calculations and understanding. The fundamental problem lies in the fact that the effects of the Coulomb interaction and of the specific material cannot be separated. Otherwise, one could calculate the interaction contributions once and for all, store them and add them every time a new material is calculated. This is an unreachable dream, but it still indicates an intriguing direction of thinking: in some model systems the effects of the Coulomb interaction can be treated exactly, or at least to a much better extent than in real systems, and attempts to use model results in order to simulate interaction effects in real materials are widespread. The most prominent example is the local density approximation (LDA) to the Kohn-Sham (KS) exchange-correlation (xc) potential $v_{\rm xc}({\bf r};[n])$ of density functional theory (DFT) Kohn and Sham (1965). The exact potential is unknown in most real materials. The LDA replaces $v_{\rm xc}({\bf r};[n])$ at a point ${\bf r}$ by its value in a homogeneous electron gas (HEG) that is calculated at the density $n({\bf r})$ of the real system in the same point ${\bf r}$. In this way, DFT profits from the existence of tabulated and interpolated Quantum Monte Carlo results Ceperley and Alder (1980). Similarly, dynamical mean field theory in the single site approximation takes the effective local self-energy from the Anderson impurity model, and although in this case results have not been tabulated, the procedure has enabled a realistic description of correlated materials Georges et al. (1996). However, in spite of numerous studies and attempts Kohn and Sham (1965); Fermi (1927); Thomas (1927); Ma and Brueckner (1968); Vosko et al. (1980); Gunnarsson et al. (1977); Alonso and Girifalco (1977, 1978); Gunnarsson et al. (1979); Gunnarsson and Jones (1980); Cuevas-Saavedra et al. (2012); Langreth and Perdew (1980); Perdew and Yue (1986); Perdew et al. (1996, 1999, 2008); Sun et al. (2015); Dobson (1998); Jung et al. (2004); García-González et al. (1996); Olsen and Thygesen (2012, 2014); Patrick and Thygesen (2015); Schmidt et al. (2017); Lu (2014); Vignale and Rasolt (1987); Zangwill and Soven (1980); Gross and Kohn (1985); Dobson et al. (1997); Ng (1989); Vignale and Kohn (1996); Vignale et al. (1997); Tokatly (2007); Gao et al. (2010); Nazarov et al. (2010); Trevisanutto et al. (2013); Panholzer et al. (2018); Giuliani and Vignale (2005); Martin (2004) it is very difficult to go beyond these simple schemes. One reason is that using results of the models is considered from the very beginning as an approximation. Breakthrough is instead often based on an exactification: a recent example is the exact factorization of the many-body wavefunction Abedi et al. (2010, 2012); Gidopoulos and Gross (2014) that includes the Born-Oppenheimer approximation as limiting case Eich and Agostini (2016) and allows developments beyond it Min et al. (2015). In a similar spirit, one may try to formulate a general, in principle exact, approach that uses model results for building a density functional for $v_{\rm xc}({\bf r};[n])$, of which the LDA would be one particular approximation.

We keep this idea as broad as possible and therefore pose the following questions: Can one exactify the idea of re-using results from one system, for example a model, to describe another system? If yes, under which conditions? And does this suggest strategies for systematic approximations? These are the questions at the heart of our work. Our answer, termed connector theory (COT), tells how one can in principle connect different systems, and how to find good approximations in practice. It is very general and designed to overcome numerous obstacles, not only the problem of interaction. It has potential for great speedup of calculations, and, as we will demonstrate, it is a powerful tool to design improved functionals, for example, but not exclusively, for DFT.

II Formalism

The idea can be summarized as follows: suppose one wishes to calculate an observable or quantity $O$ that is a function or functional of a set of parameters or function $Q$, and it can depend on additional arguments $x$, so $O=O(x;Q)$. For example, if $O$ is the xc potential $v_{\rm xc}({\bf r};[n])$, then $x$ is a spatial coordinate ${\bf r}$ and $Q$ is a function, the density $n({\bf r}^{\prime})$. Or, $O$ could be a spectral function, $x$ the frequency, and $Q$ the atomic potentials, etc. In most cases $O(x;Q)$ is difficult or impossible to calculate in a real material without approximations, or even unknown. However, it may be possible to calculate $O(x;Q)$ for some $Q$ in a restricted domain: this restricted domain defines a model. For example, when $Q$ is the electron density and we restrict it to homogeneous densities, the model is the HEG. Connector theory aims at using the model results in order to simulate systems where $Q$ lies outside the model domain. The goal is therefore to find, for a given real system where $Q\equiv Q^{R}$, another $Q\equiv Q^{c}$ that lies in the model domain, such that

[TABLE]

Note the subscript $x$ of $Q^{c}_{x}$, which indicates that $Q^{c}$ is allowed to be different for every value of the argument $x$.

In the following we suppose that, because the model is simple, it is described by only one effective parameter $\mathcal{Q}$ or, when there are more parameters, we can choose one that will be used to fulfill (1). This restriction can be dropped, but it is often useful and we keep it here for clarity. In this case, in the model the functional or multi-dimensional function can be represented by a simple function, $O(x;Q^{c}_{x})\to\mathcal{O}_{x}(\mathcal{Q}^{c}_{x})$. For example, in the HEG the one parameter $\mathcal{Q}$ is its number density $n^{h}$, so $v_{\rm xc}({\bf r};[n])\to v_{\rm xc}^{h}(n^{h})$, where $v_{\rm xc}^{h}$ corresponds to $\mathcal{O}$ and $n^{h}$ to $\mathcal{Q}$. Eq. (1) reads then

[TABLE]

Of course, even allowing $\mathcal{Q}^{c}_{x}$ to depend on $x$ it could be impossible to fulfill equality (2) if one restricts the model domain too severely. However, if the model is flexible enough such that the equation can be satisfied in principle, one can try to find the one or more $\mathcal{Q}^{c}_{x}$ for which the equality holds. In other words, for the idea to be applicable the following condition must be satisfied:

• •

[A] On its domain of definition, $\mathcal{O}_{x}(\mathcal{Q}^{c}_{x})$ must yield all values that $O(x;Q^{R})$ can take on its domain, i.e., the $Q^{R}$ of interest.

The $Q^{R}$ which define the domain of interest depend on the range of physical systems one wants to explore; this range does not necessarily include all possible physical systems. The domain of $O$, on the other hand, defines the model system. If for certain $Q^{R}$ and/or $x$ Eq. (2) cannot be fulfilled, we have to change model by changing its domain, i.e. the range of allowed $\mathcal{Q}^{c}_{x}$.

Eq. (2) is then formally solved,

[TABLE]

This operation requires inversion of $\mathcal{O}_{x}$, which brings us to a second condition. Indeed, ${\bf[A]}$ is the only necessary condition, but there is also a question of uniqueness in (3): $\mathcal{O}^{-1}$ may require boundary conditions in order to be well defined. This is not a problem of principle, but may create difficulties for the design of approximations. We therefore require:

• •

[B] When the inverse $\mathcal{O}^{-1}$ of $\mathcal{O}$ is not unique, it should at least be possible to specify a unique choice among the possibles $\mathcal{O}^{-1}_{i}$ .

Finally, the resulting $Q^{c}_{x}$ is used to get the observable in the real system using the model,

[TABLE]

Therefore, we call $\mathcal{Q}^{c}_{x}$ the connector. The scheme suggests to store $\mathcal{O}_{x}$ as numerical data, interpolated analytic expression or, in some cases, analytic result. Once a connector $\mathcal{Q}^{c}_{x}$ is given, for any real system one can then simply use these data instead of calculating $O(x;Q^{R})$, which will make calculations extremely efficient. This concept of re-using data has made the LDA a breakthrough, since numerous calculations could be performed without ever re-evaluating the interaction effects in the HEG.

In practice, however, conditions [A] and [B] settle the framework, but nothing has been gained yet: the unknown $O(x;Q^{R})$ still enters the calculation of the connector in (3). Therefore, for COT to be of practical interest we have to add:

• •

[C] The COT ansatz must suggest a strategy for approximations.

The LDA can be seen as a connector approximation, where the right hand side of (3) is replaced by the local density, $\bar{n}^{c}_{\bf r}\approx n({\bf r})$, so $v_{\rm xc}({\bf r};[n])\approx v_{\rm xc}^{h}(\bar{n}^{c}_{\bf r})$ according to (4). It was based on Kohn’s ingenious intuition of nearsightedness Kohn and Sham (1965); Kohn (1996); Prodan and Kohn (2005). A major topic of the present work is to lie out a strategy for systematic approximations, and to generalize the framework to applications well beyond DFT.

For this, we have to prove that (1) is a clever starting point. The fundamental idea is to make in the calculation of the connector (3) an equivalent approximation to $O$ and $\mathcal{O}$,

[TABLE]

The final approximate connector result is obtained as

[TABLE]

Note that the model observable $\mathcal{O}$ is supposed to be well known, but it is important to approximate it in (5) in the same way as $O$ of the real system, whereas the exact model function $\mathcal{O}$ is used in (6). In this way, the result becomes exact in two limits: when the approximation is increasingly good, even for a very restricted model, and when the domain of the model system tends towards that of the real system, even using a very poor approximation. In practice, one will not be able to do the exact calculation, but away from these limits, using the equivalent approximation for $O^{\rm approx}(x;Q^{R})$ and $\mathcal{O}^{\rm approx}_{x}(\mathcal{Q})$ still leads to error cancelling, and the limiting behaviour indicates that results can be improved in a controllable way. How far the model system can be chosen from the real system depends on the quality of the approximation, and vice versa, how rough the approximation is that one can tolerate depends on the closeness of the model and the real system. This double dependence is a source of the power of the connector approach. It implies that one can expect the approximate connector result (6) to be superior to the direct approximation $O^{\rm approx}(x;Q^{R})$, for similar computational cost. This benefit will be larger when the model contains important features of the real system, like in the examples below, the same kinetic energy, or the same Coulomb interaction. Therefore, in COT the main effort and intuition will go into the choice of a suitable model.

III A simple illustration: energy levels in a 1D potential

Let us start with a very simple case to illustrate both the generality of the idea well outside DFT-LDA, and the improvement obtained by using a connector. Take one electron in a one-dimensional potential of shape

[TABLE]

with $\omega_{0}$ and $L$ real positive parameters. We use atomic units throughout this work. For $L\to 0$ this is the oscillator potentialBalestri (2013) $V(z)=\frac{1}{2}\omega_{0}^{2}z^{2}$. The larger $\omega_{0}$, the steeper the walls, and for $\omega_{0}\to\infty$, the system turns into the infinite potential well $V_{\infty}(z)=V_{0}\theta\left(\left|z\right|-\frac{L}{2}\right)$ with $V_{0}\to\infty$. The potential is depicted in the top panel of Fig. 1 for three values of $\alpha\equiv 1/L^{2}\omega_{0}$, which is the adimensional parameter characterizing the problem, as explained below.

The $j^{th}$ energy level $E(j;\omega_{0},L)$ is our observable of interest $O(x;Q)$. Let us first look at the exact problem. For $L\neq 0$ the Schrödinger equation \bigl{[}-\frac{1}{2}\frac{d^{2}}{dz^{2}}+V(z)\bigr{]}\psi(z)=E\psi(z) can be written

[TABLE]

with rescaled coordinate $y{\equiv}z/L$, rescaled energy $\varepsilon{\equiv}L^{2}E$ and effective potential $v(y)=\left(y\pm\frac{1}{2}\right)^{2}/(2\alpha^{2})$ and $v(y)=0$, for $|y|>\frac{1}{2}$ and $|y|<\frac{1}{2}$, respectively. The parameter $\alpha=1/L^{2}\omega_{0}$ is a measure of the difference between the real system and the infinite potential well, which vanishes for $\alpha\to 0$.

To obtain the exact solution, the Schrödinger equation is solved separately in each of the three regions (I), (II) and (III), and the solutions are glued by requiring the wavefunction and its first derivative to be continuous at the zone boundaries $y=\pm\frac{1}{2}$. This has to be done numerically; details can be found in App. A. The result is shown in the bottom panel of Fig. 1 for the first energy levels $E(j;\omega_{0},L)$ as function of $\alpha(\omega_{0},L)$. For $L$ fixed and $\alpha\to 0$, one obtains the eigenvalues of the infinite potential well $\mathcal{E}_{j}(L)=\pi^{2}j^{2}/(2L^{2})$, whose separation increases linearly with $j$. For increasing $\alpha$, the eigenvalues are more and more evenly spaced, like for the harmonic oscillator.

Since the numerical solution is tedious, in practice one often uses perturbation theory to approximate the solutions of a Schrödinger equation. This would correspond to evaluating $O^{\rm approx}(x;Q^{R})$ using a perturbation approximation, taking this to be the final result, without using it subsequently in (5) and (6). We will call this in the following the direct approximation. For this perturbation expansion we choose as zero order an infinite well with width $L_{0}=L$. The expansion parameter is then $\alpha$, as defined above. The expansion to first order is carried out in App. A; the result is

[TABLE]

with $\Gamma$ the Gamma function. The first-order approximation has a very small computational cost compared to the exact numerical one. As Fig. 1 shows, and as expected from perturbation theory, the approximate energy levels are close to the exact $E(j;\omega_{0},L)$ for small $\alpha$ and deviate for larger values, becoming unphysically negative above a critical threshold.

The question is now whether COT, using the same approximation but combining it with knowledge from some model system, may yield better results. It is intuitive to suppose that the effect of a finite curvature $\omega_{0}$ can be simulated by varying the width $L$ of an infinite potential well, i.e., to suppose that each level sees a surrounding with some effective width. As model we will therefore choose all possible infinite potential wells, which means that we restrict the model domain to potentials with $\alpha\to 0$. COT tells us to search for the effective width that best represents the real system, i.e., we will search for the connector $L^{c}$ for which $E(j;\omega_{0},L)=\mathcal{E}_{j}(L^{c}_{j})$111Note that here we connect the energy levels one by one, which is a choice. One would make this choice also to connect, for example, poles of Green’s functions. Vanzini et al. (2018). Note that here the infinite potential well plays a double role: it serves as model, and as the zero order of the expansion. This possible double role of a model in the case of a perturbative approximation is further discussed in Sec. VI.1.

To use COT we have to fulfill condition [A]: this is easily achieved for any potential with positive eigenvalues if we accept to have a different connector, i.e., a different effective width, for every energy level, since the energy levels in an infinite potential well can take any positive value. The exact connector (3) would be $L^{c}_{j\omega_{0}L}=\mathcal{E}^{-1}_{j}(E(j;\omega_{0},L))$. Since the energies $\mathcal{E}_{j}(\tilde{L})$ of an infinite potential well with width $\tilde{L}$ go as $1/{\tilde{L}}^{2}$, there is no unique inverse $\mathcal{E}^{-1}_{j}(E)$. However, one of the two solutions of the resulting quadratic equation for $L^{c}_{j}$ is negative and can be discarded by using positivity as exact constraint, so our example also fulfills condition [B]. Finally, the exact connector reads:

[TABLE]

As expected, to find the exact connector would require calculation of the exact energy levels $E$, which is what one wants to avoid. Instead, the COT approximation strategy, required as condition [C], is to apply an equivalent approximation to both the real and the model system following (5). Therefore, although the infinite potential well can easily be solved exactly, we make the perturbation expansion also for the energy levels of the model. The expansion of an infinite well with width $\tilde{L}$ around width $L$ reads:

[TABLE]

Now we apply (5), i.e., we solve for the connector $L^{c,{\rm approx}}$ for which $E^{\rm approx}(j;\omega_{0},L)=\mathcal{E}_{j}^{\rm approx}(L^{c,{\rm approx}})$ using Eqs. (9) and (10), which yields

[TABLE]

Note that the model energy levels and first derivative are known as analytical functions, so there is no need to store and read data. Through $\alpha$, this approximate connector width $L^{c,{\rm approx}}$ depends on $L$ and $\omega_{0}$ but not on $j$: here, the connector is transferable, i.e., it could be calculated for one energy level and then used for all others. Such a transferability is not exact in general (see also the discussion in Subsec. V.3), but could be used as an approximation, as we will illustrate below in the next example. Eq. (11) also illustrates how the difference between $L^{c,{\rm approx}}$ and $L$ becomes smaller and smaller as $\alpha$ tends to zero, namely as the real system approaches the model system itself.

Once the connector has been obtained, the energy levels are taken from the exact solution of the model system following (6):

[TABLE]

Note again that in the last step the exact model function $\mathcal{E}_{j}$ is used, as prescribed by (6). The result (12) is different from the direct approximation (9): the latter is the first-order expansion of (12). Using the COT expression (12) does not add computational cost, but results are improved significantly, as shown in Fig. 1. They are always better than the direct approximation $E^{\rm approx}(j)$, and they remain physical over the whole parameter range, which reflects an implicit approximate resummation of perturbation theory to infinite order through the use of the exact model $\mathcal{E}_{j}$; see Sec. VI.1 for a general discussion of this point. The connector $L^{c}_{j}$ formalizes the intuitive interpretation of an effective width that determines the energy levels in the real system and shows that this effective width is approximately independent of $j$. This result indicates a promising route to obtain quick estimates for energy levels.

IV Quick estimates from a connector approximation

The toy example of the one-dimensional potential indicates that connector approximations may be used to speed up calculations or, equivalently, to improve the results of rough approximations without additional computational cost. In the following, we will examine this route by looking at a more realistic situation, using as example the calculation of a band structure. Here, the focus will not be to improve the accuracy, e.g., by simulating many-body effects, but to obtain quick results for a given effective hamiltonian.

IV.1 Connector approximation to calculate an independent-particle band structure

The perturbative connector approximation suggested by the simple example above could in principle be extended to real systems. Here, we will investigate another route, since it is important to note that COT approximations are not limited to perturbation theory. We will therefore examine energy levels for real materials, by changing the approximation: inspired by the perfectly transferable connector of above, we will use transferability as approximation.

As illustration, we will apply COT to obtain in an efficient way the KS band structure of solid cubic He Schuch and Mills (1961), with given KS potential $V({\bf r})$. For a material with a periodic potential the HEG is the simplest possible model. In the HEG, the band structure is the free electron parabola folded into the Brillouin zone, and its potential $V^{c}$ is constant in space, characterized by one number $\mathcal{V}^{c}$ $\in\mathbb{R}$ with respect to the chemical potential. This potential allows us to lower or raise the HEG energy levels at will - so condition [A] is satisfied -, and it will play the role of the connector.

For many materials it may be difficult to find a one-to-one relation between their bands and the HEG ones. Therefore, we do not connect bands but derive the band structure from the spectral function, which is the trace over the imaginary part of the KS Green’s function $G({\bf k},\omega)$. To be precise, we impose the connector equality Eq. (4) as

[TABLE]

where $G$ is the KS Green’s functions of the real material, and ${\bf K}$ are reciprocal lattice vectors. The Green’s function of the HEG is diagonal in ${\bf K}$, so we can represent it by $\mathcal{G}_{{\bf k}\omega{\bf K}}$, which is a function of ${\bf k}+{\bf K}$. It is also a function of the value $\mathcal{V}^{c}$ of the homogeneous potential,

[TABLE]

where the frequency $\omega$ includes an infinitesimal imaginary part. This expression is readily calculated as function of $(\omega-\mathcal{V}^{c})$ and $|{\bf k}+{\bf K}|$ and could be either stored or recalculated on the fly. From Eq. (13), the exact connector potential corresponding to (3) is

[TABLE]

It is important to note that for our purpose it is enough to ask for the diagonal $G_{{\bf K}{\bf K}}$: the HEG would not be able to simulate off-diagonal elements, at least not in a physically reasonable way, and therefore asking for off-diagonal elements would violate condition [A].

As one can see, for the band structure the wavevector ${\bf k}$ plays the role of $x$, and the connector potential is in principle different for each wavevector ${\bf k}$. Moreover, it should depend on the band; here, this is translated into a $\omega$-dependence.

Strictly speaking, since the spectral function is the imaginary part of the trace (here, a sum over ${\bf K}$), we could have asked for $\sum_{\bf K}{\rm Im}\,G_{{\bf K}{\bf K}}({\bf k},\omega;[V])=\sum_{\bf K}{\rm Im\,}\mathcal{G}_{{\bf k}\omega{\bf K}}(\mathcal{V}^{c}_{{\bf k}\omega})$, with a ${\bf K}$-independent connector potential $\mathcal{V}^{c}_{{\bf k}\omega}$. However, to invert this relation and select a unique answer would be much more cumbersome than to invert Eq. (13), for which condition [B] is easily satisfied.

The exact connector at a point ${\bf k}_{0}$ can be found by calculating $G_{{\bf K}{\bf K}}({\bf k}_{0},\omega;[V])$ and using (15). To get the band structure, one should do this for every ${\bf k}$, which would again be useless, as it would be equivalent to a full solution of the original problem. We therefore move on to our approximation. As anticipated, it will consist in supposing approximate transferability of the connector through cancellations in Eq. (5). Our approximation is hence to neglect the ${\bf k}$-dependence of the connector by using at all ${\bf k}$ its value at ${\bf k}_{0}$. Of course, applied directly to the band structure or spectral function, transferability would be a bad approximation, since the direct approximation would read $G_{{\bf K}{\bf K}}({\bf k},\omega;[V])\approx G_{{\bf K}{\bf K}}({\bf k}_{0},\omega;[V])$ and therefore yield only flat bands. Also in the COT scheme, we diagonalize the KS hamiltonian only once, at ${\bf k}_{0}=\Gamma$, but now this yields the exact $\mathcal{V}^{c}_{\Gamma\omega{\bf K}}$, and we then approximate

[TABLE]

according to the strategy required by condition [C]. Following (6), this yields the band structure at all ${\bf k}$-points as

[TABLE]

This corresponds to a diagonal approximation of the in principle exact Dyson equation that relates $G_{\bf k}$ and $G_{{\bf k}_{0}}$. The result of this straightforward COT approximation for the case of solid cubic He Schuch and Mills (1961) is shown in the left panel of Fig. 2 (see footnote 222The reference band structures of cubic helium and hexagonal hydrogen, with experimental lattice parametersSchuch and Mills (1961); Ishmaev et al. (1983), have been obtained from the diagonalization of the KS-LDA Hamiltonian of size 5.1 Ha and 3.4 Ha in reciprocal space, respectively. The Hamiltonian has been calculated using local Troullier-MartinsTroullier and Martins (1991) pseudopotentials in the Abinit codeGonze et al. (2005) with $\Gamma$-centered 5$\times$5$\times$5 and 3$\times$3$\times$3 $\mathbf{k}$-point grids, respectively. for computational details).

The straightforward connector result is blurred and suggests spurious energy levels, since it violates an exact constraint, namely, that the poles of the Green’s function should not depend on ${\bf K}$. Still, one can identify bands approximately. We can therefore impose the exact constraint by choosing for each group of poles the one with the largest oscillator strength: this reflects the situation where the real system and the HEG are most similar, which means, where COT should work best. Details of how this is done in practice are given in App. B. The result is given by the blue dots in the left panel of Fig. 2: the constrained connector makes the bands sharp and eliminates some spurious non-dispersing structures.

The right panel of Fig. 2 compares the ab initio and the constrained connector results. They are impressively similar, knowing that only one diagonalization, at $\Gamma$, is needed. The rest of the information is imported from the HEG. The HEG bands, aligned to the second valence band at $\Gamma$, are the white lines for comparison. Of course, the higher bands of He have some resemblance to free electrons, but the differences are significant and overall well reproduced by the connector. The lowest valence band is due to localized $s$-states and far from the additional free electron band given by the dashed white line aligned at $\Gamma$. The connector correctly predicts the almost complete absence of dispersion for this highly localized state, which is very far from a HEG or simple metal.

By construction, the further the real system is from the model, the less well will a simple approximation to COT perform. We also want to illustrate this case, so Fig. 3 displays constrained connector results for solid hydrogenIshmaev et al. (1983), obtained in the same way as for the example of helium. As the figure shows, for this material the ab initio band structure is very far from the HEG one. The main deviations from the HEG are well described by COT, for example, the existence of a double band below 5 eV between L and A, or of the second band between $\Gamma$ and M. At the same time, supplementing the information at $\Gamma$ with the HEG yields a realistic idea of the dispersion of the bands, which would be completely flat, had we applied our approximation directly to the calculation of $G$, without using the connector. The quality of the connector band structure is instead less satisfactory between M and L, as one would expect for a path far from the reference point, i.e. $\Gamma$. To improve, one could produce the exact results in a few more ${\bf k}$-points, assuming transferability of the connector over a part of the Brillouin zone only, or one may go beyond the HEG, e.g. moving to a periodic model, which is more pertinent in the spirit of COT. This will be left for future work. In its present approximation, the connector approach may not replace a full band structure calculation for all materials, but it is certainly a very fast way to already obtain much insight.

V Connector approach to the design of functionals

The potential use of connector approximations is not limited to speeding up calculations. Indeed, when the parameters $Q^{R}$ that characterize the real system constitute a function, i.e., when $O$ is functional of $Q$, (6) together with (5) stand for an approximate functional. This functional benefits from pre-calculated ingredients that are universal, in the sense that once the model is decided, they will be the same for all real systems. The design of functionals is therefore a very promising field of applications for COT, as we will illustrate in the following.

V.1 Exchange-correlation potential: density functional from a connector approximation

One of the prototype problems for the design of functionals is the search for a simple but reliable Kohn-Sham exchange-correlation potential $v_{\rm xc}({\bf r};[n])$ beyond the LDA, which is still a subject of intense research (see e.g. Gunnarsson et al. (1977); Alonso and Girifalco (1977, 1978); Gunnarsson et al. (1979); Gunnarsson and Jones (1980); Cuevas-Saavedra et al. (2012); Langreth and Perdew (1980); Perdew and Yue (1986); Perdew et al. (1996, 1999, 2008); Sun et al. (2015); Burke (2012); Cohen et al. (2012); Becke (2014); Pribram-Jones et al. (2015); Jones (2015)). Our aim here is to demonstrate that COT has the potential to suggest new density functionals, and to explain in which way they lead to an improvement with respect to the LDA.

To apply COT, for close comparison with the LDA we use again the HEG as model system. Contrary to above in the example of the band structure, now the HEG will not be defined by its potential, but by its density, in other words, the connector will be a number giving an effective homogeneous density. The exact connector condition (4) reads:

[TABLE]

where $n^{c}_{{\bf\mathbf{r}}}$ is the unknown homogeneous density that yields the exact value of the xc potential $v_{\rm xc}(\mathbf{r},[n])$ at position $\mathbf{r}$, when used to evaluate the xc potential of the HEG $v_{\rm xc}^{h}(n^{h})$. Note that $n^{c}_{\mathbf{r}}$ is a homogeneous density, but a different $n^{c}_{\mathbf{r}}$ is used for the potential in every point $\mathbf{r}$. This freedom is crucial, because otherwise the exact condition (18), which reflects condition [A], could not be fulfilled. In this way, instead, all $v_{\rm xc}$ of real systems can be reached as long as they are negative, as the HEG xc potential $v^{h}_{\rm xc}$ spans all negative values.333It should be noted that the condition will often be fulfilled straight away, but there are exceptions. In particular, if one is interested, for example, in the dissociation of linear molecules, one must be aware of the fact that the xc potential $v_{\rm xc}({\bf r})$ of these one-dimensional systems shows positive bumps, see, e.g., Perdew (1985); Buijse et al. (1989); Baerends and Gritsenko (1997); Tempel et al. (2009); Helbig et al. (2009); Wetherell et al. (2019) whereas the exchange-correlation potential in the HEG is always negative. Therefore, if one wants to connect these systems to the HEG additional action, such as an extension of the domain, would be needed to fulfill condition [A]. This, however, goes beyond the scope of the present illustration. Concerning condition [B] there is no problem, since $v_{\rm xc}^{h}$ is a monotonous function that can be inverted straightforwardly.

To illustrate that we can also fulfill [C] and obtain efficient expressions for practical use, our scope is now to find an approximate connector, for which (5) reads

[TABLE]

A first crude approximation could be to cut the Coulomb interaction $1/|{\bf r}|$ beyond a short distance $r=r_{c}$, without renormalizing the remaining short-range part of the interaction. As $v_{\rm xc}$ entirely stems from the Coulomb interaction, dropping a substantial part of the interaction decreases the exchange-correlation potential, and even makes it vanish with the cutoff radius $r_{c}\to 0$. The connector result corresponding to the very same approximation, instead, should tend to the LDA, because now we do not apply the approximation to the potential, but only to the way the potential in a given point “sees” the environment, limiting it, for $r_{c}\to 0$, to the environment close to the point ${\bf r}$ where the potential is calculated. We can therefore expect a much more meaningful result.

V.2 Benchmarking the description of non-locality

In order to illustrate this and the following discussion, we cannot compare to the exact functional $v_{\rm xc}$, which is unknown. Instead, we take an established non-local functional as target, and we will do all approximations and calculations consistently with respect to this functional. It is important to note that the aim is not to elucidate the quality of this target functional, but to use it as benchmark for the connector approach, and in particular, how well it can capture non-local effects beyond the LDA.

We choose as target functional an expression based on a weighted density approximation (WDA) of the xc hole $n_{\rm xc}$ Gunnarsson et al. (1977); Alonso and Girifalco (1977, 1978); Gunnarsson et al. (1979), with the weight function proposed in Gunnarsson and Jones (1980). The xc energy reads

[TABLE]

with $\tilde{n}({\bf r},{\bf r^{\prime}})=[n({\bf r})+n({\bf r^{\prime}})]/2$ García-González et al. (2000) and the Coulomb interaction $1/|{\bf r}-{\bf r}^{\prime}|$ appearing at two places, as factor in the integral and in the exponent to the fifth power. The functions $C$ and $\lambda$ are set by the sum rule $\int d{\bf r^{\prime}}n_{\rm xc}({\bf r},{\bf r^{\prime}-r})=-1$, and by making $E_{\rm xc}$ exact in the HEG. The functional derivative of $E_{\rm xc}^{\rm WDA}$ with respect to the density yields the exchange-correlation potential $v_{\rm xc}^{\rm WDA}$ given in the appendix [Eq. (39)], together with the functions $C$ and $\lambda$ [Eq. (38)]. This potential is our target. To be precise, we will not investigate the density that is created by the potential, but the potential itself, for a given input density. The resulting density, which requires the use of the potential in a self-consistent cycle, would be an important topic on itself Kim et al. (2013); Wasserman et al. (2017); Vuckovic et al. (2019).

The target potential $v_{\rm xc}^{\rm WDA}$ is shown in Fig. 4, for an inhomogeneous “real system” with periodic density $n(\mathbf{r})=A\cos(\mathbf{a}\cdot\mathbf{r})+B$, reciprocal lattice vector $\mathbf{a}$ and parameters $A$ and $B$. Here ${\bf a}$ and the maximum density correspond to the case of solid argon. The potential is negative for all systems examined in our work, so it can in principle be reproduced by $v_{\rm xc}^{h}$ as required by condition [A]. The figure also shows the results of the crude approximation obtained when the Coulomb interaction $1/|{\bf r}-{\bf r}^{\prime}|$ in (20) is cut beyond a very short distance $r=r_{c}$, which limits the integral boundaries. Applied directly to $v_{\rm xc}^{\rm WDA}$, this yields, as anticipated above and as one would expect, an almost vanishing xc potential, here obtained with $r_{c}=0.1\ a_{0}$ with $a_{0}$ the Bohr radius, less than 2% of the unit cell length. The connector result based on the very same approximation is now obtained by using the HEG xc potential to simulate the real system, via

[TABLE]

which is the equivalent of (6). As the figure shows, this potential does not tend to zero, but to the LDA. The reason is that the $v_{\rm xc}$ of both the real system and the HEG tend to zero in (19), but the connector construction leads to strong error cancelling. The huge improvement of the final result with respect to the direct approximation makes the LDA a prototype illustration for the fulfillment of condition [C] and for the power of the connector approach.

Still, the LDA has significant errors especially in regions of low density, because the approximation is too crude. In principle one could remedy by increasing the range of the Coulomb interaction, going systematically towards the exact result. However, this approximation would be of limited practical interest.

An alternative and more useful approximation, which also fits in our strategy required by condition [C], is a first-order expansion around a given homogeneous density $n_{0}$ Kohn and Sham (1965); Palummo et al. (1999). We will now use this approximation in the connector scheme to derive a new functional. We will then compare the result to the LDA, and benchmark its performance again by examining how well it can reproduce non-local effects, with the WDA defined as target.

Using a first-order expansion as approximation in Eq. (4) yields:

[TABLE]

where $f_{\rm xc}(|\mathbf{r}-\mathbf{r}^{\prime}|;n_{0})=\left.\delta v_{\rm xc}(\mathbf{r},[n])/\delta n(\mathbf{r}^{\prime})\right|_{n=n_{0}}$ is the static nonlocal xc kernel of the HEG with density $n_{0}$, and $f^{h}_{\rm xc}(n_{0})=\left.dv_{\rm xc}^{h}(n_{h})/dn_{h}\right|_{n_{h}=n_{0}}$ is its limit of zero wavevector, which corresponds to the case where variations are restricted such that the density remains homogeneous, i.e., we remain within the parameter space given by the model. The zeroth order term $v_{\rm xc}^{h}(n_{0})$ on both sides cancels. Finally, inverting the equation we get the approximate connector

[TABLE]

which should be used to reproduce the quantity of interest $v_{xc}$ according to the prescription of (21). For a density of the real system that varies slowly on the scale of spatial decay of $f_{\rm xc}(|\mathbf{r}-\mathbf{r}^{\prime}|)$ the connector tends to $n^{c,\rm approx}_{{\bf r}}([n])=n({\bf r})$, the LDA, which is exact in this limit. For a very quickly varying periodic density we obtain the mean density, as one would expect, and which is instead completely missed by the LDA. The approximate connector interpolates between these two limits: it displays the degree of nearsightedness of the xc potential, which is obtained from (21). For the exact $f_{\rm xc}$ of the HEG one can use parametrized data. Corradini et al. (1998); Moroni et al. (1995) This means that, similarly to the LDA, one can in principle benefit from high-level calculations performed once and for all in the model.

Up to here, the derivation and discussion of this linear-expansion connector approximated $v_{\rm xc}$ was general. In the following, we again want to test this COT approximation by reproducing $v_{\rm xc}^{\rm WDA}$. In order to simulate a realistic situation, $f_{\rm xc}$ has to be consistent with the $v_{\rm xc}$ that one is heading for. Therefore, for this test we should not use expressions for $f_{\rm xc}$ given in the literature, but we have to calculate the functional derivative of our target xc potential, and evaluate it at the constant density $n_{0}$. The resulting $f_{\rm xc}^{\rm WDA}$ can be found in the Appendix, Eq. (41).

The linear approximation is not yet well defined, because it depends on the homogeneous density $n_{0}$ around which we expand. This density may be chosen to be a different one for every point ${\bf r}$ where we calculate the potential. To focus on the improvement with respect to the LDA, we choose for $n_{0}$ the local density $n({\bf r})$, i.e. the direct approximation is a first-order expansion around the LDA.

Let us first look at the performance of the linear approximation itself, used directly on the potential without going through the connector scheme. The comparison in Fig. 4 shows that the WDA potential is well described by the direct linear approximation in regions of high density (large $|v_{\rm xc}|$), but not where the density is low. We now move on to COT, and first examine the connector density $n^{c}_{\bf r}$ given by (23) and shown in Fig. 5 for the same periodic system as above. The connector density that yields the LDA is simply the density itself, which has its minimum in the middle of the unit cell. The exact WDA connector $n^{c,WDA}_{\bf r}=(v_{\rm xc}^{h})^{-1}(v_{\rm xc}^{WDA}({\bf r};[n]))$ spans a smaller range of values than the local density, confirming that it represents an effective density that is suitably averaged over a range around the local density. It is far from trivial, even in this simple system: where the density is very low, the exact connector shows a bump, which is a feature that would not be easy to guess. One may explain it with the fact that in regions of low density the system is more far-sighted, which, in the present system, mixes in more of the higher densities, and with the non-monotonous distance dependence of effective interactions, such as $f_{\rm xc}$. The approximate connector density captures this effect very well, with a remaining discrepancy that is much smaller than the difference to the LDA.

Finally, Fig. 4 shows the exchange-correlation potential resulting from our connector. It gives by far the best result of all approximations, illustrating the fact that the connector (23) takes into account a significant amount of non-locality. It shows that, by using the very same linear expansion through COT, strong improvement is obtained, without additional cost.

Since the connector tends to the LDA for slowly varying density, the improvement with respect to LDA is particularly important when the density is quickly varying; however, some improvement is also found in the case of more slowly varying density, i.e. smaller $\mathbf{a}$, as shown in Fig. 6. The figure also illustrates how the connector result tends to the LDA, which itself tends to the exact result, when the density approaches the homogeneous limit.

V.3 A connector adapted to each object of interest

The above results confirm that COT is a promising approach to design new simple functionals capable to improve results with respect to the LDA. It was our choice to concentrate on $v_{\rm xc}$ for this illustration, because it is the unknown ingredient of the Kohn-Sham potential. Of course, other choices for illustrations would have been possible: we could have examined the exchange-correlation energy or the xc energy density $\epsilon_{\rm xc}({\bf r})$, or we could have looked at single pieces, such as $v_{\rm x}({\bf r})$ and $v_{c}({\bf r})$, or at the correlation contribution to the kinetic energy alone. It is important to understand what such a choice implies and how the different approaches are linked within the COT framework.

Let us start with the question what happens if we separate $v_{\rm xc}({\bf r})$ into two contributions, e.g., $v_{\rm x}({\bf r})+v_{\rm c}({\bf r})$. Suppose that for the systems of interest condition [A] can be fulfilled, e.g., an exact connector exists in principle, for the sum as well as for both pieces. This would mean that

[TABLE]

in principle, each contribution is calculated with its own connector, and there is no reason to suppose the connector that reproduces exactly the exchange contribution, $n^{c}_{{\bf r},{\rm x}}=(v_{\rm x}^{h})^{-1}(v_{\rm x}({\bf r};[n]))$, to be the same as the exact connector for the correlation contribution $n^{c}_{{\bf r},{\rm c}}=(v_{\rm c}^{h})^{-1}(v_{\rm c}({\bf r};[n]))$, nor to suppose that it equals $n^{c}_{{\bf r}}$. As a consequence, even when a consistent approximation is chosen for all pieces, most often also the approximate connectors will be different.

This is a general point concerning the separation of the object of interest into two or more contributions, not limited to the separation illustrated here. Such a separation may be used to improve results: first, suppose that the single contributions are of quite different nature, for example, that they show a different degree of nearsightedness, or that they suggest naturally a perturbation expansion in a different parameter, or around a different zero order. This would allow one to find a better approximation for each contribution alone and hence also to improve the final result in many cases. Even when it is not advantageous, however, the fact that connectors are in principle different for every piece of interest has to be kept in mind as soon as one is interested, for some reason, in the individual pieces: when the connector $n^{c}$ is designed to approximate the sum, it could well be that some single pieces are not well described. This would not be a failure of the COT approximation, but of the way it is used. More generally, supposing transferability of the connector is a priori an approximation. It is the approximation that we have used for the band structure of helium. That example illustrates the deviations from the exact results, including the appearance of spurious poles.

Since in the case of $v_{\rm xc}$ one is often interested in the exchange and correlation contributions separately, we illustrate this point in Fig. 7, again using the WDA, and now determining the exchange and correlation connectors separately. Here the exchange part is defined by imposing the sum rule on the exchange hole and the correct limit to the HEG, the correlation contribution being the difference to the full WDA. The exchange potential (upper panel) is an order of magnitude larger than the correlation contribution; note the change in scale. The picture for the exchange is therefore very similar to the full xc result. The correlation contribution shows a different behaviour, with the LDA being too shallow. The correlation connector well reproduces the target $v_{c}$ for the larger densities. However, for very low densities the linear expansion around the LDA yields negative connector densities, so the connector result (curve with blue filled circles) is absent for $r$ in the range $[2.43,4.32]$. This is not due to a violation of condition [A], which is indeed fulfilled, but to the fact that the approximate result of (5) falls out of the allowed domain, such that (6) cannot be evaluated. The problem can be overcome in different ways, the simplest one being a change in the approximation. Indeed, the result of the linear approximation depends on the homogeneous density $n_{0}$ around which one expands. Fig. 7 also shows the result of the direct approximation that we obtain by expanding around a much smaller density than before, here $n_{0}=n({\bf r})/6$ instead of the local density $n({\bf r})$. Contrary to before, now the direct linear expansion of $v_{\rm c}$ has a bump at low densities, but it is exaggerated, and moreover the results at larger densities are very bad. However, now the resulting connector density is always positive, and the connector $v_{c}$ can be calculated over the whole range. With respect to $n_{0}=n({\bf r})$, there is very little change at high densities: this demonstrates again the strong error cancelling inherent in the connector approach, making, as a consequence, the result also very stable. In the region of low density, instead, we now also obtain an excellent result. It should be noted, however, that in this region the result does depend on $n_{0}$ (which is what we have used to solve the problem), and further work is needed in order to define a universal prescription for the choice of $n_{0}$ for very low densities, where $v_{\rm c}$ as function of the density is very steep, and where a bad choice may therefore lead to pathologies such as the negative density above. One promising route, in particular, is the use of exacts constraints.

Finally, Fig. 7 also shows $v_{c}$ calculated as $v_{\rm c}({\bf r})\approx v_{\rm c}^{h}(n^{c}_{{\bf r}})$, i.e., using for the calculation of the correlation part the connector determined for the full exchange-correlation potential, with the original choice $n_{0}=n({\bf r})$. Clearly, the result is much worse. This illustrates our point that there is a connector for each object of interest, which implies also that one should use separate connectors if one is interested in single contributions. This is not a problem of the approach, because it does not add significant workload; one simply has to be aware of it.

The fact that the exact connectors representing single contributions to an observable are most likely different, implies that it is meaningless to discuss the relation between these contributions on the basis of a connector derived for their sum. In particular, when using only one connector one cannot require the relation between the contributions to fulfill exact constraints. Only in some limiting cases is the question well posed. For example, one can easily define the limit of extreme nearsightedness, the same for every contribution, i.e. one can unambiguously define the LDA and verify whether it fulfills the virial theoremAverill and Painter (1981); Williams and von Barth (1983); Levy and Perdew (1985) linking the correlation contribution to the total energy and the exchange-correlation energy density and potential, which it does, indeed. Averill and Painter (1981); Williams and von Barth (1983); Levy and Perdew (1985)

Another consequence of the fact that every observable has its own connector concerns the relation between energy and potential. Suppose that we have succeeded to find the connector $n^{c}_{{\bf r},e}$ that correctly describes the energy density, so

[TABLE]

We now calculate the potential, which is the functional derivative of $E_{xc}[n]$ in (25), which yields

[TABLE]

while the connector xc potential using the same connector would read

[TABLE]

which is a different expression. The connector density that correctly represents $v_{\rm xc}^{c}({\bf r})$, instead, is obtained by inverting the HEG $v_{\rm xc}^{h}$ applied to (26) and reads

[TABLE]

With these clarifications in mind, we can resume what our illustration is meant to yield and what can instead not be inferred: choosing an object of interest (here, $v_{\rm xc}$), we have shown that COT is a convenient framework to derive a simple approximate functional that reproduces the non-local dependence on the density beyond the LDA for an inhomogeneous system, using knowledge from the HEG. When evaluated in the framework of the WDA the results are highly promising, and therefore future work will be devoted to applications of the expressions to real materials, without the weighted density approximation, and to the self-consistent calculation of the density. Other observables, such as the xc energy, can also be treated within COT, but will require a different connector that will be developed separately.

VI Comparison to other model-based approximations

COT is potentially very powerful because it relies on the use of a model that can easily be calculated and/or tabulated once and for all. Other approximation schemes exist that might be used in the same way, and it is interesting to make some comparison. We will in the following concentrate on two main routes, namely, perturbation theory, and the use of a variational principle.

VI.1 Perturbation Theory

Perturbation theory (PT) is an expansion around a zero-order starting point, which can be considered to be a model. The model result as well as its derivatives could in principle be calculated and tabulated as function of some parameters, such that the results for real systems can be reconstructed efficientlyAouina et al. (2020).

To compare PT and COT we suppose that both are based on the same model. Moreover, besides the choice of a model, COT (which is in principle exact) requires the choice of an approximation, and for a straightforward comparison this will be the same PT. In other words, here we will look more in detail at the direct use of PT as compared to its use within COT.

If $\mathcal{O}$ is the observable in the model and $\mathcal{O}^{PT}$ its finite-order PT approximation, the connector result (6) is

[TABLE]

where $\Delta\mathcal{O}$ contains the higher orders neglected by the finite order perturbation theory in the model, and $\mathcal{Q}_{x}^{c{,PT}}$ is the perturbative connector. In other words, with respect to the direct PT result, the connector adds higher orders. The approximation consists in evaluating the correction not at the true $Q^{R}$, but for an effective parameter set that lies in the model domain and is expressed as $\mathcal{Q}_{x}^{c{PT}}$. Our approximation strategy defines $\mathcal{Q}_{x}^{c{PT}}$, and the improvement with respect to straightforward PT, which completely neglects the higher orders, has been illustrated in Figs. 1, 4 and 6 above. It should be stressed again that this gain comes without additional computational cost.

Finite-order PT can be optimized by choosing an appropriate zero-order starting point and by adapting the points where the derivatives are taken. The COT result $O^{c}_{x}=\mathcal{O}(\mathcal{Q}_{x}^{c{,PT}})$ can also be interpreted as an optimized zero-order PT calculation. This means that all other orders, or their sum, should be small. For the first order, this requirement reads

[TABLE]

where the integral symbolizes integration and/or summation over parameters, and where we have allowed for the possibility that the derivative is taken at some optimized parameter set $Q_{0}$ within the model domain. Solving for the optimized zero order $Q_{x}^{c}$ by supposing, as throughout this work, that the model can be represented by one effective parameter $\mathcal{Q}$, one obtains

[TABLE]

which is right the first order connector. In other words, supplementing PT with COT can be seen, alternatively, as a way to efficiently simulate higher orders, or as a way to optimize the zero-order by making some higher orders disappear. The latter point of view indicates that COT based on PT will work well if making some low order (in our example: the first order) disappear also leads to vanishing or small higher orders. This suggests one possible way to analyze a given problem in order to optimize the COT approximation, for example, through the choice of $\mathcal{Q}_{0}$.

Finally, it should again be stressed that a perturbation expansion is only one possible realization of a connector approximation, and COT can also be used to improve completely different approximations. The cutoff on the Coulomb interaction in Fig. 4 is a striking example.

VI.2 Variational search

Designing a model system such that it represents certain aspects of a real system can be seen as an optimization problem, which allows one, in some cases, to use a variational principle. In particular, when one is directly interested in energies, this can be used straightforwardly, since the model system can be seen as a trial approximation to the real system and tuned such that the energy is minimized. For other observables this is not necessarily true: one may hope that the model configuration that minimizes the energy also well describes other observables of the system, but there is no guarantee. The connector strategy, instead, applies equally to all observables.

Still, it is interesting to compare the variational and the connector approach for the calculation of an energy. We will use the example of the 1D potential of Sec. III for a simple comparison of principle and to give numerical results.

The strategy of the variational approach is to minimize the energy of the real system by varying the parameters of the wavefunctions of the model system. For our model system, the infinite potential well of width $L^{v}$, the wavefunctions are

[TABLE]

The kinetic energy in a state with even or odd $j$ is

[TABLE]

equal to the total energy of the infinite square well with width $L^{v}$. The potential energy is calculated with the potential $V(z;\omega_{0},L)$ of the real system with flat region of width $L$, so for states with odd $j$

[TABLE]

where $L^{v}>L$, and where the upper integral boundary is due to the fact that the trial wavefunction vanishes outside the region defined by $L^{v}$. An equivalent expression holds for even $j$. We now minimize $E^{\rm kin}+E^{\rm pot}$ with respect to $L^{v}$, for given $j$, $L$ and $\omega_{0}$. This has to be done numerically. The result is shown in Fig. 8, for the lowest level $j=1$ (upper panel) and for $j=4$ (lower panel). The variational method performs very well in both cases. For $j=1$, the connector outperforms the variational approach. For higher $j$, the connector does better than the variational approach for small $\alpha$ and worse for larger $\alpha$, as shown here for the case $j=4$ (lower panel), where the crossover is at $\alpha\approx 0.02$. The crossover point decreases with increasing $j$.

Both the crossover and the fact that the variational method tends to win earlier for higher levels can be understood by the fact that we have chosen a first order expansion for the connector approximation: this approximation becomes worse for larger curvature, and eventually the error cancelling in the connector approach, though impressive, is not efficient enough to repair the bad perturbation result. Also the trend with increasing $j$ can be understood: in our example we have, for all levels, expanded around a box of width $L_{0}=L$. This is a good choice for lower levels, but higher levels experience an increasingly larger effective width, as suggested by Fig. 1. This explains why the direct first order approximation becomes worse for higher levels (note the change in scale between upper and lower panel), and explains the earlier crossover of connector and variational results.

The variational approach adapts automatically to each situation. In the connector approximation based on perturbation theory, the expansion point (here, chosen to be $L_{0}=L$) introduces some flexibility. In the present example, one can use for an efficient expansion of higher levels a larger $L_{0}$, relying on the fact that improving the direct approximation also improves the connector result. Of course, for the method to be predictive such a further optimization must stem from fundamental principles, e.g., fulfillment of exact constraints, or self-consistency.

It is important to note that the computational cost for the variational approach is clearly higher than for the perturbative connector approximation. This suggests that within a significant parameter range it may be much more efficient to use COT based on a perturbation expansion, than using the variational approach. As our simple example shows, this holds even for the calculation of energies, which are naturally variational, contrary to most other observables.

Finally, we will use our simple 1-dimensional example to illustrate another aspect of COT, linked to the choice of the model system. As pointed out since the beginning, for COT to be efficient the approximation must be equivalent for the real and the model system. Now suppose we had chosen another model system, namely, the harmonic oscillator, which corresponds to $\alpha\to\infty$. This model can be used, for example, in the variational approach, with the results shown in Fig. 8: as expected, they are very good for large $\alpha$ and deviate significantly from the exact result for small $\alpha$. Note that now the performance of the variational approach is better for higher levels, since they are less influenced by the flat bottom of the potential. The harmonic oscillator can also be used as starting point for a perturbation series, where we expand the real system around $L=0$ which yields excellent results for large $\alpha$ and becomes very bad for small $\alpha$, much sooner than the variational resultBalestri (2013): in the range of $\alpha$ chosen in Fig. 8, the results would be too negative to be seen in the figure. One would now like to improve the perturbation result by using COT. To do so, one has to find an equivalent approximation in the model, i.e. the oscillator. The only parameter here is $\omega_{0}$. However, the energy levels of the harmonic oscillator are linear in $\omega_{0}$, which means, the first order expansion in the model would already be exact, contrary to the first order expansion in the real system. Therefore, in this case COT would not benefit from error cancelling, but simply yield the same result as direct perturbation theory, so we would not be able to apply the strategy required by condition [C]. In our example, if one is interested in large $\alpha$ one would therefore have to use another approximation to benefit from COT, e.g., an approximation on the kinetic energy.

Indeed, as a final note it is important to remind again that COT can in principle be combined with any approximation, not necessarily perturbation theory: one may even consider to combine it with the variational approach. Therefore, the important question is not so much how to compare other approximations to COT, but rather how to combine them with COT in the most efficient way.

VII Conclusions

In conclusion, inspired by a historical strategy to approximate density functionals using the homogeneous electron gas, we propose an in principle exact and very general approach. Its aim is to calculate once and for all, and store, a given observable or other object in a model system with high precision. These results are then used to determine the same object in real systems, via a procedure termed connector. The connector is different for every target object, and must be approximated. We suggest a strategy for a systematic connector approach, which, for a given model, makes use of an approximation in a way that leads to strong error cancelling. Indeed, we show that, and why, a given approximation is often much more powerful when used within connector theory than when directly applied to the object of interest. Of course, the choice of the model and of the approximation still determine the quality of the results and require care and physical insight. The approach opens the way for fast computational methods, which we have illustrated with a quick estimate of energy levels and bands. It can also be used to design functionals. As example, we have derived an exchange-correlation potential for DFT, using as model system the HEG. By using the weighted density approximation as benchmark functional we have demonstrated that our connector functional is able to capture non-local effects in a very efficient way. Beyond the WDA, the connector approximation to the exact functional only requires data from the HEG that are already available in the literature as interpolated expressions. More complex models may require new calculations to be done, with results that can be either interpolated with modern machine learning algorithms or stored thanks to today’s storage capacities. Such more flexible models could potentially further improve the results Aouina et al. (2020) and are envisaged for future applications. The present work sets the framework, elucidates the fundamentals and suggests directions for practical application, with a potentially huge impact on computational materials design.

Appendix A One-dimensional potential

Here we will give some details concerning the numerically exact solution of the 1-dimensional Schrödinger equation with potential (7), and concerning the first-order expansion.

The renormalized Schrödinger equation (8) can be solved analytically in each of the regions I, II and III. The condition of continuity of the wavefunction between different regions of the potentials yields an equation for an effective wavevector,Balestri (2013)

[TABLE]

where $k$ is the standard (region II) and $h$ is a generalized (region I and III) wavevector, defined by:

[TABLE]

with $\Gamma$ the Gamma function. The two solutions of Eq. (35) are $h_{\rm O}=k\tan\frac{k}{2}$ and $h_{\rm E}=-k/\tan\frac{k}{2}$, independent of $\alpha$. Using the definition of $h$ yields two sets of discrete eigenvalues $\{\varepsilon_{\rm O}\}$ and $\{\varepsilon_{\rm E}\}$. Ordering and labelling them with a single index $j\geq 1$, one obtains a series where $\{\varepsilon_{\rm O}\}$ correspond to odd and $\{\varepsilon_{\rm E}\}$ to even values of $j$.

With the definitions for $k$ and $h$ the solution of Eq. (35) for odd $j$ yields

[TABLE]

which we solve numerically for $\varepsilon_{j}(\alpha)$, and analogously for even $j$, using the Ridder method Ridders (1979) to find the roots. Fig. 9 shows the exact results; gaps in the red lines indicate regions where the root-finder algorithm was not able to converge. These gaps were filled in Fig. 1 by interpolating. This stresses the gain that one may obtain from a more efficient approach, already in this very simple example. From $\varepsilon_{j}(\alpha)$ the energy levels are obtained as $E(j;\omega_{0},L)=\varepsilon_{j}(\alpha(\omega_{0},L))/L^{2}$.

In the limit of $\alpha\to 0$ one finds $h\to\infty$. Therefore, Eq. (36) yields $\sqrt{\varepsilon}\tan\sqrt{\frac{\varepsilon}{2}}\to\infty$. Hence the eigenvalues must fulfill $\sqrt{\frac{\varepsilon}{2}}=(2n+1)\frac{\pi}{2}$ with $n$ integer, or $\sqrt{2\varepsilon}=j\pi$, $j$ odd integer, and $\varepsilon_{j}\to\varepsilon_{j}^{\rm box}\equiv\frac{\pi^{2}j^{2}}{2}$. For even $j$ one has analogous relations. To obtain the next order, we make the ansatz $\varepsilon=\varepsilon^{box}_{j}+c\sqrt{\alpha}$ and expand $\sqrt{2\varepsilon}\tan\sqrt{\frac{\varepsilon}{2}}$. This yields $c=-\frac{2\Gamma\left(\frac{1}{4}\right)}{\Gamma\left(\frac{3}{4}\right)}\varepsilon_{j}^{box}$, and finally,

[TABLE]

which is Eq. (9). This approximate solution has of course a very small computational cost, compared to the exact one, and presents no numerical difficulties.

Appendix B Choice of poles in the constrained connector approximation

In order to choose the dominant poles within the scattered connector results, we have to group poles that most probably belong to the same band. To this aim, the following procedure is applied independently at each ${\bf k}$ (and therefore, does not necessitate any hypothesis about continuation of bands): (i) For each ${\bf K}$, we determine from the spectral function the lowest energy pole, which should in principle be the same for all ${\bf K}$. In practice, one finds a distribution. (ii) We first roughly gather together those poles that lie within ten standard deviations from the average position of the first pole. These poles appear with different intensities, whereas the true spectral function for non-interacting electrons should only show one, normalized, pole for the first band. Therefore, we sum the intensities of the approximate poles until the normalization is reached; these poles are assigned to the first band. The subsequent poles are assigned to the following band. Within the first group, the pole with maximum intensity is chosen to represent the band $j=1$. (iii) The others are discarded, and the band counter moves on, from $j=1$ to $j=2$. (iv) There may be ${\bf K}$ vectors for which the lowest pole is not visible, because it shows zero intensity. If the bands are not too dense, the next pole will lie outside of the standard deviation gathering. Also in this case, the counter moves to $j=2$, but those poles will be re-considered. (v) The same procedure is applied to band $j=2$, including as candidate poles the second lowest poles or, for the ${\bf K}$ where (iv) was applied, the poles to be re-considered. The procedure then continues in the same way. Of course, it cannot restore the exact band structure and may fail for large unit cells where bands are very dense, or for materials very far from the HEG where no dominant poles may be identified, and it could certainly be refined. However, it has the advantage of being very simple and fast.

Appendix C Weighted Density Approximation

In the following, we will give more detail concerning our target functional based on the weighted density approximation (WDA) of the xc hole $n_{\rm xc}$ introduced in Gunnarsson et al. (1977); Alonso and Girifalco (1977, 1978); Gunnarsson et al. (1979), with the weight function proposed in Gunnarsson and Jones (1980). The xc energy is given in (20). The two functions $\lambda$ and $C$ are

[TABLE]

where $\varepsilon_{\rm xc}(n)$ is the exact xc energy per particle of the HEG, whose correlation part we take from the Perdew-Zunger parameterizationPerdew and Zunger (1981).

To get the target xc potential we perform the functional derivative of the xc energy, $\frac{\delta E_{\rm xc}[n]}{\delta n(\mathbf{r})}$. It yields :

[TABLE]

with $\tilde{n}({\bf r},{\bf r^{\prime}})=[n({\bf r})+n({\bf r^{\prime}})]/2$.García-González et al. (2000) Further functional derivative yields the exchange-correlation kernel, which can be expressed as the sum of six terms, $f_{\rm xc}^{\rm WDA}(|\mathbf{r}-\mathbf{r}^{\prime}|;n_{0})=\sum_{i=1}^{6}f_{i}$. Evaluated in the HEG with density $n_{0}$ these terms read:

[TABLE]

[TABLE]

acknowledgement

Stimulating discussions with many members of the Palaiseau Theoretical Spectroscopy Group, and with Ilya Tokatly and Sergio Ciuchi are gratefully acknowledged. This research was supported by a Marie Curie FP7 Integration Grant within the 7th European Union Framework Programme, the European Research Council under the EU FP7 framework program (ERC grant No. 320971), and by the Austrian science Fund FWF under Project No. J 3855-N27.