Incremental Embedding: A Density Matrix Embedding Scheme for Molecules

Incremental Embedding: A Density Matrix Embedding Scheme for Molecules

Hong-Zhou Ye Department of Chemistry, Massachusetts Institute of Technology, Cambridge, MA 02139    Matthew Welborn Department of Chemistry, Massachusetts Institute of Technology, Cambridge, MA 02139    Nathan D. Ricke Department of Chemistry, Massachusetts Institute of Technology, Cambridge, MA 02139    Troy Van Voorhis Department of Chemistry, Massachusetts Institute of Technology, Cambridge, MA 02139
July 23, 2019

The idea of using fragment embedding to circumvent the high computational scaling of accurate electronic structure methods while retaining high accuracy has been a long-standing 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, non-overlapping partition of the system or a specification of some special sites (i.e. “edge” and “center” sites), neither of which is well-defined in general for real molecules. In this work, we present a new Schmidt decomposition-based 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 atom-centered 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.

thanks: Current address: Department of Chemistry, California Institute of Technology, Pasadena, CA 91125

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 moderate-sized 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, high-level theory is only required for each individual fragment. The complicated interaction between the fragment and its large-sized 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 orbital-based embedding), subsystem density functional theorySenatore and Subbaswamy (1986); Johnson et al. (1987); Cortona (1991); Jacob and Neugebauer (2014) (density-based embedding), and dynamic mean-field theoryMetzner and Vollhardt (1989); Georges and Krauth (1992); Georges et al. (1996); Maier et al. (2005); Turkowski et al. (2012) (local Green’s function-based embedding) to name a few.

A major challenge to the development of a general fragment-based 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 low-dimensional embedding Hamiltonian is then constructed in the Schmidt space and solved accurately therein. In practice, a high-level calculation such as FCI (full configuration interactionSzabo and Ostlund (1996)), DMRG (density matrix renormalization groupWhite (1992); Verstraete et al. (2008)) or CCSD (coupled-cluster singles and doublesIII and Bartlett (1982); Ahlrichs and Scharf (1987)) is embedded in a low-level bath (usually mean-field, e.g. Hartree-FockSzabo and Ostlund (1996)), to recover the electron correlation missing at the mean-field 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, non-overlapping fragments and matches the one-particle density matrix (1PDM) between the fragment and the bath. Mathematically this is achieved by adding to the low-level bath an effective one-particle potential, which changes both the low- and high-level 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 non-overlapping 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, non-overlapping 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 full-system 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 atom-centered 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 Schmidt-space 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 one-particle 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 Foster-Boys 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


where and are the (many-body) basis states that span the fragment and the environment, respectively. Eqn (1) can be brought to the Schmidt decomposed form


by a singular value decomposition on the coefficient matrix , where ’s are the so-called 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


obtained by projecting onto the Schmidt space with operator


shares the same ground state as at the same level of theory. For example, the exactness of HF-in-HF 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


where the summation goes over fragment sites and the same number of Schmidt bath sites. If one starts with a mean-field wave function, the complicated many-body Hamiltonian in eqn (5) is truncated at the two-body level and therefore can be solved by the aforementioned accurate quantum chemical methods such as FCI, DMRG and CCSD. For this reason, nearly all Schmidt-space 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, non-overlapping 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 mean-field 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.


