# On the effective reconstruction of expectation values from ab initio quantum
embedding

Max Nusspickel  Department of Physics, King’s College London, Strand, London
WC2R 2LS, U.K.  Basil Ibrahim  Department of Physics, King’s College London,
Strand, London WC2R 2LS, U.K.  George H. Booth<EMAIL_ADDRESS>Department of Physics, King’s College London, Strand, London WC2R 2LS, U.K.

###### Abstract

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 non-
local 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 non-local 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.

##  I Introduction

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 [1, 2, 3, 4, 5, 6]. 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, whilst
remaining free from empiricism [7, 8]. 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 [9, 10, 11]. While highly related, ‘quantum embedding’
approaches from condensed matter physics are also increasingly coming to the
fore as an alternative paradigm and applied to quantum chemical and ab initio
systems [12].

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 or unoccupied mean-field character), which we will call
the ‘fragment’ space, though 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, e) optionally, a self-consistency is
performed to embed the correlated effects from the cluster model back into the
low-level full system method to update the coupling between fragment 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) [13, 14, 15], how the
interacting cluster Hamiltonian is constructed and solved [16, 17, 18, 19],
and the choice of self-consistent requirements [20, 21, 22]. 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 [23, 24, 25, 26,
27, 28, 29], while density matrix embedding theory (DMET) and its variants
work in a wave function (static) formalism [30, 31, 32, 33, 34, 35, 36, 37,
38, 39, 40, 41, 42, 43, 44], although these two approaches are not
fundamentally distinct, and can also be rigorously connected via a common
framework [45, 46, 47].

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 [48, 49, 50, 51, 52]. 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 non-local 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 non-local 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 [30, 33]. 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 and ‘black-box’ approach, and avoids the ambiguities and non-
monotonic 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 [53].

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 [21, 20], but 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 recapitulating the original DMET ‘democratic partitioning’
for expectation values, which can be computed via fragment-projection of the
reduced density matrices of each fragment [33]. 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 non-local properties (spin correlation functions) on both
molecular systems via the W4-11 test set [54], and periodic systems, finding
an efficient approach for reconstruction of non-local observables in quantum
embedding of general fragmented systems.

###  I.1 Summary of findings

Figure 1: Binding energy of the Chlorine dimer in a STO-6G basis compared to
FCI, with $E[\gamma^{{\mathbf{x}}},\Gamma^{{\mathbf{x}}}]$ corresponding to
the democratically partitioned energy expression of traditional DMET [33].
Different energies correspond to different total energy functionals,
constructed from the same solutions to the two DMET embedding problems (with
fragment spaces of the (Löwdin-orthogonalized) core and valence orbitals of
each Chlorine atom). The two cluster spaces are defined by 10 orbitals (9
fragment and 1 bath), with the full space comprising 18 orbitals. The
motivation and definition of these different energy functionals from the
embedded wave functions are given in the rest of this work, with more details
on this system (and others like it) discussed in section IV.

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 Fig. 1, the same identical fragment and DMET cluster
space is used, and the same (exact) solver and wave function are found from
each cluster [30, 31]. The only difference is the choice of energy functional
to reconstruct the total energy from the two cluster solutions, with
$E[\gamma^{\mathbf{x}},\Gamma^{\mathbf{x}}]$ representing the traditional
‘democratically partitioned’ DMET energy used primarily in the literature to
date [33].

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. 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 non-local
fashion, and minimizing the local approximations inherent to the embedding
framework.

  2. 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 [40]. Beyond this, enabling convergence with
respect to bath size of each fragment can render self-consistency entirely
unnecessary.

  3. 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. 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. 5.

If the approach above results in an intractable scaling with respect to 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 $\mathcal{O}(N^{3})$ 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
111Embedding code and documentation can be found at
https://github.com/BoothGroup/Vayesta, which interfaces with the PySCF code
[56, 57], with scripts to generate many of the results of this work (including
the input required for Fig. 1) found in the Supplementary Information (SI).

##  II Global expectation values from cluster density matrices

###  II.1 Democratic partitioning of density matrices

Expectation values in DMET derived from operators which span more than one
fragment are calculated from ‘democratically partitioned’ reduced density-
matrices (RDMs) [33]. These can be written as

$\displaystyle\gamma_{pq}$
$\displaystyle=\sum_{\mathbf{x}}^{N_{\mathrm{frag}}}(\hat{P}^{\mathbf{x}}\gamma^{\mathbf{x}})_{pq}$
(1) $\displaystyle\Gamma_{pqrs}$
$\displaystyle=\sum_{\mathbf{x}}^{N_{\mathrm{frag}}}(\hat{P}^{\mathbf{x}}\Gamma^{\mathbf{x}})_{pqrs},$
(2)

where $\gamma$ and $\Gamma$ are one- and two-body reduced density matrices
formed from the high-level solution of each cluster problem, as

$\displaystyle\gamma_{pq}^{{\mathbf{x}}}$
$\displaystyle=\braket{\Psi^{{\mathbf{x}}}}{\hat{c}_{p}^{\dagger}\hat{c}_{q}}{\Psi^{{\mathbf{x}}}},$
(3) $\displaystyle\Gamma_{pqrs}^{{\mathbf{x}}}$
$\displaystyle=\braket{\Psi^{{\mathbf{x}}}}{\hat{c}_{p}^{\dagger}\hat{c}_{r}^{\dagger}\hat{c}_{s}\hat{c}_{q}}{\Psi^{{\mathbf{x}}}},$
(4)

where the high-level cluster wave function, $|\Psi^{{\mathbf{x}}}\rangle$,
includes the contribution from the doubly occupied environmental orbitals of
each cluster. $\hat{P}^{\mathbf{x}}$ is an operator, which introduces a
projection onto the fragment subspace of the cluster ${\mathbf{x}}$ in a
symmetrically averaged fashion, e.g.

$\displaystyle(\hat{P}^{\mathbf{x}}\gamma^{\mathbf{x}})_{pq}$
$\displaystyle=\frac{1}{2}\sum_{t}^{N_{\mathrm{mo}}}(P_{pt}^{\mathbf{x}}\gamma^{\mathbf{x}}_{tq}+P_{qt}^{\mathbf{x}}\gamma^{\mathbf{x}}_{pt}),$
(5) $\displaystyle(\hat{P}^{\mathbf{x}}\Gamma^{\mathbf{x}})_{pqrs}$
$\displaystyle=\frac{1}{4}\sum_{t}^{N_{\mathrm{mo}}}(P_{pt}^{\mathbf{x}}\Gamma^{\mathbf{x}}_{tqrs}+P_{qt}^{\mathbf{x}}\Gamma^{\mathbf{x}}_{ptrs}$
(6)
$\displaystyle\qquad\qquad+P_{rt}^{\mathbf{x}}\Gamma^{\mathbf{x}}_{pqts}+P_{st}^{\mathbf{x}}\Gamma^{\mathbf{x}}_{pqrt}).$

In these, ${\mathbf{x}}$ labels both the $N_{\mathrm{frag}}$ fragments and
clusters, for which there is an unambiguous one-to-one correspondence, and
$P_{pq}^{\mathbf{x}}$ is the fragment projection matrix acting over the whole
molecular orbital (MO) space, as

$P^{{\mathbf{x}}}_{pq}=\sum_{x}^{N_{\mathrm{cl}}^{\mathbf{x}}}\left(\mathbf{C}^{T}\mathbf{S}\mathbf{C}_{\mathrm{f}}^{\mathbf{x}}\right)_{px}\left(\mathbf{C}^{T}\mathbf{S}\mathbf{C}_{\mathrm{f}}^{{\mathbf{x}}}\right)_{qx},$
(7)

