Self-Consistent Density-Functional Embedding: A Novel Approach for Density-Functional Approximations.

In the present work, we introduce a self-consistent density-functional embedding technique, which leaves the realm of standard energy-functional approaches in density functional theory and targets directly the density-to-potential mapping that lies at its heart. Inspired by the density matrix embedding theory, we project the full system onto a set of small interacting fragments that can be solved accurately. Based on the rigorous relation of density and potential in density functional theory, we then invert the fragment densities to local potentials. Combining these results in a continuous manner provides an update for the Kohn-Sham potential of the full system, which is then used to update the projection. We benchmark our approach for molecular bond stretching in one and two dimensions and show that, in these cases, the scheme converges to accurate approximations for densities and Kohn-Sham potentials. We demonstrate that the known steps and peaks of the exact exchange-correlation potential are reproduced by our method with remarkable accuracy.

Introduction

Over the past decades, density functional theory (DFT) has become a well-established and successful method that is able to accurately describe molecular and condensed matter systems. One reason for its success can be attributed to its computational efficiency as all physical observables of interest are functionals of the ground-state density n(r).[1] The most popular technique to find the density of the system accurately is the Kohn–Sham (KS) DFT, where the density of the full interacting system is computed via an auxiliary noninteracting system.[2] All interactions and correlations of the interacting system are mimicked by the so-called exchange-correlation (xc) potential, which is usually determined as the derivative of the xc energy functional Exc[n] that is unknown and has to be approximated in practice.[2−6] A remaining challenge is to find functional approximations describing other wanted observables O[n].

Another issue with DFT is that, although significant progress in functional development over the years has been achieved, approximate DFT functionals usually still struggle to describe systems with strongly correlated electrons.[7] The dissociation limit of the H2 molecule is a good example for a simple system that is not easy to describe with commonly used approximate DFT functionals. Those functionals that are optimized to be able to mimic the dissociation of H2[8−13] do not perform equally good on other problems.[14]

There are alternative methods that are able to describe strongly correlated electrons accurately. One big group are wave-function methods, such as full configuration interaction (FCI) methods[15] and density-matrix-renormalization group (DMRG).[16] These methods, although becoming more and more efficient, still have high computational cost and thus are only able to describe relatively small systems.

A pathway to use accurate methods on a larger scale is provided by embedding theories. The general idea behind embedding consists of dividing a system into one or more fragments of interest and an environment, which is then considered only indirectly. With this partition, the need of performing an expensive calculation on the full system is circumvented. An established group of embedding theories are various density-functional embedding methods[17−22] that have been successfully applied to a large range of complex systems, such as molecule adsorption on metallic surfaces,[23] proton transfer reactions in solution,[24] and photophysical properties of natural light-harvesting complexes,[25] to name a few. They provide ways of calculating a system that is weakly bounded to an environment by representing the environment by an external field. Opposed to that, embedding methods such as dynamical-mean-field theory (DMFT),[26−28] density-matrix-embedding theory (DMET),[29−31] and density-embedding theory (DET)[32,33] consider correlations between the system and the environment more explicitly and thus are successful in describing systems with strongly correlated electrons. This is achieved by mapping the full system onto a fragment that is embedded into a, in some cases, correlated, bath. In the latter two methods, only the fragment is described accurately, while the rest of the system is described with a lower-level calculation. Here, the challenge is the connection between the high-level and low-level calculations.[29,31−36]

All mentioned embedding methods are tailored to describe the behavior of the fragments accurately. Opposed to that, we use in the present work the embedding idea to improve our large-scale description of the full system by including insights from small fragments. To this end, we introduce a feedback algorithm, which combines DMET with density inversions based on the one-to-one correspondence of density and potential in exact DFT. Our approach is closely related to the DET[32] realization of DMET with the main difference being that we partition the system into strongly overlapping fragments. Although the idea of overlapping fragments was considered before in bootstrap embedding,[34] the procedure itself is different to ours. This results in a self-consistent density-functional embedding (SDE) technique, which allows to explicitly construct approximations to the xc potential with increasing accuracy. Here, no optimized-effective potential (OEP)[37,38] procedure needs to be employed since it is not the energy that is approximated but directly the local xc potential. In our context, we are not using an explicit expression of the xc potential in terms of the density but rather employ a direct numerical construction. We use the embedding to find numerically local approximations to the density-potential mappings that give direct access to the xc potential. Once the optimal KS system is obtained, we gain information about observables from those interacting fragment wave functions, which serve as an approximation to the full interacting wave function. To put it differently, we approximate the involved potential-density maps of standard DFT as a combination of local maps. In the limit that the fragment describes the full system, we find also the exact KS potential.