In other words, this limits the size of fragments one can use in a Schmidt decomposition. If the system in concern is half-filled (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 large-sized basis set is often needed in order to include as much dynamic correlation as possible. This will make the system far from half-filling 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 high-level theory in a HF bath. To optimize the embedding, a one-particle effective potential is added to the HF bath:


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


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 least-squares senseWouters et al. (2016). Nevertheless, once the matching is achieved, the total energy can be expressed as a sum of fragment energies


where the summation goes over all (non-overlapping) fragments ; partitions


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 alBulik et al. (2014a) simplified the matching condition in eqn (8) by requiring only the diagonal of the 1PDM to be matched, i.e.


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 non-overlapping 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 one-site embedding scheme:


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 non-overlapping 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


which can be transformed to an eigenvalue problem of a dressed Hamiltonian


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 semi-definite 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


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 two-body interactions due to the mean-field 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 Schmidt-space matrix representation. For instance, can be elegantly represented by a -by- matrix in the complete one-site Schmidt basis:


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 one-site 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


Then a Schmidt decomposition of on site gives the mean-field approximation to


which is merely the matrix representation of in the one-site Schmidt basis of site derived from the mean-field bath. A better approximation can be obtained by the following SR


where could be any other site. According to the exact nature of SR, contains the correlation between and at the level of two-site 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 double-counting 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 two-site fragment that is most relevant to it, and then assemble them to construct an approximate that is free of double-counting.

Specifically, consider the following partition of Hamiltonian


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 two-site fragments anchored by . With this partition in hand, we define the incremental Hamiltonian from two sites to one site as


where is the matrix representation of in the one-site Schmidt basis of site obtained by the following SR


and is its mean-field counterpart


The physical meaning of is clear: it accumulates the correlations between site and all other sites that are missing at the mean-field level. Adding this correction to , we have a better approximation to ,


One can expect to be of the quality of two-site embeddings, since each piece of the Hamiltonian is improved by the embedding calculation involving the most relevant two-site fragment. This is also schematically illustrated in FIG. 1.

Figure 1: Schematic illustration of IE in a four-site lattice model. The mean-field approximation (brown) is improved by incremental Hamiltonians from three two-site embedding calculations: (red), (blue), and (yellow). The site density and pair-density derived from (green) can in return be used to constrain the embedding calculations.

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


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 site-specific chemical potentials to tune the population for each site. This is useful when IE is performed in a non-self-consistent manner (see Sec. III.6).

iii.3 Expectation Values in Ie

In last section, we discussed how a one-site effective Hamiltonian can be constructed from successively improving the mean-field description by incorporating the correlation with every other site. Once done, the ground state for each site is approximated by a four-dimensional vector, , in the one-site 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


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 two-site fragments in a way similar to eqn (20)


With this partition, a matrix representation of in the Schmidt basis of IE can be constructed by the same procedure described in eqns (2124)


which gives the site energy for at the level of IE


Finally, the total energy for IE is simply a sum of all site energies


Eqn (30) can be viewed as an unambiguous generalization of the one-site DMET energy in eqn (12) to two-site 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 (2628)], 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 (pair-densities, ).

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


where is an appropriate constant that ensures the series of equations terminate appropriately (vide infra); is the incremental Hamiltonian from sites to one site,


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


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

  1. exact when the system is at half-filling (i.e. ), and

  2. an average of all -site embedding calculations otherwise,

the following expressions for can be derived


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 ,


where is the same set of coefficients given by eqn (34). With this in hand, the site energy


and the total energy


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 one-site 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 mean-field 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 pair-densities 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 pair-densities for a two-site fragment . We can achieve this by introducing the following Lagrangian


where is the embedding Hamiltonian; and are the target values; is short for . Making stationary leads to the following eigenvalue equation


where is dressed by a constraint potential


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 pair-densities. 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 self-consistency is reached, making the theory a closed loop (FIG. 1).

The discussion above immediately suggests an algorithm to optimize the densities in IE self-consistently. Here, we state it for IE for the sake of simplicity, and the generalization to larger fragments should be straightforward.

  1. Solve the HF wave function for the whole system; obtain the mean-field Hamiltonians for all sites.

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

  3. Perform embedding calculations for all two-site fragments; in each calculation, constrain the site densities and pair-densities of the fragment sites to match and , respectively.

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

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

  6. 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 self-consistent version, we note that IE can also be formulated as a non-self-consistent 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 site-specific 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 one-shot 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 self-consistency 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 one-site Hamiltonian and hence need to be evaluated only once. For each fragment, there are two potential rate-limiting steps: (i) the basis transformation of integrals (including the partitioned Hamiltonian ) from the site basis to the Schmidt basis, and (ii) the high-level 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 mean-field bath not being optimized in IE. This feature makes the basis transformation step usually negligible, especially in the self-consistent 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 Foster-Boys styleBoys (1960) are used in all other cases. The necessary atomic integrals are generated by Psi4Parrish et al. (2017). The Foster-Boys localization is performed in Q-ChemShao et al. (2015). Molecular geometries are also optimized in Q-Chem at the B3LYPBecke (1993)/cc-pVTZJr. (1989) level and can be found in Supporting Information. All embedding calculations, including IE, BE and one-site DMET, are performed using the electronic structure program, frankenstein,Ye (2017) developed by one of the authors. Spin-restricted HF (RHF) and FCI are used as bath and high-level solvers, respectively. In self-consistent IE, both site densities and pair-densities are self-consistently determined based on an RHF guess. In the non-self-consistent version, only the site densities are constrained to the RHF values due to the bad quality of the mean-field pair-densities. For one-site DMET, we abandon the self-consistency 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 Head-Gordon (2002); Chan (2004); Ghosh et al. (2008); Sharma and Chan (2012); Olivares-Amaya et al. (2015).

