Max Nusspickel,Basil Ibrahim,George H. Booth
Effective Reconstruction of Expectation Values from Ab Initio Quantum Embedding
Quantum embedding is an appealing route to fragment a large interacting quantum system into several smaller auxiliary “cluster” problems to exploit the locality of the correlated physics. In this work, we critically review approaches to recombine these fragmented solutions in order to compute nonlocal expectation values, including the total energy. Starting from the democratic partitioning of expectation values used in density matrix embedding theory, we motivate and develop a number of alternative approaches, numerically demonstrating their efficiency and improved accuracy as a function of increasing cluster size for both energetics and nonlocal two-body observables in molecular and solid state systems. These approaches consider the N-representability of the resulting expectation values via an implicit global wave function across the clusters, as well as the importance of including contributions to expectation values spanning multiple fragments simultaneously, thereby alleviating the fundamental locality approximation of the embedding. We clearly demonstrate the value of these introduced functionals for reliable extraction of observables and robust and systematic convergence as the cluster size increases, allowing for significantly smaller clusters to be used for a desired accuracy compared to traditional approaches in ab initio wave function quantum embedding.
Effective Reconstruction of Expectation Values from Ab Initio Quantum Embedding

Quantum embedding is an appealing route to fragment a large interacting quantum system into several smaller auxiliary “cluster” problems to exploit the locality of the correlated physics. In this work, we critically review approaches to recombine these fragmented solutions in order to compute nonlocal expectation values, including the total energy. Starting from the democratic partitioning of expectation values used in density matrix embedding theory, we motivate and develop a number of alternative approaches, numerically demonstrating their efficiency and improved accuracy as a function of increasing cluster size for both energetics and nonlocal two-body observables in molecular and solid state systems. These approaches consider the N-representability of the resulting expectation values via an implicit global wave function across the clusters, as well as the importance of including contributions to expectation values spanning multiple fragments simultaneously, thereby alleviating the fundamental locality approximation of the embedding. We clearly demonstrate the value of these introduced functionals for reliable extraction of observables and robust and systematic convergence as the cluster size increases, allowing for significantly smaller clusters to be used for a desired accuracy compared to traditional approaches in ab initio wave function quantum embedding.

Quantum chemical methods to describe explicit correlations in an ab initio many-electron system can be highly accurate, though their applicability is often stymied by a steep computational scaling with respect to system size, which (despite significant recent progress) limits their use for extended systems. To combat this, the locality of this correlated physics is increasingly exploited, enabling a reduction in scaling to be competitive compared to mean-field or density functional approaches, while remaining free from empiricism. The field of “local correlation” methods in quantum chemistry generally build these locality constraints in the particle-hole excitation picture of the system, localizing each of these spaces separately. While highly related, “quantum embedding” approaches from condensed matter physics are also increasingly coming to the fore as an alternative paradigm and being applied to quantum chemical and ab initio systems. A loose (and necessarily imperfect) characterization of a key difference in these approaches could be that quantum embedding does not build this locality from a particle-hole picture—rather, a fully local set of “atomic-orbital-like” degrees of freedom are chosen initially (which will in general have neither fully occupied nor unoccupied mean-field character), which we will call the “fragment” space, though it is also often called the “impurity” space for historical reasons in traditional quantum embedding literature. A larger space is then constructed by augmenting these fragment orbitals with additional orbitals (often called “bath” orbitals). These are designed to reproduce the quantum fluctuations, entanglement, and/or hybridization between the fragment and the rest of the system, as characterized by some tractable (generally mean-field) level of theory which can be performed on the full system. These individual local quantum problems of the fragment and bath orbitals define a “cluster”, which is then solved to provide the correlated properties of the original fragment space, potentially with a subsequent self-consistency then applied to update the original mean-field/low-level theory on the full system. The general algorithm in most quantum embeddings is therefore summarized as: a) fragment the system; b) for each fragment, construct a bath space describing the coupling to the wider system; c) solve an interacting problem in the cluster space of each fragment via a “high-level” correlated method; d) extract properties of the system; and e) optionally, perform a self-consistency to embed the correlated effects from the cluster model back into the low-level full system method to update the coupling between the fragments and environment. There are a large number of choices and variations within this general framework, including (but not limited to) how the bath space is defined (including the choice of “low-level” theory), how the interacting cluster Hamiltonian is constructed and solved, and the choice of self-consistent requirements. Furthermore, the fundamental quantum variables by which these quantities over the different spaces are characterized can vary, with dynamical mean-field theory and its variants working in a Green’s function (dynamical) formalism, while density matrix embedding theory (DMET) and its variants work in a wave function (static) formalism. However, these two approaches are not fundamentally distinct and can also be rigorously connected via a common framework. Much recent progress has been made in these various quantum embeddings and their application to ab initio systems, including the use of quantum computation as a high-level solver. The key point in all of these embedding approaches, however, is that the scaling with respect to the full system size is defined by only the scaling of the low-level (often mean-field) method, given the local nature of the auxiliary cluster problems. Furthermore, the choice of high-level solver and arbitrary atomic-orbital-like fragmentation allows for spaces which capture strong (albeit still local) correlation effects, beyond the traditional constraints of the particle-hole picture of most local quantum chemistry. In this work, we focus on a critical aspect of quantum embedding which we believe has received less attention, but which has substantial ramifications for its accuracy and applicability. This concerns how nonlocal properties (including the total energy) of the full system can be reconstructed from the independent cluster solutions of each fragment. We will assess the effect of the inherent locality approximation of quantum embedding on the convergence of different functionals of these nonlocal expectation values, and motivate and demonstrate new approaches which substantially accelerate convergence with respect to fragment and bath size in the embedding. While this is quite a technical work, the outcome is that general design principles by which different functionals can be devised become clear, including the N-representability of these estimates. Here, we focus on wave-function-based quantum embedding (we believe that the ability and approach for constructing appropriate functionals in a Green’s function perspective is clearer). We start from the density matrix embedding theory as the parent wave function approach to quantum embedding. We always consider fragments consisting of a single atom only and, where we seek to systematically improve expectation values, enlarge the cluster spaces by adding additional interacting bath orbitals. We believe that this is an efficient, “black-box” approach and avoids the ambiguities and nonmonotonic improvements in the alternative of defining a systematic expansion in the fragment sizes for ab initio systems, which also can suffer from reducing the symmetry of the problem. The expansion of the bath space is defined from the approximate interacting density matrix (or instantaneous hybridization) between the fragment and environment at a simple approximate second-order perturbation theory (MP2) and is controlled by a single cutoff parameter as detailed in ref (53). This provides cluster-specific bath natural orbitals (BNOs) as a controllable, reliable, and well-defined expansion of the bath space (and hence overall cluster space) of a given fragment. Furthermore, in this work, we will neglect considerations of self-consistency of the original mean-field (beyond where self-consistency is required for meaningful extraction and comparison of expectation values, e.g. to control total electron number). More extensive self-consistency to qualitatively change the original full-system mean-field will be considered in future work, but it is unlikely to change the conclusions of this work, especially as convergence with cluster size (either via expansion of the fragment or, as in this work, an interacting bath) obviates the effect of self-consistency. We start the paper by recapitulating the original DMET “democratic partitioning” for expectation values, which can be computed via fragment-projection of the reduced density matrices of each fragment. We then describe an improved approach for the total energy based around a cumulant functional for the two-body effects. We move on to an approach based around direct projection of wave function amplitudes, rather than density matrices, formally satisfying N-representability conditions (not satisfied in the aforementioned density-matrix approaches) and substantially improving expectation values. Finally, practical approaches and further approximations to these will be described to define an efficient protocol for arbitrary expectation values and high-level wave function descriptions. These different approaches are benchmarked for energetics and other nonlocal properties (spin correlation functions) on both molecular systems via the W4–11 test set and periodic systems, finding an efficient approach for reconstruction of nonlocal observables in quantum embedding of general fragmented systems.