The paper is organized as follows. In Section , we introduce the proposed SDE method step by step. In Section , we present the Hamiltonian for two electrons in a heteroatomic model system in one and two dimensions, which we use to benchmark our approach. The results for the energy, density, and KS potential of the introduced systems are shown in Section . Finally, a summary of the SDE method and an outlook toward more general applications are given in Section .

Theory

Density Functional Theory

In this section, we introduce key aspects of DFT and issues of standard DFT approximations that we wish to address with our approach. Based on the Hohenberg–Kohn (HK) theorem,[1] in KS DFT,[2] the ground-state density n(r) of a target interacting many-body system is obtained through a set of auxiliary one-body (KS) equations with an effective local (KS) potential vKS[n](r) (atomic units are used throughout the paper)The difference between the KS potential and the external potential vext(r) of the interacting system is the so-called Hartree-exchange-correlation (Hxc) potential vHxc[n](r) that accounts for all the interactions and kinetic correlations of the interacting system. This potential is usually obtained by approximating the corresponding Hxc energy functional EHxc[n] and then taking the functional derivative of the latter with respect to the density.

Although highly successful, there are several issues with this approach. From a formal perspective, it is shown that the exact functionals as defined by Lieb are not functionally differentiable.[39] So, to provide the main ingredient, regularizations need to be done.[40] Further, it is very hard to systematically increase the accuracy of known approximate functionals.[41] Also, even if we had an accurate approximate functional, it would usually be given in terms of KS orbitals, and a numerically demanding OEP procedure[37,38] would be needed to obtain the KS potential. Furthermore, there is the often overlooked but important issue of how to construct other observables from the KS Slater determinant as any observable that cannot be expressed directly in terms of the density needs to be approximated in terms of the latter.