V Results

v.1 Radial Expansion of the Hydrogen Ring Model

We select the minimal-basis 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 Schmidt-space embedding with a HF bath according to our discussion in Sec. II.1 because the system is at half-filling. 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.

Figure 2: Total energy error per atom (in kcal/mol) of the radial expansion of STO-3G \ceH10.

In FIG. 2 the energy errors per atom in the radial expansion of STO-3G \ceH10 are plotted for IE and BE, respectively. The non-self-consistent 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 pair-densities of edge sites are made to match that of the center site. Due to the half-filled 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 spin-restricted. At intermediate geometries ( ), however, the convergence with fragment size is not monotonic for IE: the level suffers from severe over-correlation and is worse than ; this over-correlation 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 Hubbard-like and therefore optimal for the latter.

v.2 Single Bond Breaking

Figure 3: Potential energy surfaces of (a) stretching one \ceC-H bond of a methane molecule and (b) symmetrically dissociating an ethane molecule into two methyl radicals, predicted by the non-self-consistent IE and one-site DMET. The STO-3G basis set is used in both cases.

Now we consider two real molecules, \ceCH4 and \ceC2H6 in their minimal basis (STO-3G), to which BE is no longer applicable. Both molecules show merely a small deviation from the ideal half-filled configuration. Therefore, we can still expect good performance from the embedding calculations. Specifically we are interested in the energetics of the following two single-bond breaking processes:


First we consider the PESs obtained by non-self-consistent IE as shown in FIG. 3, along with non-self-consistent, one-site 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 density-driven, we should expect good accuracy at equilibrium geometries as well as increasing error along the dissociation processes. This is indeed the case for one-site 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.

Figure 4: The same calculation as shown in FIG. 3, with site densities and pair-densities self-consistently determined by IE.

The performance of non-self-consistent IE – especially the under-correlation at equilibrium geometry at the level – is an indication that the error of this method is not purely density-driven. To confirm this inference, we repeat the IE calculations above but self-consistently determine the site densities and pair-densities. The results are presented in FIG. 4. By comparing it to FIG. 3, one can clearly see that imposing self-consistency 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 self-consistency has opposite effects: the PESs are shifted upwards slightly at both levels (though has a much smaller amplitude) compared to the non-self-consistent results, making the under-correlation problem even more severe therein. A scrutiny on the comparison of the site densities and pair-densities obtained by all these methods (Figure S1 in Supporting Information) shows that the change from FIG. 3 to FIG. 4 is not density-driven, as the self-consistent densities and pair-densities are consistently worse than the non-self-consistent counterparts. Nevertheless, the results of IE seem to be stable, especially near equilibrium geometries. In those cases, one can safely abandon the self-consistency condition without losing much accuracy.

v.3 Correlation Energies at Equilibrium Geometry

Figure 5: Equilibrium geometry total energy errors (in kcal/mol) obtained by the non-self-consistent IE for several molecules in different basis sets: (a) STO-3G, (b) 3-21G and (c) 6-311G. One-site DMET results are also included for comparison.

As a final example, we investigate the effect of larger basis sets. As mentioned in Sec. II.1, any deviation from half-filling deteriorates the performance of Schmidt-space 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 non-self-consistent IE, along with one-site DMET for comparison. Three basis sets of increasing size are used: STO-3G, 3-21G and 6-311G. For the minimal basis, all molecules are close to being half-filled. 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. One-site DMET, on the other hand, gives good results for \ceCH4 and \ceC2H6 (as already seen in FIG. 3) but under-correlates 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 Schmidt-space embedding methods are not density-driven: compared to one-site 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 half-filling renders the under-correlation problem of one-site 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 one-site 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 double-zeta 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 triple-zeta 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 3-21G to 6-311G, which happens to cancel the over-correlation errors in the double-zeta basis accidentally (Figure S2 in Supporting Information). This phenomenon, observed in both one-site 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

