Incremental Embedding: A Density Matrix Embedding Scheme for Molecules
Abstract
The idea of using fragment embedding to circumvent the high computational scaling of accurate electronic structure methods while retaining high accuracy has been a longstanding goal for quantum chemists. Traditional fragment embedding methods mainly focus on systems composed of weakly correlated parts and are insufficient when division across chemical bonds is unavoidable. Recently, density matrix embedding theory (DMET) and other methods based on the Schmidt decomposition have emerged as a fresh approach to this problem. Despite their success on model systems, these methods can prove difficult for realistic systems because they rely on either a rigid, nonoverlapping partition of the system or a specification of some special sites (i.e. “edge” and “center” sites), neither of which is welldefined in general for real molecules. In this work, we present a new Schmidt decompositionbased embedding scheme called Incremental Embedding that allows the combination of arbitrary overlapping fragments without the knowledge of edge sites. This method forms a convergent hierarchy in the sense that higher accuracy can be obtained by using fragments involving more sites. The computational scaling for the first few levels is lower than that of most correlated wave function methods. We present results for several small molecules in atomcentered Gaussian basis sets and demonstrate that Incremental Embedding converges quickly with fragment size and recovers most static correlation in small basis sets even when truncated at the second lowest level.
I Introduction
The fast and accurate calculation of quantum mechanical properties of molecules and materials is one of the major unsolved problems in quantum chemistry. The computational cost of most accurate electronic structure methods rises sharply with system size, limiting their applications to small systems and/or moderatesized basis sets.Bartlett and Musiał (2007); Hachmann et al. (2007); Chattopadhyay et al. (2015); Taffet and Scholes (2017); Chien et al. (); Stemmle et al. (2018) This scaling challenge can be potentially circumvented by fragment embedding, where the system is divided into smaller fragments, and the computationally involved, highlevel theory is only required for each individual fragment. The complicated interaction between the fragment and its largesized surroundings is then approximated by the interaction with an effective bath that mimics the rest of the system. The idea of fragment embedding serves as the basis for many methods, including fragment molecular orbital theoryKitaura et al. (1999); Fedorov and Kitaura (2004); Khaliullin et al. (2006); Fedorov and Kitaura (2007) (localized molecular orbitalbased embedding), subsystem density functional theorySenatore and Subbaswamy (1986); Johnson et al. (1987); Cortona (1991); Jacob and Neugebauer (2014) (densitybased embedding), and dynamic meanfield theoryMetzner and Vollhardt (1989); Georges and Krauth (1992); Georges et al. (1996); Maier et al. (2005); Turkowski et al. (2012) (local Green’s functionbased embedding) to name a few.
A major challenge to the development of a general fragmentbased method is the treatment of chemical bonds between fragments. Recently, Schmidt decompositionKlich (2006); Peschel and Eisler (2009); Peschel (2012) has been used for embedding fragments that are strongly correlated to a bath, which occurs when embedding fragments across chemical bonds. For each fragment, the Schmidt decomposition transform the rest of the system into an entangled, effective bath which is of same dimension as the fragment. A lowdimensional embedding Hamiltonian is then constructed in the Schmidt space and solved accurately therein. In practice, a highlevel calculation such as FCI (full configuration interactionSzabo and Ostlund (1996)), DMRG (density matrix renormalization groupWhite (1992); Verstraete et al. (2008)) or CCSD (coupledcluster singles and doublesIII and Bartlett (1982); Ahlrichs and Scharf (1987)) is embedded in a lowlevel bath (usually meanfield, e.g. HartreeFockSzabo and Ostlund (1996)), to recover the electron correlation missing at the meanfield level.
In order to optimize the embedding, some matching conditions are usually imposed. So far there have been two main classes. In the first class, DMET (density matrix embedding theoryKnizia and Chan (2012, 2013)) and DET (density embedding theoryBulik et al. (2014a, b)), one uses rigid, nonoverlapping fragments and matches the oneparticle density matrix (1PDM) between the fragment and the bath. Mathematically this is achieved by adding to the lowlevel bath an effective oneparticle potential, which changes both the low and highlevel 1PDMs. This effective potential is then tuned to satisfy the matching condition. This approach has shown good performance on model systems such as the Hubbard model and atomic rings/chains, even in the strong correlation domainKnizia and Chan (2012, 2013); Wouters et al. (2016); Zheng and Chan (2016); Bulik et al. (2014a, b). As with many fragment embedding methods, however, the restriction to nonoverlapping fragments results in persistent edge effects and slow convergence with fragment sizeWouters et al. (2016); Zheng et al. (2017). In BE (Bootstrap EmbeddingWelborn et al. (2016)), one instead uses overlapping fragments and requires in the overlapping region the “edge” sites from one fragment agree on density matrix elements with the “center” sites from another fragment. As long as one can make a clear distinction between the edge sites (usually on the boundary of a fragment) and the center sites (usually the most embedded part), this scheme helps to get rid of the edge effects and leads to faster convergence as demonstrated on the Hubbard modelWelborn et al. (2016).
Though successful for model systems, these methods encounter difficulties with real molecules. On the one hand, the need for a rigid, nonoverlapping partition of the system makes DMET/DET ambiguous when high symmetries such as translational invariance are lost. As a consequence, applications to realistic systems have been so far restricted to atomic ringsKnizia and Chan (2013); Wouters et al. (2016), chainsKnizia and Chan (2013) or simple polymersBulik et al. (2014b) in small basis sets where one atom or monomer with several basis functions can be considered as a fragment. One attempt at modeling real molecules was made by Wouters et al. Wouters et al. (2016) who performed DMET calculations for the potential energy surface of a symmetric reaction, \ceC12H25F\bond…F^ ¿ C12H25F\bond…F^. However, the DMET result was less accurate than that of the fullsystem CCSD even with a fragment as large as four \ceCH2 groups. BE, on the other hand, can do overlapping fragments but requires a clear definition of edge and center sites, which is also ambiguous in systems lacking certain symmetries. This was demonstrated by Ricke et al. Ricke et al. (2017) who performed a number of BE calculations on 2D Hubbard model with fragments of different shapes and choices of center sites. According to their report, the combination that gives the best energetics is not always intuitive.
We present here a scheme, Incremental Embedding, that allows the combination of arbitrary overlapping fragments without the knowledge of edge sites. As a proof of concept, we test this new fragment embedding scheme on several molecules in atomcentered Gaussian basis sets. Numerical results suggest that this method converges quickly with fragment size at equilibrium geometry and recovers most static correlation in minimal basis sets, but the performance deteriorates in either bond dissociation limit or larger basis sets. We show that this arises due to the nature of the HF Schmidt bath and point out some possible solutions.
This article is organized as follows. In Sec. II, we give the theoretical background through briefly reviewing the Schmidt decomposition and existing Schmidtspace fragment embedding methods. In Sec. III, we formulate our theory of Incremental Embedding based on a new concept, Schmidt reduction, introduced therein. In Sec. IV, we give the computational details. Then in Sec. V, we present numerical results on several molecular systems as a proof of concept. In Sec. VI, we discuss a potential problem of using HF as the bath. Finally in Sec. VII, we conclude this work by pointing out several future directions.
Ii Background
In the following, terminologies from lattice model are used for the formal discussion. All results can be adapted to realistic systems by replacing the “site basis” with the appropriate oneparticle basis in the corresponding scenario. For example, for molecules this could be either symmetrically orthogonalized atomic orbitalsSzabo and Ostlund (1996) (SOAO) or localized molecular orbitals (LMO) given by some flavor of orbital localization methods (e.g. the FosterBoys schemeBoys (1960)).
ii.1 Schmidt Decomposition
Suppose the system consists of two parts, the fragment (which we assume to be the minority) and the environment, such that the Hilbert space for the whole system is a direct product of the two subsystems, i.e. . Any state therefore has the following tensor product decomposition
(1) 
where and are the (manybody) basis states that span the fragment and the environment, respectively. Eqn (1) can be brought to the Schmidt decomposed form
(2) 
by a singular value decomposition on the coefficient matrix , where ’s are the socalled Schmidt entangled bath states and ’s are the singular values. The benefit of the Schmidt decomposition is that the length of the expansion is limited by the number of linearly independent fragment states (those with nonzero ). Thus the decomposed form has a manageable length no matter how large the original system is, as long as the fragment is not too large. Eqn (2) is exact in the sense that if is a ground state of at some level of theory, the embedding Hamiltonian
(3) 
obtained by projecting onto the Schmidt space with operator
(4) 
shares the same ground state as at the same level of theory. For example, the exactness of HFinHF embedding has been verified explicitly in Ref. 29.
In electronic structure theory, the Schmidt decomposition is usually performed in the site basis. In general, for a fragment composed of sites can be written as
(5) 
where the summation goes over fragment sites and the same number of Schmidt bath sites. If one starts with a meanfield wave function, the complicated manybody Hamiltonian in eqn (5) is truncated at the twobody level and therefore can be solved by the aforementioned accurate quantum chemical methods such as FCI, DMRG and CCSD. For this reason, nearly all Schmidtspace fragment embeddings – including this work – use a HF bath, with only a few exceptionsTsuchimochi et al. (2015). Then in order to optimize the embedding, one often needs to impose some flavor of matching conditions. So far there have been two main classes: DMET/DET match the fragment and bath using rigid, nonoverlapping fragments, and BE matches fragment to fragment when they overlap. We will review both classes in the following two subsections.
Before that, let us put special emphasis on two potential problems underlying the HF bath. First, if some site is unentangled with all other sites in the bath wave function, it gives vanishing singular value in eqn (2) and is therefore redundant for the embedding. This basis set degradation problem is in fact not rare for realistic systems with a meanfield bath, to which we will get back in Sec. VI. Second, for a site, electron system the maximum number of sites that can be entangled with each other is upper bounded by the number of electrons or holes in the system, i.e.
(6) 
In other words, this limits the size of fragments one can use in a Schmidt decomposition. If the system in concern is halffilled (i.e. ), one can expect to approach the exact solution by resorting to fragments of larger size. If this is not the case, however, such convergence is vague. In practice, a largesized basis set is often needed in order to include as much dynamic correlation as possible. This will make the system far from halffilling and therefore deteriorate the performance of embedding with a HF bath. We will also get back to this point later in Sec. V.3.
ii.2 DMET and DET
In 2012, Knizia and ChanKnizia and Chan (2012) proposed the idea to embed a highlevel theory in a HF bath. To optimize the embedding, a oneparticle effective potential is added to the HF bath:
(7) 
where is restricted to be a single Slater determinant. The bath 1PDM is then made to match the fragment 1PDM by tuning this effective potential
(8) 
where denotes the (correlated) embedding wave function, and indices and go over fragment sites only. Note that in generating [eqn (3)], is involved only in the bath part. Therefore, gains its dependency on only through the bath. As demonstrated by Tsuchimochi et al., the matching condition in eqn (8) is not always exactly satisfiableTsuchimochi et al. (2015) and therefore is optimized in a leastsquares senseWouters et al. (2016). Nevertheless, once the matching is achieved, the total energy can be expressed as a sum of fragment energies
(9) 
where the summation goes over all (nonoverlapping) fragments ; partitions
(10) 
such that each of them only involves terms that belong to the fragment as well as half the interactions between the fragment and the bath to avoid double counting. DMET is extremely powerful for strongly correlated model systems such as the Hubbard model Knizia and Chan (2012); Zheng and Chan (2016) and atomic ringsKnizia and Chan (2013); Wouters et al. (2016) where unique partitions of the system are obvious due to the high symmetry. Two years later, Bulik et al. Bulik et al. (2014a) simplified the matching condition in eqn (8) by requiring only the diagonal of the 1PDM to be matched, i.e.
(11) 
where goes over fragment sites only. This variant is named DET and has shown performance that are comparable with the original DMETBulik et al. (2014a, b).
When applying DMET/DET to realistic systems, the absence of high symmetry makes an unambiguous nonoverlapping partition of the system difficult or even impossible. Different partitions often lead to different results and there is no apparent way to evaluate the quality of those different choices. In the most general scenario, perhaps the best one can do is to adopt the following onesite embedding scheme:
(12) 
where index goes over all sites and is simply the “site energy”. Eqn (12) is free of the ambiguity problem by construction, but the generalization to fragments of larger size is not obvious in the context of DMET. We will see in Sec. III how this scheme could be improved systematically by Incremental Embedding.
ii.3 Bootstrap Embedding
In addition to the ambiguity of fragment partition, the restriction to nonoverlapping fragments also results in slow convergence with fragment size due to the persistent edge effects. Recently, Welborn et al. Welborn et al. (2016) proposed the BE scheme in order to eliminate the edge effects in certain situations. The motivation is that when several fragments overlap, the edge sites in one fragment might be the center of another. Therefore by matching properties such as 1PDM and/or 2PDM of the former to those of the latter, one can expect improving the description of the edge sites without deteriorating the center sites. Mathematically this is formulated as a constrained optimization. Suppose fragment overlaps partially with fragment , and ’s center sites are edge sites of . Then the matching condition between and is satisfied by making the following Lagrangian stationary
(13) 
which can be transformed to an eigenvalue problem of a dressed Hamiltonian
(14) 
In other words, the Lagrange multipliers play the role of a constraint potential to satisfy the matching conditions. Similar to the direct optimization method in DFTWu and Yang (2003); Wu and Van Voorhis (2005), Ricke et al. proved that has a negative semidefinite Hessian, making its stationary point a maximum and rendering eqn (14) numerically favorableRicke et al. (2017). Results on model systems have suggested that this method indeed leads to faster convergence with fragment sizeWelborn et al. (2016).
Despite its success on model systems, the generalization of BE to a general molecular system is challenging. This is because the distinction between the edge and the center sites is often vague unless certain symmetries such as translational invariance, are present. Nevertheless, the idea of matching among overlapping fragments is of significant importance. We will see that it provides – after being combined with the technique we are going to introduce in the next section – a path towards realistic systems.
Iii Theory
iii.1 Schmidt Reduction
We first introduce a tool that enables us to encode the information from a larger embedding space to a smaller one. Suppose we start with a wave function , perform the Schmidt decomposition with an site fragment, and obtain the embedding Hamiltonian in the resulting electron, site space (assuming all sites are entangled). Then we solve for its ground state , perform a second Schmidt decomposition involving sites, and obtain a new embedding Hamiltonian in the resulting electron, site space. Overall, the process can be summarized as
(15) 
We call this process a Schmidt reduction (SR) from sites to sites, or SR for short. Due to the exact nature of the Schmidt decomposition, and share the same ground state even though the latter could be of much smaller dimension.
If one starts with a HF wave function , has a simple form involving only one and twobody interactions due to the meanfield nature of its bath. However, is correlated and so is the Schmidt bath derived from it. This renders complicated and awkward to deal with in practice [cf. eqn (5)]. The workaround here is to go to the Schmidtspace matrix representation. For instance, can be elegantly represented by a by matrix in the complete onesite Schmidt basis:
(16) 
One of the important consequences of eqn (15) is that it suggests an obvious way to combine overlapping fragments: reduce each of them onto the sites where they all overlap. It is then natural to require that the properties of all fragments agree on those common sites. This provides a powerful set of matching conditions that are not based on the discrimination of edge and center sites. In the following, we introduce one realization of this idea, Incremental Embedding.
iii.2 Incremental Embedding from Two Sites to One Site
One way to exploit the strength of Schmidt reduction is what we introduce in this work, Incremental Embedding (henceforth abbreviated as IE). The goal is to construct a onesite effective Hamiltonian for any given site such that it contains correlation at the level of site embeddings (). In this section, we focus on the lowest order, . The resulting theory is named IE from two sites to one site, or IE for short. The generalization to will be presented in Sec. III.4.
Suppose we have already solved the HF wave function for a system described by the following site Hamiltonian
(17) 
Then a Schmidt decomposition of on site gives the meanfield approximation to
(18) 
which is merely the matrix representation of in the onesite Schmidt basis of site derived from the meanfield bath. A better approximation can be obtained by the following SR
(19) 
where could be any other site. According to the exact nature of SR, contains the correlation between and at the level of twosite embedding, and thus is an improvement over . However, approximating by is problematic: different choices of in general give different ’s and it is hard to determine which choice is better than others.
These observations motivate one to consider how to appropriately accumulate contributions from multiple sites. A simple sum of ’s over several sites leads to severe doublecounting problem, since each matrix alone is a representation of the full Hamiltonian . To this end, we propose to first divide the Hamiltonian into pieces, obtain the SR matrix representation for each piece using the twosite fragment that is most relevant to it, and then assemble them to construct an approximate that is free of doublecounting.
Specifically, consider the following partition of Hamiltonian
(20) 
where we require each include terms that (i) belong to the fragment or (ii) are interactions between and other sites. If some term is shared by fragments, we simply attach a factor of to average it over all relevant fragments. In other words, eqn (20) distributes evenly to all twosite fragments anchored by . With this partition in hand, we define the incremental Hamiltonian from two sites to one site as
(21) 
where is the matrix representation of in the onesite Schmidt basis of site obtained by the following SR
(22) 
and is its meanfield counterpart
(23) 
The physical meaning of is clear: it accumulates the correlations between site and all other sites that are missing at the meanfield level. Adding this correction to , we have a better approximation to ,
(24) 
One can expect to be of the quality of twosite embeddings, since each piece of the Hamiltonian is improved by the embedding calculation involving the most relevant twosite fragment. This is also schematically illustrated in FIG. 1.
Once we obtain the effective Hamiltonians for all sites, we can readily determine their ground states (as the eigenvectors of the lowest eigenvalue) and compute the site densities (vide infra). In general, they do not add up to the correct number of electrons, because each is generated from multiple fragments of different chemical potentials. This violation in the conservation of particle number can be fixed by introducing a global chemical potential , which is determined by solving
(25) 
where the dependent site densities are obtained from solving (as opposed to the bare Hamiltonians); in the Schmidt basis shown in eqn (16). We note that it is also possible to apply a set of sitespecific chemical potentials to tune the population for each site. This is useful when IE is performed in a nonselfconsistent manner (see Sec. III.6).
iii.3 Expectation Values in Ie
In last section, we discussed how a onesite effective Hamiltonian can be constructed from successively improving the meanfield description by incorporating the correlation with every other site. Once done, the ground state for each site is approximated by a fourdimensional vector, , in the onesite Schmidt basis [eqn (16)]. The ground state expectation value of any given operator can then be evaluated by summing contributions from each site. Let us take the total energy as an example, whose corresponding operator is the Hamiltonian . First we obtain a partition of over all sites
(26) 
which is a special case of eqn (10) with each fragment involving only one site. For each site , we further partition by distributing it evenly to all relevant twosite fragments in a way similar to eqn (20)
(27) 
With this partition, a matrix representation of in the Schmidt basis of IE can be constructed by the same procedure described in eqns (21 – 24)
(28) 
which gives the site energy for at the level of IE
(29) 
Finally, the total energy for IE is simply a sum of all site energies
(30) 
Eqn (30) can be viewed as an unambiguous generalization of the onesite DMET energy in eqn (12) to twosite fragments. The generalization to an arbitrary number of fragment sites is presented in the following section.
Note that the process above becomes extremely simple, if the operator involves only one site, say . In that scenario, one can bypass the partition and summation steps [eqn (26 – 28)], and obtain the matrix representation of that operator using any fragments involving . Examples of this type include the diagonal elements of 1PDM (site densities, ) and 2PDM (pairdensities, ).
iii.4 Generalization to Fragments of Arbitrary Size
In this section, we generalize IE to fragments of arbitrary size, in a way that is similar to the method of increments commonly used in local correlation methodsPulay (1983); Pulay and Saebø (1986); Saebø and Pulay (1987); Hampel and Werner (1996). The end results can be summarized in the following recursive formula for
(31) 
where is an appropriate constant that ensures the series of equations terminate appropriately (vide infra); is the incremental Hamiltonian from sites to one site,
(32) 
where terms like and are matrix representations of the sum of relevant pieces of as defined in eqn (20). For example, can be obtained by the following SR
(33) 
The physical meaning of is also straightforward: it is the correction from site embedding calculations that are not included in any site embedding calculations.
The constant coefficients arise due to the difference between traditional incremental methods (such as the aforementioned local correlation methods) and IE. In local correlation methods, a hierarchy similar to eqn (31), but with for all , can be derived, which terminates when all sites are involved (i.e. ). In that situation, the local correlation method is exact in the sense that it is equivalent to applying the same correlation method to all sites. In IE, on the other hand, the highest level one can go with eqn (31) is limited by the maximum number of entangled sites, [eqn (6)]. If one requires that the highest level of IE be