Here, we avoid these issues by following a different path, which involves no explicit approximate expression for EHxc[n] or vHxc[n]. Instead, we first introduce a formal approach that employs density-potential mappings of DFT directly (see, e.g., ref (42)) and then make this approach practical by applying approximations to it. Following the HK theorem, for a given density n(, there is a interacting system with the external potential v[n(] that produces this density. Also, exactly the same density can be reproduced by the noninteracting system with the potential vS[n(]. Hence, an interacting density n((r) can be uniquely inverted to both an interacting potential v[n(](r) and a noninteracting potential vS[n(](r). The Hxc potential is then defined by the difference of those two potentialsSolving the single-particle eigenvalue equations (eq ) for vKS[vext, n(], we obtain the updated density n(. Starting with some initial density n(0), this scheme converges at the true ground-state density n that is produced by the external potential vext = v[n], and we have also found the noninteracting potential vKS[n] = vS[n] to reproduce this density.

Note that the fixed-point iteration scheme introduced above does not need any explicit expression of an energy functional. However, it is obvious that the scheme itself is not practical at all. To avoid solving the exact Schrödinger equation (SE) for one interacting system with vext, we ended up performing inversions not only to obtain the noninteracting vS[n(], which in principle, is feasible,[43−45] but also to obtain the interacting v[n(], which would involve solving the interacting SE multiple times at each step and hence increase the numerical complexity of the problem instead of decreasing it.

The method we present in this paper targets directly at approximating the fixed-point iteration scheme in a way that no inversion for v[n] is necessary. Within our approach, the connection between v[n] and n is given by a projection (that we introduce in Section ), and the exact SE is solved in smaller subsystems.

Self-Consistent Density-Functional Embedding Method

The fundamental idea of the SDE approach is to replace the mapping between the global KS potential and the corresponding density by dividing the system into a set of fragments {i} and mapping those onto a set of auxiliary interacting systems with a corresponding set of external potentials {v}, interacting wave functions {| Ψ⟩}, and densities {n}. Here, no interacting inversion is needed, and we also get an approximated mapping between the KS Slater determinant |Φ0⟩ and the ground-state wave function of the system |Ψ0⟩.

The SDE method is depicted schematically in Figure . It consists of the following parts, to each of which we assign a distinct subsection:

Figure 1

General SDE idea: properties of an interacting electronic system with an external potential vext and a ground-state wave function |Ψ0⟩ are fully determined by its electronic density n(r), which can be uniquely reproduced by a noninteracting system (KS system). The interacting system is divided into fragments. For each fragment (orange), the system is projected onto a smaller auxiliary interacting (embedded) system. The embedded system consists of the fragment, which remains unchanged by the projection, and the part of the system that includes interaction and correlation with the fragment (depicted in violet). Each of the embedded systems is then solved on a wave-function level, yielding an accurate density that then can be uniquely mapped onto an auxiliary noninteracting system with the same density. These accurate local potentials are then used to improve the global KS description of the full system. The whole process is repeated self-consistently until convergence of the global KS potential is reached.

The full system is described in terms of its ground-state density n(r) by means of KS DFT, as we have discussed in Section .

The system is divided into fragments. Our proposed partition differs significantly from partition DFT[46,47] or DMET,[31] and we will introduce our “continuous partition” in Section .

For each fragment, the full system is projected onto an embedded system, where the fragment is embedded into an effective bath. In this paper, the choice for the projector is inspired by the DMET approach, which we explain in detail in Section .

For each fragment, an accurate calculation is performed with a wave-function method. The fragment wave functions are then used to calculate accurate fragment densities. These wave functions also serve as a local approximation to the mapping between the KS Slater determinant and the ground-state wave function |Φ0[n]⟩→|Ψ0[n]⟩, from which we can directly calculate correlated observables via O[n] = ⟨Ψ[n]|Ô|Ψ[n]⟩. We explain how this calculation is performed in practice in Section .

Finally, for each fragment i, an auxiliary noninteracting system is found that reproduces the density n, and the set of obtained potentials {vS[n]} is then used to update the global KS potential. How this is done in practice is explained in Section . The SDE scheme is applied self-consistently, and the algorithm is also explained in Section .

As we divide our system into fragments in real space, we will, for the sake of convenience, consider only systems that are discretized on a real-space grid throughout the paper.

Continuous Partition

We continue by considering the problem of dividing the full problem into fragments. Generally, the fragments have to cover the full system and should be selected small enough to be calculated with required accuracy.

In embedding approaches such as subsystem DFT[22] and also in the framework of partition DFT,[46] the system is divided into non-overlapping fragments, which are weakly bounded to one another. In other words, the partition is dictated by density distribution and correlations within the system and cannot be chosen arbitrarily. Therefore, those approaches are not applicable when connections along fragments become important.

In DMET,[29−31] the system is also divided into non-overlapping fragments. The partition itself can be chosen arbitrarily as particle transfer between the fragment and the rest of the system is possible within this approach. The size of the fragments is dictated mostly by the correlation length in the system.[29] Hence, the amount of correlation, which is captured with the DMET method is constrained by the size of the fragment. Thus, by increasing the fragment size, a convergence toward the exact solution is feasible, which makes the method systematically improvable. Dividing the system into non-overlapping fragments, however, causes artificial discontinuities in local observables such as density1,[34] which sometimes also leads to convergence problems.[36] This is one reason why DMET can have convergence problems when applied to inhomogeneous systems.[36] For such systems, a simple single-shot embedding is usually performed,[31] which still provides very good results for the energies, which is, after all, the target of the DMET method.

In SDE, we employ the same type of projection as in DMET (see Section ), but since we are particularly targeting the density, we further introduce a partition that guarantees that all fragments connect smoothly to one another. Specifically, we define a continuous partition, where the system is covered by overlapping fragments, as depicted in Figure . The idea of using overlaps to remove edges of a fragment was first introduced in bootstrap embedding,[34] where densities of overlapping fragments are matched to one another through additional self-consistency loops. Here, we use a much simpler scheme. We sweep through the system by just going one grid point further for each fragment calculation, and when computing local observables such as the density, we only take into account the grid point in the middle of each fragment. Hence, our partition is constructed such that the local observables are continuous on the real-space grid. The accuracy can be improved by selecting the grid spacing appropriately. In practice, this has to be balanced with the computational cost as for any real-space implementation.

Figure 2

Visualization of the partition procedure: To obtain a continuous density, we sweep through the system by just going one site forward for each fragment calculation. Then, only the physical properties of the centering site are taken into account when considering local observables. The upper image (a) shows the partition in 1D, whereas the lower image (b) illustrates the partition in 2D. Projections P onto embedded systems as well as effective bases depicted by different kinds of crosses are explained in Section . This partition procedure can be extended to 3D in a straight-forward manner (not shown).

Projection onto the Embedded System

Having decided on how to divide the system into fragments, we now treat each fragment separately and find an effective description for the corresponding embedded system (see Figure ). We want the embedded system to be such that it describes the physics on the fragment as accurately as possible. As depicted in Figure , we have to project the full system onto an embedded system for each fragment.

Out of a manifold of possible projections,[21,26,29] we adopt here the projection used in DMET[29−31] as it provides an efficient way of including static correlations between the fragment and the rest of the system, which we call bath from now on.

The DMET method can be understood as a complete active space (CAS) calculation under the assumption that the fragment basis functions are always in the active space. What then remains to be found are the orbitals that build up the remaining part of the active space, which we here call the correlated bath. It is constructed such that it has the same dimensionality Nfrag as the basis that spans the fragment.

Since the construction of the DMET projection has already been introduced in the literature multiple times,[29−31] we leave the step-by-step instructions on how we do it in practice to the Appendix, and we give a visualization of the projection in Figure . By solving a mean-field Hamiltonian for the full system, we obtain a new smaller set of orbitals in which we then express the interacting Hamiltonian Ĥ of the full system to obtain the Hamiltonian Ĥemb for the embedded system.

Figure 3

Visualization of the decomposition of the system into the fragment and bath and the projection onto the embedding part (CAS) and environment. The dots depict the sites, which correspond to our chosen initial basis set, and the crosses depict the orbitals after projecting. To describe the physics of the fragment, only the embedding part is considered.

Fragment Calculation

For each fragment i, we obtain the embedding Hamiltonian Ĥembas described in Section and then diagonalize it to obtain the embedding wave function |Ψemb⟩ and the corresponding density n of this embedded system.

In the present work, we use exact diagonalization (ED) to solve for the ground-state wave function of the embedded system. We emphasize that also other solvers, such as DMRG,[16,48,49] coupled cluster,[50−52] selective CI approaches,[53] or Monte Carlo methods,[54,55] can be used for the fragment calculation.

The correlated embedding wave functions can then be used to calculate the energy of the full system E0 or any other correlated observable. As described in ref (30), the energy of the full system E0 can be approximated as a sum of fragment energies, which are calculated by taking a partial trace of the corresponding embedding density matrix2 ρ̂emb = |Ψemb⟩⟨Ψemb|. In the SDE approach, for each fragment i, only the middle site α is considered for obtaining properties of the full system (see Section ). The site value of an observable is calculated by tracing out all CAS orbitals of the fragment i except for the site α (see Appendix for a detailed definition of CAS). An observable of the full system is then obtained as a sum over values for every site. Hence, taking for example the energy, we adopt the formula from ref (30) towhere N denotes the number of fragments, which is equal to the number of grid points. Here, we have approximated the full wave function by a set of fragment wave functions. The correlation length that can be captured within this approximation is limited by the fragment size.

The formula above can be applied to any other observable. Thus, we circumvent the usual problem in DFT of finding explicit functional dependence O[n] between an observable of interest O and the density n by simply using the embedding wave functions instead of the density.

Before moving on to improving the KS description of the full system, we have to add an additional constrain to the fragment calculations. As in DMET or partition DFT, we have to make sure that, when patching the system back together, we retain the correct particle number in the full systemFollowing ref (31), we achieve this by adding and self-consistently optimizing a chemical potential μ to the embedding Hamiltonian of each fragmentwhere n̂ denotes the density operator on site α, and the index α runs over all fragment sites. The constant μ in eq is added only to the fragment part of the embedding Hamiltonian to achieve a correct particle distribution between the fragment and the environment. In other words, the chemical potential is a Lagrange multiplier, which assures that the constraint in eq is fulfilled.

Self-Consistency

So far, we have discussed how, starting from an initial guess for the KS potential (we usually start with vKS = vext), we project the full system onto a set of interacting embedded systems with {Hemb↔|Ψemb⟩↔nemb}. We now want to use this set of quantities to update the KS potential of the full system.

For each fragment i, the Hamiltonian contains a one-body part ĥemb and a two-body part Ŵemb

Following the KS construction, the corresponding density nemb can be reproduced by an auxiliary noninteracting system withwhere the correlations are mimicked by the Hxc potential v̂Hxc, emb, which is defined as the difference of one-body terms of the interacting and noninteracting systems. In practice, this potential is obtained either by analytical[56] or numerical inversion[43−45] or by a robust minimization routine as usually employed in DMET.[31] The specific inversion scheme that is used to compute the results presented later in this paper will be introduced in Section , together with the model Hamiltonians we use for Section . The idea of using the density rather than the one-body reduced density matrix (1RDM) to obtain the embedding potential has already been introduced by Bulik et al. within the DET approach.[32] This method, however, does not target a KS system that, in principle, reproduces exactly the density of the interacting system. A major difference between DET and our approach is, besides the connection to unique mappings in DFT, the continuous partition introduced in Section that results in an accurate approximation of the global KS potential.

To this aim, we approximate the Hxc potential of the full system vHxc on each site α by the corresponding value of v̂Hxc, emb on the same site.

The KS potential is then updated according to eq as

This yields the new KS Hamiltonian ĤKS = T̂ + V̂KS, which is then used to calculate a new set of projections P. This is done until convergence (see algorithm in Figure ). Eventually, we obtain an accurate density and KS potential from which also correlated observables can be calculated as described in eq . The SDE algorithm can be improved by increasing the fragment size, and it converges to the exact solution. Note that the choice of reproducing accurately the density of the interacting embedded system by a noninteracting one, which has also been used in DET,[32] is crucial as it is based on rigorous one-to-one relations between densities and potentials in DFT and gives us a well-defined target for the inversion. This would not be the case with any other quantity such as, for example, the 1RDM (which is used in DMET), since the 1RDM of an interacting system cannot be reproduced exactly by a noninteracting one, as for example, comprehensively discussed in ref (32).

Figure 4

Visualization of the SDE algorithm: The full system can be uniquely mapped onto a noninteracting KS system. The system is divided into overlapping fragments such that a continuous reconstruction of the full system is possible. An initial guess for the global KS system is made, from which a projection is built for each fragment. Then, for each fragment, the embedding Hamiltonian is calculated, and the corresponding ground-state wave function and density are computed. A self-consistency cycle is added to maintain the correct particle number. As soon as the correct particle number is ensured in the full system, the density of every fragment is inverted and yields an updated vHxc on each site independently. This potential is then used to update the KS system. The procedure is repeated until self-consistency. In pink, we mark those parts of the algorithm that are close to the DMET approach.

As in SDE, we adopt the projection from DMET, to illustrate the distinction between the two methods, that we mark in Figure in pink, which are parts of the algorithm that SDE shares with DMET. Both methods coincide only for fragment size Nfrag = 1 as only then there is no difference in partition (single-site fragments cannot overlap) and also between the density and 1RDM on the fragment (as there are no off-diagonal elements).

To complete the introduction of the SDE method, we now turn to its numerical cost. The cost of fragment calculations in SDE grows exponentially with the fragment size Nfrag, and the cost for the underlying calculation of the noninteracting system grows quadratically with the total number of grid points N. This has to be multiplied by the number of fragments, which is also N, and the needed self-consistency iterations η, yielding a total scaling of 42 ·  · N3 · η. This is, of course, more expensive than a usual DFT calculation (that is, N2 · η in local density approximation (LDA)), but cheaper than the exponentially growing cost of an FCI calculation.

Diatomic Molecule Model in One and Two Dimensions

In this section, we introduce the model Hamiltonians, which we use to validate our approach (see Section ), and also the inversion scheme used for all results.

The SDE approach so far is valid for all closed systems that can be represented by a time-independent Schrödinger equation. To benchmark our method and to show its efficiency, we describe the two-electron bond stretching of a heteroatomic molecule in one and two dimensions (see Figure ).

Figure 5

Visualization of the 1D H2 molecule. The real space is discretized on a grid with N sites. The two atoms are modeled through a symmetric double-well potential vext1D.

We model this system with the following Hamiltonian3 on a 1D/2D real-space grid[57]where ĉ† and ĉ are the usual creation and annihilation operators of an electron, respectively, with spin σ on site i, and n̂ = ĉ†ĉ is the corresponding density operator. In 2D, the index i becomes a double index withThe spacing Δx is determined by the box size L in direction x and the number of grid points N, and as the external potential, we employ a double-well potential vext.

The first part of the Hamiltonian takes into account the kinetic energy of the molecule by means of a next-neighbor hopping term. The second term in eq is the external potential, which mimics the ions of the molecule and depends on the considered dimension. In the one-dimensional case, the external potential on each point is given bywith . The numbers z1 and z2 determine the depth of each well. In our case, they take values between 0 and 2, and we will characterize the potential by their difference Δz = z1 – z2. In the two-dimensional case, the external potential takes the formaccounting for both the charge distributions of the ions in x and y directions.

The third term of the Hamiltonian takes into account the interaction of the electrons. We model the electronic interaction as well as the core potentials by the soft-Coulomb interaction, which avoids the singularity at zero distance. To do so, we include a softening parameter α = 1.

One reason for choosing a problem that only includes two electrons is that, for this example, we can analytically invert the density n of the interacting problem to yield the potential vS[n] of the auxiliary noninteracting system that has the same density. As the ground state of a two-electron problem is always a singlet, it is valid thatInserting this property into the one-body equations (eq ) yields[56]where v[n] is the external potential of the interacting system, which yields the same density n. The constant ε0 can be chosen arbitrarily as it only fixes the gauge. We choose it such that v̂Hxc(r) vanishes at the boundaries. The formula above is given in the real-space domain, but it can be applied to any quantum lattice system4 as there is a one-to-one correspondence between tdensity and potential for those systems.[58] The exact inversion formula can therefore be applied to every embedded system with two electrons, hence, to every embedded system resulting from our model. Note that, although the exact inversion formula can only be used for the special case of two electrons, there are different ways to expand this toward the treatment of more particles. The analytic inversion can either be replaced by numerical inversion schemes[43−45] or by the robust minimization scheme used in DMET.[31]

Results

To demonstrate the feasibility of our approach, we calculate densities, KS potentials, and total energies of model Hamiltonians introduced in Section . Although our numerical results are limited to 1D and 2D model systems, we still discuss cases that are notoriously difficult to capture for standard KS DFT.

Dissociation of the One-Dimensional H2 Molecule

Common DFT functionals such as the local density approximation (LDA[2]) or generalized gradient approximations (GGA[3,4]) fail to describe the dissociation limit of the H2 molecule. This failure is attributed to the so-called static correlation error, which is related to fractional spin states.[14] Common approximate functionals, however, violate this condition and predict wrong energies for fractional spin states, resulting in the wrong dissociation limit.

Although there are methods such as the strictly correlated electron functional,[9] functionals based on the random phase approximation (RPA)[8,59] and on GW combined with RPA,[60] or the exchange-correlation potential by Baerends et al.,[61,62] which were designed to overcome these issues, modeling the bond stretching of H2 remains a challenging test for any new functional.

In Figure , we show how the SDE method performs in this test case. We plot the ground-state energy of the Hamiltonian in eq with Δz = 0 as a function of interatomic distance calculated with FCI, one-dimensional LDA-DFT,[63] one-site DMET (DMET(1)) that is equivalent to one-site SDE (SDE(1))5, single-shot DMET with five fragment sites (DMET(5s-s))6, and five-site SDE (SDE(5)). The initial guess for the projection for both SDE and DMET is built from the one-body part of the Hamiltonian in eq . The exact (FCI) energy curve shows the following well-known behavior: when varying the distance of the two core potentials d, the curve has a minimum corresponding to a stable molecule. For smaller core distances, the energy grows due to the repulsion of the two cores. Increasing the distance d→∞ leads to the vanishing of the binding energy, resulting in two separate atoms.

Figure 6

Ground-state energy of the 1D H2 molecule, calculated with FCI (black dashed line), one-dimensional LDA (green dashed line), five-site single-shot DMET (turquoise circles), five-site SDE (blue stars), and single-site DMET/SDE (red dash-dotted line). While LDA and DMET(1)/SDE(1) fail to describe the correct long-distance behavior, both DMET(5s-s) and SDE(5) show excellent agreement with the exact result. The following set of parameters has been used (see Section ): number of real-space grid sites N = 120, box size L = 20, potential well difference Δz = 0, and softening parameter α = 1. Atomic units (a.u.) are used throughout the paper.

As discussed above, LDA does not predict the correct dissociation behavior of H2 due to the static correlation error, and the energy of the two separated atoms is overestimated. One-site embedding methods DMET(1)/SDE(1) also fail to describe this behavior correctly as static correlation cannot be captured with such small fragment sizes. They perform even worse than LDA for large distances.

In contrast, both SDE and single-shot DMET show excellent agreement with FCI for Nfrag = 5. Both curves are on top of the FCI result. DMET even results in slightly better energies for intermediate distances. This might seem surprising at first glance, but the SDE algorithm is optimized to provide good densities and potentials and, as widely discussed in the literature,[41] this does not necessarily go hand in hand with more accurate energies. The difference in energy between SDE and DMET is, however, negligible, and in the next section, we show that SDE indeed does provide excellent densities and KS potentials.

Peaks and Steps in the KS Potential

For the H2 model, the KS system needs to describe the repulsion of the two electrons. As the system does not include an actual interaction term, this repulsion needs to be mimicked by the KS potential. As has been investigated in various works,[56,61,62] we expect to see a peak that prevents the two electrons from being at the same atom. In Figure , we plot the density and KS potential obtained with SDE for fragment sizes of 1, 5, and 9 sites and compare them with the exact density and exact KS potential.

Figure 7

Density distribution n(x) and KS potential vKS(x) with SDE(5) (blue solid line), SDE(9) (orange solid line), SDE(1) (red dash-dotted line), LDA (green dashed line), and FCI (black dashed line). The exact and SDE solutions for fragments sizes larger than one agree quantitatively. The SDE KS potential in these cases shows the expected peak in the center, which mimics the electron–electron interaction. For Nfrag = 5, this peak is slightly overestimated but converges quickly to a quantitatively exact result for Nfrag = 9. The SDE(1) and LDA results, on the other hand, differ significantly from the exact solution. The peak in the KS potential is missing completely. The following set of parameters has been used: N = 120, L = 20, and d = 10.

The density from the SDE calculations for the two larger fragment sizes agrees quantitatively with the exact density. We also see a peak at position x = 0 in the KS potential for both SDE(5) and SDE(9) calculations. This peak is slightly overestimated for Nfrag = 5 but agrees quantitatively with the exact solution as the fragments get bigger (Nfrag = 9). The SDE(1)/DMET(1) results are also plotted. As already discussed in the case of the energy, both density and potential deviate strongly from the exact solution. The peak in the KS potential accounting for strong correlations in the system is missing completely, and hence, also the density distribution deviates strongly from the exact solution. The same applies to results obtained with LDA.

Further, we compare SDE densities to the ones from our real-space implementation of single-shot DMET that showed good results for ground-state energies of the model in the previous section. In Figure , we plot the deviation of the approximate densities Δn from the exact ones (FCI) for both methods for Nfrag = 5. We see that the DMET density deviates stronger from the exact solution than the SDE density. Furthermore, in DMET, we clearly see a peculiarly shaped density, especially at fragment boundaries. This behavior is caused by the fact that there is no smooth connection between the fragments. This comparison reveals the need of our type of partitioning to have accurate densities.

Figure 8

Deviation of densities Δn from FCI reference results for five-site SDE (blue solid line) and five-site single-shot DMET (turquoise solid line with circles). SDE density exhibits smooth behavior, while DMET density shows discontinuities at fragment boundaries. The following set of parameters has been used: N = 120, L = 20, and d = 10.

As the next challenge, we consider more general situations such as bond stretching of heteroatomic molecules, such as LiH, that can also be modeled by the Hamiltonian of eq by considering an asymmetric external potential. The SDE results are plotted in Figure , and also here, we observe excellent agreement with exact results for both density and potential. We observe an asymmetric density distribution, which is mimicked by a KS potential that, in addition to the peak observed in the symmetric case in Figure , has a step between the two wells. The appearance of the step and its importance in KS DFT is, to this day, a widely discussed issue in the literature.[64−67]

Figure 9

Density distribution n(x) and KS potential vKS(x) for an asymmetric external potential with SDE(5) (blue solid line), SDE(9) (orange solid line), and FCI (black dashed line). Both SDE results agree with the exact solution and show expected peak and step in the KS potential. The following set of parameters has been used: N = 120, L = 20, and d = 10.

Even though approximate functionals, for example, those based on the exact-exchange approximation, do reproduce the step in the KS potential,[68] to the best of our knowledge, so far, there does not exist any approximate energy functional that can reproduce both peaks and steps[69] at the same time. Within the SDE approach, we achieve both claims, and that is why we believe that, with SDE, we provide a new path toward accurate KS potentials even for strongly correlated systems.

Convergence Behavior

In contrast to conventional DFT approaches, the SDE method can be improved systematically simply by increasing the size of the fragments, and from our numerical evidence, we deduce that this improvement is systematic. In Figures and 11, we see the deviation of our results from the exact solution for different properties Q of the system, integrated over the whole spacewhere Δx is the grid constant.

Figure 10

Integrated deviation of the density (upper graph) and KS potential (lower graph) of the SDE calculation from the exact solution for weakly static correlated (d = 0) and strongly static correlated electrons (d = 10). In both cases, we observe a decrease in the error between the two calculations. While, in the weakly correlated case, the error estimate is higher for small fragments and decreases faster, in the strongly correlated case already, the calculations for small fragments are very good and decrease slower. Already for (Nfrag = 3), the error is of the order of Δn ≤ 10–4. Parameters for d = 0: N = 120, L = 10, Δz = 0, and α = 1; parameters for d = 10: N = 120, L = 20, Δz = 0, and α = 1.

Figure 11

Difference of the total energy between the SDE and the exact solution ΔE0 with and without rescaling with respect to the particle number. We consider two different core distances (d = 0, upper graph, and d = 10, lower graph), which correspond to weak and strong correlation between the electrons. For the weakly static correlated system, already for Nfrag = 9, the error between the two calculations is below our selected accuracy limit. For strongly static correlated electrons, d = 10, we observe that the energy estimate of the SDE calculations for Nfrag ≥ 9 is too low compared to the exact solution. The deviation in energy is very low for small fragment sizes (ΔE0 ≤ 10–5). Parameters for d = 0: N = 120, L = 10, Δz = 0, and α = 1; parameters for d = 10: N = 120, L = 20, Δz = 0, and α = 1.

In Figure , we plot the deviation of the density Δn and KS potential ΔvKS between the SDE calculation and the exact result. We consider two different core distances (d = 0 and d = 10), which correspond to weak and strong static correlation between the electrons. In both cases and for both chosen properties, we observe a monotonous decrease in ΔQ with increasing fragment size up to a quantitative agreement of the two solutions. Already for the smallest considered fragment size Nfrag = 3, the deviations are relatively small, which is of the order of the fourth digit for the density Δn ≤ 10–4 and of the order of the first digit for the KS potential ΔvKS = 10–1.

In Figure , we show the deviation of the total energy E0 of the SDE method from the exact calculation. Again, we consider one example with weakly static correlated electrons and one example with strongly static correlated electrons. For weakly correlated electrons, the difference in energy decreases, and already for an fragment size of Nfrag = 7, the deviation from the exact solution is below a chemical accuracy of 1.6 mhartree.

For strongly (static) correlated electrons, we observe that the SDE energy becomes smaller than the exact energy for a range of fragments between Nfrag = 9 and Nfrag = 20. This is because the SDE method is not variational and the estimate for the energy therefore can also be lower than the exact energy. Also, for this observable though, already for small fragments, our estimate is of order ΔE0 ≤ 10–5, which is far below the chemical accuracy.

Since we approximate the wave function of the full system by a set of fragment wave functions, the total particle number calculated with fragment wave functions is not necessarily correct. The employed optimization of the chemical potential leads to the correct number for up to a desired accuracy (). As the energy difference is of the same order of magnitude, we further rescale the energy with respect to the particle numberto see if we achieve a better convergence behavior. We indeed do as we can also see in Figure . Nonetheless, the calculated energy can still be lower than the exact energy, meaning that we still observe the nonvariational nature of our approximation.

Application to Systems in 2D

To demonstrate that the SDE method can be applied to higher-dimensional models, we here discuss the H2 molecule and a model heteroatomic molecule in two dimensions.

In Figure , we plot the density n, KS potential vKS, external potential vext, Hartree-exchange-correlation potential vHxc, and deviations from the exact solution Δn and ΔvHxc for the two-dimensional H2 model. We observe a homogeneous density distribution around the two core potentials that is consistent with the external potential. The Hartree-exchange-correlation potential, which mimics the interactions of the electrons as well the kinetic correlations in the interacting case, shows a peak in the middle of the molecule. Our observations are consistent with the exact solution of this problem.

Figure 12

H2 molecule in two dimensions. Plotted are the density n, the Hartree-exchange-correlation potential vHxc, as well as their difference from the exact reference Δn and ΔvHxc, respectively, the KS potential vKS, and the external potential vext with SDE(4 × 4). We observe a homogeneous density consistent with the external potential. vHxc shows the peak accounting for the interactions of the two electrons. We observe good agreement with the exact reference. The following set of parameters has been used: N = 40, N = 20, L = 20, L = 10, d = 10, and Δz = 0.

For a model heteroatomic molecule, we plot in Figure the same properties as for H2. The density for the heteroatomic molecule in the two-dimensional case is asymmetrically distributed between the two cores, again consistent with the external potential. In the Hartree-exchange-correlation potential, additional to the peak accounting for the interaction of the electrons, we also observe a step that accounts for the asymmetric distribution of the density.

Figure 13

Heteroatomic molecule in two dimensions. Plotted are the density n, the Hartree-exchange-correlation potential vHxc, as well as their difference from the exact reference Δn and ΔvHxc, respectively, the KS potential vKS, and the external potential vext with SDE(4 × 4). We observe an asymmetric density consistent with the external potential. vHxc again shows the peak accounting for the interactions of the two electrons. Additionally, a step accounting for the asymmetric distribution of the density can be observed. Again, we observe good agreement with the exact reference. The following set of parameters has been used: N = 40, N = 20, L = 20, L = 10, d = 10, and Δz = 0.5.

Conclusions and Outlook

We present a self-consistent density-functional embedding (SDE) approach, which is a way to apply KS DFT without any explicit functional expressions but approximating the density to potential mapping. To this end, we employ a DMET-inspired algorithm, which reconstructs the KS potential of the global system from its local fragments. To achieve this goal, we use fragments that strongly overlap with each other. Observables O are calculated through a set of fragment wave functions that avoid the need of explicit functionals O[n]. SDE yields accurate results for two-electron systems in one and two dimensions for moderate fragment sizes. Not only we can very accurately reproduce the exact potential energy surfaces of these systems but also the peaks and steps in the KS potential predicted by the exact solution. Additionally, from our numerical evidence, the SDE method appears to be systematically improvable by increasing the size of the fragment, and it converges to the exact solution.

To calculate larger fragment sizes and particle numbers with SDE, a wide range of solvers based on DMRG,[16,48,49] coupled cluster,[50−52] selective CI,[53] or quantum Monte Carlo[54,55] can be included into the algorithm. Further, to treat larger particle numbers, the analytic inversion scheme in eq has to be substituted by a numeric one, as for example, proposed in refs (43−45), or simply be replaced by robust optimization schemes, as in conventional DMET.[31] We expect to face one challenge with respect to the treatment of larger systems, and that is the storage and projection of the electron–electron interaction term, which numerically is stored in a large tensor of fourth order (and thus also grows by fourth order with respect to the system size). To treat larger systems, either we have to find an efficient way of storing the interaction tensor of the original system and then project it to the embedded system or we could employ the noninteracting bath picture from DMET,[31] which circumvents the treatment of the interaction tensor for the full system altogether.

In this work, we provide a promising group of methods that combine functional methods with embedding schemes, yielding systematically improvable results. Work to extend the method to larger systems is underway.