Rigorous screened interactions for realistic correlated electron systems

We derive a widely-applicable first principles approach for determining two-body, static effective interactions for low-energy Hamiltonians with quantitative accuracy. The algebraic construction rigorously conserves all instantaneous two-point correlation functions in a chosen model space at the level of the random phase approximation, improving upon the traditional uncontrolled static approximations. Applied to screened interactions within a quantum embedding framework, we demonstrate these faithfully describe the relaxation of local subspaces via downfolding high-energy physics in molecular systems, as well as enabling a systematically improvable description of the long-range plasmonic contributions in extended graphene.

Building effective models that retain the physics of interest but strip away extraneous complexity is central to progress in understanding physical mechanisms and emergent behaviour in complex systems.Nowhere is this more true than interacting electron systems, where effective Hamiltonians occupy a central position in both condensed matter and quantum chemistry, from the Hubbard model to ligand field theory [1][2][3][4].A push is now underway to erode the boundary between phenomenological effective models with empirical parameters, and ab initio modelling with material specificity, where the act of constructing a material-specific effective Hamiltonian is increasingly the first step in a larger workflow involving its subsequent solution within a multi-method approach.This is critical to extend the scope of accurate yet computationally demanding many-body methods, placing significant urgency on approaches in which the relevant physics outside a low-energy model space is rigorously downfolded or renormalized into the effective Hamiltonian [5][6][7].In this work we propose a simple and efficient approach to effective interactions, demonstrating quantitative accuracy and specificity resulting from integrating out ab initio high-energy and long-range physics, motivated by the exact conservation of instantaneous expectation values.
A salient example in the need for accurate renormalized interactions is quantum embedding, describing correlated many-body phenomena within a local subspace [8][9][10][11].The missing interactions with states outside the subspace should renormalize the effective subspace interactions in real materials, generally necessitating the use of effective interactions that are often parameterized empirically via Hubbard or Hund terms [12][13][14], or downfolded from other theories [11,[15][16][17].In quantum chemistry, analogous subspaces are often described by a 'complete active space' (CAS), whereby a small number of low-energy mean-field orbitals are chosen for an accurate treatment of the strong correlation in 'multi-reference' approaches [18].More broadly, a wide range of both qualitative and quantitative studies into correlated materials rely on first obtaining appropriately screened effective interactions of a simplified model from an ab initio starting point [5,19].
It is often argued that the random phase approximation (RPA) contains much of the physics missing from low-energy models, and is the appropriate theory in which to construct effective interactions [11,13,20,21].As an infinite resummation of the bubble diagrams, it correctly describes the long-range charge fluctuations, plasmons, high-energy scattering and many-body dispersion that predominates in high-density metallic and semiconducting extended systems [22][23][24].While its traditional formulation lacks exchange or ladder diagrams, their contributions generally decay more rapidly and can often be captured within the model space, with approaches to screen beyond RPA still an active research area [15,16,25].
The constrained RPA (cRPA) method has therefore become a widespread choice for deriving low-energy effective interactions from first principles, applied from molecules to Mott insulators and routinely as a precursor to quantum embedding methods [5,8,13,16,[25][26][27][28][29][30][31][32][33][34][35][36][37][38][39][40][41].Since the RPA is a well-defined diagrammatic theory, the bubbles corresponding to the polarizability within the model space (Π m (ω)) can be removed from the total polarization (Π ext (ω) = Π(ω) − Π m (ω)).The screened interactions, U(ω), are then found from the infinite RPA resummation coupling this space to the external degrees of freedom via this scattering channel, as where v denotes the bare Coulomb interaction.In this way, double counting of correlated effects are avoided once the resulting effective model is solved, yet direct scattering events to all orders between the two spaces are included [26].The cRPA approach is also internally consistent: an RPA calculation on the resulting effective Hamiltonian with interactions U(ω), will give exactly the results of the RPA on the full system.More specifically, the cRPA ensures that the RPA density-density (dd) response function, χ(ω), of the model space with effective interactions U(ω), is identical to the projection of the RPA full system dd-response into the model space.In this way, the effective interactions can also be seen to be 'state universal', providing the correctly renormalized interactions for the entire RPA spectrum of the model subspace.
While this puts the cRPA approach on a solid footing, few methods are computationally tractable with the resulting dynamical interactions of Eq. 1. Practical necessity therefore forces the widespread approximation of taking the static, ω → 0 limit of the screened interaction of Eq. 1 (static-cRPA).This uncontrolled approximation is qualitatively justified in capturing the relevant long-wavelength behaviour where the RPA is accurate [25], however forces us to entirely give up on rigorous conservation of any expectation values from the RPA.There is significant scepticism of its accuracy, with the cRPA interactions often over-screening the physics [15,16,25,28,42].Far from purely a small quantitative shift, this can result in wrong energy ordering of states and phases in the resulting model, or missing spectral weight transfer to plasmon satellites [14,35,43].
A further serious technical limitation of the cRPA approach is that the model space must be chosen from a selection of low-energy bands, rather than directly as local degrees of freedom.Specifically, it requires the irreducible polarizability (defined by the reference meanfield) to have no coupling between excitations within the model space and ones in the rest of the system, otherwise U(ω) alone is not sufficient to reproduce χ(ω).Effective local interactions are therefore found via localization after determining the screening of these lowenergy bands.However, this is problematic, missing out on screening within the low-energy bands themselves if the resulting local interactions are subsequently truncated (e.g. to Hubbard/Hund form), and causing further difficulties when it is not possible to fully represent the relevant subspace (e.g.atomic d-shells) in the low energy bands due to significant hybridization with other states [29,44].In addition, the resulting cRPA screened interactions modify the reference mean-field state, resulting in unintended changes to the subspace density and bandstructure [28,42,45].
To motivate our approach, we ask an alternative pertinent question: What physical quantities can we rigorously conserve at the RPA level in a chosen model space, under the constraint that the resulting effective interactions must remain static and two-body?This leads to the development of the 'moment-constrained' RPA (mRPA) approach, which exactly conserves the instantaneous part of the two-point dd-response in the model subspace, as well as the reference state.We argue that conserving physical expectation values in the construction of effective static interactions provides a more rigorous foundation than the widespread static-cRPA approximations.
Moment-constrained RPA:-The RPA can be formulated as a quasi-bosonic eigenvalue problem in the space of particle-hole excitations and de-excitations of a refer-ence state [46], with all blocks being of dimension given by the product of the number of hole and particle states.We define A = ∆ +v, where ∆ is a diagonal matrix of particle-hole excitation energies (provided by a mean-field), defining the poles of Π(ω), with B = v providing the coupling between the excitations and de-excitations via the bare Coulomb interaction, v, in the particle-hole (ph) channel.
The diagonal matrix Ω provides the poles of the RPA dd-response function, χ(ω), with the residues defined by the amplitudes of the excitations and de-excitations, X and Y respectively.These form a biorthogonal set of eigenvectors, providing the relations (X + Y) −1 = (X − Y) T and (X − Y) −1 = (X + Y) T .We can expand χ(ω) as a Laurent series, with its dynamics fully characterized by the moments of its spectral distribution [40,47], as The zeroth moment of the distribution, η (0) , characterizes the instantaneous part of the correlated dd-response, as ⟨(ĉ summed over the same-spin particle-hole (de)excitations, denoted by indices (a, i) and (b, j).This describes the correlated contribution to the two-body reduced density matrix, and all resulting static expectation values at the RPA level [48][49][50].It is this quantity (through to first order) which we aim to rigorously conserve within a chosen model space with our effective interactions.Crucially, this can be achieved while maintaining static renormalized effective interactions in the model space, preserving symmetries, and without changing the reference meanfield in the model space.
The structure of the RPA equations imposes a relation between the first two dd-moments [47], Inserting the definitions of A and B, we find an equation linear in the interaction.This can be analytically inverted to find an interaction kernel that under the RPA gives rise to a desired η (0) and η (1) , Therefore, by substituting all quantities on the RHS of Eq. 6 for their projection of the full RPA into the chosen model space, an effective static model space interaction can be found, U mRPA .This ensures that the model space RPA with v → U mRPA rigorously conserves all full system RPA expectation values derived from η (0) and η (1)  by construction, which includes all instantaneous correlators in the model space.A similar approach to conserve higher-order moments results in equations non-linear in the resulting interaction, and subsequently necessitates a dynamical component to the resulting interaction, as expected to recover the dynamical cRPA limit of the full dd-response [40].Focusing on conservation of just the first two dd-moments enables a fully static model space screened interaction.The model space with U mRPA can also reproduce the subspace contribution to the RPA correlation energy due to its dependence on η (0) , as E RPA corr = 1 2 Tr[η (0) (A + B) − A], exploited in App.D to compute energetic corrections to the model space.Beyond this, the first two spectral moments of the GW self-energy in the subspace are also exactly described with mRPA interactions, indicating a formal conservation of certain correlated one-body properties [47].Symmetries, including spin-independence of the resulting effective interactions are also exactly preserved, derived in App.B. These rigorous properties provide a robust footing for use of U mRPA in subsequent correlated treatments.
Building U mRPA via Eq.6 requires the component of η (0) and ∆ in the subspace.Construction directly from Eqs. 2-3 entails a prohibitive O N 6 scaling, but this is reduced to an efficient O N 4 following the approach of App. A. The only constraint on the choice of subspace that the cluster excitation space is an orthogonal projection of the full system excitation space, i.e. that particleand hole-like character of the subspace orbitals as defined by the reference state is preserved.While this is trivially true where cRPA is valid via a subspace selected from mean-field bands, this is a far looser requirement allowing for a subspace where a non-diagonal Π(ω) couples the subspace to its environment to be considered.This admits direct mRPA screening of arbitrary (e.g.local atomic-like) subspaces by at most doubling their size by conserving the reference mean-field density matrix over the subspace [51,52], or local spaces formed by localizing hole and virtual bands separately before screening.This direct screening of local subspaces therefore includes screening via long-range low-energy band transitions precluded in traditional cRPA, enabling direct application to ab initio quantum embedding clusters.
Interestingly, the resulting U mRPA only screens interactions in the ph channel (which is expected to be the dominant long-range contribution to screening), rather than the full four-point interaction.Other approximations (e.g.T -matrix for pp diagrams [53][54][55][56]) would screen other interaction channels in an analogous formulation, and future work can consider the effect on mRPA interactions from these other channels.It is justified that only the ph channel interaction is screened in mRPA unlike cRPA, since RPA itself is fully determined by this component of the interaction, and mRPA is constructed to describe this RPA physics as opposed to the use of the screening equation (Eq. 1) in cRPA.This feature of only screening the ph interactions, along with the conservation of the reference density, ensures that the occupied reference bandstructure (Hartree-Fock or Kohn-Sham) well as the widely used static limit, and the static-by-construction mRPA interactions.Inset: Error in the RPA moments as an expansion of the dynamical subspace dd-response.The first two mRPA dd-moments are exact by construction, yet also exhibit marginally smaller errors for higher moments compared to the static cRPA interaction.Moment error is computed as the mean squared error of η (n) over the subspace, normalized by the exact ||η (n) || 2 at each order.
is unchanged with mRPA screened interactions.
Results:-We begin by benchmarking the screening of the low-energy Hartree-Fock orbitals of Benzene (an active space of five hole and five particle states).Benzene has previously been considered a paradigmatic example in this context [25], while also small enough to enable comparison to high-level reference results.In Fig. 1 we show the bare Coulomb, dynamical and static-cRPA, and mRPA subspace interactions traced over all channels.The U mRPA is generally less screened than static-cRPA, noting the recent evidence that static-cRPA overscreens interactions [15,28,42].We can consider the accuracy of resulting subspace observables by comparing to coupledcluster (CCSD) [57][58][59].CCSD encodes the wavefunction in T -amplitudes, and we can consider the projection of the full-system T -amplitudes into the active space as an accurate subspace description.This is used to benchmark subspace-only CCSD calculations with the different static interactions [60].We find the mean squared error of the subspace CCSD T 2 -amplitudes to be 1.10 for bare subspace interactions, reduced to 0.25 for static-cRPA and only 0.20 for mRPA, indicating that all ground-state expectation values within the subspace can expected to be more faithfully reproduced with mRPA interactions than the alternatives.
While we expect ground states to be particularly faithful given the mRPA conservation of instantaneous cor- relators, we can also compare the ability of these effective interactions to reproduce the full subspace excitation spectrum.Fig. 1 quantifies this via the errors in the subspace RPA moment expansion of Eq. 3 which fully characterizes the dd-response.While the first two moments are exact for mRPA by construction (and η (1) for static-cRPA), the errors in higher moments are also marginally reduced in their relative error compared to static-cRPA, indicating that the fidelity of the subspace dd-response over all frequencies is at least as accurate as static-cRPA.We now consider the diverse W4-11 test set of 150 molecules, exhibiting a wide range of bonding, radical and correlated physics [61].Moving towards direct application of mRPA to quantum embedding methodologies, we consider fragmenting each molecule into individual atoms comprising a minimal set of their intrinsic valence atomic orbitals (IAOs) [62].These IAO fragments are augmented with an interacting bath of at most the same dimension as the IAO fragment, following the static density matrix embedding (DMET) approach [51,52,[63][64][65][66].This defines multiple small atomic subspaces for each molecule, which can be individually solved via CCSD with either bare or mRPA interactions.The resulting subspace CCSD states are recombined to provide a total energy estimator over the subspaces for each molecule (see Refs. 66 and 67 for more details).These are compared to the energy of the exact full-system CCSD projected into these subspaces.This discrepancy quantifies the ability of the mRPA to account for relaxation of these embedded atomic descriptions due to the neglected interactions with the rest of the molecule.
Figure 2 shows this error aggregated over the test set in increasingly large basis sets, where screening of these lowenergy subspaces by higher-energy scattering dominates.The mRPA interactions screen the fragments defined by the embedding, resulting in a substantial ∼ 66% reduction in both the mean absolute error and standard deviation of the energy error across the data set.This consistent reduction in error across larger basis sets points for these diverse systems attests to the broad applicability of the mRPA interactions for correlated systems.We note that cRPA cannot be easily compared in this context, due to the difficulties in directly screening subspaces that are not mean-field orbitals, as discussed previously.
Finally, we consider the application of mRPA to extended systems where long-range collective plasmons strongly renormalize local properties, which are difficult to describe in a local subspace [25,68].Furthermore, we demonstrate systematic improvability of this screened local subspace.This is achieved in a consistent framework via extending the bath space, formally including additional states which exactly minimize the error of the subspace RPA η (0) with bare interactions, thereby spanning physics at longer length-scales in the model subspaces.The algebraic construction of these additional bath states relies on evaluation of the same quantities as the mRPA interactions, and systematically enlarges the local correlated subspaces in an optimal way to completeness, with details in App. C.This approach is inspired by the unscreened perturbative bath expansion of Ref. 66, but here adapted for a screened embedding.
We converge semi-metallic graphene sheets with a fully ab initio CCSD description, embedding atomic fragments and systematically enlarging the bath space for each fragment.On enlarging, the bath interactions provide increased screening of the fragments explicitly, and the mRPA static screening is correspondingly reduced in a fashion that precludes double-counting of the fragment screening and naturally converges to an in-method (CCSD) exact limit.Figure 3 shows fully ab initio twobody instantaneous charge and spin fluctuations in the local atomic 2p z space, on top of a symmetry-preserving Hartree-Fock reference [67,69].The correlated treatment quantitatively changes these two-body correlators, but bare subspace interactions overestimate the magnitude of the changes from mean-field.Importantly, the bias in these correlators with the bare interactions appears unable to be compensated for with increasing bath, indicating the importance of coupling to truly long-range q → 0 plasmonic modes for an appropriate relaxation of these local properties, impractical in a local embedding with bare interactions.
In contrast, the mRPA screened interactions fold in this coupling, resulting in a rapid and stable convergence with increasing bath size, corroborated with smaller 4×4 k-point meshes where comparison to the full system CCSD result is possible and the same qualitative behaviour is observed.While U mRPA is static, it nevertheless integrates over all RPA diagrams, including plasmonic contributions required to appropriately relax the subspace.This even impacts on the convergence of magnetic correlators in the subspace, despite not being described in the η (0) description of RPA density fluctuations.The results indicate that the plasmonic coupling therefore suppresses the tendency towards formation of atomic magnetic moments in graphene.In App.D we also demonstrate improvements in non-local magnetic fluctuations as well as energetics in this system, and the insensitivity of these results to choice of underlying reference mean-field theory.
In summary, the developed 'mRPA' efficiently inte-grates over all external RPA diagrams to provide a manifestly static effective interaction for low energy or local models, conserving all instantaneous RPA correlators in a chosen subspace by construction.For multi-method workflows, we show this provides a systematic and principled route for the inclusion of long-range and high-energy screening in a local correlation framework.This also opens opportunities in the push for fully improvable and non-empirical quantum embedding, with reliable convergence in extended systems coupled to long-range coherent quasiparticles.moments We consider an efficient approach for the construction of the subspace projection of the RPA quantities η (0) and η (1) .These are defined by their relation in Eq. 5, as required for the subsequent computation of the subspace screened interaction U mRPA via Eq. 6.This is similar to GW theory, or from an iterative approach to the direct construction and inversion of the (I + Q(z)) matrix with local projectors, which will be considered in future work.Furthermore, since the η (1) quantity is diagonal, finding the subspace component of this quantity is trivial.
Finally, we note that this O N 4 formal scaling is unchanged if you want to find the mRPA screened interactions of a single subspace (even one whose size grows linearly with system size) or multiple local subspaces simultaneously, since the full system Q(z) matrix construction and inversion is still performed once, and other contractions do not scale worse than this.This feature is exploited in the main text in the results of Fig. 2, where each atom is considered part of its own screened subspace interactions, and thus the number of subspaces which you need to find the screened interactions for increases linearly with system size.Nevertheless, the formal scaling of finding all mRPA screened subspace interactions is still O N 4 overall.

Appendix B: Spin independence of mRPA interactions
In the case of a spin-conserving reference state where the irreducible RPA polarizability is the same for both spin-channels, and where the interaction kernel is spinindependent (e.g. in the absence of magnetic fields), the mRPA interaction resulting from the above procedure is also spin-independent.This allows for straightforward compatibility with methods to solve the effective Hamiltonian designed for use with the bare Coulomb interaction in 'restricted' formulations, and furthermore allows for efficiency improvements in the construction of this effective interaction.
In such a case, VV T is the same for all spin-conserving transitions, so the Coulomb coupling between two spinconserving excitations depends only on the spatial orbitals.This combined with the fact that ∆ is spinblocked means that when combining Eqs.A8 and A15 all terms containing the expansion of the Coulomb interaction over an auxiliary space intermediate are spinindependent.The only spin-dependent contribution is thus the product of the two non-Coulombic terms in these equations, ∆∆ −1 = I, and all other terms must be spinindependent.
When separated into blocks of single-particle excitations in each (α, β) spin channel the dRPA zeroth ddmoment must as such take the form where η (0) s is a spin-independent, correlation-induced part of the zeroth dd-moment, depending only on the spatial component of the orbital basis.In this case, we expand the inverse of η (0) in the spin-orbital basis via a Neumann series, as We can therefore write the inverse of the zeroth ddmoment compactly as Insertion of this definition into Eq.6 gives the spinindependent mRPA interactions for restricted references as Given that this only requires an inversion of a matrix in Eq.B6 which is half the size of the analogous spinned quantity in Eq. 6, this can result in a more computationally efficient implementation with respect to model space size, where applicable.

Appendix C: RPA interacting bath orbitals
We seek to define a bath expansion for a quantum embedding to systematically expand the single particle space, in an optimally efficient way to describe the physics of the zeroth moment of the dd-response at the RPA level (η (0) ia,jb ).From a viewpoint of the static ex-pectation values, this can be considered as the optimal single-particle bath expansion to span the fragment correlations of the two-body reduced density matrix at the RPA level.While the use of mRPA interactions in the cluster ensures that the subspace η (0) is exactly reproduced at the level of RPA, this is not true when the cluster is solved at a higher level of theory, and thus a bath expansion can provide a route to systematically control the remaining error in this beyond-RPA description of the cluster.
We define a bath expansion which will optimally minimize the error in the cluster η (0) with bare interactions, ensuring that the bath spans these most relevant degrees of freedom at the level of the resulting solver, and in the process minimize the difference between the mRPA and bare interactions throughout the cluster for a given number of bath orbitals.This approach is related to a similar (interacting) bath expansion which was introduced in Ref. [66] which optimally captured the one-body physics from a bare (Møller-Plesset) second-order perturbative level of theory.We expect that a screened perturbative bath expansion based on RPA to be more suitable for extended and polarizable systems, while noting that its construction relies on knowledge of the zeroth moment of the dd-response which is already required to construct the mRPA interactions in the resulting cluster, making the approach particularly straightforward and consistent in this framework.
We start from a cluster formed from an arbitrary set of (generally local) 'fragment' orbitals.These are augmented with the bath space from density matrix embedding theory (DMET), which ensures that the reference single determinant state is exactly spanned by a single state in the cluster, and that the mean-field density matrix of the fragment can be faithfully reproduced [51,52,[63][64][65][66].We note that this construction also ensures the applicability of mRPA directly to this cluster.To obtain the RPA bath expansion orbitals to further augment this initial cluster, we form one-particle/hole pseudo-density matrices representing the correlated hole and particle fluctuations at the level of RPA which couple this initial cluster to its wider environment, Here, ĩ, (ã) indices denote the local hole (particle) states associated with the initial cluster, x, which can be efficiently implemented via projection operators, while other indices represent environmental hole (particle) orbitals in the complementary space [66].Note that all indices represent spin-orbitals, which in spin colinear systems will require that all indices have the same spin label.As such, from this point on all indices will be assumed to represent the same spin channel.The eigenvectors of these matrices represent orbitals of a bath expansion of the correlated environment in occupied and unoccupied states, with their corresponding eigenvalues quantifying their importance in describing the full RPA physics of η (0) in the cluster.These eigenvalues can then be used as a systematic truncation via a threshold in order to define the size of the bath space, analogous to the approach in Ref. 66.The mRPA (or bare) interactions can then be found in this subspace for the high-level solver in a quantum embedding framework.The required η (0) iãjã components can be calculated in O N 4 computational cost based on the approach in App.A, with additional techniques (not pursued here) expected to be able to reduce this to O N 3 via restriction of the full excitation space in the construction of η (0) iãjã , or working with a doubly factorized bare interaction [47].
To provide some further insight into the nature of these bath orbitals, we can explicitly consider the relation between the two-body density matrix, Γ pqrs , and the zeroth moment of the four-point dd-response in an arbitrary basis as an expansion over all neutral excitations, pqrs is the zeroth moment of the dd-response, which is related to the η (0) quantity central to this work via summation over the excitations and de-excitation manifold, We find that the RPA pseudo-density matrices for each cluster x used to define the RPA bath expansion in Eqs.C1-C2 are related to a hypothetical RPA one-body density matrix, γ, as While this does not guarantee positive definiteness of the quantities in Eqs.C1-C2, since the RPA density matrices are not well-defined [75] (hence the use of pseudo-density matrix), it gives some physical justification for such a property.It also demonstrates that the expansion captures the one-body coupling to the environment, despite originating from a purely two-body correlated quantity.

Appendix D: Non-local and energetic expectation values in graphene
In this section, we expand on the results of mRPA and quantum embedding calculations on the graphene system, considering the effect of different reference states, as  [57,58,65,76].However, no double-counting results from the use of Kohn-Sham references in order to define the orbitals and their energies in the irreducible polarizability instead of Hartree-Fock.There may be advantages in this in terms of convergence of the mean-field solution in metallic systems and accuracy of the resulting mRPA interactions in extended systems.However, too much dependence on this starting reference would also be a cause of concern and an indication of also inheriting the empiricism in the choice of mean-field theory.
In Fig. 4 we present analogous calculations to those of Fig. 3 in the main text, showing the local charge and spin correlators with bare and mRPA interactions in the quantum embedding of increasing cluster sizes, but now on a Kohn-Sham reference state with the PBE functional.This can change both the construction of the RPA bath orbitals of App.C defining the single-particle space of the cluster, as well as the mRPA interactions which are used in its subsequent solution.While the mean-field gap is quite different to the Hartree-Fock reference, the convergence of the CCSD local correlators from the clusters for both mRPA and bare interactions is almost identical.This is reassuring that it is genuinely describing the effect of the post-mean-field physics, and inherits the largely reference-independent property of the CCSD solver in the full system limit.Future work will consider a selfconsistent approach to the mRPA interactions, formally 20  removing the mean-field dependence.
We also consider the accuracy of non-local expectation values in a quantum embedding framework with mRPA screened interactions, where the operators span more than one fragment (impurity).One approach is to simply enlarge the fragment space to encompass the range of the operator of interest.However, this is inefficient and will not be suitable for full system expectation values such as total energies.Reference 67 considered approaches to couple the solutions of different embedded quantum clusters for non-local expectation values.We therefore follow this prescription to compute atomic 2p z spin correlators to quantify the fluctuations between two nearest-neighbour atoms (labelled A and B), as ⟨ ŜA z ŜB z ⟩ in Fig. 5, via construction of the two-body reduced density matrix (denoted Γ[(γ, K * )[Ψ x ]] in Ref. 67).We highlight that the atomic fragmentation means that these correlators span different embedded cluster solutions.
These results again show the impact of the mRPA screening in the cluster solutions for this system, even for non-local expectation values which are not explicitly conserved in the mRPA construction.Plasmonic modes over the longest length-scales are folded in to the local interactions, reducing the nearest-neighbour magnetic fluctuations, and enabling high-accuracy and systematic convergence compared to bare interacting clusters.
When extending this analysis to total energies, some more care is required.The magnetic correlators considered extend beyond the length of any single fragment, but the union of all fragments do span the considered operators.This is not always the case however, and in our embedding protocol where fragments are intrinsic atomic orbitals of minimal size, the space of high-energy unoccupied states in large basis sets is not contained in any fragment.While the mRPA interactions will relax and improve the description of the correlated low-energy frag- ment spaces when solved in isolation, they will not account for the incomplete trace in the evaluation of nonlocal expectation values over these high energy degrees of freedom not spanned by any of the solved subspaces.
When considering total energies with variational solvers of fragmented (yet incomplete) subspaces, this incompleteness in the space spanned by the Hamiltonian in the union of the subspaces will serve to increase the total energy of the fragmented system.Similarly, the use of mRPA interactions in the subspaces will serve to also increase the energy of the individual subspace descriptions.Despite the fact that the description of these subspaces with mRPA interactions is closer to the projection of the full system description into this subspace, the use of mRPA interactions will therefore necessarily serve to increase the total energy error.This is shown in Fig. 6, where the correlation energy error of graphene increases with mRPA interactions.This energy over the embedded fragments is defined via the 'projected' energy estimate, introduced in Ref. 66, and expanded on in Ref. 67.It should be noted that this argument of energy error increase with mRPA interactions is not necessarily true when considering energy differences.
We therefore seek a way to improve the convergence of energy estimators in keeping with the spirit of mRPA, to include the energetic contributions from high-energy and long-range contributions that are not contained within any subspace description.This would therefore directly include the energetic effect of e.g.plasmon modes, rather than just the effect of their coupling on the low-energy descriptions.The brute-force approach via an (interacting) bath expansion can be seen in Fig. 6 to be poorly convergent, while previous DMET approaches have considered direct enlargement of fragment spaces to include these energetic contributions, but substantially increases the cost of the subspace solutions [65].
As an alternative approach, we formulate a full-system energetic correction to these incomplete fragmented subspace expectation values via the mRPA construction.We begin by noting that we can construct a global two-body density matrix contribution from the RPA via the relation given in Eq.C3. χ (0) pqrs then relates to the cumulant (or connected) component of the two-body density matrix, K pqrs , as pqsr − γ ps (δ rq − γ rq ). (D1) This represents the correlated, two-body correlations not decomposable as a combination of one-body correlations.The contribution of these irreducible two-body fluctuations to the total energy is given by We therefore obtain an overall correlation energy contribution from these irreducible two-body fluctuations as To use this as a correction to our correlated energy estimate obtained via quantum embedding of our incomplete fragmented wavefunction, we must remove local contributions from each cluster to avoid double counting.This can be achieved simply in our quantum embedding approach by computing a contribution to these two-body energetic fluctuations in each subspace, by projecting one index into the fragment space, and the remaining indices restricted to the subspace of consideration [67].This provides a decomposition of the total RPA two-body energy for a cluster x as where P x is defined as symmetrically projecting one of the occupied indices of the full-system η (0) into the fragment space, with all other indices projected into the cluster space of x (as in Eq. 25 of [67]).
The final non-local energy correction is therefore obtained by deducting the local subspace contributions from the global expectation, giving the result This is equivalent to only including contributions in Eq.D6 where all indices are not contained within any of their cluster space.While this choice of decomposition into subspace contributions to avoid double-counting is specific for our approach for computing expectation values, the general approach could be easily adapted for other embedding or subspace methodologies in order to define a non-local RPA energy correction.
Application of this RPA energetic two-body correction for the incompleteness of the cluster spaces in the atomic IAO quantum embedding of graphene is shown in Fig. 6.This is found to substantially improve the rate of convergence with respect to the size of the explicit CCSD clusters towards the full CCSD limit, with the energy correction accurately compensating for the long-range bubble contributions not spanned by the clusters.Other energy contributions are expected to decay more rapidly, and therefore more easily described by the increasing local cluster spaces.

FIG. 1 :
FIG.1: Trace of the screened and bare interaction, U pq,pq in a (10,10) CAS space of Benzene in a cc-pVDZ basis.Shown are the dynamical cRPA interactions, as well as the widely used static limit, and the static-by-construction mRPA interactions.Inset: Error in the RPA moments as an expansion of the dynamical subspace dd-response.The first two mRPA dd-moments are exact by construction, yet also exhibit marginally smaller errors for higher moments compared to the static cRPA interaction.Moment error is computed as the mean squared error of η (n) over the subspace, normalized by the exact ||η (n) || 2 at each order.

FIG. 2 :
FIG. 2: Distribution of energy errors in atomically fragmented DMET-CCSD with bare or mRPA interactions across the W4-11 test set of 150 molecules, in increasingly large basis sets.Errors are computed comparing to the full-system CCSD of each molecule projected into each subspace.
× 8 Graphene with increasing size of embedded cluster, for both bare Coulomb and mRPA cluster interactions.Inset: Error in these two-body correlators in 4 × 4 k−point meshes, where comparison to exact CCSD (on the PBE reference) is possible.well as the impact of the mRPA screening on non-local expectation values.The results of the main text were performed for a 4 × 4 and 8 × 8 Γ-centered k-point mesh, in a cc-pVDZ basis with a restricted Hartree-Fock reference state, via the Vayesta quantum embedding package, which interfaces with the pyscf codebase