exact when the system is at halffilling (i.e. ), and

an average of all site embedding calculations otherwise,
the following expressions for can be derived
(34) 
Note that when , eqn (34) gives and hence formally reduces to traditional local correlation methods. Moreover, eqn (34) gives for , which is also consistent with IE [eqn (24)].
In order to generalize the energy evaluation scheme, we need to generalize eqn (28) to multiple sites. To that end, a recursive formula similar to eqn (31) can be derived for ,
(35) 
where is the same set of coefficients given by eqn (34). With this in hand, the site energy
(36) 
and the total energy
(37) 
for IE can be straightforwardly evaluated, where is the lowest eigenvector of (with a proper chemical potential). Eqn (37) can be viewed as an unambiguous generalization of the onesite DMET energy in eqn (12) to fragments composed of an arbitrary number of sites.
iii.5 Matching Conditions
So far we have not touched one of the most powerful ingredients in embedding calculations – the matching condition. From the discussion above, constructing requires embedding calculations for all site fragments involving site . Without any constraints, these overlapping fragments in general will not agree with one another on their common site, with the only exception where the exact bath is used as opposed to the meanfield approximation. This observation indicates that one can optimize these embedding calculations by forcing the match to happen.
Suppose we have obtained for all sites. Solving them under an appropriate chemical potential, we can compute and as described in Sec. III.3. These values are our current best estimation of the diagonal elements of the exact 1PDM and 2PDM. Naturally, we can require the site densities and pairdensities of all fragments match them. Mathematically, the problem of constraining certain density matrix elements to given values can be formulated as a constrained optimization, and has already been addressed in BE [eqn (13) and (14)].Welborn et al. (2016); Ricke et al. (2017) Here, we adapt the method to IE. Suppose we want to constrain both site densities and pairdensities for a twosite fragment . We can achieve this by introducing the following Lagrangian
(38) 
where is the embedding Hamiltonian; and are the target values; is short for . Making stationary leads to the following eigenvalue equation
(39) 
where is dressed by a constraint potential
(40) 
Eqns (39) and (40) enable us to apply the desired constraints to fragment calculations readily in a ground state formalism. Note that the constraint potential only exists in obtaining , and should not be included in other steps of SR.
iii.6 Density Optimization
Once the matching conditions are imposed to each fragment in IE, one can construct a new set of and recompute the site densities and pairdensities. In general, these values are of better quality compared to the old estimation, due to the embedding being optimized by the matching conditions. In return, these new densities can be used to constrain further embedding calculations which will generate of even better quality. This process can be repeated until selfconsistency is reached, making the theory a closed loop (FIG. 1).
The discussion above immediately suggests an algorithm to optimize the densities in IE selfconsistently. Here, we state it for IE for the sake of simplicity, and the generalization to larger fragments should be straightforward.