In Figure 1, we show results from the different total energy functionals described in this work as a representative illustration of the significant difference that this choice can make in a simple system (chlorine dimer, minimal basis, two atomic fragments). For all the different total energy functionals shown in Figure 1, the same identical fragment and DMET cluster space is used, and the same (exact) solver and wave function are found from each cluster. The only difference is the choice of energy functional to reconstruct the total energy from the two cluster solutions, with E[γx, Γx] representing the traditional “democratically partitioned” DMET energy used primarily in the literature to date. As a summary regarding the reconstruction of expectation values from cluster solutions in this work, we find the main conclusions to be the following: 1. It is advantageous to separate factorizable products of lower-rank contributions to expectation values where possible. In this way, we can construct e.g. the two-electron energy from the one-body density matrix and two-body cumulant, rather than the two-body density matrix directly. We demonstrate that this improves estimators via the inclusion of implicit cross-cluster contributions to the expectation values. Many of the improvements in this work arise from implicitly building in contributions to expectation values (or the wave function itself) from products of different cluster contributions, coupling the individual cluster solutions in a nonlocal fashion, and minimizing the local approximations inherent to the embedding framework. 2. N-representability of estimators can be used as an effective guiding principle, i.e. that they can be derived from a valid “global” wave function from the combined cluster wave function solutions. We show how this can be achieved exactly. Defining this “global” wave function ensures conservation of many quantum numbers, e.g. total electron number, which removes the necessity for costly (and sometimes ill-defined) self-consistency conditions and can allow for variational estimators. Beyond this, enabling convergence with respect to bath size of each fragment can render self-consistency entirely unnecessary. 3. It is possible to construct estimators from this well-defined global wave function by casting the expectation value as a functional of the wave function amplitudes of each cluster, rather than combining density matrices or observables directly. Each factor of the wave function amplitudes in the expectation value has their occupied indices (symmetrically) projected onto the fragment spaces of each cluster to avoid double-counting, and a sum over all cluster solutions for each factor of the wave function parameters is introduced. 4. Along with dramatically improved estimates, this wave function approach furthermore avoids the requirement for optimizations of chemical potentials, as the global electron number is strictly conserved, and the condition of the union of the fragment spaces is only that they span the occupied, rather than the full, orbital space of the original system in order to converge to exact results (e.g., as the (interacting) bath is expanded to completeness for each cluster), vastly reducing the burden of large fragment spaces spanning virtual orbitals for calculations in realistic basis sets. 5. If the approach above results in an intractable scaling with respect to the number of fragments in the system, then a principled approximation can be made, which we show can fortuitously lead to effective cancellation of errors and an even faster convergence of the ground state energy to the exact global expectation values. We motivate and evidence these conclusions in the following sections, resulting in our final recommendation for an approach to one- and two-body expectation values from DMET and related wave function quantum embeddings in ab initio systems, while retaining at most an scaling with system size in the evaluation of these observables. We detail the practical implementation of these schemes for both an exact high-level solver, and an (arbitrary order) coupled-cluster framework. All results can be reproduced from our recently released Vayesta codebase for quantum embedding (Note: Embedding code and documentation can be found at https://github.com/BoothGroup/Vayesta), which interfaces with the PySCF code, with scripts to generate many of the results of this work (including the input required for Figure 1) found in the Supporting Information (SI).

Expectation values in DMET derived from operators which span more than one fragment are calculated from ‘democratically partitioned’ reduced density-matrices (RDMs). These can be written aswhere γ and Γ are one- and two-body reduced density matrices formed from the high-level solution of each cluster problem, aswhere the high-level cluster wave function, |Ψx⟩, includes the contribution from the doubly occupied environmental orbitals of each cluster. is an operator, which introduces a projection onto the fragment subspace of the cluster x in a symmetrically averaged fashion, e.g.In these, x labels both the Nfrag fragments and clusters, for which there is an unambiguous one-to-one correspondence, and Ppqx is the fragment projection matrix acting over the whole molecular orbital (MO) space, aswith C representing the MO coefficients, S the atomic orbital (AO) overlap matrix, and Cfx the coefficients of the fragment orbitals in cluster x. Note that a post mean-field cluster solver will only modify the 1-RDM in the cluster–cluster part, whereas the 2-RDM also acquires contributions in off-diagonal cluster–environment parts. As a result, the projection in eq 1 can be performed in the respective cluster spaces, whereas in eq 2 it has to be performed in the full system space, in order to take these changes in the off-diagonal parts into account. The projection onto the fragment space of each cluster is required since the bath spaces overlap significantly between different clusters, and we must project out any double counting arising from contributions from overlapping bath spaces. In contrast, the partitioning of the system into fragment spaces is considered to be a disjoint fragmentation, where the set of all fragment orbitals are an orthonormal set with no overlapping fragment spaces in different clusters. From the democratically partitioned density-matrices, the total energy can be calculated asNote that, in practice, one can avoid forming the full system density-matrices to calculate the total energy and instead calculate energy contributions directly from the individual cluster density-matrices over purely the cluster degrees of freedom. The contribution from the doubly occupied environmental orbitals can be integrated out by forming the Coulomb and exchange potential of the unentangled occupied orbitals of each DMET cluster via an effective one-electron potential. This leads to the more common expression for the DMET energy, equivalent to eq 8 via construction of democratically partitioned density matrices of eq 5, aswhere γx and Γx refer to the cluster reduced density matrices, with the projector purely acting in this cluster space and Nclx denoting the number of orbitals in the cluster x. The energetic effect of these states (static Coulomb and exchange contributions) is then included via the construction of the , which includes the potential from these unentangled states to the one-body Hamiltonian as where γcorex is the density matrix from these core states. We can exploit fragments that are (by symmetry) in identical chemical environments by only computing the cluster solutions and energy contributions of symmetry-unique fragments. In the rest of this work, the expression E[γx, Γx] will denote this democratically partitioned energy functional, shorthand for E[{γx}, {Γx}], denoting that the energy is computed from the set of individually constructed cluster one- and two-body RDMs. This energy expression is exact when two conditions are satisfied. First, it requires that the fragmentation of the full system is complete, i.e. that the union of the fragment spaces spans all degrees of freedom of the system. This condition ensures that the trace of the sum of the different fragment projectors is exactly equal to the total number of orbitals in the system, or alternatively, thatWhile it is a relatively mild condition to ensure that the combined fragment spaces span the generally localizable occupied space, ensuring that they span the (generally much larger) high-energy virtual space is harder to achieve and leads to much larger fragment spaces. This has required DMET simulations in realistic basis sets to augment fragment spaces with projected atomic orbitals (PAOs) to ensure this condition is fulfilled for reasonable results. The second criteria which must be fulfilled is that the individual cluster density matrices must be exact, which in general will require the clusters of each fragment themselves to be enlarged to completeness, either by increasing the size of the fragment or (interacting) bath space. This ensures that |Ψx⟩ → |Ψ⟩ and the projected density matrices of the clusters (eq 5) are equivalent to the projections of the exact density matrix over the whole system. Combined with the completeness of the projector (eq 10), this will lead to the exact energy from eq 9. Away from this exact limit, there are a number of drawbacks to this approach to compute properties from the DMET solutions. First, the reconstructed full-system density matrices (eqs 1 and 2) from the projected cluster solutions are not N-representable, meaning that they cannot be derived from a valid wave function. This can be seen as the democratically partitioned 1-RDM of eq 1,cannot be rewritten as a simple expectation value of some wave function . As a specific consequence, this can result in eigenvalues (occupation numbers) becoming negative or greater than two (in a restricted basis), violating the Pauli principle, and removing any variational guiding principle in the method. Furthermore, conserved quantities and good quantum numbers in the individual cluster solutions such as electron number (N), spin and its z-projection (Sz, S2), and other symmetries are not maintained in these composite full system descriptions. To mitigate some of these effects, a global chemical potential (or potentially a fragment-specific chemical potential) is almost universally optimized in DMET to ensure that at least an exact, integer number of electrons is recovered in these democratically partitioned density matrices. This can move electrons between the fragment and bath of each correlated cluster solution, such that the known global electron number is maintained as a constraint. While this corrects one known global quantum number in the density matrices, it does not correct others and does not in general restore N-representability of the full system RDMs. Furthermore, this necessitates costly additional self-consistent loops over the high-level calculations. Furthermore, this requirement of a chemical potential optimization to a known total electron number further underlines the importance of the constraint of eq 10, ensuring that the full (occupied and virtual) space is spanned by the fragments, as all can be partially occupied in the correlated state and the democratically partitioned density matrices must trace to the correct electron number. Numerical demonstration of the breaking of these N-representability constraints in the democratically partitioned RDMs will be given in Sec. 4 (with and without a global chemical potential optimization), also showing the deleterious effect on properties and computed energetics of the system that result.

In the following, we propose a simple alternative for the construction of democratically partitioned two-body density-matrices from DMET clusters, from which expectation values such as the energy can be calculated. Instead of partitioning the 2-RDM directly as in eq 2, we partition the two-body cumulant, K̃, defined (in a restricted basis) (Note: the two-body cumulant would often be denoted by λ, but we define by K̃ to avoid potential confusions with the lambda-amplitudes of coupled-cluster theory used later in the text) viaThe noncumulant (disconnected) contributions to the 2-RDM can then be reconstructed from the democratically partitioned one-body density-matrix, given by eq 1, such that eq 2 is replaced byThe difference between eqs 2 and 13 lies purely in the noncumulant contribution to the two-body density matrix, which can be written as the product of democratically partitioned 1-RDMs. In the standard DMET partitioning of eq 2, these are taken from a single embedding problem at a time, whereas the partitioning of eq 13 contains “cross-cluster” contributions, which can be seen by inserting eq 1 into the first term of eq 13:In this way, the nonlocal (correlated) one-body physics of two distinct clusters, x ≠ y, contribute to global two-body expectation values; the same is not possible in eq 2. This is expected to be important in cases where the orbitals p and q are far from the orbitals r and s, and are not spanned together in any single DMET cluster. In this case, the conventional DMET partitioning of eq 2 will not account for the relaxation of the external Coulomb and exchange potential of a fragment due to the (potentially correlation-induced) density changes within the other orbital set. In contrast, this will be implicitly included in the partitioning according to eq 13. This is the key physics where we expect a partitioning of cumulants to be superior for two-body physics to the traditional democratic partitioning approach of Sec. 2.1 for nonlocal expectation values. We will denote any total energies resulting from the democratically partitioned cumulant approach described in this section as E[γx, Kx] in the rest of this work. We present a simple example showing the difference between the partitioned density matrices and the partitioned cumulant approach in Figure 2 for a representative system of N2 in a minimal basis set. The fragment space consists of the five symmetrically (Löwdin) orthogonalized atomic orbitals of a single atom (the 1s, 2s, 2px, 2py, 2pz spaces) with the DMET bath consisting of an additional three bath orbitals consistent with the bond order of the dimer in a minimal basis. The DMET cluster in this example therefore contains 10 electrons in 8 orbitals, which is compared to the full system of 14 electrons and 10 orbitals and is solved with exact diagonalization (FCI). It is found that regardless of whether the fragment chemical potential is optimized or not, the democratically partitioned cumulant results in a substantial improvement in the energy, with the partitioned cumulant being almost unaffected by this chemical potential optimization. This is further consistent with the Cl2 dimer shown in Figure 1 and analyzed in more detail in Sec. 4, while the improvement in nonenergetic properties (the spin–spin correlation function) from the cumulant partitioning will be shown in Sec. 6.3. A minor technical detail to mention is that we use a slightly different partitioning in practice than that of eq 12 throughout this work. The correlated cluster 1-RDM can be decomposed as γ = γ0 + Δγ, where γ0 is the reference mean-field density-matrix. We then work with a slightly modified cumulant definition, K, asThe approximate cumulant K of eq 15 therefore contains the (Δγ)2-terms, which are not present in the true cumulant K̃, with the relation between the two given asThe final expression for the 2-RDM is thenwhere (in contrast to the democratically partitioned 2-RDM of eq 2), the projection of the cumulants can be performed in the cluster space only, since the correlated cluster solver will only lead to a nonzero cumulant in this space. This allows the total energy functional in terms of γx and Kx to be written aswhere EHF[γ0] and F[γ0] are the Hartree–Fock energy and Fock matrix corresponding to the reference mean-field density-matrix γ0, and Δγ is the correlated part of the democratically partitioned 1-RDM, formed from the cluster density matrices, γx, as shown in eq 1. Note that eq 18, in contrast to the conventional DMET energy functional (eq 9) which it replaces, does not involve a cluster-specific effective core-Hamiltonian. We find that democratically partitioning the approximate cumulant instead of the true cumulant gives almost identical results for all systems tested in this paper. However, in the case of the true cumulant, the Fock matrix of the correlated 1-RDM, F[γ0 + Δγ], would need to be constructed to calculate the energy, a step that scales as with respect to the system size N. For this reason, we use the democratically partitioned approximate cumulant and calculate the total energy according to eq 18, referred to as E[γx, Kx] in the results of this paper.

While the democratically partitioned cumulant described in the previous section can dramatically improve two-body properties (as will be further shown later), it still suffers from the same fundamental problem discussed in Sec. 2.1. This is that the resulting density matrices are not generally N-representable, i.e. they do not correspond to a physical Fermionic wave function over the system. In this section, we go a step further and present an alternative paradigm of directly combining the cluster wave functions rather than RDMs to overcome this representability issue. This was first proposed in a correlation energy functional in ref (53) (described in Sec. 3.3) but is now generalized and expanded. The basic idea is to consider an implicit full-system wave function reconstructed from the cluster wave functions themselves, asFrom this global state, expectation values can be defined, which resolves N-representability issues and results in improved estimators. For instance, the 1-RDM can be written asAs long as all projected cluster wave functions are physical, Fermionic wave functions in their own right, then the linear superposition of eq 19 is a physical, Fermionic wave function as well. As a result, the 1-RDM of eq 20 is trivially N-representable, as it is simply with the “global” wave function |Ψ⟩ given by eq 19. Higher order density-matrices can be constructed in the same way and are thus also N-representable. Furthermore, the full system energy computed from RDMs constructed in this fashion will be variational if the cluster solvers are themselves variational methods. It is not difficult to rationalize that expectation values would also be improved from the formulation as sketched for the 1-RDM in eq 20 compared to the democratically partitioned RDMs. This is because a large number of “cross-cluster” contributions to these expectation values are included, where wave function amplitudes from different clusters are combined beyond their mean-field contributions found in eq 11.

So far we have not specified what we mean by the projection of a cluster wave function as in practice. Since the DMET bath orbitals guarantee that the mean-field reference determinant |Φ⟩ is represented exactly in each cluster (ignoring the unentangled environment orbitals), it is convenient to define in terms of its action on the correlated part of the wave function only; i.e. our full system wave function can be written aswhere |ΔΨx⟩ = |Ψx⟩ – |Φ⟩, and intermediate normalization is assumed. We then choose to represent the correlated cluster wave function part |ΔΨ⟩ in the basis of particle–hole excitations around the reference determinant in a linear wave function ansatzwhere we use i, j, ... (a, b, ...) to represent general occupied (virtual) orbitals in this work. We choose the fragment projection of this wave function to be defined by its symmetric action on the occupied coefficient indices, for examplein the case of the single and double excitation coefficients in the cluster, C1x and C2x, with generalization to higher-order excitation levels straightforward. Indices in these expressions correspond to the cluster canonicalized orbitals, with a critical consequence of the DMET (and BNO) bath construction being that the n-fold excitations of the cluster space are entirely spanned by the n-fold excitations of the full system due to the coincidence of their reference states. The fragment projection operators acting in the occupied-only space of the cluster can be formulated aswhere Cfx represents the AO coefficients of the fragment orbitals (x) of fragment x, and Cx represents the coefficients of the occupied cluster orbitals (in contrast to the projector defined in the full-system MO space of eq 7), with i, k representing occupied orbitals of cluster x. We note that this projection operator is not diagonal in the occupied orbital basis of the cluster. Furthermore, we have crucially chosen to only apply the projection to the occupied dimensions of the wave function coefficients in eqs 23 and 24. This is an important difference compared to the projection of density matrices according to eq 5, where the projector always needs to act on general indices, which enumerate both occupied and virtual orbitals. The significance of this is that the initial fragmentation of the system now only needs to ensure that the entire occupied space is spanned by the union of all fragments, in order for the projector to be complete. If this is satisfied, the only remaining approximation results from the deviation of each cluster wave function from exactness (which can be systematically resolved via increasing the cluster/bath space). This is a significant advantage, ensuring that the expectation values become exact as the bath space of each cluster is increased, without requiring the fragment spaces to span the virtual space of the system. The requirement for completeness of the fragment projectors is now relaxed from eq 10 and can be written aswhere Ne is the number of electrons. This quality is especially impactful in larger basis sets required for quantitative accuracy, as fragment sizes are now independent of the size of this basis, with the required virtual space captured instead by the bath expansion. A fragmentation spanning the occupied space is, for example, guaranteed by the choice of intrinsic atomic orbitals (IAOs) for the fragment spaces, used in the main results of this work unless otherwise specified. We note that neither the representation of the wave function amplitudes in a particle–hole basis nor the choice of projection of just the occupied space in itself precludes the solver from capturing strong correlations, since we do not assume a truncation of the excitation levels which are represented in the wave function. The ability to capture strong correlation effects (for an exact high-level solver) is determined by the suitability of the cluster space and is invariant to the choice of representation of the resulting wave function. It could, however, be possible that a more efficient projection operator could be formulated for, e.g., strongly correlated lattice models, where the requirement of the fragmentation spanning the virtual space is not difficult to fulfill. However, for ab initio quantum embeddings, we do not anticipate that this would be beneficial.

While all wave functions can be written in the linear form of eq 22, we now show that there are benefits in casting the wave function to an exponential parametrization for partitioning in this space. This form (common to coupled-cluster methods) can be written asThe amplitudes can be converted between the linear and exponential forms more generally via the relationsFor ground state wave functions, the norm of the T-amplitudes of the exponential ansatz generally decays more quickly with respect to the excitation rank than the C-amplitudes of the linear ansatz. The exponential wave function ansatz thus allows for an alternative way to partition the wave function by projecting the T-amplitudes, rather than C-amplitudes, which by analogy with eqs 23 and 24, can be written as To illustrate the difference between the partitioning of the global wave function via its C- and T-amplitudes, we can compare the wave function ansätze truncated at the double-excitation level. For the partitioning of the linear C-amplitude representation of the wave function (Sec. 3.1), we obtainwhile for the exponential form of the wave function, we achieveWhile the first two terms in eqs 33 and 34 are equivalent and involve only a single summation over clusters, the last term of eq 34 (representing disconnected double excitations) has a double loop over pairs of clusters. As a result, a single excitation in cluster x and a different single excitation in cluster y can contribute together to form a double excitation of the global wave function, even if the two clusters are far apart and do not overlap. Similarly, for wave functions truncated at the triple excitation level, two independent clusters (e.g., in T2xT1y) or three independent clusters (e.g., in T1xT1yT1z) can combine to contribute disconnected triple excitations to the global wave function. The same is not possible for the linear wave function partitioning, for which both connected and disconnected contributions must always come from a single cluster only. This partitioning of the global wave function at the level of the T-amplitudes therefore results in nonlocal “cross-cluster” information being built into the solution from the quantum embedding—this time on the level of the reconstructed global wave function itself, rather than similar cross-cluster information being built into the 2-RDM (see Sec. 2.2) or other expectation values (see eq 20). The beneficial impact of these cross-cluster contributions is a key tenet of this work and serves to mitigate the impact of the local approximations inherent to quantum embedding by not treating the cluster solutions as contributing independently to final expectation values. Finally, we should stress that while this framework of partitioning the global T-amplitudes of the wave function naturally fits with the use of coupled-cluster as a high-level solver, we are not inherently constrained to this. If we avoid truncating the rank of the T-operator, then any wave function can be cast into the exponential form of eq 27 for the high-level solution of the cluster problems and benefit in this way. In particular, a linear wave function obtained from an FCI cluster solver can be converted to the exponential representation via eqs 28–30, projected, and finally recombined into a global exponential wave function to benefit from these cross-cluster contributions. Nevertheless, the partitioning of the global T-amplitudes of the wave function is also particularly suited to coupled-cluster high-level solvers, which are increasingly being used as a compact and accurate approach which enables access to larger fragment and bath spaces and have recently been shown to be accurate in a variety of embedding contexts.

In this section, we consider the use of this implicit “global” wave function to compute expectation values of operators which commute with the Hamiltonian and can therefore be calculated from functionals which are linear in the wave function. A trivial example of this is the total energy, as its associated operator is the Hamiltonian itself. Projecting the time-independent Schrödinger equation from the left with the Hartree–Fock determinant ⟨Φ| and assuming intermediate normalization ⟨Φ|Ψ⟩ = 1, it follows thatwhich is linear in Ψ. This is the traditional energy functional of, e.g., coupled-cluster theory, but holds for any state. Since the Hamiltonian operator of an electronic structure problem involves up to two-body interactions, only the double excitations contribute to the correlation energywhere the single excitations do not contribute for a Hartree–Fock reference state, due to Brillouin’s theorem (though a single-body contribution can easily be included if using a non-Hartree–Fock reference). Projecting the C2-amplitudes in eq 36 according to eq 24 enables us to partition the correlation energy into cluster contributions as These contributions can be formulated from the individual cluster solutions, where the orbital indices run over the occupied or virtual states of each cluster, and where the projection operator of eq 25 can be constructed in this cluster orbital space, to evaluate the energy efficiently in a cost linear in the number of clusters. This is the energy functional introduced and used in ref (53), and it will be denoted by E[Ψx] throughout this work to indicate that it is constructed directly from the cluster wave function amplitudes without requiring intermediate cluster density matrices. The use of this functional for energies is common to quantum chemical local correlation methods, for example within the PNO or cluster-in-molecules approaches. However, the “fragments” in these methods are then generally defined in terms of (localized) purely occupied orbitals, instead of the general “atomic-orbital-like” fragment spaces which are used in the context of quantum embedding and which in this work require the use of explicit (nondiagonal) projection operators as described. This functional can also be used within a fragment projection of the T-amplitudes of an exponential representation of the high-level cluster wave function (see Sec. 3.2). In this case, intercluster contributions arise from terms. This physics is neglected in the C-amplitude projection of above in the case of no overlap between clusters x and y. However, in the examples in this work, the projection of the C- and T-amplitudes gave almost indistinguishable results—unsurprising given the small contribution of the overall T12 contribution to the energy and the fact that the local portion of this is still included. For this reason, and to reduce the computational cost of this energy functional and keep consistency with ref (53), all results in this work denoted E[Ψx] will correspond to the projection of the C-amplitudes of the cluster state (even for coupled-cluster high-level solvers) as shown in eq 37, though the projection of T-amplitudes may be preferred in the future.

General expectation values, with operators that do not commute with the Hamiltonian, require a quadratic functional of the wave function, of the form . In the case of one- and two-body expectation values, O1 and O2, these can also be expressed in terms of the one- and two-body reduced density-matrices (which in turn depend quadratically on the wave function) asThe one- and two-body RDMs are thus important quantities to be able to compute for expectation values beyond the energy. Of course, with knowledge of the density matrices, one can also calculate the total energy via eq 8. For most quantum chemical methods, where a globally stationary state is found, the result will be identical to the projected energy of eq 35. However, within the quantum embedding framework, the wave function is partitioned into fragment contributions which are then calculated in an approximate fashion, and we do not expect this equivalency to hold. Therefore, while the individual cluster wave functions are invariant to this choice of energy functional, the reconstructed full system total energy may not be. Indeed, there is likely a benefit in calculating the total energy as a symmetric (quadratic) expectation value (or equivalently from the RDMs), since the error in the energy will reduce quadratically with the error in the wave function, rather than linearly as in the energy functional of eq 35. In sections 3.1 and 3.2, we showed how global wave functions can be defined in terms of contributions of projected cluster wave functions. In principle, these global wave functions then fully and uniquely determine their corresponding full system reduced density-matrices. For example, the 1-RDM of a partitioned linear wave function can be written aswhere λ and μ enumerate the excitations (to all included orders) within the respective cluster in the particle-hole basis. Note that in contrast to the democratically partitioned 1-RDM of eq 11, the 1-RDM of a partitioned wave function contains a summation over cluster pairs, thus incorporating simultaneous excitations within two different clusters, as indicated above and in eq 20. For a partitioned exponential wave function, the global state itself already contains nested summations over projected cluster wave function T-amplitude contributions (see eq 34). The resulting 1-RDM will thus have even higher order contributions than the C-amplitude projection (going up to three simultaneous cluster contributions for a high-level wave function represented up to T2, and up to five simultaneous clusters for the 2-RDM). Both approaches result in explicitly N-representable RDMs, but again the T-amplitude projection can be argued to have a larger number of physical nonlocal “cross-cluster” contributions and therefore represents a better approximation. In Sec. 5, we will detail the technical aspects of how all of these cross-cluster contributions can be efficiently computed while retaining no more than DFT scaling in the full system size (or number of clusters). Total energies computed from these T-amplitude projected RDMs as described in this section will be denoted E[(γ, Γ)[Ψx]] throughout this work, to signify that these are density matrices derived from the global partitioned wave function.

We consider a simple example to corroborate the claims in the work so far, before describing specific implementational details in Sec. 5 and considering a more extensive set of test systems in Sec. 6. We consider the binding of Cl2 in a minimal basis. The choice of a minimal basis is intentional, since this means that we can easily compare to exact FCI results and that the same (Löwdin orthogonalized) atomic orbital space can be used for the fragments in all approaches. For larger basis sets, the democratically partitioned density matrix and cumulant approaches require all orbitals to be spanned by the fragment spaces. However, the projected wave function approaches require only (at least) the occupied space to be assigned to fragments for the projector to be complete. Therefore, for more realistic basis sets, the fragment spaces would generally be chosen according to different criteria, with the fragment spaces of the former needing to grow with basis size. Using a minimal basis therefore avoids this issue and allows the different methods to be compared on the same footing. Consequently, the DMET cluster space in this example consists of 18 electrons in 9 fragment orbitals plus a single DMET bath orbital (due to the minimal basis) in each of the two symmetrically equivalent clusters, while the full space consists of 34 electrons in 18 orbitals. The interacting cluster Hamiltonian is solved exactly via FCI, and we compare the four different energy expressions outlined to far in Figure 3: the democratically partitioned RDM energy common to DMET approaches to date (E[γx, Γx] from Sec. 2.1), the democratically partitioned cumulant energy (E[γx, Kx] from Sec. 2.2), the linear energy functional (E[Ψx] from Sec. 3.3), and the energy from the RDMs of the T-amplitude projected global wave function (E[(γ, Γ)[Ψx]] from Sec. 3.4). We include calculations with and without a global chemical potential optimization to ensure that the right number of electrons are conserved for the full system 1-RDM, but no further optimization of a correlation potential is considered here. When used (right column of Figure 3), this fragment chemical potential, μfrag, is added to the one-electron cluster Hamiltonian within the fragment space only and optimized such that the density-matrix traces to the correct number of electrons (these μfrag-optimized results are the same as shown in Figure 1). All results are subject to identical locality approximation of DMET, whereby the full system is represented only by the two wave functions in the cluster subspaces, allowing for a faithful comparison of the quality of the reconstructed expectation values. An example Python script, showing how these calculations can be performed using PySCF and Vayesta, can be found in the SI. We find in Figure 3 that the democratically partitioned RDM energy is changed significantly by the optimization of the fragment chemical potential, but in neither case are the results reasonable, with unphysical results without μfrag-optimization and substantial overbinding of the dimer when this is included. These results are also found in many other systems (barring hydrogen dimers and chains, where we find the democratically partitioned density matrix energy accurate, as has been noted elsewhere). The democratically partitioned cumulant approach is also changed by the μfrag-optimization, but is much less sensitive to this, with both results already significantly improved over the democratically partitioned RDM energy. The chemical potential optimization further drops the nonparallelity error of the democratically partitioned cumulant approach from ∼6 mEH to just over 1 mEH in this system. We also show the energies derived from partitioned wave functions over the system, calculating the energy from either a linear or quadratic functional (via the RDMs) of the wave function probability amplitudes. It is shown that these are also very accurate, which we can ascribe to two properties of these functionals. The first is that they are intrinsically N-representable, ensuring that they fulfill the physical constraints of being derivable from a wave function. As a consequence, they necessarily correspond to the correct, integer number of electrons, and therefore obviate the necessity of chemical potential optimization which has no effect (therefore also substantially reducing the cost of the calculations). The second rationalization for the particularly good performance of the partitioned wave function RDM energy is due to the introduction of contributions from products of wave function amplitudes from different clusters, coupling the cluster solutions together. Furthermore, for the E[(γ, Γ)[Ψx]] energy, these contributions result in an error in the energy which is rigorously quadratic rather than linear in the wave function error, and introduces a variationality (for a variational method) into the results, which results in the smallest nonparallelity error of 0.5 mEH. In order to analyze the different approaches in more detail, in Figure 4, we compare the 1- and 2-RDMs resulting from the democratically partitioned RDMs (γ[γx] for the 1-RDM, and Γ[Γx] or Γ[γx, Kx] for the 2-RDM from the direct democratic partitioning and the cumulant partitioning, respectively). We also compare to the RDMs constructed from the partitioned wave function expansion (γ[Ψx] and Γ[Ψx]). We consider the norm of the errors compared to the exact full system 1- and 2-RDMs, as well as the electron number error from the trace of the 1-RDMs. We also introduce a measure of the N-representability error of the 1-RDM, Ω, which sums the deviation from allowable occupations of the system orbitals aswhere ni are the eigenvalues of γ. While the overall norm error of the 1-RDMs seem relatively similar between the approaches (especially after chemical potential optimization), it can be seen that there is still a significant electron number error (without chemical potential optimization) and N-representability error in this quantity. While the chemical potential optimization fixes the electron number error, it is actually found to increase the N-representability error. In contrast, the partitioned wave function approach ensures that the electron numbers are exactly fulfilled even without a chemical potential and that the N-representability error is strictly eliminated. For the 2-RDM errors, it is clear that significant improvements are found by going to a democratically partitioned cumulant compared to the 2-RDM directly, further indicating that it is the noncumulant error in the 2-RDM which is contributing to the significant energy errors in the standard DMET energy functional.

While the previous sections describe the principle and numerical advantages behind the reconstruction of a global wave function, in practice we will often want this to be implicit as the reconstruction of the global correlated wave function will in general be prohibitive in cost. Instead, we want to directly compute the RDMs or expectation values of interest from this state in low-polynomial time with respect to system size, and without combining the wave function amplitudes into an explicit global state. This allows for the cluster wave functions to remain distributed in their fragmented cluster representations. We focus in this section on the efficient computation of the 1- and 2-RDMs, from which all two-body properties can be computed, with the principle of an implicit global wave function demonstrated in Sec. 3.4 underlying this construction. We restrict ourselves here to a projection of the T-amplitudes of the global state, represented in an exponential form. This is consistent with the approach in the previous section, maximizes the number of nonlocal combinations of cluster solutions, and is a natural choice for coupled-cluster high-level solvers (although as explained in Sec. 3.2, other solvers can be cast in this form). Furthermore, we focus on an efficient implementation of the cluster solver truncated at the T2-amplitudes, i.e. CCSD (or MP2), which is a cost-effective high-level cluster solver for ab initio systems, yet is sufficiently computationally inexpensive for us to rigorously demonstrate the convergence of these functionals with respect to bath size in Sec. 6. Extensions to other high-level cluster solvers (e.g., FCI) can proceed via eq 40 and will be described in future work. We stress here that this approach to compute expectation values from a global wave function with alternative (i.e., strongly correlated) solvers would not lead to an increase in scaling with respect to full system size compared to the algorithm shown. Indeed, the fact that coupled-cluster RDMs involve high-degree polynomials of the cluster solution variables (i.e., Tn amplitudes) results in the formal scaling with system size being generally larger than solvers such as FCI (where the RDMs are just quadratic in the cluster Cn variables). Explicit reconstruction of the global CCSD T-amplitudes (which we call the ‘global T2 algorithm’) would require a memory overhead scaling of and computational scaling of and for the 1- and 2-RDMs, respectively (where N is a measure of full system size) via the usual CCSD equations. We reduce this scaling down to cost in memory and time for the 1-RDM construction via direct use of the cluster amplitudes and construction of appropriate intermediates without introducing any additional approximations. This is described in Sec. 5.1, which we denote the “distributed T2 algorithm”, and is crucial to ensure that this step is a subleading scaling compared to the initial mean-field calculation and for applicability to large systems. For the construction of the 2-RDM and properties derived from it, the requirement of a similar quadratic scaling with system size necessitates the introduction of a further approximation, which is described in Sec. 5.2, with rigorous validation of this further approximation demonstrated in Sec. 6.

In order to construct density-matrices at the CCSD level, both wave function T-amplitudes and the Lagrange multipliers (or Λ-amplitudes) are required. These are optimized within each cluster after the T-amplitudes are found, giving rise to a set of Λx amplitudes for each cluster, x. The principle of reconstruction of a “global” set of Λ-amplitudes follows symmetric projection of the occupied indices of each cluster Λx-amplitudes onto the fragment space (defined analogously to eqs 31 and 32), and a summation over the cluster solutions, as All expectation values in coupled cluster (including the 1-RDM considered here) involve sums of expressions which are polynomials in the T-amplitudes and up to linear in the Λ-amplitudes. The central idea of the efficient “distributed” approach is to avoid the explicit reconstruction of the global T- and Λ-amplitudes entirely, and instead iterate over tuples of cluster solutions, directly summing in their corresponding contribution to the 1-RDM. The number of clusters which need to be looped over in order to include all cross-cluster contributions to the desired expectation value (or 1-RDM in this instance) depends on the maximum order of the polynomial in the expression, i.e. the total number of T- or Λ-amplitudes which are contracted together. For the 1-RDM, this maximum order is formally three, corresponding to contractions of the type Λ1T12. This corresponds to the maximum rank of simultaneous “cross-cluster” information in the construction of this quantity (i.e., triplets of cluster solutions are required). However, we can reduce this by first explicitly constructing the global T1- and Λ1-amplitudes from the cluster solutions according to eqs 31 and 43, since this will only require a scaling of and in memory and time, respectively, which is an acceptable overhead compared with the initial mean-field calculation. This is in contrast to the analogous explicit construction of the global T2-amplitudes, which we have already argued is prohibitive. The fragment projection for T1 and Λ1 is therefore performed in each cluster and then projected back into the full system MO space and summed over all clusters. Furthermore, we store these global T1- and Λ1-amplitudes on every memory partition of a distributed memory parallel calculation. The key result is that we now only need to loop over cluster tuples up to the maximum rank of distributed (i.e., T2 or Λ2) amplitudes in the 1-RDM expression, which is only quadratic, resulting from an L2T2 contraction. This requires only an double summation over pairs of cluster amplitudes (which can be further reduced asymptotically down to as described later). If certain contractions require a T1 to be contracted with a cluster-distributed T2, then the global T1 amplitudes can be projected into the required cluster space of the T2. We demonstrate this distributed amplitude algorithm for the example of the CCSD contribution to the occupied–occupied block of the 1-RDM, given asThe first term of this contribution can readily be calculated from the explicitly combined global singles amplitudes with scaling. For the second term, we can calculate the contribution arising from a specific cluster pair (x, y), by using the Λ2- and T2-amplitudes corresponding to these clusters, aswhere we use the notation ix (ax) to indicate an occupied (virtual) orbital of cluster x, and Sxy to represents the overlap matrix between the occupied (or virtual, as indicated by the subscript indices) orbitals of cluster x with cluster y, given byThese Sxy matrices can be precomputed in for all cluster pairs, where the columns of Cx are the relevant orbitals of cluster x in the AO representation, and S is the AO overlap matrix. We then sum the 1-DM contributions of eq 46 over all cluster pairs, according towhere R represents the overlap between cluster MOs and MOs of the full system, i.e. To analyze the computational scaling of the distributed amplitude algorithm, we note that the number of clusters, AOs, and full system MOs grows linearly with the system size, whereas the number of cluster orbitals per cluster remains constant. As a result, the calculation of all cluster pair contributions (eq 46) scales as , while eqs 47–49 can be evaluated in time, if intermediates are formed as indicated by the brackets. An analogous loop over cluster pairs can also be used to calculate the virtual–virtual and mixed occupied-virtual blocks of the 1-RDM, resulting in an overall scaling algorithm. In Figure 5 we compare the scaling of the global amplitude and distributed amplitude algorithms for a series of alkanes between C6H14 and C22H46 in a cc-pVDZ basis and atomic IAO fragmentation. To extend beyond a minimal DMET bath size, we augment the bath expansion with “bath natural orbitals” (BNO), which can be systematically enlarged in order to describe the beyond-mean-field coupling of the fragment to the environment at the level of approximate MP2 theory. This bath construction is similar to pair natural orbitals (PNOs)—though traditional PNOs are applied only for the virtual coupling of a localized electron pair. The completeness of this bath space is controlled by a threshold η, where the bath becomes complete as η → 0 and is limited to just the traditional DMET bath as η becomes large. More details on the construction and motivation of this bath expansion can be found in ref (53). From the timing data for the 1-RDM construction in Figure 5 for two different bath space sizes, we fit the polynomial t(N) = aNb to the last six data points of each curve to determine the overall computational scaling. For the global amplitude algorithm, this is close to the expected scaling with b ∼ 5, whereas for the distributed amplitude algorithm, we find an exponent close to 2, indicating a quadratically scaling algorithm, instead of the expected exponent of 3. This indicates the small prefactor of the efficient operations of eqs 47–49, whereas the evaluation of eq 46, although scaling as , has a significant larger prefactor and dominates the computational costs for the tested system sizes. Despite this favorable scaling, we employ a number of additional techniques to further improve the efficiency of this algorithm. These include a singular value decomposition of each cluster overlap space, in order to find the most compact domain for describing the intercluster physics. In the large system limit, this also naturally leads to a scaling in the 1-RDM construction, due to the fact that each cluster will only have appreciable overlap with other clusters. In addition, we consider efficiency gains that can be made from k-point sampled periodic systems, effective distributed memory parallelism, and a compressed representation of the projected amplitudes. These technical improvements in the 1-RDM construction are detailed in Appendix A.

In principle, an analogous approach to the above can be used to recover any expectation value from the implicit partitioned wave function, with the scaling determined by the maximum number of products of Tn or Λn in any of the terms, where n > 1. For the 2-RDM, this increases to three (from two for the 1-RDM), which would necessitate looping over triplets of cluster solutions (noting that a similar SVD approach to screen the overlap would still asymptotically reduce to linear scaling). Note that including the T1 and Λ1 contributions, the 2-RDM would include simultaneous contributions of up to five cluster amplitudes, but by explicitly reconstructing these T1 and L1 amplitudes globally, only looping over three clusters (x, y, z) would be required. However, including the projections back to the full space and noting the storage of the global 2-RDM, these scalings are likely prohibitive if we aim for the quantum embedding to maintain a scaling with full system size which is no more than mean-field theory. To return to (at most) a scaling approach, an additional approximation is thus required. As a first step, we split the 2-RDM into products of 1-RDMs and an (approximate) cumulant, K, defined according to eq 15. The noncumulant part can be treated without further approximations, using the efficient calculation of the 1-RDM presented in the previous section. This ensures that all the cross-cluster contributions in the noncumulant part are exactly included from an implicit valid wave function. The cumulant part, on the other hand, will be treated in an “in-cluster approximation”, meaning that all products of wave function amplitudes are taken within a single cluster at a time, rather than explicitly including triplets of clusters. In order to illustrate this, let us consider a contribution to the cumulant, δK, which involves a triple product of double (de)excitation amplitudes. If evaluated exactly according to the implicit partitioned wave function, this would require contributions from all triplets (x, y, z) of clusters, asIn the in-cluster approximation, this contribution is replaced bywhich only contains amplitudes from cluster x. Note that eq 51 only contains a single fragment-space projector, which is applied to the Λ-amplitude, compared to eq 50, where every amplitude is projected. In general, the number of projectors in a term has to agree with the number of nested loops over clusters the given term is summed over. This ensures preservation of a core principle of the embedding: that we guarantee an exact result, free from double-counting, as the bath space increases to completeness, as long as the sum of the fragment spaces spans the entire occupied space (eq 26). Instead of projecting the Λ-amplitude, we could have also chosen to project the first or second T-amplitude in eq 51, resulting in a different in-cluster approximation, which nevertheless becomes exact in the full-bath limit. However, projecting the Λ-amplitudes yields a simple recipe, since all contributions to the 2-RDM cumulant are at most linear in Λ1 or Λ2. Some limited experimentation with other options for the projected amplitudes did not seem to significantly change the results. The only term this prescription does not work for is the Λ0 terms of coupled-cluster expressions, which results in an MP2-like contributionin which case we use projected T-amplitudes, according toWe note here that the “global” T1- and Λ1-amplitudes as defined in the last section are not used in this “in-cluster” cumulant approximation, ensuring that only a single cluster summation is used for all terms. Relaxation of this constraint will be investigated in the future. This “in-cluster” approximation to the two-body cumulant reduces its computation to sums over single cluster wave function contributions, avoiding the computational effort associated with the fully nonlocal approach (as used for the 1-RDM in Sec. 5.1). The computational effort to compute the two-body cumulant contributions is now negligible compared to the 1-RDM construction, which requires contributions from all pairs of clusters. However, we do not want to build the whole full-system cumulant (which would require a prohibitive scaling in memory), but rather calculate expectation values directly from this sum over the independent cluster wave functions. We can therefore directly construct the energy functional with this “in-cluster” approximated cumulant, which we denote E[(γ, K*)[Ψx]) in this work (the asterisk denoting the “in-cluster” approximation to the two-body cumulant), aswhere indicates that each term in the expression for the (approximate) cumulant of eq 15 is computed from products of contributions from the same cluster solution (the “in-cluster approximation”). This allows the two-body contributions from each cluster to be computed entirely independently, as opposed to Δγ[Ψx], which is assumed to contain all relevant products of different cluster solutions (see Sec. 5.1). Other static two-body expectation values can be derived using this “in-cluster” two-body cumulant analogously. In Sec. 6.3 and Appendix B, we go beyond energetics to consider the two-point instantaneous spin correlation function, which can be computed in this same framework via reconstruction of a global 1-RDM, with cluster-local two-body contributions via the “in-cluster” approximated cumulant. While this approach bears similarities to the democratically partitioned cumulant for two-body properties of Sec. 2.2 (specifically that the disconnected part of the 2-RDM is treated separately to include cross-cluster contributions), it has significant advantages. In particular, the 1-RDM part is treated in a fully N-representable way, with contributions from all triplets of different cluster wave functions, as well as the fact that the projections acting on the occupied spaces mean that fragments only need to span the complete occupied space to ensure convergence to “in-method” exactness as the bath spaces of each fragment are enlarged. We consider this a critical aspect to allow for practical extension of the embedding to realistic basis sets. Although this “in-cluster” approximation is justified primarily by numerical expediency, it is still a controllable approximation in the same fashion as the embedding approximation itself, via an expansion of the bath space, guaranteeing systematic improvability of the approximation to exactness. However, it is also clear that the “in-cluster” approximation to the two-body cumulant does result in a loss of strict N-representability of the resulting 2-RDM. This is because it is no longer derivable from a full system wave function, since there is no single set of “global” T- and Λ-amplitudes used for all expressions. This can be seen in the first set of results in Figure 1, which presents all energy functionals used in this work. While the error in the energy functional which builds the RDMs from the explicit full system partitioned wave function (E[(γ, Γ)[Ψx]]) is variational with respect to the true solution, this is lost when the low-scaling “in-cluster” approximation is used for the two-body cumulant (E[(γ, K*)[Ψx]]). In the next section, we will rigorously benchmark these different approaches to computing energies and nonenergetic spin–spin correlation functions, and analyze their convergence in realistic molecular and extended systems with systematic expansion of the bath space of each cluster. Surprisingly, we find that energies derived with the “in-cluster” cumulant approximation exhibit an improved convergence to the exact complete bath limit of the embedding theory, due to more favorable cancellation of errors.

We now turn to a larger selection of systems to consider the convergence of these different energy estimators to their “in-method” exact limit, as the size of the interacting bath expansion of each cluster is enlarged. For this expansion, we consider the ‘cluster-specific bath natural orbitals’ (BNO) introduced in ref (53), which can be used as a systematic way to increase the orbital space controlled by a single threshold parameter η, and we will use CCSD as a cluster solver throughout. As η is reduced toward zero, the bath space increases in size toward completeness, capturing longer-ranged and higher-energy coupling of the fragment to its environment. The use of CCSD as a solver allows us to access bath sizes sufficient to converge to “in-method” exactness, and we will move to consider stronger-correlation solvers in future work. We will initially compare the energies from the three introduced energy functionals which are derived from reconstructed wave function forms. As a recap, these are E[Ψx]: Introduced in ref (53) and Sec. 3.3, this is the cheapest approach to energy computation of eq 37, requiring no communication between independent additive cluster contributions. E[(γ, Γ)[Ψx]]: Described in Sec. 3.4, this computes the energy from the exact density matrices from the reconstructed wave function partitioned over the clusters. This scales worse than the mean-field for the exact reconstructed Γ, and rapidly becomes the dominant cost in the calculation. It will primarily be used for reference. E[(γ, K*)[Ψx]]: Described in Sec. 5.2, this approximates the cumulant contribution to the 2-RDM via the “in-cluster” approximation, which ensures low-scaling with system size for this nonfactorizable two-body energy contribution. However, the 2-RDM loses strict N-representability. The 1-RDM contributions are still N-representable from the partitioned wave function according to the efficient algorithm of Sec. 5.1. We do not compare the energies from the density matrix derived expectation values of Sec. 2 due to the different constraints on the fragment spaces imposed. For the wave function derived expectation values above, all of these will converge to in-method exactness as the bath space increases as long as the combined fragment spaces span the occupied space, while the density matrix expectation values require the combined fragment spaces to span all degrees of freedom. This allows us to consider simple minimal-size atomic intrinsic atomic orbital fragments, with the only convergence criteria now being the bath size controlled by η.

We first consider the comparison of these energy functionals for the W4–11 test set of 152 small molecules consisting of first- and second-row elements, which covers a broad range of bonding character, with static and dynamic correlated physics present, as well as a mix of high-spin, low-spin, and open-shell systems. In a cc-pVDZ basis, canonical CCSD calculations on the full test set is tractable to obtain a “ground-truth” total energy for all molecules, and we use ROHF/UHF and UCCSD for open-shell cases. Using the information from the corresponding Hartree–Fock calculation, we then fragment each molecule into a minimal basis set of IAOs on each atom, which we use as the fragments of each clusters, which ensures that the total occupied space of the full system is spanned by the fragments (we omit the beryllium dimer test system due to a lack of virtual orbitals in the default IAO basis). The choice of single-atom fragments simplifies their selection, but future work to consider larger (or smaller) disjoint fragmentation causes no conceptual or practical difficulties. To complete the cluster space of each atom, we augment this with the DMET bath space (which must be at most equal in size to the fragment space), and then with the BNOs defined from a given threshold, η. Smaller η values correspond to larger bath spaces, with η > 1 just returning to the smallest (DMET) bath size. Each cluster is then solved independently at the level of CCSD (restricted or unrestricted), and the total energies reconstructed from the three wave function-based energy functionals provided above, with no self-consistency or chemical potential optimization. These results are shown in Figure 6, with their absolute per electron energy errors aggregated for different values of η in Figure 7. It is clear that the bath expansion of these IAO atomic fragments results in a rapid, monotonic, and systematic convergence of all energy estimates to in-method exactness, obtaining agreement with canonical CCSD calculations to within tight errors. The largest of these errors, in both MAE and standard deviation over the different molecules in the test set, arises from the E[Ψx] functional, which is to be largely expected, since this functional does not have any “cross-cluster” contributions to the energy, relying on a sum over independent energy contributions from the fragments. A significant improvement in the rate of convergence can be found from the E[(γ, Γ)[Ψx]] functional, which contains many products of contributions between different cluster solutions in its construction. Furthermore, this energy is rigorously Ψ-derivable, and is therefore expected to exhibit an error which is quadratic in the implicitly reconstructed wave function error. Though the energy is not strictly variational due to the nonvariationality of coupled-cluster theory, we also find it to be fully variational compared to the canonical CCSD in all cases. However, somewhat surprisingly, the E[(γ, K*)[Ψx]] energy functional converges to in-method exactness at the fastest rate, obtaining the least energy error across the test set for each bath truncation. This approach was motivated as an efficient approximation to the E[(γ, Γ)[Ψx]] energy, but due to favorable cancellation of errors between the approximated cumulant energy contribution and the inherent locality approximation in the bath truncation, appears to be the preferred approach (as well as being particularly efficient). For intermediate cluster sizes, the reduction in per-electron MAE can be nearly an order of magnitude compared to the E[Ψx] energy. Due to the approximation in this cumulant, nonvariational total energies can result (compared to the canonical result); however the appearance of these nonvariational results are very rare, and only observed for already tightly converged energies. This is because the bath truncation itself is a variational approximation (since this defines a subspace of the full variational freedom of the cluster), which explains the favorable cancellation of errors and maintenance of generally variational results, even in this approximation which breaks strict N-representability in the 2-RDM.

A key feature of the embedding approach proposed is the applicability to both molecular and periodic systems, enabled by ensuring all steps scale no worse than the initial Hartree–Fock calculation. We therefore also consider the convergence of structural properties of crystalline fcc diamond in an all-electron cc-pVTZ basis and a 5 × 5 × 5 k-point sampling. At 7500 orbitals and 1500 electrons, this is well beyond the capabilities of canonical CCSD for a direct comparison, but we can observe the relative convergence of the different energy estimators and the resulting equation of state. The supercell is split into 250 individual embedded cluster problems, comprising the minimal basis IAOs of a single carbon atom as the fragment (6 orbitals), the DMET bath space (also 6 orbitals), and the BNOs (whose size depends on the η parameter). The size of the resulting clusters varies between 50 and ∼250 orbitals. Due to the translational and rotational symmetry, all embedded problems are symmetry-equivalent, and therefore only a single cluster needs to be solved, with the linear energy functional trivially exploiting this symmetry, while a slightly more involved exploitation of the symmetry to avoid explicitly solving multiple cluster problems in the accumulation of the long-range RDMs for other expectation values is described in Appendix A. We note that exploitation of this symmetry information is not required and only used here for numerical efficiency, with identical results achieved by computing all cluster contributions independently. Results are presented in Figure 8 for the linear energy functional E[Ψx] and the “in-cluster” cumulant functional, E[(γ, K*)[Ψx]], for a range of lattice parameters of the unit cell. The E[(γ, Γ)[Ψx]] functional is omitted due to its unfavorable scaling to these system sizes. The input to our Vayesta code to generate these results is given in the SI, where periodic integrals for the system are generated from PySCF. Changing the lattice vectors leads to an effective change in the basis set coverage at different cell sizes (since the same number of basis functions span different volumes of the supercell, changing the effective completeness of the overall basis). To compensate for this, we correct all energies for this ‘basis set superposition error’ (BSSE) via an adaptation of the counterpoise correction method as described in ref (53). The BSSE per atom is estimated by computing two independent molecular systems, the first for a single carbon atom (30 basis functions) and the second for a single atom in addition to the basis functions of all carbon atoms within a 3 × 3 × 3 supercell (1620 basis functions in total), sufficient to describe the changing basis set coverage for the carbon atom under consideration. The difference between the two energies is the estimated BSSE error, with this correction changing across the shown lattice range by ∼7.2 mEH, which slightly lowers the energy at large lattice parameters relative to the more compressed cells, leading to a small expansion of the equilibrium volume. The equation of state for the two different energy functionals and BNO bath threshold (η) values are shown in Figure 8. These demonstrate the systematic improvement in the total energy as η decreases. However, while the total energy is not fully converged, we can fit these equations of state to a Birch–Murnaghan form about approximately ±6% of the equilibrium value in order to estimate structural properties for each of these curves and gauge their convergence. We consider the convergence of the equilibrium lattice parameter and bulk modulus in Figure 9, showing convergence to experimental values, corrected for zero-point vibrational effects. We find a small yet consistently improved convergence for both of these properties from the in-cluster approximated cumulant energy functional at each cluster size compared to the linear energy functional derived properties, in keeping with the conclusions from the convergence of the energy in the molecular test set of Figure 7. Our most accurate result (η = 10–7.5, resulting from a single CCSD calculation on ∼330 orbitals), gives a bulk modulus which agrees with experiment by less than 10 GPa, while the equilibrium lattice parameter is in error by only ∼0.002 Å. By way of comparison, literature values from commonly used exchange-correlation functionals of DFT for structural properties give equilibrium lattice constants in error by between +0.02 Å (PBE) and −0.015 Å (HSESol) compared to experimental values, an order of magnitude greater than these results.

While accurate energetics are important, they are far from the only nonlocal expectation values of interest, and we should also gauge the validity of the cumulant-approximated 2-RDM approach for other two-body expectation values. We therefore evaluate instantaneous spin–spin correlation functions between pairs of atoms (A, B), given bywhere PA is a projector onto an appropriately chosen subspace representing the atom A. We define these atomic projectors via symmetrically (Löwdin) orthogonalized atomic orbitals (SAO) for each atom, as given in Appendix B. Note that while the total spin–spin correlations of the full system, , are necessarily zero in a spin-restricted formalism (though not in unrestricted), the local contribution from a given atom pair is generally nonzero, indicating the conditional spin-density between different points in the system. We adapt the fragmentation and BNO bath expansion procedure as detailed in Sec. 6.1 to use a spin-broken UHF reference, spin-dependent bath orbitals, a UCCSD solver, and generalizations of all expectation value accumulation in a unrestricted formalism. We can then consider the convergence of the spin–spin correlation function with bath size for the n-propyl radical (shown in Figure 10) in the cc-pVTZ basis between the radical C1 position (atom A) and other atoms in the system (atom B). In order to compare results from the traditional democratically partitioned 2-RDM of Sec. 2.1 (Γ[Γx]), the democratically partitioned cumulant of Sec. 2.2 (Γ[γx, Kx]) and the in-cluster cumulant approximated 2-RDM derived from the partitioned wave function of Sec. 5.2 (Γ[(γ, K*)[Ψx]]), we choose a complete fragmentation of the system via atomic IAO ⊕ PAO fragments, ensuring that the full virtual space of the system is also spanned by the fragmentation. This is a necessary condition for the democratically partitioned approaches to approach exactness as the bath space becomes complete, but it is a sufficient but not necessary condition for Γ[(γ, K*)[Ψx]], which only requires fragments to span the occupied space. However, choosing this common fragmentation will allow for results to be compared on the same footing, despite requiring a larger fragment space than would be most efficient for the wave function partitioned RDMs. In Appendix A, we detail how an efficient in-cluster approximated form for Γ[(γ, K*)[Ψx]] can be utilized for this spin–spin correlation function, ensuring that all four-index contractions are performed within each cluster, to maintain cubic scaling with system size. Figure 11 plots the convergence of these spin correlation functions from the various functionals as the average cluster (bath) size increases. Beyond just convergence of any single expectation value, these results also allow us to consider the accuracy of these two-body spin-correlators as a function of their range, indicating the ability of each to compensate for the locality approximation inherent in the embedding for smaller cluster sizes, and converge highly nonlocal two-point observables. While there is only a small difference between the traditional partitioning of the 2-RDM and cumulant partitioning schemes for the on-site correlation function (top left plot), this changes for the longer-ranged correlations in the lower rows of the figure, where the partitioned cumulant approach of Sec. 2.2 converges significantly quicker and more smoothly to the full system UCCSD limit, where the direct democratic partitioning of the 2-RDM can become highly erratic. It is particularly noteworthy that the democratically partitioned 2-RDM results can be worse than then uncorrelated UHF values, even when going to average cluster sizes as large as ≈100, at which point the partitioned cumulant results are almost indistinguishable from full-system CCSD values. However, as expected, the wave function derived expectation values are by far the most reliable across the range of cluster sizes, even with the “in-cluster” approximated cumulant, with an improved and more systematic convergence even for the fully local spin correlator. As such, our results support the use of this approach for both energetics and other one- and two-body properties where possible in this wave function embedding context.

In this work, we have critically assessed approaches to reconstruct expectation values from wave function quantum embedding methods. In doing so, we have motivated and compared a number of approaches for calculating total energies across the fragmented system, as well as considering the two-point spin correlation function as an example nonenergetic and nonlocal quantity via the reconstructed reduced density matrices. A fundamental difference exists between the two classes of approaches, where expectation values are computed from either partitioned density matrices, or partitioned wave functions, particularly due to the difference these approaches impose on the choice of fragmentation if a bath expansion of the cluster is to tend to exactness. In the former case, it was found that partitioning the cumulant rather than the 2-RDM directly gave rise to significantly improved two-body properties and energies, and it is likely to find applications in solvers where access to wave function amplitudes is difficult (e.g., the increasing use of quantum computing algorithms as high-level solvers). In the latter case where an implicit partitioned wave function is considered, total energy expressions are motivated and derived from a simple linear energy functional, and from density matrices constructed directly from this wave function. If obtained via these RDMs, the total energy from the embedding is variational (for a variational solver) and N-representable, but incurs significant overhead in the 2-RDM construction. A rigorous and efficient quadratically scaling algorithm is developed for this N-representable 1-RDM, and an “in-cluster” approximation for the 2-RDM is motivated, neglecting many of the cross-cluster contributions in the 2-cumulant. Taken together, these two developments are found to provide the most efficient convergence of energetics and nonlocal expectation values across a range of systems via a systematic bath expansion of the cluster, while only requiring an atomic IAO fragmentation of the system. We believe that this approach will prove important going forward in the development of robust and accurate wave function embedding techniques for both chemical and extended systems. There are a number of further directions we aim to explore as a result of these insights. In systems where the electronic state is required to qualitatively change (e.g., quantum phase transitions), it is possible that a brute-force bath expansion (while necessarily still improvable to exactness) is not going to be the most efficient route in order to capture this changing physics, and feedback from the correlated state to the underlying mean-field reference is required. The identification of an implicit wave function, as well as a formally N-representable 1-RDM opens new avenues for robust and accurate self-consistency conditions. For instance, the self-consistency could maximize the overlap between this implicit partitioned wave function and a mean-field state by minimizing the partitioned T1 amplitudes. Alternatively, coupling between wave function amplitudes of different cluster are also being formulated. Additionally, the approach of “projected”-DMET could be reconsidered with a full-system 1-RDM which contains more of the nonlocal physics and is fully N-representable compared to the standard democratic partitioning. Finally, the application of a broader range of solvers which can treat significantly strongly correlated systems, as well as excited states, spectra, and analytic gradients or forces with these new perspectives of an implicit embedded wave function are ongoing avenues of investigation.