with $\mathbf{C}$ representing the MO coefficients, $\mathbf{S}$ the atomic
orbital (AO) overlap matrix, and $\mathbf{C}_{\mathrm{f}}^{\mathbf{x}}$ the
coefficients of the fragment orbitals in cluster ${\mathbf{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 as

$E_{\mathrm{tot}}=E_{\mathrm{nuc}}+\sum_{pq}^{N_{\mathrm{mo}}}h_{pq}\gamma_{qp}+\frac{1}{2}\sum_{pqrs}^{N_{\mathrm{mo}}}(pq|rs)\Gamma_{pqrs}.$
(8)

Note 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 [33], equivalent to Eq. (8) via construction of democratically
partitioned density matrices of Eq. (5), as

$\begin{split}E[\gamma^{\mathbf{x}},\Gamma^{\mathbf{x}}]=E_{\mathrm{nuc}}+&\sum_{{\mathbf{x}}}^{N_{\mathrm{frag}}}\left[\sum_{pq}^{N_{\mathrm{cl}}^{\mathbf{x}}}\tilde{h}^{{\mathbf{x}}}_{pq}(\hat{P}^{\mathbf{x}}\gamma^{\mathbf{x}})_{qp}\right.\\\
&+\left.\frac{1}{2}\sum_{pqrs}^{N_{\mathrm{cl}}^{\mathbf{x}}}(pq|rs)(\hat{P}^{\mathbf{x}}\Gamma^{\mathbf{x}})_{pqrs}\right],\end{split}$
(9)

where $\gamma^{\mathbf{x}}$ and $\Gamma^{\mathbf{x}}$ refer to the cluster
reduced density matrices, with the projector purely acting in this cluster
space, and $N_{\mathrm{cl}}^{{\mathbf{x}}}$ denoting the number of orbitals in
the cluster ${\mathbf{x}}$. The energetic effect of these states (static
Coulomb and exchange contributions) is then included via the construction of
the $\tilde{h}^{{\mathbf{x}}}$, which includes the potential from these
unentangled states to the one-body hamiltonian as
$\frac{1}{2}\sum_{mn}(pq||mn)\gamma^{{\mathbf{x}}}_{\mathrm{core},mn}$ where
$\gamma^{{\mathbf{x}}}_{\mathrm{core}}$ 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[\gamma^{\mathbf{x}},\Gamma^{\mathbf{x}}]$ will denote this democratically
partitioned energy functional, shorthand for
$E[\\{\gamma^{\mathbf{x}}\\},\\{\Gamma^{\mathbf{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, that

$\sum_{{\mathbf{x}}}^{N_{\textrm{frag}}}{\hat{P}}^{{\mathbf{x}}}=\hat{\mathds{1}}.$
(10)

While 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)[58] to ensure this condition is fulfilled for reasonable results [35].
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
$|\Psi^{\mathbf{x}}\rangle\rightarrow|\Psi\rangle$ 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. Firstly, 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),

$\gamma_{pq}=\sum_{\mathbf{x}}^{N_{\mathrm{frag}}}\hat{P}^{\mathbf{x}}\braket{\Psi^{\mathbf{x}}}{c_{p}^{\dagger}c_{q}}{\Psi^{\mathbf{x}}},$
(11)

cannot be rewritten as a simple expectation value
$\braket{\tilde{\Psi}}{c_{p}^{\dagger}c_{q}}{\tilde{\Psi}}$ of some wave
function $\ket{\tilde{\Psi}}$. 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 [40]. Furthermore, conserved quantities and
good quantum numbers in the individual cluster solutions such as electron
number ($N$), spin and its $z$-projection ($S_{z}$, $S^{2}$) 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 [33, 59]. 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. IV (with and without a global chemical potential
optimization), also showing the deleterious effect on properties and computed
energetics of the system that result.

###  II.2 Democratic partitioning of cumulants

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, ${\tilde{K}}$, defined (in a restricted basis)222We note that the
two-body cumulant would often be denoted by $\lambda$, but define by
${\tilde{K}}$ to avoid potential confusions with the lambda-amplitudes of
coupled-cluster theory used later in the text. via

$\Gamma_{pqrs}=\gamma_{pq}\gamma_{rs}-\frac{1}{2}\gamma_{pr}\gamma_{sq}+{\tilde{K}}_{pqrs}.$
(12)

The non-cumulant (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 by

$\Gamma_{pqrs}=\gamma_{pq}\gamma_{rs}-\frac{1}{2}\gamma_{pr}\gamma_{sq}+\sum_{\mathbf{x}}^{N_{\mathrm{frag}}}(\hat{P}^{\mathbf{x}}{\tilde{K}}^{\mathbf{x}})_{pqrs}.$
(13)

The difference between Eq. (2) and (13) lies purely in the non-cumulant
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):

$\gamma_{pq}\gamma_{rs}=\sum_{\mathbf{x}}^{N_{\mathrm{frag}}}\sum_{\mathbf{y}}^{N_{\mathrm{frag}}}(\hat{P}^{\mathbf{x}}\gamma^{\mathbf{x}})_{pq}(\hat{P}^{\mathbf{y}}\gamma^{\mathbf{y}})_{rs}.$
(14)

In this way, the non-local (correlated) one-body physics of two distinct
clusters, ${\mathbf{x}}\neq{\mathbf{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. II.1 for non-local expectation values. We will
denote any total energies resulting from the democratically partitioned
cumulant approach described in this section as
$E[\gamma^{{\mathbf{x}}},K^{{\mathbf{x}}}]$ in the rest of this work.

Figure 2: Dissociation curve of N2 in the STO-6G basis, calculated from two
atomic DMET embedding problems. The same DMET-FCI cluster solution is used in
both the partitioned density matrix (Sec. II.1) and partitioned cumulant (Sec.
II.2) energy expressions, $E[\gamma^{{\mathbf{x}}},\Gamma^{{\mathbf{x}}}]$ and
$E[\gamma^{{\mathbf{x}}},K^{{\mathbf{x}}}]$ respectively, defined as given by
Eqs. (9) and (18). Plots of the right have optimized a global chemical
potential to ensure that the democratically partitioned 1RDM traces to the
correct number of electrons, while the plots on the left omit this chemical
potential optimization.

We present a simple example showing the difference between the partitioned
density matrices and the partitioned cumulant approach in Fig. 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 Fig. 1 and analyzed in more detail in
Sec. IV, while the improvement in non-energetic properties (the spin–spin
correlation function) from the cumulant partitioning will be shown in Sec.
VI.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
$\gamma=\gamma^{0}+\Delta\gamma$, where $\gamma^{0}$ is the reference mean-
field density-matrix. We then work with a slightly modified cumulant
definition, ${K}$, as

$\begin{split}\Gamma_{pqrs}&=\gamma_{pq}^{0}\gamma_{rs}^{0}+\gamma_{pq}^{0}\Delta\gamma_{rs}+\Delta\gamma_{pq}\gamma_{rs}^{0}\\\
&-\frac{1}{2}(\gamma_{pr}^{0}\gamma_{sq}^{0}+\gamma_{pr}^{0}\Delta\gamma_{sq}+\Delta\gamma_{pr}\gamma_{sq}^{0})+{K}_{pqrs}.\end{split}$
(15)

The approximate cumulant ${K}$ of Eq. (15) therefore contains the
$(\Delta\gamma)^{2}$-terms which are not present in the true cumulant
${\tilde{K}}$, with the relation between the two given as

${K}_{pqrs}={\tilde{K}}_{pqrs}+\Delta\gamma_{pq}\Delta\gamma_{rs}-\frac{1}{2}\Delta\gamma_{pr}\Delta\gamma_{sq}.$
(16)

The final expression for the 2-RDM is then

$\begin{split}\Gamma_{pqrs}&=\gamma_{pq}^{0}\gamma_{rs}^{0}+\gamma_{pq}^{0}\Delta\gamma_{rs}+\Delta\gamma_{pq}\gamma_{rs}^{0}\\\
&-\frac{1}{2}(\gamma_{pr}^{0}\gamma_{sq}^{0}+\gamma_{pr}^{0}\Delta\gamma_{sq}+\Delta\gamma_{pr}\gamma_{sq}^{0})\\\
&+\sum_{\mathbf{x}}^{N_{\mathrm{frag}}}(\hat{P}^{\mathbf{x}}{K}^{\mathbf{x}})_{pqrs},\end{split}$
(17)

where (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 non-zero cumulant in this
space. This allows the total energy functional in terms of
$\gamma^{\mathbf{x}}$ and ${K}^{\mathbf{x}}$ to be written as

$\begin{split}E[\gamma^{\mathbf{x}},{K}^{\mathbf{x}}]&=E_{\mathrm{HF}}[\gamma^{0}]+\sum_{pq}F_{pq}[\gamma^{0}]\Delta\gamma_{qp}\\\
&+\frac{1}{2}\sum_{{\mathbf{x}}}^{N_{\mathrm{frag}}}\sum_{pqrs}^{N_{\mathrm{cl}}^{\mathbf{x}}}(pq|rs)(\hat{P}^{\mathbf{x}}{{K}}^{\mathbf{x}})_{pqrs},\end{split}$
(18)

where $E_{\mathrm{HF}}[\gamma^{0}]$ and $F[\gamma^{0}]$ are the Hartree–Fock
energy and Fock matrix corresponding to the reference mean-field density-
matrix $\gamma^{0}$, and $\Delta\gamma$ is the correlated part of the
democratically-partitioned 1-RDM, formed from the cluster density matrices,
$\gamma^{{\mathbf{x}}}$, as shown in Eq. (1). Note that Eq. (18), in contrast
to the conventional DMET energy functional (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[\gamma^{0}+\Delta\gamma]$, would need to be constructed
to calculate the energy, a step that scales as $\mathcal{O}(N^{4})$ 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[\gamma^{{\mathbf{x}}},K^{{\mathbf{x}}}]$ in the
results of this paper.

##  III Global expectation values from cluster wave functions

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.
II.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. III.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, as

$\ket{\Psi}=\sum_{\mathbf{x}}^{N_{\mathrm{frag}}}\hat{P}^{\mathbf{x}}\ket{\Psi^{\mathbf{x}}}.$
(19)

From 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 as

$\gamma^{(\text{part.
WF})}_{pq}=\sum_{\mathbf{xy}}^{N_{\mathrm{frag}}}\braket{\hat{P}^{\mathbf{x}}\Psi^{\mathbf{x}}}{\hat{c}_{p}^{\dagger}\hat{c}_{q}}{\hat{P}^{\mathbf{y}}\Psi^{\mathbf{y}}}$
(20)

As long as all projected cluster wave functions
$\hat{P}^{\mathbf{x}}\ket{\Psi^{\mathbf{x}}}$ are physical, fermionic wave
function 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
$\braket{\Psi}{c_{p}^{\dagger}c_{q}}{\Psi}$ with the “global” wave function
$\ket{\Psi}$ 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).

###  III.1 Partitioning of general wave functions

So far we have not specified what we mean by the projection of a cluster wave
function as $\hat{P}^{\mathbf{x}}\ket{\Psi^{\mathbf{x}}}$ in practice. Since
the DMET bath orbitals guarantee that the mean-field reference determinant
$\ket{\Phi}$ is represented exactly in each cluster (ignoring the unentangled
environment orbitals) [30], it is convenient to define $\hat{P}^{\mathbf{x}}$
in terms of its action on the correlated part of the wave function only, i.e.
our full system wave function can be written as

$\ket{\Psi}=\ket{\Phi}+\sum_{\mathbf{x}}^{N_{\mathrm{frag}}}\hat{P}^{\mathbf{x}}\ket{\Delta\Psi^{\mathbf{x}}},$
(21)

where $\ket{\Delta\Psi^{\mathbf{x}}}=\ket{\Psi^{\mathbf{x}}}-\ket{\Phi}$ and
intermediate normalization is assumed.

We then choose to represent the correlated cluster wave function part
$\ket{\Delta\Psi}$ in the basis of particle–hole excitations around the
reference determinant in a linear wave function ansatz

$\begin{split}\ket{\Delta\Psi}&=(\hat{C}_{1}+\hat{C}_{2}+\dots)\ket{\Phi}\\\
&=\sum_{i}\sum_{a}c_{i}^{a}\ket{\Phi_{i}^{a}}+\sum_{ij}\sum_{ab}c_{ij}^{ab}\ket{\Phi_{ij}^{ab}}+\dots,\end{split}$
(22)

where we use $i,j,\dots$ ($a,b,\dots$) 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
example

$\displaystyle(\hat{P}^{\mathbf{x}}C_{1}^{\mathbf{x}})_{i}^{a}=\sum_{k}P^{\mathbf{x}}_{ik}(C^{\mathbf{x}}_{1})_{k}^{a},$
(23)
$\displaystyle(\hat{P}^{\mathbf{x}}C_{2}^{\mathbf{x}})_{ij}^{ab}=\frac{1}{2}\sum_{k}\left[P^{\mathbf{x}}_{ik}(C^{\mathbf{x}}_{2})_{kj}^{ab}+P^{\mathbf{x}}_{jk}(C^{\mathbf{x}}_{2})_{ik}^{ab}\right],$
(24)

in the case of the single and double excitation coefficients in the cluster,
$C^{{\mathbf{x}}}_{1}$ and $C^{{\mathbf{x}}}_{2}$, 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 as

$P^{{\mathbf{x}}}_{ik}=\sum_{x\in{\mathbf{x}}}\left[(\mathbf{C}^{{\mathbf{x}}})^{T}\mathbf{S}\mathbf{C}_{\mathrm{f}}^{\mathbf{x}}\right]_{ix}\left[(\mathbf{C}^{{\mathbf{x}}})^{T}\mathbf{S}\mathbf{C}_{\mathrm{f}}^{{\mathbf{x}}}\right]_{kx}$
(25)

where $\mathbf{C}_{\mathrm{f}}^{{\mathbf{x}}}$ represents the AO coefficients
of the fragment orbitals ($x$) of fragment ${\mathbf{x}}$, and
$\mathbf{C}^{{\mathbf{x}}}$ 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
${\mathbf{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, 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 as

$\mathrm{Tr}\left[\sum_{{\mathbf{x}}}^{N_{\textrm{frag}}}{\hat{P}}_{{\mathbf{x}}}\gamma^{0}{\hat{P}}_{{\mathbf{x}}}\right]=N_{e},$
(26)

where $N_{e}$ 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[61], 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 correlation, 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 fulfil. However,
for ab initio quantum embeddings, we do not anticipate that this would be
beneficial.

###  III.2 Partitioning of exponential wave function forms

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
parameterization for partitioning in this space. This form (common to coupled-
cluster methods[62]) can be written as

$\begin{split}\ket{\Delta\Psi}&=(\mathrm{e}^{\hat{T}_{1}+\hat{T}_{2}+\dots}-1)\ket{\Phi}=\sum_{i}^{N_{\mathrm{occ}}}\sum_{a}^{N_{\mathrm{vir}}}t_{i}^{a}\ket{\Phi_{i}^{a}}\\\
&+\sum_{ij}^{N_{\mathrm{occ}}}\sum_{ab}^{N_{\mathrm{vir}}}\left(t_{ij}^{ab}+t_{i}^{a}t_{j}^{b}\right)\ket{\Phi_{ij}^{ab}}+\dots\end{split}$
(27)

The amplitudes can be converted between the linear and exponential forms more
generally [63] via the relations

$\displaystyle\hat{C}_{1}$ $\displaystyle=\hat{T}_{1},$ (28)
$\displaystyle\hat{C}_{2}$
$\displaystyle=\hat{T}_{2}+\frac{(\hat{T}_{1})^{2}}{2},$ (29)
$\displaystyle\hat{C}_{3}$
$\displaystyle=\hat{T}_{3}+\frac{(\hat{T}_{2}\hat{T}_{1}+\hat{T}_{1}\hat{T}_{2})}{2}+\frac{(\hat{T}_{1})^{3}}{3!},$
(30) $\displaystyle\dots.$

For 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, 24) can be written as

$\displaystyle(\hat{P}^{\mathbf{x}}T_{1}^{\mathbf{x}})_{i}^{a}=\sum_{k}^{N_{\mathrm{occ}}^{\mathbf{x}}}P^{\mathbf{x}}_{ik}(T^{\mathbf{x}}_{1})_{k}^{a},$
(31)
$\displaystyle(\hat{P}^{\mathbf{x}}T_{2}^{\mathbf{x}})_{ij}^{ab}=\frac{1}{2}\sum_{k}^{N_{\mathrm{occ}}^{\mathbf{x}}}\left[P^{\mathbf{x}}_{ik}(T^{\mathbf{x}}_{2})_{kj}^{ab}+P^{\mathbf{x}}_{jk}(T^{\mathbf{x}}_{2})_{ik}^{ab}\right],$
(32)

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. III.1), we
obtain

$\begin{split}\ket{\Delta\Psi}=&\sum_{\mathbf{x}}^{N_{\mathrm{frag}}}\sum_{i}^{N_{\mathrm{occ}}^{\mathbf{x}}}\sum_{a}^{N_{\mathrm{vir}}^{\mathbf{x}}}(\hat{P}^{\mathbf{x}}C^{\mathbf{x}}_{1})_{i}^{a}\ket{\Phi_{i}^{a}}\\\
+&\sum_{\mathbf{x}}^{N_{\mathrm{frag}}}\sum_{ij}^{N_{\mathrm{occ}}^{\mathbf{x}}}\sum_{ab}^{N_{\mathrm{vir}}^{\mathbf{x}}}(\hat{P}^{\mathbf{x}}C^{\mathbf{x}}_{2})_{ij}^{ab}\ket{\Phi_{ij}^{ab}}\end{split}$
(33)

while for the exponential form of the wave function, we achieve

$\begin{split}\ket{\Delta\Psi}=&\sum_{\mathbf{x}}^{N_{\mathrm{frag}}}\sum_{i}^{N_{\mathrm{occ}}^{\mathbf{x}}}\sum_{a}^{N_{\mathrm{vir}}^{\mathbf{x}}}(\hat{P}^{\mathbf{x}}T^{\mathbf{x}}_{1})_{i}^{a}\ket{\Phi_{i}^{a}}\\\
+&\sum_{\mathbf{x}}^{N_{\mathrm{frag}}}\sum_{ij}^{N_{\mathrm{occ}}^{\mathbf{x}}}\sum_{ab}^{N_{\mathrm{occ}}^{\mathbf{x}}}(\hat{P}^{\mathbf{x}}T^{\mathbf{x}}_{2})_{ij}^{ab}\ket{\Phi_{ij}^{ab}}\\\
+&\sum_{{\mathbf{x}}{\mathbf{y}}}^{N_{\mathrm{frag}}}\sum_{i}^{N_{\mathrm{occ}}^{\mathbf{x}}}\sum_{a}^{N_{\mathrm{vir}}^{\mathbf{x}}}\sum_{j}^{N_{\mathrm{occ}}^{\mathbf{y}}}\sum_{b}^{N_{\mathrm{vir}}^{\mathbf{y}}}(\hat{P}^{\mathbf{x}}T^{\mathbf{x}}_{1})_{i}^{a}(\hat{P}^{\mathbf{y}}T^{\mathbf{y}}_{1})_{j}^{b}\ket{\Phi_{ij}^{ab}}.\end{split}$
(34)

While the first two terms in Eqs. (33, 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 ${\mathbf{x}}$ and a different single
excitation in cluster ${\mathbf{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
$T_{2}^{\mathbf{x}}T_{1}^{\mathbf{y}}$), or three independent clusters (e.g.
in $T_{1}^{\mathbf{x}}T_{1}^{\mathbf{y}}T_{1}^{\mathbf{z}}$) 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 non-local ‘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. II.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 a 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 also particularly suited to coupled-cluster high-level
solvers, which are increasingly being used as a compact and accurate approach
which enable access to larger fragment and bath spaces, and have recently been
shown to be accurate in a variety of embedding contexts [64, 65, 66, 67, 53].

###  III.3 Expectation values from linear functionals

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 which 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
$\langle\Phi|$ and assuming intermediate normalization
$\braket{\Phi}{\Psi}=1$, it follows that

$E[\Psi]=\braket{\Phi}{\hat{H}}{\Psi},$ (35)

which is linear in $\Psi$. 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 energy

$E_{\mathrm{corr}}=E-E_{\mathrm{HF}}=\sum_{ij}^{N_{\mathrm{occ}}}\sum_{ab}^{N_{\mathrm{vir}}}c_{ij}^{ab}\left[2(ia|jb)-(ib|ja)\right],$
(36)

where 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
$C_{2}$-amplitudes in Eq. (36) according to Eq. (24) enables us to partition
the correlation energy into cluster contributions as

$\begin{split}E[\Psi^{{\mathbf{x}}}]&=E_{\mathrm{HF}}+\sum_{\mathbf{x}}^{N_{\mathrm{frag}}}E_{\mathrm{corr}}^{\mathbf{x}}\\\
&=E_{\mathrm{HF}}+\sum_{{\mathbf{x}}}^{N_{\mathrm{frag}}}\sum_{ij}^{N_{\mathrm{occ}}^{\mathbf{x}}}\sum_{ab}^{N_{\mathrm{vir}}^{\mathbf{x}}}(\hat{P}^{\mathbf{x}}C^{\mathbf{x}}_{2})_{ij}^{ab}\left[2(ia|jb)-(ib|ja)\right].\end{split}$
(37)

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 will be denoted by $E[\Psi^{{\mathbf{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 [68] or cluster-in-molecules
approaches [69]. 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 (non-
diagonal) 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. III.2). In this case, inter-cluster contributions arise
from
$(\hat{P}^{{\mathbf{x}}}T_{1}^{{\mathbf{x}}})(\hat{P}^{{\mathbf{y}}}T_{1}^{{\mathbf{y}}})$
terms. This physics is neglected in the $C$-amplitude projection of above in
the case of no overlap between clusters ${\mathbf{x}}$ and ${\mathbf{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 $T_{1}^{2}$ 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[\Psi^{{\mathbf{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.

###  III.4 Partitioned wave function density matrices

General expectation values, with operators that do not commute with the
Hamiltonian, require a quadratic functional of the wave function, of the form
$\braket{\Psi}{\hat{O}}{\Psi}$. In the case of one- and two-body expectation
values, $O_{1}$ and $O_{2}$, 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) as

$\displaystyle O_{1}[\Psi]$
$\displaystyle=\braket{\Psi}{\hat{O}_{1}}{\Psi}=\mathrm{Tr}\>[\gamma_{pq}\braket{p}{\hat{O}_{1}}{q}],$
(38) $\displaystyle O_{2}[\Psi]$
$\displaystyle=\braket{\Psi}{\hat{O}_{2}}{\Psi}=\mathrm{Tr}\>[\Gamma_{pqrs}\braket{pr}{\hat{O}_{2}}{qs}].$
(39)

The 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
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 III.1 and III.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 as

$\displaystyle\gamma_{pq}$ $\displaystyle=$
$\displaystyle\braket{\Psi}{\hat{c}_{p}^{\dagger}\hat{c}_{q}}{\Psi}=\braket{\Phi+\Delta\Psi}{\hat{c}_{p}^{\dagger}\hat{c}_{q}}{\Phi+\Delta\Psi}$
$\displaystyle=$
$\displaystyle\gamma_{pq}^{0}+\sum_{\mathbf{x}}^{N_{\mathrm{frag}}}\sum_{\lambda\in{\mathbf{x}}}(\hat{P}^{\mathbf{x}}C^{\mathbf{x}})_{\lambda}\braket{\Phi}{\hat{c}_{p}^{\dagger}\hat{c}_{q}}{\Phi_{\lambda}}$
$\displaystyle+$
$\displaystyle\sum_{\mathbf{x}}^{N_{\mathrm{frag}}}\sum_{\lambda\in{\mathbf{x}}}(\hat{P}^{\mathbf{x}}C^{\mathbf{x}})_{\lambda}\braket{\Phi_{\lambda}}{\hat{c}_{p}^{\dagger}\hat{c}_{q}}{\Phi}$
$\displaystyle+$
$\displaystyle\sum_{{\mathbf{x}},{\mathbf{y}}}^{N_{\mathrm{frag}}}\sum_{\lambda\in{\mathbf{x}}}\sum_{\mu\in{\mathbf{y}}}(\hat{P}^{\mathbf{x}}C^{\mathbf{x}})_{\lambda}(\hat{P}^{\mathbf{y}}C^{\mathbf{y}})_{\mu}\braket{\Phi_{\lambda}}{\hat{c}_{p}^{\dagger}\hat{c}_{q}}{\Phi_{\mu}},$
(41)

where $\lambda$ and $\mu$ 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 $T_{2}$, 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
non-local ‘cross-cluster’ contributions and therefore represents a better
approximation. In Sec. V 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[(\gamma,\Gamma)[\Psi^{{\mathbf{x}}}]]$
throughout this work, to signify that these are density matrices derived from
the global partitioned wave function.

##  IV A simple example: Chlorine dimer

Before describing specific implementational details in Sec. V, we consider a
simple example to corroborate the claims in the work so far, before a more
extensive set of test systems is considered in Sec. VI. 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. The DMET cluster space in this example
therefore 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.

Figure 3: Comparison of different proposed DMET energy expressions for the
binding of Cl2 in a STO-6G basis. Each cluster consists of the atomic orbitals
of one atom, with a single DMET bath orbital on the other, with all clusters
solved exactly. Right column results include a fragment chemical potential
optimization to constrain the total number of electrons, while this is omitted
in the left column.

The interacting cluster Hamiltonian is solved exactly via FCI [70], and we
compare the four different energy expressions outlined to far in Fig. 3: the
democratically-partitioned RDM energy common to DMET approaches to date
($E[\gamma^{{\mathbf{x}}},\Gamma^{{\mathbf{x}}}]$ from Sec. II.1), the
democratically partitioned cumulant energy
($E[\gamma^{{\mathbf{x}}},K^{{\mathbf{x}}}]$ from Sec. II.2), the linear
energy functional ($E[\Psi^{{\mathbf{x}}}]$ from Sec. III.3) and the energy
from the RDMs of the $T$-amplitude projected global wave function
($E[(\gamma,\Gamma)[\Psi^{{\mathbf{x}}}]]$ from Sec. III.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 Fig. 3), this fragment chemical potential,
$\mu_{\mathrm{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 [33] (these $\mu_{\mathrm{frag}}$-optimized
results are the same as shown in Fig. 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 Fig. 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
$\mu_{\mathrm{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 [31,
33, 46, 21]). The democratically partitioned cumulant approach is also changed
by the $\mu_{\mathrm{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
non-parallelity error of the democratically partitioned cumulant approach from
$\sim 6\,\mathrm{m}E_{\mathrm{H}}$ to just over $1\,\mathrm{m}E_{\mathrm{H}}$
in this system.

We also show the energies derived from partitioned wave functions over the
system, either calculating the energy from 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 fulfil 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 obviates 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[(\gamma,\Gamma)[\Psi^{{\mathbf{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 non-
parallelity error of 0.5m$E_{H}$.

Figure 4: Analysis of the 1- and 2-RDMs from direct democratic partitioning
(dotted lines) compared to RDMs constructed from partitioning of the wave
function (solid line), for the same Cl2 system as Fig. 3. The results are
shown without (left) and with (right) chemical potential optimization over the
DMET clusters. Top to bottom rows show the errors in the norm, trace and
$N$-representability conditions (Eq. 42) of the 1-RDM for the different
approaches, with the bottom row showing the norm error in the 2-RDM including
the democratically partitioned cumulant approach.

In order to analyze the different approaches in more detail, in Fig. 4 we
compare the 1- and 2-RDMs resulting from the democratically partitioned RDMs
($\gamma[\gamma^{{\mathbf{x}}}]$ for the 1-RDM, and
$\Gamma[\Gamma^{{\mathbf{x}}}]$ or
$\Gamma[\gamma^{{\mathbf{x}}},K^{{\mathbf{x}}}]$ 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
($\gamma[\Psi^{{\mathbf{x}}}]$ and $\Gamma[\Psi^{{\mathbf{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, $\Omega$,
which sums the deviation from allowable occupations of the system orbitals as

$\Omega=\sum_{i:n_{i}<0}|n_{i}|+\sum_{i:n_{i}>2}|n_{i}-2|,$ (42)

where $n_{i}$ are the eigenvalues of $\gamma$.

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 non-cumulant error in the 2-RDM which is contributing to the
significant energy errors in the standard DMET energy functional.

##  V From concepts to practice

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. III.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, will maximize the number of non-local
combinations of cluster solutions, and is a natural choice for coupled-cluster
high-level solvers (although as explained in Sec. III.2, other solvers can be
cast in this form). Furthermore, we focus on an efficient implementation of
the cluster solver truncated at the $T_{2}$-amplitudes, i.e. CCSD (or MP2),
which is a cost-effective high-level cluster solver for ab initio systems [64,
65, 66, 67, 53], yet is sufficiently computationally inexpensive for us to
rigorously demonstrate the convergence of these functionals with respect to
bath size in Sec. VI. Extensions to other high-level cluster solvers (e.g.
FCI) can proceed via Eq. III.4, 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. $T_{n}$ 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 $C_{n}$
variables).

Explicit reconstruction of the global CCSD $T$-amplitudes (which we call the
‘global $T_{2}$ algorithm’) would require a memory overhead scaling as
$\mathcal{O}(N^{4})$ and computational scaling of $\mathcal{O}(N^{5})$ and
$\mathcal{O}(N^{6})$ for the 1- and 2-RDMs respectively (where $N$ is a
measure of full system size) via the usual CCSD equations [71]. We reduce this
scaling down to $\mathcal{O}(N^{2})$ cost in memory and $\mathcal{O}(N^{3})$
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. V.1, which we denote the
‘distributed $T_{2}$ algorithm’, and is crucial to ensure that this step is a
sub-leading 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 necessitate the introduction of a further approximation,
which is described in Sec. V.2, with rigorous validation of this further
approximation demonstrated in Sec. VI.

###  V.1 One-RDM from cluster wave functions

In order to construct density-matrices at the CCSD level, both wave function
$T$-amplitudes and the Lagrange multipliers (or $\Lambda$-amplitudes) are
required. These are optimized within each cluster after the $T$-amplitudes are
found, giving rise to a set of $\Lambda^{{\mathbf{x}}}$ amplitudes for each
cluster, ${\mathbf{x}}$. The principle of reconstruction of a ‘global’ set of
$\Lambda$-amplitudes follows symmetric projection of the occupied indices of
each cluster $\Lambda^{{\mathbf{x}}}$-amplitudes onto the fragment space
(defined analogously to Eqs. 31 and 32), and a summation over the cluster
solutions, as

$\displaystyle\lambda_{a}^{i}=\sum_{\mathbf{x}}^{N_{\mathrm{frag}}}(\hat{P}^{\mathbf{x}}\Lambda_{1}^{\mathbf{x}})_{a}^{i},$
(43)
$\displaystyle\lambda_{ab}^{ij}=\sum_{\mathbf{x}}^{N_{\mathrm{frag}}}(\hat{P}^{\mathbf{x}}\Lambda_{2}^{\mathbf{x}})_{ab}^{ij}.$
(44)

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 $\Lambda$-amplitudes. The central idea of the
efficient ‘distributed’ approach is to avoid the explicit reconstruction of
the global $T$\- and $\Lambda$-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 needs 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
$\Lambda$-amplitudes which are contracted together.

For the 1-RDM, this maximum order is formally three, corresponding to
contractions of the type $\Lambda_{1}T_{1}^{2}$. 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 $T_{1}$\- and
$\Lambda_{1}$-amplitudes from the cluster solutions according to Eqs. (31,
43), since this will only require a scaling of $\mathcal{O}(N^{2})$ and
$\mathcal{O}(N^{3})$ 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 $T_{2}$-amplitudes, which
we have already argued is prohibitive. The fragment projection for $T_{1}$ and
$\Lambda_{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 $T_{1}$\- and $\Lambda_{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. $T_{2}$ or $\Lambda_{2}$) amplitudes in the 1-RDM
expression, which is only quadratic, resulting from an $L_{2}T_{2}$
contraction. This requires only an $\mathcal{O}(N^{2})$ double summation over
pairs of cluster amplitudes (which can be further reduced asymptotically down
to $\mathcal{O}(N)$ as described later). If certain contractions require a
$T_{1}$ to be contracted with a cluster distributed $T_{2}$, then the global
$T_{1}$ amplitudes can be projected into the required cluster space of the
$T_{2}$.

We demonstrate this distributed amplitude algorithm for the example of the
CCSD contribution to the occupied–occupied block of the 1-RDM, given as

$\begin{split}\Delta\gamma_{ij}&=\Delta\gamma_{ij}^{(\Lambda_{1},T_{1})}+\Delta\gamma_{ij}^{(\Lambda_{2},T_{2})}\\\
&=\sum_{a}^{N_{\mathrm{vir}}}\lambda_{a}^{i}t_{j}^{a}+\sum_{k}^{N_{\mathrm{occ}}}\sum_{ab}^{N_{\mathrm{vir}}}\lambda_{ab}^{ik}\left(2t_{jk}^{ab}-t_{jk}^{ba}\right).\end{split}$
(45)

The first term of this contribution can readily be calculated from the
explicitly combined global singles amplitudes with $\mathcal{O}(N^{3})$
scaling. For the second term, we can calculate the contribution arising from a
specific cluster pair (${\mathbf{x}},{\mathbf{y}}$), by using the
$\Lambda_{2}$\- and $T_{2}$-amplitudes corresponding to these clusters, as

$\begin{split}\Delta\gamma^{\mathbf{xy}}_{i_{x}j_{y}}=&\sum_{k_{x}}^{N_{\mathrm{occ}}^{{\mathbf{x}}}}\sum_{a_{x}b_{x}}^{N_{\mathrm{vir}}^{{\mathbf{x}}}}\sum_{k_{y}}^{N_{\mathrm{occ}}^{{\mathbf{y}}}}\sum_{a_{y}b_{y}}^{N_{\mathrm{vir}}^{{\mathbf{y}}}}S_{a_{x}a_{y}}^{\mathbf{xy}}S_{b_{x}b_{y}}^{\mathbf{xy}}S_{k_{x}k_{y}}^{\mathbf{xy}}\\\
&[\hat{P}^{\mathbf{x}}\Lambda^{\mathbf{x}}_{2}]_{a_{x}b_{x}}^{i_{x}k_{x}}\left(2[\hat{P}^{\mathbf{y}}T^{\mathbf{y}}_{2}]_{j_{y}k_{y}}^{a_{y}b_{y}}-[\hat{P}^{\mathbf{y}}T^{\mathbf{y}}_{2}]_{j_{y}k_{y}}^{b_{y}a_{y}}\right),\end{split}$
(46)

where we use the notation $i_{x}$ ($a_{x}$) to indicate an occupied (virtual)
orbital of cluster ${\mathbf{x}}$, and $S^{\mathbf{xy}}$ to represents the
overlap matrix between the occupied (or virtual, as indicated by the subscript
indices) orbitals of cluster ${\mathbf{x}}$ with cluster ${\mathbf{y}}$, given
by

$S_{p_{x}p_{y}}^{\mathbf{xy}}=\sum_{\alpha}^{N_{\mathrm{ao}}}C_{\alpha
p_{x}}^{\mathbf{x}}\left[\sum_{\beta}^{N_{\mathrm{ao}}}S_{\alpha\beta}C_{\beta
p_{y}}^{\mathbf{y}}\right].$ (47)

These $S^{\mathbf{xy}}$ matrices can be precomputed in $\mathcal{O}(N^{3})$
for all cluster pairs, where the columns of $C^{\mathbf{x}}$ are the relevant
orbitals of cluster ${\mathbf{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 to

$\begin{split}\Delta\gamma_{ij}^{(\Lambda_{2},T_{2})}=\sum_{\mathbf{x}}^{N_{\mathrm{frag}}}\sum_{i_{x}}^{N_{\mathrm{occ}}^{\mathbf{x}}}R_{ii_{x}}^{\mathbf{x}}\left[\sum_{\mathbf{y}}^{N_{\mathrm{frag}}}\sum_{j_{y}}^{N_{\mathrm{occ}}^{\mathbf{y}}}R_{jj_{y}}^{\mathbf{y}}\Delta\gamma^{{\mathbf{xy}}}_{i_{x}j_{y}}\right],\end{split}$
(48)

where $R$ represents the overlap between cluster MOs and MOs of the full
system, i.e.

$R_{pp_{x}}^{{\mathbf{x}}}=\sum_{\alpha}^{N_{\mathrm{ao}}}C_{\alpha
p}\left[\sum_{\beta}^{N_{\mathrm{ao}}}S_{\alpha\beta}C_{\beta
p_{x}}^{\mathbf{x}}\right].$ (49)

To analyse 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 (46) scales as $\mathcal{O}(N^{2})$, while Eqs. (47–49) can be
evaluated in $\mathcal{O}(N^{3})$ 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 $\mathcal{O}(N^{3})$ scaling algorithm.

In Fig. 5 we compare the scaling of the global amplitude and distributed
amplitude algorithms, for a series of alkanes between hexane (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 [53]. This bath construction is similar to pair natural
orbitals (PNO)—though traditional PNOs are applied only for the virtual
coupling of a localized electron pair [72, 73, 68]. The completeness of this
bath space is controlled by a threshold $\eta$, where the bath becomes
complete as $\eta\rightarrow 0$, and is limited to just the traditional DMET
bath as $\eta$ becomes large. More details on the construction and motivation
of this bath expansion can be found in Ref. 53.

Figure 5: Timings for the construction of the full system 1-RDM from the
partitioned wave function for two different bath truncation thresholds,
$\eta=10^{-4}$ and $\eta=10^{-6}$, for the ‘global’ and ‘distributed’ $T_{2}$
algorithms. Increasing the bath size has little effect on the global $T_{2}$
algorithm, whereas it leads to an increase in the prefactor (but not the
exponent) of the distributed $T_{2}$ algorithm.

From the timing data for the 1-RDM construction in Fig. 5 for two different
bath space sizes, we fit the polynomial $t(N)=aN^{b}$ 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\sim
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
$\mathcal{O}(N^{3})$ operations of Eqs. (47–49), whereas the evaluation of Eq.
(46), although scaling as $\mathcal{O}(N^{2})$, 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 inter-cluster physics. In the large system
limit, this also naturally leads to a $\mathcal{O}(N)$ scaling in the 1-RDM
construction, due to the fact that each cluster will only have appreciable
overlap with $\mathcal{O}(1)$ other clusters. In addition, we consider
efficiency gains that can be made from $\mathbf{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.

###  V.2 ‘In-cluster’ approximation to the two-RDM

In principle, an analogous approach to 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 $T_{n}$ or
$\Lambda_{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 $T_{1}$ and $\Lambda_{1}$ contributions, the 2-RDM would include
simultaneous contributions of up to five cluster amplitudes, but by explicitly
reconstructing these $T_{1}$ and $L_{1}$ amplitudes globally, only looping
over three clusters (${\mathbf{x}},{\mathbf{y}},{\mathbf{z}}$) would be
required. However, including the projections back to the full space and noting
the $\mathcal{O}(N^{4})$ 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 $\mathcal{O}(N^{3})$ 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 non-cumulant 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 non-
cumulant 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,
$\delta{K}$, which involves a $(\Lambda_{2}T_{2}^{2})$ triple product of
double (de)excitation amplitudes. If evaluated exactly according to the
implicit partitioned wave function, this would require contributions from all
triplets $({\mathbf{x}},{\mathbf{y}},{\mathbf{z}})$ of clusters, as

$\begin{split}\delta{K}^{{\mathbf{xyz}}}_{i_{y}a_{z}j_{y}b_{z}}=&\sum_{k_{x}l_{x}}\sum_{c_{x}d_{x}}\sum_{c_{y}d_{y}}\sum_{k_{z}k_{z}}S^{\mathbf{xy}}_{c_{x}c_{y}}S^{\mathbf{xy}}_{d_{x}d_{y}}S^{\mathbf{xz}}_{k_{x}k_{z}}S^{\mathbf{xz}}_{l_{x}l_{z}}\\\
&[\hat{P}^{\mathbf{x}}\Lambda^{\mathbf{x}}_{2}]_{c_{x}d_{x}}^{k_{x}l_{x}}[\hat{P}^{\mathbf{y}}T^{\mathbf{y}}_{2}]_{i_{y}j_{y}}^{c_{y}d_{y}}[\hat{P}^{\mathbf{z}}T^{\mathbf{z}}_{2}]_{k_{z}l_{z}}^{a_{z}b_{z}}.\end{split}$
(50)

In the in-cluster approximation, this contribution is replaced by

$\delta{K}^{{\mathbf{x}}}_{i_{x}a_{x}j_{x}b_{x}}=\sum_{k_{x}l_{x}}\sum_{c_{x}d_{x}}[\hat{P}^{\mathbf{x}}\Lambda^{\mathbf{x}}_{2}]_{c_{x}d_{x}}^{k_{x}l_{x}}[T^{\mathbf{x}}_{2}]_{i_{x}j_{x}}^{c_{x}d_{x}}[T^{\mathbf{x}}_{2}]_{k_{x}l_{x}}^{a_{x}b_{x}},$
(51)

which only contains amplitudes from cluster ${\mathbf{x}}$. Note that Eq. (51)
only contains a single fragment-space projector, which is applied to the
$\Lambda$-amplitude, compared to Eq. (50), where every amplitude is projected.
In general, the number of projectors in a term has to agree to 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 $\Lambda$-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
$\Lambda$-amplitudes yields a simple recipe, since all contributions to the
2-RDM cumulant are at most linear in $\Lambda_{1}$ or $\Lambda_{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 $\Lambda^{0}$ terms of coupled-cluster expressions,
which results in a MP2-like contribution

$\delta{K}_{iajb}=t_{ij}^{ab}+t_{i}^{a}t_{j}^{b}-\frac{1}{2}\left[t_{ij}^{ba}+t_{i}^{b}t_{j}^{a}\right],$
(52)

in which case we use projected $T$-amplitudes, according to

$\begin{split}\delta{K}_{i_{x}a_{x}j_{x}b_{x}}^{\mathbf{x}}&=[\hat{P}^{\mathbf{x}}T^{\mathbf{x}}_{2}]_{i_{x}j_{x}}^{a_{x}b_{x}}+[\hat{P}^{\mathbf{x}}T_{1}^{\mathbf{x}}]_{i_{x}}^{a_{x}}[T_{1}^{\mathbf{x}}]_{j_{x}}^{b_{x}}\\\
&-\frac{1}{2}\left([\hat{P}^{\mathbf{x}}T^{\mathbf{x}}_{2}]_{i_{x}j_{x}}^{b_{x}a_{x}}+[\hat{P}^{\mathbf{x}}T_{1}^{\mathbf{x}}]_{i_{x}}^{b_{x}}[T_{1}^{\mathbf{x}}]_{j_{x}}^{a_{x}}\right).\end{split}$
(53)

We note here that the ‘global’ $T_{1}$\- and $\Lambda_{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 non-local approach (as used
for the 1-RDM in Sec. V.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 $\mathcal{O}(N^{4})$ 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[(\gamma,K^{*})[\Psi^{{\mathbf{x}}}])$ in this work (the asterisk denoting
the ‘in-cluster’ approximation to the two-body cumulant), as

$\begin{split}E[(\gamma,{K}^{*})[\Psi^{\mathbf{x}}]]&=E_{\mathrm{HF}}[\gamma^{0}]+\sum_{pq}F_{pq}[\gamma^{0}]\Delta\gamma_{qp}[\Psi^{\mathbf{x}}]\\\
&+\frac{1}{2}\sum_{{\mathbf{x}}}^{N_{\mathrm{frag}}}\sum_{pqrs}^{N_{\mathrm{cl}}^{\mathbf{x}}}(pq|rs){{K}^{*}}[\hat{P}^{\mathbf{x}}\Psi^{\mathbf{x}}])_{pqrs},\end{split}$
(54)

where ${{K}^{*}}[\hat{P}^{\mathbf{x}}\Psi^{\mathbf{x}}]$ 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
$\Delta\gamma[\Psi^{{\mathbf{x}}}]$ which is assumed to contain all relevant
products of different cluster solutions (see Sec. V.1). Other static two-body
expectation values can be derived using this ‘in-cluster’ two-body cumulant
analogously. In Sec. VI.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. II.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
$\Lambda$-amplitudes used for all expressions. This can be seen in the first
set of results in Fig. 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[(\gamma,\Gamma)[\Psi^{{\mathbf{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[(\gamma,K^{*})[\Psi^{{\mathbf{x}}}]]$).

In the next section, we will rigorously benchmark these different approaches
to computing energies and non-energetic 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.

##  VI Results

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 $\eta$, and we will
use CCSD as a cluster solver throughout. As $\eta$ is reduced towards zero,
the bath space increases in size towards 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[\Psi^{{\mathbf{x}}}]$: Introduced in Ref. 53 and Sec. III.3, this is the
cheapest approach to energy computation of Eq. (37), requiring no
communication between independent additive cluster contributions.

  * •

$E[(\gamma,\Gamma)[\Psi^{{\mathbf{x}}}]]$: Described in Sec. III.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 $\Gamma$, and rapidly becomes the dominant
cost in the calculation. It will primarily be used for reference.

  * •

$E[(\gamma,K^{*})[\Psi^{{\mathbf{x}}}]]$: Described in Sec. V.2, this
approximates the cumulant contribution to the 2-RDM via the ‘in-cluster’
approximation, which ensures low-scaling with system size for this non-
factorizable 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.
V.1.

We do not compare the energies from the density matrix derived expectation
values of Sec. II, 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 $\eta$.

###  VI.1 Molecular results: W4-11 test set

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 as well as
dynamic correlated physics present, as well as a mix of high-spin, low-spin
and open-shell systems [54]. 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 [61] (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,
$\eta$. Smaller $\eta$ values correspond to larger bath spaces, with $\eta>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.

Figure 6: Convergence of the reconstructed system energy from individual
cluster solutions to the canonical CCSD with respect to the maximum cluster
size for any atomic fragment in any molecule in the W4-11 test set. Cluster
sizes are controlled by the bath threshold ($\eta$), starting from just the
DMET bath, for the wave function derived energy functionals. Points indicate
the mean absolute error (MAE) per electron in the total energy, while error
bars show the standard deviation over the test set. Inset shows a magnified
version of the larger bath threshold results. Figure 7: Violin plots showing
the distribution of total energy errors per electron over the W4-11 test set
of molecules compared to canonical CCSD calculations, for the different
partitioned wave function energy estimators. This is performed with a bath
expansion of each IAO atomic fragment for different values of $\eta$, which
corresponds to a maximum cluster size as indicated in Fig. 6. Total
distribution ranges show the minimum and maximum errors over the test set.

These results are shown in Fig. 6, with their absolute per electron energy
errors aggregated for different values of $\eta$ in Fig. 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[\Psi^{{\mathbf{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[(\gamma,\Gamma)[\Psi^{{\mathbf{x}}}]]$ functional, which contains many
products of contributions between different cluster solutions in its
construction. Furthermore, this energy is rigorously $\Psi$-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 non-variationality 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[(\gamma,K^{*})[\Psi^{{\mathbf{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[(\gamma,\Gamma)[\Psi^{{\mathbf{x}}}]]$ energy, but due to favourable
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[\Psi^{{\mathbf{x}}}]$ energy. Due to the
approximation in this cumulant non-variational total energies can result
(compared to the canonical result), however the appearance of these non-
variational 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 favourable
cancellation of errors, and maintenance of generally variational results, even
in this approximation which breaks strict $N$-representability in the 2-RDM.

###  VI.2 Solid state results: Diamond

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$\times$5$\times$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 of 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 $\eta$ parameter. The size of the resulting clusters vary
between 50 and $\sim 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 Fig. 8 for the linear energy functional
$E[\Psi^{{\mathbf{x}}}]$ and the ‘in-cluster’ cumulant functional,
$E[(\gamma,K^{*})[\Psi^{{\mathbf{x}}}]]$, for a range of lattice parameters of
the unit cell. The $E[(\gamma,\Gamma)[\Psi^{{\mathbf{x}}}]]$ functional is
omitted due to its unfavourable 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[5, 74]. 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[75] as
described in Ref. 53. The BSSE per atom is estimated by computing two
independent molecular sytems, the first for 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\times 3\times 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 $\sim 7.2$ m$E_{\textrm{H}}$, 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 ($\eta$) values are shown in Fig. 8. These demonstrate the
systematic improvement in the total energy as $\eta$ decreases. However, while
the total energy is not fully converged, we can fit these equations of state
to a Birch-Murnaghan form[76] about approximately $\pm 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 Fig. 9, showing convergence to
experimental values, corrected for zero-point vibrational effects [77, 5]. 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 Fig. 7. Our most accurate result ($\eta=10^{-7.5}$,
resulting from a single CCSD calculation on $\sim 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 $\sim 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[77], an order of magnitude greater than these results.

Figure 8: Equation of state for diamond, using both the linear energy
functional ($E[\Psi^{{\mathbf{x}}}]$, purple) and the RDM energy functional
with in-cluster approximated cumulant
($E[(\gamma,K^{*})[\Psi^{{\mathbf{x}}}]]$, red). Data from different lattice
parameters are fit to a Birch-Murnaghan equation of state for different values
of $\eta$, systematically expanding the size of the cluster space. $5\times
5\times 5$ ${\bf k}$-points are used in the supercell, with an all-electron
cc-pVTZ basis, while small energy corrections for the basis set superposition
error are also included. The calculated energies are shown as circles, the
fits are solid lines and the equilibrium lattice parameter from the fit is
indicated by pentagons. Figure 9: Convergence of equilibrium lattice
parameters (circles, left $y$-axis) and bulk modulus (diamonds, right
$y$-axis) for a $5\times 5\times 5$ diamond supercell as the total number of
cluster orbitals increases (controlled by the $\eta$ value). Values are
estimated from a Birch-Murnaghan fit to the equation of state shown in Fig. 8.
Results are shown for both the linear energy functional
($E[\Psi^{{\mathbf{x}}}]$) and the RDM energy functional with in-cluster
approximated cumulant ($E[(\gamma,K^{*})[\Psi^{{\mathbf{x}}}]]$), with
experimental values shown as horizontal lines [77, 5].

###  VI.3 Beyond energetics: Non-local spin correlation functions

While accurate energetics are important, they are far from the only non-local
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 by

$\begin{split}\braket{\hat{S}_{z}^{A}\hat{S}_{z}^{B}}&=\frac{1}{4}\sum_{ijkl}P_{ij}^{A}P_{kl}^{B}\left(\Gamma_{ijkl}^{\alpha\alpha}-\Gamma_{ijkl}^{\alpha\beta}-\Gamma_{ijkl}^{\beta\alpha}+\Gamma_{ijkl}^{\beta\beta}\right)\\\
&+\frac{1}{4}\sum_{ijk}P_{ik}^{A}P_{jk}^{B}\left(\gamma_{ij}^{\alpha}+\gamma_{ij}^{\beta}\right),\end{split}$
(55)

where $P^{A}$ 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, $\braket{S_{z}^{2}}=\sum_{AB}\braket{\hat{S}_{z}^{A}\hat{S}_{z}^{B}}$,
are necessarily zero in a spin-restricted formalism (though not in
unrestricted), the local contribution from a given atom pair is generally non-
zero, indicating the conditional spin-density between different points in the
system.

We adapt the fragmentation and BNO bath expansion procedure as detailed in
Sec. VI.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 Fig.
10) in the cc-pVTZ basis between the radical C1 position (atom $A$) and other
atoms in the system (atom $B$).

Figure 10: Labels of the atoms in the $n$-propyl radical used in this work.

In order to compare results from the traditional democratically partitioned
2-RDM of Sec. II.1 ($\Gamma[\Gamma^{{\mathbf{x}}}]$), the democratically
partitioned cumulant of Sec. II.2
($\Gamma[\gamma^{{\mathbf{x}}},K^{{\mathbf{x}}}]$) and the in-cluster cumulant
approximated 2-RDM derived from the partitioned wave function of Sec. V.2
($\Gamma[(\gamma,K^{*})[\Psi^{{\mathbf{x}}}]]$), we choose a complete
fragmentation of the system via atomic IAO $\oplus$ 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 is a sufficient but not necessary condition for the
$\Gamma[(\gamma,K^{*})[\Psi^{{\mathbf{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 $\Gamma[(\gamma,K^{*})[\Psi^{{\mathbf{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:  $S_{z}^{2}$-correlation functions in the $n$-propyl radical
between the radical position (C1) and other carbon (first column) and hydrogen
atoms (second column) in the molecule (see Fig. 10). Spin correlators are
derived from the 2-RDM computed in different ways; traditional democratic
partitioning of Sec. II.1 (denoted $\Gamma[\Gamma^{{\mathbf{x}}}]$, green),
democratic partitioning of the cumulant of Sec. II.2 (denoted
$\Gamma[\gamma^{{\mathbf{x}}},K^{{\mathbf{x}}}]$, yellow) and the wave
function derived 2-RDM with in-cluster approximated cumulant of Sec. V
(denoted $\Gamma[(\gamma,K^{*})[\Psi^{{\mathbf{x}}}]]$, red). The calculation
was performed in the cc-pVTZ basis set (188 orbitals) with SAO projectors to
evaluate the spin correlation functions, and IAO $\oplus$ PAO atomic
fragmentation.

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 non-local 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. II.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 $\approx$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.

##  VII Conclusions and Outlook

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 non-energetic and non-local quantity
via the reconstructed reduced density matrices. A fundamental difference
exists between 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 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[48, 49, 50, 51, 52]).

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 non-local 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 $T_{1}$ amplitudes [78].
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 non-
local physics and is fully $N$-representable compared to the standard
democratic partitioning [20]. 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.

## Acknowledgements

G.H.B. gratefully acknowledges support from the Royal Society via a University
Research Fellowship, as well as funding from the European Union’s Horizon 2020
research and innovation programme under grant agreement No. 759063. We are
grateful to the UK Materials and Molecular Modelling Hub for computational
resources, which is partially funded by EPSRC (EP/P020194/1 and EP/T022213/1).

## References

  * Booth  _et al._ [2013] G. H. Booth, A. Grüneis, G. Kresse, and A. Alavi, “Towards an exact description of electronic wavefunctions in real solids,” Nature 493, 365–370 (2013).
  * Pisani  _et al._ [2012] C. Pisani, M. Schütz, S. Casassa, D. Usvyat, L. Maschio, M. Lorenz, and A. Erba, “Cryscor: a program for the post-hartree–fock treatment of periodic systems,” Phys. Chem. Chem. Phys. 14, 7615–7628 (2012).
  * Zhang and Grüneis [2019] I. Y. Zhang and A. Grüneis, “Coupled Cluster Theory in Materials Science,” Frontiers in Materials 6 (2019).
  * Gruber  _et al._ [2018] T. Gruber, K. Liao, T. Tsatsoulis, F. Hummel, and A. Grüneis, “Applying the Coupled-Cluster Ansatz to Solids and Surfaces in the Thermodynamic Limit,” Physical Review X 8, 021043 (2018).
  * McClain  _et al._ [2017] J. McClain, Q. Sun, G. K.-L. Chan, and T. C. Berkelbach, “Gaussian-Based Coupled-Cluster Theory for the Ground-State and Band Structure of Solids,” Journal of Chemical Theory and Computation 13, 1209–1218 (2017).
  * Xing and Lin [2022] X. Xing and L. Lin, “Staggered mesh method for correlation energy calculations of solids: Random phase approximation in direct ring coupled cluster doubles and adiabatic connection formalisms,” Journal of Chemical Theory and Computation 18, 763–775 (2022), pMID: 34989566, https://doi.org/10.1021/acs.jctc.1c00985 .
  * Herbert [2019] J. M. Herbert, “Fantasy versus reality in fragment-based quantum chemistry,” The Journal of Chemical Physics 151, 170901 (2019), https://doi.org/10.1063/1.5126216 .
  * Gordon  _et al._ [2012] M. S. Gordon, D. G. Fedorov, S. R. Pruitt, and L. V. Slipchenko, “Fragmentation methods: A route to accurate calculations on large systems,” Chemical Reviews 112, 632–672 (2012), pMID: 21866983, https://doi.org/10.1021/cr200093j .
  * Saebø [2002] S. Saebø, “Low-scaling methods for electron correlation,” in  _Computational Chemistry: Reviews of Current Trends_ (World Scientific, 2002) pp. 63–87.
  * Usvyat, Maschio, and Schütz [2018] D. Usvyat, L. Maschio, and M. Schütz, “Periodic and fragment models based on the local correlation approach,” WIREs Computational Molecular Science 8, 1–27 (2018).
  * Schäfer  _et al._ [2021] T. Schäfer, F. Libisch, G. Kresse, and A. Grüneis, “Local embedding of coupled cluster theory into the random phase approximation using plane waves,” The Journal of Chemical Physics 154, 011101 (2021).
  * Sun and Chan [2016] Q. Sun and G. K.-l. Chan, “Quantum Embedding Theories,” Accounts of Chemical Research 49, 2705–2712 (2016).
  * Fertitta and Booth [2018] E. Fertitta and G. H. Booth, “Rigorous wave function embedding with dynamical fluctuations,” Physical Review B 98, 235132 (2018).
  * Nusspickel and Booth [2020a] M. Nusspickel and G. H. Booth, “Efficient compression of the environment of an open quantum system,” Phys. Rev. B 102, 165107 (2020a).
  * Nusspickel and Booth [2020b] M. Nusspickel and G. H. Booth, “Frequency-dependent and algebraic bath states for a dynamical mean-field theory with compact support,” Phys. Rev. B 101, 045126 (2020b).
  * Potthoff [2001] M. Potthoff, “Two-site dynamical mean-field theory,” Phys. Rev. B 64, 165114 (2001).
  * Lu  _et al._ [2019] Y. Lu, X. Cao, P. Hansmann, and M. W. Haverkort, “Natural-orbital impurity solver and projection approach for green’s functions,” Phys. Rev. B 100, 115134 (2019).
  * Ganahl  _et al._ [2015] M. Ganahl, M. Aichhorn, H. G. Evertz, P. Thunström, K. Held, and F. Verstraete, “Efficient dmft impurity solver using real-time dynamics with matrix product states,” Phys. Rev. B 92, 155132 (2015).
  * Scott and Booth [2021] C. J. C. Scott and G. H. Booth, “Extending density matrix embedding: A static two-particle theory,” Phys. Rev. B 104, 245114 (2021).
  * Wu  _et al._ [2019] X. Wu, Z.-H. Cui, Y. Tong, M. Lindsey, G. K.-L. Chan, and L. Lin, “Projected density matrix embedding theory with applications to the two-dimensional hubbard model,” The Journal of Chemical Physics 151, 064108 (2019), https://doi.org/10.1063/1.5108818 .
  * Wu  _et al._ [2020] X. Wu, M. Lindsey, T. Zhou, Y. Tong, and L. Lin, “Enhancing robustness and efficiency of density matrix embedding theory via semidefinite programming and local correlation potential fitting,” Phys. Rev. B 102, 085123 (2020).
  * Faulstich  _et al._ [2022] F. M. Faulstich, R. Kim, Z.-H. Cui, Z. Wen, G. Kin-Lic Chan, and L. Lin, “Pure state v-representability of density matrix embedding theory,” Journal of Chemical Theory and Computation 18, 851–864 (2022), pMID: 35084855, https://doi.org/10.1021/acs.jctc.1c01061 .
  * Metzner and Vollhardt [1989] W. Metzner and D. Vollhardt, “Correlated Lattice Fermions in $\infty$ Dimensions,” Physical Review Letters 62, 324–327 (1989).
  * Georges  _et al._ [1996] A. Georges, G. Kotliar, W. Krauth, and M. J. Rozenberg, “Dynamical mean-field theory of strongly correlated fermion systems and the limit of infinite dimensions,” Reviews of Modern Physics 68, 13–125 (1996).
  * Kotliar  _et al._ [2006] G. Kotliar, S. Y. Savrasov, K. Haule, V. S. Oudovenko, O. Parcollet, and C. A. Marianetti, “Electronic structure calculations with dynamical mean-field theory,” Reviews of Modern Physics 78, 865–951 (2006).
  * Karolak  _et al._ [2010] M. Karolak, G. Ulm, T. Wehling, V. Mazurenko, A. Poteryaev, and A. Lichtenstein, “Double counting in LDA+DMFT—The example of NiO,” Journal of Electron Spectroscopy and Related Phenomena 181, 11–15 (2010), proceedings of International Workshop on Strong Correlations and Angle-Resolved Photoemission Spectroscopy 2009\.
  * Zhu, Cui, and Chan [2020] T. Zhu, Z.-H. Cui, and G. K.-L. Chan, “Efficient formulation of ab initio quantum embedding in periodic systems: Dynamical mean-field theory,” Journal of Chemical Theory and Computation 16, 141–153 (2020).
  * Rusakov  _et al._ [2019] A. A. Rusakov, S. Iskakov, L. N. Tran, and D. Zgid, “Self-energy embedding theory (seet) for periodic systems,” Journal of Chemical Theory and Computation 15, 229–240 (2019).
  * Zhu and Chan [2021] T. Zhu and G. K.-L. Chan, “Ab Initio Full Cell $GW$ \+ DMFT for Correlated Materials,” Physical Review X 11, 021006 (2021).
  * Knizia and Chan [2012] G. Knizia and G. K.-L. Chan, “Density Matrix Embedding: A Simple Alternative to Dynamical Mean-Field Theory,” Physical Review Letters 109, 186404 (2012).
  * Knizia and Chan [2013] G. Knizia and G. K.-l. Chan, “Density Matrix Embedding: A Strong-Coupling Quantum Embedding Theory,” Journal of Chemical Theory and Computation 9, 1428–1432 (2013).
  * Bulik, Scuseria, and Dukelsky [2014] I. W. Bulik, G. E. Scuseria, and J. Dukelsky, “Density matrix embedding from broken symmetry lattice mean fields,” Phys. Rev. B 89, 035140 (2014).
  * Wouters  _et al._ [2016] S. Wouters, C. A. Jiménez-Hoyos, Q. Sun, and G. K. Chan, “A Practical Guide to Density Matrix Embedding Theory in Quantum Chemistry,” Journal of Chemical Theory and Computation 12, 2706–2719 (2016).
  * Welborn, Tsuchimochi, and Van Voorhis [2016] M. Welborn, T. Tsuchimochi, and T. Van Voorhis, “Bootstrap embedding: An internally consistent fragment-based method,” J. Chem. Phys. 145, 074102 (2016).
  * Cui, Zhu, and Chan [2020] Z.-H. Cui, T. Zhu, and G. K.-L. Chan, “Efficient Implementation of Ab Initio Quantum Embedding in Periodic Systems: Density Matrix Embedding Theory,” Journal of Chemical Theory and Computation 16, 119–129 (2020).
  * Pham, Hermes, and Gagliardi [2020] H. Q. Pham, M. R. Hermes, and L. Gagliardi, “Periodic electronic structure calculations with the density matrix embedding theory,” Journal of Chemical Theory and Computation 16, 130–140 (2020).
  * Lau, Knizia, and Berkelbach [2021] B. T. Lau, G. Knizia, and T. C. Berkelbach, “Regional Embedding Enables High-Level Quantum Chemistry for Surface Science,” Journal of Physical Chemistry Letters 12, 1104–1109 (2021).
  * Sekaran  _et al._ [2021] S. Sekaran, M. Tsuchiizu, M. Saubanère, and E. Fromager, “Householder-transformed density matrix functional embedding theory,” Phys. Rev. B 104, 035121 (2021).
  * Hermes and Gagliardi [2019] M. R. Hermes and L. Gagliardi, “Multiconfigurational self-consistent field theory with density matrix embedding: The localized active space self-consistent field method,” Journal of Chemical Theory and Computation 15, 972–986 (2019), pMID: 30620876, https://doi.org/10.1021/acs.jctc.8b01009 .
  * Lin and Lindsey [2022] L. Lin and M. Lindsey, “Variational embedding for quantum many-body problems,” Communications on Pure and Applied Mathematics 75, 2033–2068 (2022), https://onlinelibrary.wiley.com/doi/pdf/10.1002/cpa.21984 .
  * Cui  _et al._ [2022] Z.-H. Cui, H. Zhai, X. Zhang, and G. K.-L. Chan, “Systematic electronic structure in the cuprate parent state from quantum many-body simulations,” Science 377, 1192–1198 (2022), https://www.science.org/doi/pdf/10.1126/science.abm2295 .
  * Cui  _et al._ [2020] Z.-H. Cui, C. Sun, U. Ray, B.-X. Zheng, Q. Sun, and G. K.-L. Chan, “Ground-state phase diagram of the three-band hubbard model from density matrix embedding theory,” Phys. Rev. Research 2, 043259 (2020).
  * Mitra  _et al._ [2022] A. Mitra, M. R. Hermes, M. Cho, V. Agarawal, and L. Gagliardi, “Periodic density matrix embedding for co adsorption on the mgo(001) surface,” The Journal of Physical Chemistry Letters 13, 7483–7489 (2022), pMID: 35939641, https://doi.org/10.1021/acs.jpclett.2c01915 .
  * Cao  _et al._ [2022] C. Cao, J. Sun, X. Yuan, H.-S. Hu, H. Q. Pham, and D. Lv, “Ab initio quantum simulation of strongly correlated materials with quantum embedding,”  (2022).
  * Lee  _et al._ [2019] T.-H. Lee, T. Ayral, Y.-X. Yao, N. Lanata, and G. Kotliar, “Rotationally invariant slave-boson and density matrix embedding theory: Unified framework and comparative study on the one-dimensional and two-dimensional hubbard model,” Phys. Rev. B 99, 115129 (2019).
  * Fertitta and Booth [2019] E. Fertitta and G. H. Booth, “Energy-weighted density matrix embedding of open correlated chemical fragments,” The Journal of Chemical Physics 151, 014115 (2019).
  * Sriluckshmy  _et al._ [2021] P. V. Sriluckshmy, M. Nusspickel, E. Fertitta, and G. H. Booth, “Fully algebraic and self-consistent effective dynamics in a static quantum embedding,” Physical Review B 103, 085131 (2021).
  * Mineh and Montanaro [2022] L. Mineh and A. Montanaro, “Solving the hubbard model using density matrix embedding theory and the variational quantum eigensolver,” Phys. Rev. B 105, 125117 (2022).
  * [49] J. J. M. Kirsopp, C. Di Paola, D. Z. Manrique, M. Krompiec, G. Greene-Diniz, W. Guba, A. Meyder, D. Wolf, M. Strahm, and D. Muñoz Ramo, “Quantum computational quantification of protein–ligand interactions,” International Journal of Quantum Chemistry n/a, e26975, https://onlinelibrary.wiley.com/doi/pdf/10.1002/qua.26975 .
  * Li  _et al._ [2022] W. Li, Z. Huang, C. Cao, Y. Huang, Z. Shuai, X. Sun, J. Sun, X. Yuan, and D. Lv, “Toward practical quantum embedding simulation of realistic chemical systems on near-term quantum computers,” Chem. Sci. 13, 8953–8962 (2022).
  * Kawashima  _et al._ [2021] Y. Kawashima, E. Lloyd, M. P. Coons, Y. Nam, S. Matsuura, A. J. Garza, S. Johri, L. Huntington, V. Senicourt, A. O. Maksymov, J. H. V. Nguyen, J. Kim, N. Alidoust, A. Zaribafiyan, and T. Yamazaki, “Optimizing electronic structure simulations on a trapped-ion quantum computer using problem decomposition,” Communications Physics 4, 245 (2021).
  * Tilly  _et al._ [2021] J. Tilly, P. V. Sriluckshmy, A. Patel, E. Fontana, I. Rungger, E. Grant, R. Anderson, J. Tennyson, and G. H. Booth, “Reduced density matrix sampling: Self-consistent embedding and multiscale electronic structure on current generation quantum computers,” Phys. Rev. Research 3, 033230 (2021).
  * Nusspickel and Booth [2022] M. Nusspickel and G. H. Booth, “Systematic improvability in quantum embedding for real materials,” Phys. Rev. X 12, 011046 (2022).
  * Karton, Daon, and Martin [2011] A. Karton, S. Daon, and J. M. Martin, “W4-11: A high-confidence benchmark dataset for computational thermochemistry derived from first-principles w4 data,” Chemical Physics Letters 510, 165–178 (2011).
  * Note [1] Embedding code and documentation can be found at https://github.com/BoothGroup/Vayesta.
  * Sun  _et al._ [2020] Q. Sun, X. Zhang, S. Banerjee, P. Bao, M. Barbry, N. S. Blunt, N. A. Bogdanov, G. H. Booth, J. Chen, Z.-H. Cui, J. J. Eriksen, Y. Gao, S. Guo, J. Hermann, M. R. Hermes, K. Koh, P. Koval, S. Lehtola, Z. Li, J. Liu, N. Mardirossian, J. D. McClain, M. Motta, B. Mussard, H. Q. Pham, A. Pulkin, W. Purwanto, P. J. Robinson, E. Ronca, E. R. Sayfutyarova, M. Scheurer, H. F. Schurkus, J. E. T. Smith, C. Sun, S.-N. Sun, S. Upadhyay, L. K. Wagner, X. Wang, A. White, J. D. Whitfield, M. J. Williamson, S. Wouters, J. Yang, J. M. Yu, T. Zhu, T. C. Berkelbach, S. Sharma, A. Y. Sokolov, and G. K.-L. Chan, “Recent developments in the PySCF program package,” The Journal of Chemical Physics 153, 024109 (2020).
  * Sun  _et al._ [2018] Q. Sun, T. C. Berkelbach, N. S. Blunt, G. H. Booth, S. Guo, Z. Li, J. Liu, J. D. McClain, E. R. Sayfutyarova, S. Sharma, S. Wouters, and G. K. Chan, “PySCF: the Python‐based simulations of chemistry framework,” WIREs Computational Molecular Science 8 (2018), 10.1002/wcms.1340.
  * Pulay [1983] P. Pulay, “Localizability of dynamic electron correlation,” Chemical Physics Letters 100, 151–154 (1983).
  * Bulik, Chen, and Scuseria [2014] I. W. Bulik, W. Chen, and G. E. Scuseria, “Electron correlation in solids via density embedding theory,” The Journal of Chemical Physics 141, 054113 (2014), https://doi.org/10.1063/1.4891861 .
  * Note [2] We note that the two-body cumulant would often be denoted by $\lambda$, but define by ${\tilde{K}}$ to avoid potential confusions with the lambda-amplitudes of coupled-cluster theory used later in the text.
  * Knizia [2013] G. Knizia, “Intrinsic Atomic Orbitals: An Unbiased Bridge between Quantum Theory and Chemical Concepts,” Journal of Chemical Theory and Computation 9, 4834–4843 (2013).
  * Bartlett and Musiał [2007] R. J. Bartlett and M. Musiał, “Coupled-cluster theory in quantum chemistry,” Rev. Mod. Phys. 79, 291–352 (2007).
  * Lehtola  _et al._ [2017] S. Lehtola, N. M. Tubman, K. B. Whaley, and M. Head-Gordon, “Cluster decomposition of full configuration interaction wave functions: A tool for chemical interpretation of systems with strong correlation,” The Journal of Chemical Physics 147, 154105 (2017), https://doi.org/10.1063/1.4996044 .
  * Shee and Zgid [2019] A. Shee and D. Zgid, “Coupled Cluster as an Impurity Solver for Green’s Function Embedding Methods,” Journal of Chemical Theory and Computation 15, 6010–6024 (2019).
  * Zhu  _et al._ [2019] T. Zhu, C. A. Jiménez-Hoyos, J. McClain, T. C. Berkelbach, and G. K.-L. Chan, “Coupled-cluster impurity solvers for dynamical mean-field theory,” Physical Review B 100, 115154 (2019).
  * Shee, Yeh, and Zgid [2022] A. Shee, C.-N. Yeh, and D. Zgid, “Exploring coupled cluster green’s function as a method for treating system and environment in green’s function embedding methods,” Journal of Chemical Theory and Computation 18, 664–676 (2022), pMID: 34989565, https://doi.org/10.1021/acs.jctc.1c00712 .