Figure 6: A re-plot of FIG. 2 for the total energy of \ceH_10 near equilibrium geometry. Data labelled with an asterisk are computed in the new basis STO-3G.

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 STO-3G, obtained by adding one orbital to each STO-3G 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 one-electron atomic integrals between the new orbitals and the original orbitals vanish by symmetry. The only non-vanishing parts are the two-electron integrals involving an even number of orbitals such as . Unfortunately, these non-trivial interactions are not captured by the mean-field wave function due to it being a one-electron 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 one-site 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 mean-field bath could be common in realistic systems, especially when large-sized bases are used.

Vii Conclusion

In conclusion, we have introduced Incremental Embedding, a new Schmidt-space 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 one-site effective Hamiltonians for all sites by hierarchically incorporating corrections from embeddings involving two-site fragments, three-site fragments, and etc to the mean-field 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 many-site generalization of one-site DMET. It can be made either self-consistent or non-self-consistent. 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 atom-centered 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 self-consistency in site densities and pair-densities improves the performance of IE considerably near equilibrium geometry, through an approach that is not density-driven. For larger bases, both IE and one-site 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 Schmidt-space 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 high-level 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 Foster-Boys 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 Pipek-Mezey schemePipek and Mezey (1989); Boughton and Pulay (1993) and the Edmiston-Ruedenberg schemeEdmiston and Ruedenberg (1963). Last but not least, the conflict between the half-filled embedding space and the maximum number of entangled sites in a HF bath calls for a better bath wave function. In this regard, the Hartree-Fock-BogoliubovBogoliubov et al. (1958); Bogoliubov (1959) (HFB) wave function might be a good candidate because it is (i) always half-filled in the quasi-particle space, and (ii) still a mean-field theory and therefore retains the simplicity of embedding Hamiltonians in a mean-field bath.

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. CHE-1464804). TV is a David and Lucille Packard Foundation Fellow.


  • Bartlett and Musiał (2007) R. J. Bartlett and M. Musiał, Rev. Mod. Phys. 79, 291 (2007), URL
  • Hachmann et al. (2007) J. Hachmann, J. J. Dorando, M. Avilés, and G. K.-L. Chan, J. Chem. Phys. 127, 134309 (2007), URL
  • Chattopadhyay et al. (2015) S. Chattopadhyay, R. K. Chaudhuri, and U. S. Mahapatra, J. Comput. Chem. 36, 907 (2015), ISSN 1096-987X, URL
  • Taffet and Scholes (2017) E. J. Taffet and G. D. Scholes, Chem. Phys. (2017), ISSN 0301-0104, URL
  • 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
  • Stemmle et al. (2018) C. Stemmle, B. Paulus, and O. Legeza, Phys. Rev. A 97, 022505 (2018), URL
  • 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. Head-Gordon, 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 1759-0884.
  • 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 1678-4448.
  • 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énez-Hoyos, 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 0009-2614.
  • Pulay and Saebø (1986) P. Pulay and S. Saebø, Theor. Chim. Acta 69, 357 (1986), ISSN 1432-2234.
  • 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
  • Godby et al. (1988) R. W. Godby, M. Schlüter, and L. J. Sham, Phys. Rev. B 37, 10159 (1988), URL
  • van Schilfgaarde et al. (2006) M. van Schilfgaarde, T. Kotani, and S. Faleev, Phys. Rev. Lett. 96, 226402 (2006), URL
  • 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, (2017).
  • Chan and Head-Gordon (2002) G. K.-L. Chan and M. Head-Gordon, 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).
  • Olivares-Amaya et al. (2015) R. Olivares-Amaya, 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 1096-987X.
  • 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 1521-3979.
  • Bogoliubov (1959) N. N. Bogoliubov, Sov. Phys. Usp. 67, 236 (1959).
Comments 0
Request Comment
You are adding the first comment!
How to quickly get a good reply:
  • Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
  • Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
  • Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minimum 40 characters and the title a minimum of 5 characters
Add comment
Loading ...
This is a comment super asjknd jkasnjk adsnkj
The feedback must be of minumum 40 characters
The feedback must be of minumum 40 characters

You are asking your first question!
How to quickly get a good answer:
  • Keep your question short and to the point
  • Check for grammar or spelling errors.
  • Phrase it like a question
Test description