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 E xc [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 H 2 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 H 2 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 dynamicalmean-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 2, we introduce the proposed SDE method step by step. In Section 3, 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 4. Finally, a summary of the SDE method and an outlook toward more general applications are given in Section 5.

THEORY
2.1. 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 v KS [n](r) (atomic units are used throughout the paper) The difference between the KS potential and the external potential v ext (r) of the interacting system is the so-called Hartree-exchange-correlation (Hxc) potential v Hxc [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 E Hxc [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 E Hxc [n] or v Hxc [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 (i) , there is a interacting system with the external potential v[n (i) ] that produces this density. Also, exactly the same density can be reproduced by the noninteracting system with the potential v S [n (i) ]. Hence, an interacting density n (i) (r) can be uniquely inverted to both an interacting potential v[n (i) ](r) and a noninteracting potential v S [n (i) ](r). The Hxc potential is then defined by the difference of those two potentials Solving the single-particle eigenvalue equations (eq 1) for v KS [v ext , n (i) ], we obtain the updated density n (i + 1) . Starting with some initial density n (0) , this scheme converges at the true ground-state density n that is produced by the external potential v ext = v[n], and we have also found the noninteracting potential v KS [n] = v S [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 Schrodinger equation (SE) for one interacting system with v ext , we ended up performing inversions not only to obtain the noninteracting v S [n (i) ], which in principle, is feasible, 43−45 but also to obtain the interacting v[n (i) ], 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 2.4), and the exact SE is solved in smaller subsystems.
2.2. 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 i }, interacting wave functions {| Ψ⟩ i }, and densities {n i }. 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 ⟩.

Article
The SDE method is depicted schematically in Figure 1. It consists of the following parts, to each of which we assign a distinct subsection: 1. 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 2.1. 2. 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 2.3. 3. 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 2. We explain how this calculation is performed in practice in Section 2.5. 5. Finally, for each fragment i, an auxiliary noninteracting system is found that reproduces the density n i , and the set of obtained potentials {v S [n i ]} is then used to update the global KS potential. How this is done in practice is explained in Section 2.6. The SDE scheme is applied selfconsistently, and the algorithm is also explained in Section 2.6. 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.
2.3. 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 nonoverlapping 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 nonoverlapping fragments, however, causes artificial discontinuities in local observables such as density 1 , 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 2.4), 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 2. 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 selfconsistency 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 1. General SDE idea: properties of an interacting electronic system with an external potential v ext 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.

Journal of Chemical Theory and Computation
Article 2.4. 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 1). We want the embedded system to be such that it describes the physics on the fragment as accurately as possible. As depicted in Figure 1, 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 N frag 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  3. 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 Ĥe mb for the embedded system.
2.5. Fragment Calculation. For each fragment i, we obtain the embedding Hamiltonian Ĥe mb i as described in Section 2.4 and then diagonalize it to obtain the embedding wave function |Ψ emb ⟩ i and the corresponding density n i 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 E 0 or any other correlated observable. As described in ref 30, the energy of the full system E 0 can be approximated as a sum of fragment energies, which are calculated by taking a partial trace of the corresponding embedding density matrix 2 ρê mb i = |Ψ emb ⟩ i ⟨Ψ emb | i . In the SDE approach, for each fragment i, only the middle site α i is considered for obtaining properties of the full system (see Section 2.3). The site value of an observable is calculated by tracing out all CAS orbitals of the fragment i except for the site α i (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 to where 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 system Following ref 31, we achieve this by adding and self-consistently optimizing a chemical potential μ to the embedding Hamiltonian of each fragment where nα denotes the density operator on site α, and the index α runs over all fragment sites. The constant μ in eq 5 is added only to the fragment part of the embedding Hamiltonian to achieve a correct particle distribution between the fragment and the  . 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.

Journal of Chemical Theory and Computation
Article environment. In other words, the chemical potential is a Lagrange multiplier, which assures that the constraint in eq 4 is fulfilled.
2.6. Self-Consistency. So far, we have discussed how, starting from an initial guess for the KS potential (we usually start with v KS = v ext ), we project the full system onto a set of interacting embedded systems with {H emb i ↔|Ψ emb ⟩ i ↔n emb i }. 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 ĥe mb i and a two-body part Ŵe mb Following the KS construction, the corresponding density n emb i can be reproduced by an auxiliary noninteracting system witĥ=̂+̂[ where the correlations are mimicked by the Hxc potential v̂H xc, emb i , 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 3, together with the model Hamiltonians we use for Section 4. 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 2.3 that results in an accurate approximation of the global KS potential.
To this aim, we approximate the Hxc potential of the full system v Hxc on each site α i by the corresponding value of vĤ xc, emb i on the same site.
The KS potential is then updated according to eq 2 as This yields the new KS Hamiltonian ĤK S = T̂+ V̂K S , which is then used to calculate a new set of projections P i . This is done until convergence (see algorithm in Figure 4). Eventually, we obtain an accurate density and KS potential from which also correlated observables can be calculated as described in eq 3. 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-toone 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.
As in SDE, we adopt the projection from DMET, to illustrate the distinction between the two methods, that we mark in Figure  4 in pink, which are parts of the algorithm that SDE shares with DMET. Both methods coincide only for fragment size N frag = 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 N frag , 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 4 2 · N frag · N 3 · η. This is, of course, more expensive than a usual DFT calculation (that is, N 2 · η 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 4), 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 Schrodinger equation. To benchmark our method and to show its efficiency, we describe . 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 v Hxc on each site independently. This potential is then used to update the KS system. The procedure is repeated until selfconsistency. In pink, we mark those parts of the algorithm that are close to the DMET approach.

Journal of Chemical Theory and Computation
Article the two-electron bond stretching of a heteroatomic molecule in one and two dimensions (see Figure 5). We model this system with the following Hamiltonian 3 on a 1D/2D real-space grid 57 where cî , σ † and cî , σ are the usual creation and annihilation operators of an electron, respectively, with spin σ on site i, and nî , σ = cî , σ † cî , σ is the corresponding density operator. In 2D, the index i becomes a double index with The 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 v ext . 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 10 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 by . The numbers z 1 and z 2 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 = z 1 − z 2 . In the two-dimensional case, the external potential takes the form ,ext 1D 2 1 2 2 x y (13) accounting 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 v S [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 that φ = | | n r r ( ) 2 ( ) 0 2 (14) Inserting this property into the one-body equations (eq 1) 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Ĥ xc (r) vanishes at the boundaries. The formula above is given in the real-space domain, but it can be applied to any quantum lattice system 4 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 3. 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. 4.1. Dissociation of the One-Dimensional H 2 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 H 2 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 H 2 remains a challenging test for any new functional.
In Figure 6, we show how the SDE method performs in this test case. We plot the ground-state energy of the Hamiltonian in eq 10 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 10. The exact (FCI) energy curve shows the following well-known behavior:

Journal of Chemical Theory and Computation
Article 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.
As discussed above, LDA does not predict the correct dissociation behavior of H 2 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 N frag = 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 H 2 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 7, 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.
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 N frag = 5 but agrees quantitatively with the exact solution as the fragments get bigger (N frag = 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 realspace implementation of single-shot DMET that showed good results for ground-state energies of the model in the previous section. In Figure 8, we plot the deviation of the approximate densities Δn from the exact ones (FCI) for both methods for N frag = 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.
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 10 by considering an asymmetric external potential. The SDE results are plotted in Figure 9, and also here, we observe excellent  . 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 N frag = 5, this peak is slightly overestimated but converges quickly to a quantitatively exact result for N frag = 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.

Journal of Chemical Theory and Computation
Article 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 7, 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 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.
4.3. 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 10 and 11, we see the deviation of our results from the exact solution for different properties Q of the system, integrated over the whole space where Δx is the grid constant. In Figure 10, we plot the deviation of the density Δn and KS potential Δv KS 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 N frag = 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 Δv KS = 10 −1 .
In Figure 10, we show the deviation of the total energy E 0 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 N frag = 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 N frag = 9 and N frag = 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   Figure 11. Difference of the total energy between the SDE and the exact solution ΔE 0 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 N frag = 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 N frag ≥ 9 is too low compared to the exact solution. The deviation in energy is very low for small fragment sizes (ΔE 0 ≤ 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.

Journal of Chemical Theory and Computation
Article observable though, already for small fragments, our estimate is of order ΔE 0 ≤ 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 ( |⟨̂⟩ −̂| < − 10 5 ). As the energy difference is of the same order of magnitude, we further rescale the energy with respect to the particle number 17) to see if we achieve a better convergence behavior. We indeed do as we can also see in Figure 11. 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 H 2 molecule and a model heteroatomic molecule in two dimensions.
In Figure 12, we plot the density n, KS potential v KS , external potential v ext , Hartree-exchange-correlation potential v Hxc , and deviations from the exact solution Δn and Δv Hxc for the twodimensional H 2 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.
For a model heteroatomic molecule, we plot in Figure 13 the same properties as for H 2 . 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.

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 15 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 Figure 12. H 2 molecule in two dimensions. Plotted are the density n, the Hartree-exchange-correlation potential v Hxc , as well as their difference from the exact reference Δn and Δv Hxc , respectively, the KS potential v KS , and the external potential v ext with SDE(4 × 4). We observe a homogeneous density consistent with the external potential. v Hxc 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 x = 40, N y = 20, L x = 20, L y = 10, d = 10, and Δz = 0. Figure 13. Heteroatomic molecule in two dimensions. Plotted are the density n, the Hartree-exchange-correlation potential v Hxc , as well as their difference from the exact reference Δn and Δv Hxc , respectively, the KS potential v KS , and the external potential v ext with SDE(4 × 4). We observe an asymmetric density consistent with the external potential. v Hxc 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 x = 40, N y = 20, L x = 20, L y = 10, d = 10, and Δz = 0.5.

Journal of Chemical Theory and Computation
Article systematically improvable results. Work to extend the method to larger systems is underway.

■ APPENDIX The Construction of the Projection
Here, we give step-by-step instructions on how the projection in DMET is constructed and how it is modified in SDE to account for different particle numbers in the system.
As discussed in Section 2.4, the projection is nothing but a single-particle basis transformation optimized to describe the physics of the fragment. The new basis is found as follows: 1. The Hamiltonian of the full system with M electrons is approximated by a noninteracting single-particle Hamiltonian hm f 7 with corresponding single-particle eigenvalue equation hm f φ j (r) = ε j φ j (r). From this, we calculate the spin-summed 1RDM in the grid basis. It is a N × N matrix that reads The fact that only the lowest M/2 eigenvectors contribute is a direct consequence of the fact that this 1RDM is built from a noninteracting wave function. In the case of an interacting one, all N eigenvectors φ j (r) would contribute with some occupation number λ i , which lies between 0 and 2. 2. Having set up the 1RDM matrix γ, we separate it in different submatrix blocks, namely, one that would correspond to those grid points that belong purely to the fragment, two blocks that contain one grid point on the fragment and one on the bath, and one block that has only bath grid points. ∏ ∂ μ i k j j j j j j j j j j j j j j j j j j j j j j j j j j j j j j j j j j j j y { z z z z z z z z z z z z z z z z z z z z z z z z z z z z z z z z z z z z i k j j j j j j j j y { z z z z z z z z Here, without loss of generality, we have assumed that the fragment contains the first N frag grid points of the system. 3. We then diagonalize the submatrix γ bath . Its eigenvalues λj will be all between 0 and 2 containing up to N frag eigenvalues with 0<λj <2. 70 4. From the eigenvectors φj , bath (r) with 0 <λj<2, we build the correlated bath orbitals of our CAS in real-space basis as 5. Having obtained in this way the correlated bath orbitals, we also use a set of orbitals to describe the fragment. The fragment orbitals will be as many as the size of the fragment, and each of them will have a coefficient 1 at a specific fragment point and 0 elsewhere. The embedding Hamiltonian Ĥe mb is constructed by projecting the full Hamiltonian Ĥin the subspace that is spanned by the set of the aforementioned orbitals. In other words, Ĥe mb is obtained by a basis transformation of the original Hamiltonian.
The number of correlated bath orbitals in the CAS is equal to N frag as long as 2N frag <M<2(N − N frag ) holds; 70 otherwise, the number of correlated bath orbitals is smaller. As DMET was constructed for Hubbard-type lattice systems, for which the condition above mostly holds, in DMET, the orbital construction that we just described is used without modifications.
In SDE, we now modify the orbital construction of DMET to get N frag correlated bath orbitals regardless of the particle number M. For the low particle numbers that are considered in this paper, we achieve this by artificially including correlations into the 1RDM of the full system by including higher-energy single-particle orbitals. To do so, we adjust the formula of eq 18 to  (20) with some small value η and then continue with the orbital construction from 2. The actual value of η is not of great importance as it is only used to include higher-lying orbitals into the 1RDM, and the same CAS would be obtained for different values of η. In our implementation, η = 0.01 is chosen.
Note that eq 20 is valid only for M<2N frag . For large particle numbers M>2(N − N frag ), the procedure can be adapted in a straight-forward manner due to particle−hole symmetry. 4 Of course, the formula has to be adapted to the Hamiltonian of this system. 5 DMET(1) is the only version of DMET that we could apply to the model systems studied here self-consistently. It is also the only case in which DMET and SDE results coincide (for details, see Section 2.6 ). 6 Note that results of the single-shot approach in DMET are independent of the type of matching as no update of the meanfield embedding potential is performed. Hence, the results could also have been labeled DET(5s-s) if the same initial embedding potential would have been used. 7 In DMET, the choice of hm f is not fixed, and usually, the Fock operator is used. In SDE, hm f is the single-particle operator of the KS system. ■ REFERENCES