Solve the HF wave function for the whole system; obtain the meanfield Hamiltonians for all sites.

Obtain some guess densities and (e.g. HF).

Perform embedding calculations for all twosite fragments; in each calculation, constrain the site densities and pairdensities of the fragment sites to match and , respectively.

For each site , Schmidt reduce all involved fragments to , and compute according to eqn (24).

Diagonalize under an appropriate global chemical potential ; obtain the ground states , and recompute and .

If the new densities do not match and , go back to step 2 with the new guess densities; otherwise, the density optimization is converged, and the ground state expectation values of desired operators can be computed using .
In addition to the selfconsistent version, we note that IE can also be formulated as a nonselfconsistent theory. In terms of the algorithm, the main difference lies in step 5: instead of solving for a global chemical potential, one determines a set of sitespecific chemical potentials , such that for each site the population matches the guess densities. In other words, the densities are not optimized and IE is used in a oneshot style, similar to the method.Hybertsen and Louie (1986); Godby et al. (1988); van Schilfgaarde et al. (2006) This approximation would be useful when (i) the system is large and hence full selfconsistency is expensive, and (ii) the quality of the guess densities is reasonable. We will examine the performance of this approximation in Sec. V.
iii.7 Computational Scaling
We end this section by briefly discussing the computational scaling. For IE, the total work is dominated by the embedding calculations of site fragments. Symmetries can effectively reduce this number by a constant factor, as fragments related by symmetry operations will give same onesite Hamiltonian and hence need to be evaluated only once. For each fragment, there are two potential ratelimiting steps: (i) the basis transformation of integrals (including the partitioned Hamiltonian ) from the site basis to the Schmidt basis, and (ii) the highlevel calculation (in this paper, FCI). For small fragments (), the former dominates with a scaling, which will exceed the aforementioned scaling if . Fortunately, the basis transformation needs to be performed only once, due to the meanfield bath not being optimized in IE. This feature makes the basis transformation step usually negligible, especially in the selfconsistent version. Under these conditions, IE has a scaling of , which is lower than most accurate quantum chemical methods [] if the incremental expansion [eqn (31)] can be truncated at a small .
Iv Computational Details
In the following computational work, we will examine the performance of IE using several small molecules. The symmetrically orthogonalized atomic orbitalsSzabo and Ostlund (1996) (SOAOs) are used as site basis for the radial expansion of the hydrogen ring, while localized molecular orbitals (LMOs) of the FosterBoys styleBoys (1960) are used in all other cases. The necessary atomic integrals are generated by Psi4Parrish et al. (2017). The FosterBoys localization is performed in QChemShao et al. (2015). Molecular geometries are also optimized in QChem at the B3LYPBecke (1993)/ccpVTZJr. (1989) level and can be found in Supporting Information. All embedding calculations, including IE, BE and onesite DMET, are performed using the electronic structure program, frankenstein,Ye (2017) developed by one of the authors. Spinrestricted HF (RHF) and FCI are used as bath and highlevel solvers, respectively. In selfconsistent IE, both site densities and pairdensities are selfconsistently determined based on an RHF guess. In the nonselfconsistent version, only the site densities are constrained to the RHF values due to the bad quality of the meanfield pairdensities. For onesite DMET, we abandon the selfconsistency and also constrain the site densities to the RHF values. For all systems tested in this work, the exact solutions are accessible and obtained by the Block DMRG codeChan and HeadGordon (2002); Chan (2004); Ghosh et al. (2008); Sharma and Chan (2012); OlivaresAmaya et al. (2015).
V Results
v.1 Radial Expansion of the Hydrogen Ring Model
We select the minimalbasis hydrogen ring model as our first example for several reasons. First it can be viewed as the simplest generalization of the Hubbard model towards real molecules, covering both weakly correlated domain (near equilibrium geometry) and strongly correlated domain (dissociation limit). Second, it is an “easy” case for Schmidtspace embedding with a HF bath according to our discussion in Sec. II.1 because the system is at halffilling. Last but not least, the cyclic symmetry makes all sites (which are symmetrically orthogonalized orbitals in this case) equivalent. This not only renders BE applicable for comparison, but also tremendously reduces the computational work for IE so that the trend of convergence with fragment size can be examined thoroughly.
In FIG. 2 the energy errors per atom in the radial expansion of STO3G \ceH10 are plotted for IE and BE, respectively. The nonselfconsistent version of IE is used since the site densities are completely determined by the cyclic symmetry. For BE, fragments involving two and three adjacent sites are used, but only in the latter is there the distinction of center and edge sites. In that case, the pairdensities of edge sites are made to match that of the center site. Due to the halffilled configuration, both methods become exact when using fragments composed of five sites. Overall, both methods are very accurate even at the site level. The error is consistently small ( kcal/mol per atom) for all geometries tested here. Near the equilibrium position ( ), systematic improvements are observed for both methods when fragments of larger size are used. In the dissociation limit, correct asymptotic behavior is recovered in all cases even though the HF bath is spinrestricted. At intermediate geometries ( ), however, the convergence with fragment size is not monotonic for IE: the level suffers from severe overcorrelation and is worse than ; this overcorrelation is only ameliorated by corrections from the level. Quite the contrary, BE continues to reduce the error by using a larger fragment and shows better accuracy compared to IE with the same fragment size. Nevertheless, the performance loss of IE compared to BE is somewhat expected in this specific example, since the model is Hubbardlike and therefore optimal for the latter.
v.2 Single Bond Breaking
Now we consider two real molecules, \ceCH4 and \ceC2H6 in their minimal basis (STO3G), to which BE is no longer applicable. Both molecules show merely a small deviation from the ideal halffilled configuration. Therefore, we can still expect good performance from the embedding calculations. Specifically we are interested in the energetics of the following two singlebond breaking processes:
(41) 
First we consider the PESs obtained by nonselfconsistent IE as shown in FIG. 3, along with nonselfconsistent, onesite DMET results for comparison. In both methods, the site densities are constrained to the RHF values, whose quality is high near equilibrium geometry but deteriorates quickly as the bond is stretched (see Figure S1 in Supporting Information). If the error is mainly densitydriven, we should expect good accuracy at equilibrium geometries as well as increasing error along the dissociation processes. This is indeed the case for onesite DMET (green), as can be seen from the growing gap between the embedding solution and the exact one in FIG. 3. IE (red) shows similar trend when approaching the dissociation limit, but recovers only a limited amount of correlation energies at equilibrium position. On the contrary, IE (blue) predicts equilibrium energies of very high accuracy (almost overlap the exact solution) and also improves the asymptotic behaviors significantly. This indicates that IE could effectively mitigate the sensitivity to the quality of the underlying approximate densities by using fragments of larger size.
The performance of nonselfconsistent IE – especially the undercorrelation at equilibrium geometry at the level – is an indication that the error of this method is not purely densitydriven. To confirm this inference, we repeat the IE calculations above but selfconsistently determine the site densities and pairdensities. The results are presented in FIG. 4. By comparing it to FIG. 3, one can clearly see that imposing selfconsistency significantly improves the results at equilibrium geometry at the level, and keeps the high accuracy of IE at the same time. In the dissociation limit, however, imposing selfconsistency has opposite effects: the PESs are shifted upwards slightly at both levels (though has a much smaller amplitude) compared to the nonselfconsistent results, making the undercorrelation problem even more severe therein. A scrutiny on the comparison of the site densities and pairdensities obtained by all these methods (Figure S1 in Supporting Information) shows that the change from FIG. 3 to FIG. 4 is not densitydriven, as the selfconsistent densities and pairdensities are consistently worse than the nonselfconsistent counterparts. Nevertheless, the results of IE seem to be stable, especially near equilibrium geometries. In those cases, one can safely abandon the selfconsistency condition without losing much accuracy.
v.3 Correlation Energies at Equilibrium Geometry
As a final example, we investigate the effect of larger basis sets. As mentioned in Sec. II.1, any deviation from halffilling deteriorates the performance of Schmidtspace embeddings using HF bath wave functions. In practice, however, large basis sets are often essential to recover the dynamic part of electron correlation. It is thus of significant importance to examine how IE behaves in large basis sets. In FIG. 5, we present in terms of bar plot the error of total energies for four molecules at equilibrium geometries predicted by nonselfconsistent IE, along with onesite DMET for comparison. Three basis sets of increasing size are used: STO3G, 321G and 6311G. For the minimal basis, all molecules are close to being halffilled. IE shows a consistent improvement with fragment size in all cases except for \ceC2H4, and the errors of the energies are within several kcal/mol’s. Onesite DMET, on the other hand, gives good results for \ceCH4 and \ceC2H6 (as already seen in FIG. 3) but undercorrelates badly for the other two unsaturated molecules. Since the densities are the same for all three cases, these results again confirm our conjecture that the errors of Schmidtspace embedding methods are not densitydriven: compared to onesite DMET, IE successfully captures the entanglement among different sites, especially when fragments of larger size are used.
When the basis size is increased, the deviation from halffilling renders the undercorrelation problem of onesite DMET even worse, as can be seen from the increasing heights of the green bars in FIG. 5b and c. In terms of absolute values, one can see clearly that this is because the correlation energies recovered by onesite DMET remain unchanged or even drop slightly as the basis size gets larger, while the exact correlation energies always go up (Figure S2 in Supporting Information). The same phenomenon appears in IE, but the trend there is more complicated. At first, switching from minimal to doublezeta basis leads to overestimation of correlation energies for all molecules. Moreover, the convergence with fragment size gets reversed: predicts more negative numbers than does, making it worse by going to a higher level of theory. Moving further to the triplezeta basis, however, sets IE back on track: both and energies are of high accuracy. This occurs as a consequence of error cancellation: the correlation energies recovered by IE do not increase when going from 321G to 6311G, which happens to cancel the overcorrelation errors in the doublezeta basis accidentally (Figure S2 in Supporting Information). This phenomenon, observed in both onesite DMET and IE, is in fact related to the basis unentanglement problem of the RHF bath. In next section, we will discuss this problem more thoroughly using a specific example.
Vi Discussion
In Sec. II.1 we briefly mentioned that the unentanglement in bath wave function (in this work, RHF) leads to linear dependency in the Schmidt decomposition [eqn (2)] and effectively degrades the basis. Here we present a specific example that embodies this problem. We repeat our calculations in Sec. V.1 using a homemade basis set STO3G, obtained by adding one orbital to each STO3G hydrogen (assuming the atomic ring lies in the plane). The PESs obtained in the new basis are presented in FIG. 6 along with the original results for comparison. The first thing to notice is that the RHF solution remains unchanged. Population analysis suggests that the set of orbitals are not occupied at all. This is expected by our construction of the new basis since all of the oneelectron atomic integrals between the new orbitals and the original orbitals vanish by symmetry. The only nonvanishing parts are the twoelectron integrals involving an even number of orbitals such as . Unfortunately, these nontrivial interactions are not captured by the meanfield wave function due to it being a oneelectron theory. As a result, the Schmidt space in the new basis is exactly the same as in the original one, which explains the concurrence of two IE curves in FIG. 6. For the exact solution, however, those interactions do contribute to the total energy, and one can see a lower energy in the new basis.
This specifically designed example helps us to understand the trend shown in FIG. 5. As mentioned in Sec. V.3, in terms of absolute values, the correlation energies recovered by both onesite DMET and IE cease to increase once the basis set reaches a certain size (see also Figure S2 in Supporting Information). This observation indicates that the basis unentanglement problem of the meanfield bath could be common in realistic systems, especially when largesized bases are used.
Vii Conclusion
In conclusion, we have introduced Incremental Embedding, a new Schmidtspace fragment embedding scheme that allows the combination of arbitrary overlapping fragments without the knowledge of edge sites. Underlying this new method is one of the key concept introduced in this work, Schmidt reduction, which allows information to be encoded from a large embedding space to a smaller one. Based on this technique, IE constructs onesite effective Hamiltonians for all sites by hierarchically incorporating corrections from embeddings involving twosite fragments, threesite fragments, and etc to the meanfield approximation. The potential double counting problem is avoided by an elaborate application of the method of increments. This method can be viewed as an unambiguous manysite generalization of onesite DMET. It can be made either selfconsistent or nonselfconsistent. The computational scaling is for the lowest few levels in the hierarchy, which are much lower than most correlated wave function theories.
Numerical simulations on small molecules in atomcentered Gaussian bases suggest that the convergence with fragment size is quick in small bases; most of the electron correlation is recovered for all molecules tested, even when truncated at the level. Imposing selfconsistency in site densities and pairdensities improves the performance of IE considerably near equilibrium geometry, through an approach that is not densitydriven. For larger bases, both IE and onesite DMET recover only a fraction of the correlation energy, which can be attributed to the more general basis unentanglement problem of the RHF bath. In summary, this work marks the first attempt of applying Schmidtspace embedding methods to realistic molecular systems using overlapping fragments.
In the future, IE can be extended in a number of directions. First of all, up to this point we restrict ourselves to FCI for the embedding Hamiltonian, which is computationally expensive and can only be applied to fragments of limited size. One can, of course, pursue other highlevel solvers such as DMRG and CCSD, which have better computational scaling and can therefore be extended to fragments of larger size. Second, as for the site basis, we restrict ourselves in this work to SOAOs or LMOs by the FosterBoys scheme for the sake of simplicity, but other choices do exist. In analogy to local correlation methods such as local MP2Pulay (1983); Pulay and Saebø (1986); Saebø and Pulay (1987) and local CCSDHampel and Werner (1996), perhaps the most straightforward way is to explore the possibilities of using other LMOs such as those given by the PipekMezey schemePipek and Mezey (1989); Boughton and Pulay (1993) and the EdmistonRuedenberg schemeEdmiston and Ruedenberg (1963). Last but not least, the conflict between the halffilled embedding space and the maximum number of entangled sites in a HF bath calls for a better bath wave function. In this regard, the HartreeFockBogoliubovBogoliubov et al. (1958); Bogoliubov (1959) (HFB) wave function might be a good candidate because it is (i) always halffilled in the quasiparticle space, and (ii) still a meanfield theory and therefore retains the simplicity of embedding Hamiltonians in a meanfield bath.
Acknowledgements.
HY thanks Dr. Tianyu Zhu for the discussion on the method of increments. This work was funded by a grant from the NSF (Grant No. CHE1464804). TV is a David and Lucille Packard Foundation Fellow.References
 Bartlett and Musiał (2007) R. J. Bartlett and M. Musiał, Rev. Mod. Phys. 79, 291 (2007), URL https://link.aps.org/doi/10.1103/RevModPhys.79.291.
 Hachmann et al. (2007) J. Hachmann, J. J. Dorando, M. Avilés, and G. K.L. Chan, J. Chem. Phys. 127, 134309 (2007), URL https://doi.org/10.1063/1.2768362.
 Chattopadhyay et al. (2015) S. Chattopadhyay, R. K. Chaudhuri, and U. S. Mahapatra, J. Comput. Chem. 36, 907 (2015), ISSN 1096987X, URL http://dx.doi.org/10.1002/jcc.23873.
 Taffet and Scholes (2017) E. J. Taffet and G. D. Scholes, Chem. Phys. (2017), ISSN 03010104, URL http://www.sciencedirect.com/science/article/pii/S0301010417307826.
 Chien et al. (0) A. D. Chien, A. A. Holmes, M. Otten, C. J. Umrigar, S. Sharma, and P. M. Zimmerman, J. Phys. Chem. A 0, null (0), URL https://doi.org/10.1021/acs.jpca.8b01554.
 Stemmle et al. (2018) C. Stemmle, B. Paulus, and O. Legeza, Phys. Rev. A 97, 022505 (2018), URL https://link.aps.org/doi/10.1103/PhysRevA.97.022505.
 Kitaura et al. (1999) K. Kitaura, E. Ikeo, T. Asada, T. Nakano, and M. Uebayasi, Chem. Phys. Lett. 313, 701 (1999).
 Fedorov and Kitaura (2004) D. G. Fedorov and K. Kitaura, J. Chem. Phys. 120, 6832 (2004).
 Khaliullin et al. (2006) R. Z. Khaliullin, M. HeadGordon, and A. T. Bell, J. Chem. Phys. 124, 204105 (2006).
 Fedorov and Kitaura (2007) D. G. Fedorov and K. Kitaura, J. Phys. Chem. A 111, 6904 (2007).
 Senatore and Subbaswamy (1986) G. Senatore and K. R. Subbaswamy, Phys. Rev. B 34, 5754 (1986).
 Johnson et al. (1987) M. D. Johnson, K. R. Subbaswamy, and G. Senatore, Phys. Rev. B 36, 9202 (1987).
 Cortona (1991) P. Cortona, Phys. Rev. B 44, 8454 (1991).
 Jacob and Neugebauer (2014) C. R. Jacob and J. Neugebauer, WIREs. Comput. Mol. Sci. 4, 325 (2014), ISSN 17590884.
 Metzner and Vollhardt (1989) W. Metzner and D. Vollhardt, Phys. Rev. Lett. 62, 324 (1989).
 Georges and Krauth (1992) A. Georges and W. Krauth, Phys. Rev. Lett. 69, 1240 (1992).
 Georges et al. (1996) A. Georges, G. Kotliar, W. Krauth, and M. J. Rozenberg, Rev. Mod. Phys. 68, 13 (1996).
 Maier et al. (2005) T. Maier, M. Jarrell, T. Pruschke, and M. H. Hettler, Rev. Mod. Phys. 77, 1027 (2005).
 Turkowski et al. (2012) V. Turkowski, A. Kabir, N. Nayyar, and T. S. Rahman, J. Chem. Phys. 136, 114108 (2012).
 Klich (2006) I. Klich, J. Phys. A: Math. Gen. 39, L85 (2006).
 Peschel and Eisler (2009) I. Peschel and V. Eisler, J. Phys. A: Math. Theor. 42, 504003 (2009).
 Peschel (2012) I. Peschel, Braz. J. Phys. 42, 267 (2012), ISSN 16784448.
 Szabo and Ostlund (1996) A. Szabo and N. S. Ostlund, Modern Quantum Chemistry: Introduction to Advanced Electronic Structure Theory (Dover Publications Inc., Mineola, New York, 1996).
 White (1992) S. R. White, Phys. Rev. Lett. 69, 2863 (1992).
 Verstraete et al. (2008) F. Verstraete, V. Murg, and J. I. Cirac, Adv. Phys. 57, 143 (2008).
 III and Bartlett (1982) G. D. P. III and R. J. Bartlett, J. Chem. Phys. 76, 1910 (1982).
 Ahlrichs and Scharf (1987) R. Ahlrichs and P. Scharf, The Coupled Pair Approximation (John Wiley & Sons, Inc., Hoboken, NJ, USA, 1987), vol. 67, ISBN 9780470142936.
 Knizia and Chan (2012) G. Knizia and G. K.L. Chan, Phys. Rev. Lett. 109, 186404 (2012).
 Knizia and Chan (2013) G. Knizia and G. K.L. Chan, J. Chem. Theory Comput. 9, 1428 (2013).
 Bulik et al. (2014a) I. W. Bulik, G. E. Scuseria, and J. Dukelsky, Phys. Rev. B 89, 035140 (2014a).
 Bulik et al. (2014b) I. W. Bulik, W. Chen, and G. E. Scuseria, J. Chem. Phys. 141, 054113 (2014b).
 Wouters et al. (2016) S. Wouters, C. A. JiménezHoyos, Q. Sun, and G. K.L. Chan, J. Chem. Theory Comput. 12, 2706 (2016), pMID: 27159268.
 Zheng and Chan (2016) B.X. Zheng and G. K.L. Chan, Phys. Rev. B 93, 035126 (2016).
 Zheng et al. (2017) B.X. Zheng, J. S. Kretchmer, H. Shi, S. Zhang, and G. K.L. Chan, Phys. Rev. B 95, 045103 (2017).
 Welborn et al. (2016) M. Welborn, T. Tsuchimochi, and T. Van Voorhis, J. Chem. Phys. 145, 074102 (2016).
 Ricke et al. (2017) N. Ricke, M. Welborn, H.Z. Ye, and T. Van Voorhis, Mol. Phys. 115, 2242 (2017).
 Boys (1960) S. F. Boys, Rev. Mod. Phys. 32, 296 (1960).
 Tsuchimochi et al. (2015) T. Tsuchimochi, M. Welborn, and T. Van Voorhis, J. Chem. Phys. 143, 024107 (2015).
 Wu and Yang (2003) Q. Wu and W. Yang, J. Chem. Phys. 118, 2498 (2003).
 Wu and Van Voorhis (2005) Q. Wu and T. Van Voorhis, Phys. Rev. A 72, 024502 (2005).
 Pulay (1983) P. Pulay, Chem. Phys. Lett. 100, 151 (1983), ISSN 00092614.
 Pulay and Saebø (1986) P. Pulay and S. Saebø, Theor. Chim. Acta 69, 357 (1986), ISSN 14322234.
 Saebø and Pulay (1987) S. Saebø and P. Pulay, J. Chem. Phys. 86, 914 (1987).
 Hampel and Werner (1996) C. Hampel and H. Werner, J. Chem. Phys. 104, 6286 (1996).
 Hybertsen and Louie (1986) M. S. Hybertsen and S. G. Louie, Phys. Rev. B 34, 5390 (1986), URL https://link.aps.org/doi/10.1103/PhysRevB.34.5390.
 Godby et al. (1988) R. W. Godby, M. Schlüter, and L. J. Sham, Phys. Rev. B 37, 10159 (1988), URL https://link.aps.org/doi/10.1103/PhysRevB.37.10159.
 van Schilfgaarde et al. (2006) M. van Schilfgaarde, T. Kotani, and S. Faleev, Phys. Rev. Lett. 96, 226402 (2006), URL https://link.aps.org/doi/10.1103/PhysRevLett.96.226402.
 Parrish et al. (2017) R. M. Parrish, L. A. Burns, D. G. A. Smith, A. C. Simmonett, A. E. DePrince, E. G. Hohenstein, U. Bozkaya, A. Y. Sokolov, R. Di Remigio, R. M. Richard, et al., J. Chem. Theory Comput. 13, 3185 (2017), pMID: 28489372.
 Shao et al. (2015) Y. Shao, Z. Gan, E. Epifanovsky, A. T. Gilbert, M. Wormit, J. Kussmann, A. W. Lange, A. Behn, J. Deng, X. Feng, et al., Mol. Phys. 113, 184 (2015).
 Becke (1993) A. D. Becke, J. Chem. Phys. 98, 1372 (1993).
 Jr. (1989) T. H. D. Jr., J. Chem. Phys. 90, 1007 (1989).
 Ye (2017) H.Z. Ye, Frankenstein embedding: A density matrix embedding theory for molecules, https://github.com/hongzhouye/frankenstein (2017).
 Chan and HeadGordon (2002) G. K.L. Chan and M. HeadGordon, J. Chem. Phys. 116, 4462 (2002).
 Chan (2004) G. K.L. Chan, J. Chem. Phys. 120, 3172 (2004).
 Ghosh et al. (2008) D. Ghosh, J. Hachmann, T. Yanai, and G. K.L. Chan, J. Chem. Phys. 128, 144117 (2008).
 Sharma and Chan (2012) S. Sharma and G. K.L. Chan, J. Chem. Phys. 136, 124121 (2012).
 OlivaresAmaya et al. (2015) R. OlivaresAmaya, W. Hu, N. Nakatani, S. Sharma, J. Yang, and G. K.L. Chan, J. Chem. Phys. 142, 034102 (2015).
 Pipek and Mezey (1989) J. Pipek and P. G. Mezey, J. Chem. Phys. 90, 4916 (1989).
 Boughton and Pulay (1993) J. W. Boughton and P. Pulay, J. Comput. Chem. 14, 736 (1993), ISSN 1096987X.
 Edmiston and Ruedenberg (1963) C. Edmiston and K. Ruedenberg, Rev. Mod. Phys. 35, 457 (1963).
 Bogoliubov et al. (1958) N. N. Bogoliubov, V. V. Tolmachov, and D. V. Å irkov, Fortschr. Phys. 6, 605 (1958), ISSN 15213979.
 Bogoliubov (1959) N. N. Bogoliubov, Sov. Phys. Usp. 67, 236 (1959).