Post-Kohn–Sham Random-Phase Approximation and Correction Terms in the Expectation-Value Coupled-Cluster Formulation

Using expectation-value coupled-cluster theory and many-body perturbation theory (MBPT), we formulate a series of corrections to the post-Kohn–Sham (post-KS) random-phase approximation (RPA) energy. The beyond-RPA terms are of two types: those accounting for the non-Hartree–Fock reference and those introducing the coupled-cluster doubles non-ring contractions. The contributions of the former type, introduced via the semicanonical orbital basis, drastically reduce the binding strength in noncovalent systems. The good accuracy is recovered by the attractive third-order doubles correction referred to as Ec2g. The existing RPA approaches based on KS orbitals neglect most of the proposed corrections but can perform well thanks to error cancellation. The proposed method accounts for every contribution in the state-of-the-art renormalized second-order perturbation theory (rPT2) approach but adds additional terms which initially contribute in the third order of MBPT. The cost of energy evaluation scales as noniterative in the implementation with low-rank tensor decomposition. The numerical tests of the proposed approach demonstrate accurate results for noncovalent dimers of polar molecules and for the challenging many-body noncovalent cluster of CH4···(H2O)20.


INTRODUCTION
The random-phase approximation (RPA), or the ring-diagram approximation, originated in electron gas physics as a method developed by, among others, Gell-Mann and Brueckner to sum up the most divergent correlation energy terms in the highdensity regime. 1,2−6 The RPA response is the central quantity used to formulate the correlation energy approximation via the adiabaticconnection fluctuation−dissipation formula where the integration is over the coupling constant λ and frequency u.−9 At long intermolecular distances, the RPA correlation energy correctly reduces to the Casimir−Polder dispersion energy with correlated response functions of the monomers. 5,10χ RPA accounts for the electrodynamic screening in large polarizable noncovalent complexes, which is a neglected effect in, e.g., second-order Møller−Plesset theory. 11,12−17 Still, there are some obvious limitations of RPA, e.g., the lack of single excitations and the presence of exclusion-principle violating terms.Within DFT, a way to improve RPA is to use an ab initio Kohn−Sham Hamiltonian 18−20 and to combine it with a more sophisticated model of the exchange−correlation kernel in the screening equation. 21,22An alternative, which we consider in this work, is to apply many-body perturbation theory (MBPT) and coupled-cluster techniques to account for the beyond-RPA contributions.The state-of-the-art approach of this type is the renormalized second-order perturbation theory (rPT2) method of Ren et al., 23 where E c RPA is supplemented with the renormalized singles energy, E c rSE , and the second-order screened exchange energy, 24 E c SOSEX .Both E c rSE and E c SOSEX are well-justified in MBPT.The singles energy accounts for a subset of terms that arise for the non-Hartree−Fock reference.The second-order screened exchange (SOSEX) energy removes the second-order Pauli-exclusion violating contributions in E c RPA . 25However, using both of those corrections at the same time leads to a systematic overbinding in hydrogen-bonded systems. 23−28 The goal of our work is to develop a systematic approach where a number of third-and higher-order terms beyond the rPT2 energy are included.The work is organized as follows.We start from the correlation energy expressed as the expectation value of the normal-ordered Hamiltonian, H N , with contributions evaluated directly via the one-electron reduced density matrix (1-RDM) and the two-electron reduced density matrix (2-RDM) cumulant.The density matrices are then parametrized using the RPA doubles amplitudes, T 2 , and the mean-field approximation of 1-RDM.This leads to an expansion where the leading term is the conventional RPA correlation energy, followed by the singles and non-ring doubles corrections.Among the derived terms, we select a subset of contributions with the N ( ) 4 scaling.In the numerical section, we test the proposed approach for noncovalent systems: the dispersion-dominated and hydrogenbonded dimers of the S66 × 8 data set 29 and the many-body cluster of the methane molecule in a dodecahedral water cage. 13 THEORY 2.1.Definitions.We assume a closed-shell system with real orbitals, but an open-shell generalization of the presented equations is possible.The naming of the indices is as follows.
• Occupied (spin)orbitals: i, j, k, l • Virtual (spin)orbitals: a, b, c, d • General (spin)orbitals: p, q, r, s • Eigenvectors of the RPA doubles amplitudes: κ, μ, ν • Tensor hypercontraction vectors 30 (grid points): γ, δ We distinguish between spin orbital and orbital indices by the label above the summation symbol, ∑ spin-orb and ∑ orb .The two electron integrals (pq|rs) are defined as Where applicable, the RPA amplitudes, T ij ab , are treated as matrices with compound indices ai and bj.The κth eigenvector of this matrix is denoted as U ai,κ .The matrix product involving T 2 and, for example, the Coulomb integral matrix V, is defined as The direct-ring approximation of an algebraic formula is understood as the explicit part of which can be expressed with the matrix−matrix multiplication defined in eq 3 or, equivalently, can be graphically represented by ring diagrams with no exchange interactions.The MBPT and coupled-cluster diagrams follow the left−right convention described in ref 31.
The total electronic Hamiltonian, H, is used in the normalordered form.Thus, H is the sum where E HF is the reference mean-field energy and H N is the normal-ordered part of The above decomposition of H N , with the zeroth-order part h 0 and perturbation δh + V N , defines the MBPT contributions.The zeroth-order mean-field Hamiltonian = { } † h h p q pq pq 0 spin orb 0 (7)   has the single-determinant ground state |Ψ 0 ⟩, which defines the split into the occupied and virtual orbital subspaces.The normal ordering of operators with respect to |Ψ 0 ⟩ is indicated by curly braces.h 0 can be any mean-field model.Here, we will consider the Kohn−Sham Hamiltonian in the PBE exchange− correlation approximation 32 and the generalized many-body theory (GMBPT) Hamiltonian 33 based on PBE.The perturbation in eq 6 is composed of the one-electron part and the two-electron part in the form of the normal-ordered Coulomb interaction operator The one-electron perturbation originates from the difference between the matrix elements of the Fock matrix i.e., the one-electron part of H N , and the elements of the reference mean-field Hamiltonian, h 0 .The Fock matrix h is computed using the 1-RDM, ρ 0 , corresponding to |Ψ 0 ⟩.

Decomposition of the Correlation Energy.
Let us consider the correlation energy relative to the reference determinant |Ψ 0 ⟩, defined as the expectation value of the normal-ordered Hamiltonian with the interacting wave function A direct evaluation of eq 12 requires a model of 1-RDM and 2-RDM p r sq qs pr (14)   To arrive at a useful truncation scheme, we split Γ into the antisymmetrized product of 1-RDMs and the size-extensive remainder, Λ, referred to as the 2-RDM cumulant 34−36 = + qs pr pq rs ps rq qs pr (15)   Using the above definitions and Wick's theorem, the complete matrix element of H N on the rhs of eq 12 can be built from the one-particle contributions

Journal of Chemical Theory and Computation
|{ }| | = † p q pq (16)   and the two-particle contributions p r sq pq rs ps rq qs pr (17)   where δρ is the difference between the fully interacting 1-RDM and the diagonal, noninteracting density matrix, ρ 0 , corresponding to the single determinant Combining eqs 16, 17, and 12 yields the decomposition of the total correlation energy into The total 1-RDM contribution, , is composed of the linear component, which originates from the single-particle part of and the quadratic part, which originates from the trace of V with the antisymmetrized product of 1-RDMs  (21)   The cumulant part of the correlation energy, E c Λ , is the trace of V with the cumulant matrix ( ) pqrs qs pr c spin orb (22)   The working equations for E c 1RDM,lin , E c 1RDM,quad , and E c Λ require parametrization of the wave function.In what follows, we will assume that the reduced density matrices originate from the coupled-cluster doubles (CCD) wave function (23)   in which the cluster operator, T 2 , accounts only for the ring (RPA) terms. 37Infinite-order ring formulas for ρ and Λ follow from the expectation-value coupled-cluster theory. 38In addition to the CCD contribution, the effect of single excitations, which is the leading-order contribution to E c 1RDM,lin and E c 1RDM,quad , will be accounted for by using mean-field 1-RDMs built from the eigenvectors of the Fock matrix.We will use the terms "singles correction", "mean-field change energy", and "1-RDM energy" interchangeably for E c 1RDM .

2.3.
Choice of h 0 .In traditional post-KS RPA, the noninteracting system is described by the KS Hamiltonian h KS ; the one-particle perturbation is then where h is the Fock matrix, as in the Hartree−Fock theory, but computed with KS orbitals.However, as noted by Bartlett et al., 33,39,40 δh KS includes large diagonal matrix elements which make finite-order perturbation theory poorly behaved.GMBPT, employed throughout this work, addresses this issue. 39n GMBPT, 39 the occupied−occupied (oo) and virtual− virtual (vv) blocks of the one-particle perturbation are combined into the zeroth-order Hamiltonian defined as 18,33 The GMBPT perturbation, δh = h − h 0 , is limited to the remaining ov and vo blocks of the Fock matrix.
The occupied and virtual orbital subspaces are still defined by the underlying KS Hamiltonian, but the KS orbitals and orbital energies in all MBPT formulas are replaced by their semicanonical counterparts.The semicanonical orbitals, C, and orbital energies, ε, are the eigenvectors and eigenvalues of the oo and vv blocks, respectively, of the Fock matrix In the semicanonical basis, h 0 is diagonal and the one-electron perturbation δh has the form The semicanonical basis has been originally applied by Verma and Bartlett in their implementation of self-consistent RPA. 18rom the MBPT perspective, the justification of using the semicanonical basis for the non-Hartree−Fock mean-field reference is clear.Take for instance the second-order RPA correlation term The semicanonical-orbital E c RPA, (2) is equivalent to the infinite sum of the KS-orbital terms with perturbation δh KS , analogous to diagram V in Table 1 (see also Section 4 of the Supporting Information).Only the first term in the series is included in the traditional RPA approach based on the KS orbitals.From now on, we assume that all of the orbitals and orbital energies are semicanonical.The KS orbitals will be used only in numerical comparisons to obtain the results of the existing RPA approaches.
2.4.RPA Doubles Amplitudes.The RPA doubleexcitation amplitude matrix, T 2 , is a key intermediate replacing the infinite sums of ring contributions in ρ and Λ.An efficient evaluation of T 2 is critical in our approach to avoid additional cost relative to the conventional RPA in the dielectric matrix formulation.
We express the amplitudes in terms of the singlet excitation operators 41

Journal of Chemical Theory and Computation
The closed-shell direct-ring double excitation amplitudes ai bj 2 orb (32)   are the solution of the ring coupled-cluster doubles equation 37 Remarkably, the exact solution of eq 33 can be computed in a single shot with the matrices already available in the Cholesky/frequency-integral RPA implementations.Let us denote the matrix of the vovo Coulomb integrals as As usual, V will be used in a Cholesky-decomposed form 4 where R ai,ξ denotes the Cholesky vectors and N Chol is proportional to the orbital basis-set dimension.Additionally, let us represent the orbital energy differences as a diagonal matrix With those definitions, the matrix form of the doubles equation is The solution of eq 37 assumes the form where the auxiliary matrix is determined from the RPA response matrix, χ RPA , i.e., the solution of the RPA screening equation at full coupling strength The computational form of eq 39 employs the Cholesky decomposition of V and is given by where χ(u) is the closed-shell noninteracting density-response function at frequency u and Π(u) is the auxiliary matrix defined in ref 15: In our implementation, the RPA amplitudes are available only in the form of a set of eigenvalues a κ and eigenvectors U ai,κ : The numerical rank of T 2 is a multiple of the number of the Cholesky vectors of V and, consequently, scales linearly with the system size.The eigenvalues of T 2 satisfy which prevents division by 0 in the matrix formulas for δρ and Λ (see Section 4 of the Supporting Information).We note that the direct-ring amplitudes T 2 lack the positive part of the eigenvalue spectrum, cf.evaluated with the RPA doubles, T 2 , of eq 38.We will derive the CCD correction to the RPA correlation energy in the form of The RPA doubles are easy to compute but accurate only through first order.Fortunately, Wigner's 2n + 1 rule 38,43 applies to eq 46, which means that the first-order T 2 is sufficiently accurate to get the full set of pure-doubles thirdorder terms.Using ρ and Λ in the quadratic CCD approximation and in the infinite-order ring approximation (Section 4 in the Supporting Information), one can show that E c CCD expands into non ring 4th and higher order terms where the first two terms on the rhs are the direct-ring RPA correlation energy, eq 1, or equivalently in the coupled-cluster formulation 37 and the SOSEX correction 24 = | E ai bj T ( ) The remaining explicit contributions, E c 2b , E c 2c , ..., E c 2l , originate from the CCD cumulant and, together with E c SOSEX , form the complete set of third-order beyond-RPA energy corrections of the pure doubles type.The orbital-level formulas are listed in Table 2.The remainder of eq 48 contains non-ring contractions of T 2 and V, where T 2 is at least in the third power.The details of the derivation of eq 48 are described in Section 5 of the Supporting Information.
In the workable variant of ΔE c CCD , we will retain only a small subset of terms from eq 48 to stay within the cost of a conventional RPA calculation.As proposed by Masur    Numerical prefactors related to the permutational symmetry of Λ qs pr are indicated in the first column in front of each summation symbol.The quadratic CCD approximation of the cumulant matrix is denoted as Λ CCD [2] .

Journal of Chemical Theory and Computation
the context of local coupled-cluster calculations, 44 the criterion for the selection of the CCD diagrams can be their decay rate in remote subsystems.It has been concluded in ref 44 that the non-ring diagrams analogous to those in Figure 1 are particularly relevant in noncovalent interactions due to the slow 1/R 6 decay.To account for this contribution in the third order, we include the E c 2g term which also happens to be one of the two least-expensive beyond-RPA terms in eq 48.The other inexpensive term which we include is the SOSEX correction, 24 which decays fast with R but is required to have the complete second order of MBPT.E c SOSEX is well-established in the literature, but here we apply it in the semicanonical basis, as is the case for all correlation energy components.To summarize, is the final variant of the doubles correction to be paired with the 1-RDM (singles) correction in the implemented variant of the beyond-RPA approach.The N ( ) 4 implementation of eq 52 is described in Section 2.7.
2.6.1-RDM Correction.At this point, the proposed approximation of E c still misses the contributions from single excitations.A subset of those contributions in the second and third orders is presented in Table 1 as diagrams I, II, III, and IV.To obtain the complete set of third-order single-excitation terms, one should also add eight additional diagrams with two V operators and a single perturbation δh (see diagrams 5−12 in Figure 14b of ref 31.).We will focus on diagrams I−IV because those can be accounted for simply by using the meanfield approximation of E c 1RDM,lin and E c 1RDM,quad , as follows.Consider the eigenproblem of the full one-electron Hamiltonian where the zeroth order 1-RDM, ρ 0 , in square brackets indicates that the Fock operator is computed using the Kohn−Sham density matrix.The mean-field 1-RDM (54)   includes pure δh terms up to infinite order.Now, inserting the mean-field 1-RDM difference As shown in Section 3 of the Supporting Information, eqs 56 and 57 expand into terms I−IV with the exact MBPT prefactors.The mean-field E c 1RDM,lin is exactly equivalent to the Hartree−Fock singles correction of Klimešet al. 27 Equation 56 also agrees through the third order with the singles correction in the rPT2 method of Ren et al. 23 In the fourth order, E c 1RDM,lin includes additional pure δh diagrams not accounted for in the rPT2 energy, although this difference appears to be numerically insignificant.Contributions III and IV which originate from E c 1RDM,quad are not accounted for in any of the existing RPA singles corrections. 23,27Higher-order correlation effects involving singles can be added by using the coupledcluster expectation-value techniques with some approximation of T 1 , but this is beyond the scope of this work.See Section 3 of the Supporting Information for the details of the MBPT expansion of E c 1RDM,lin and E c 1RDM,quad .

Tensor Hypercontraction of Coulomb Integrals. The cost of E c
SOSEX and E c 2g in their base form is N ( ) 5 but can be reduced to N ( ) 4 with low-rank tensor decomposition.The eigen decomposition of eq 44 decouples the orbital pair indices ai and bj in T ij ab .The tensor hypercontraction (THC) decomposition 30,30,45−47 decouples the indices in the Coulomb integral The matrix Z γδ is the THC grid representation of the Coulomb operator.W pγ is the vector of the pth orbital values on the THC molecular grid. 30,45The THC molecular grid size, N THC , is on the order of a few hundred nodes per atom. 48The computational form of E c SOSEX and E c 2g in terms of the amplitude eigenvectors and THC vectors is where [YXU] γδμ is an intermediate defined as with Y aγ and X iγ representing the W pγ vectors transformed to the molecular orbital basis C pa and C pi are the molecular-orbital coefficients of the ath virtual orbital and the ith occupied orbital, respectively.The matrix elements [YXU] γδμ are formed on the fly and discarded immediately after being used.Equations 59 and 60 contain only a single loop over the eigenvector index μ, which means that the eigenvectors of T 2 can be generated, used to update E c SOSEX and E c 2g and immediately discarded to decrease memory usage.

NUMERICAL RESULTS
The numerical section is designed to assess the differences in noncovalent interaction energies between the proposed approach and rPT2, 23 the state of the art in the quartic-scaling category.The total energy in the implemented variant of our approach gathers the terms defined in eqs 1, 56, 57, 50, and 51 and the reference mean-field energy E HF of eq 5 The total rPT2 energy is 23 The terms enclosed in parentheses are evaluated in the molecular-orbital basis indicated by the lower label.The beyond-RPA MBPT terms accounted for in eq 64 are summarized in Table 1 and in Figure 1.The rPT2 completely lacks third-order diagrams III, IV, V, and VI in Table 1 as well as the diagrams in Figure 1.The renormalized singles energy, E c rSE in eq 65, is numerically almost indistinguishable from E c 1RDM,lin , but formally those terms differ in the fourth order (see Section 3 of the Supporting Information).All beyond-RPA calculations were preceded by the PBE selfconsistent field 32 on a dense molecular grid with 150 radial and 590 spherical points.The correlation energy components E c RPA , E c SOSEX , and E c 2g were obtained within the frozen-core approximation and extrapolated to the basis set limit using the formula of Halkier et al. 50(aug-cc-pVTZ → aug-cc-pVQZ).The frequency integration grid applied in eq 38 was the same as the direct RPA grid in refs 15 and 16.The Cholesky decomposition of the Coulomb integrals with the diagonal elements threshold τ = 10 −7 , defined as given in ref 51, was applied both in the SCF step and in the RPA correlation energy calculation.The tensor hypercontraction decomposition of atomic-orbital Coulomb integrals was implemented according to ref 48 with the THC cutoff threshold ε = 10 −3 .The THC decomposition was applied only in the post-SCF part of the energy.A recalculation of the results for CH 4 •••H 2 O and CH 4 ...(H 2 O) 2 with a tighter THC threshold ϵ = 10 −5 changed the individual dimer and trimer interaction energy contributions by less than 10 −5 kcal/mol.All noncovalent interaction energies were computed with molecular geometries fixed at their interacting cluster values, i.e., without the effect of geometry relaxation.All energies contributing to the n-body (nonadditive) interaction energy were computed in the basis set of the corresponding n-body cluster.

Interactions of Dimers.
We test the interaction energies on a set of small-molecule dimers for which the dispersion-to-electrostatic ratio (disp/elst) ranges from disp/ elst = 0.4 for the methanol-methylamine complex up to disp/ elst = 4.5 for the ethene−pentane dimer.The disp/elst ratio is defined as in ref 29.The systems are a subset of the S66 × 8 database. 29The reference CCSD(T) and CCSD interaction energy data are taken from ref 49.
The disp/elst ratio is the key descriptor that determines the behavior of the RPA methods.In the low-dispersion (polar) systems, the errors of our approach are small, ca.0.2 kcal/mol for methanol-methylamine and only a few hundredths of a kcal/mol for methylamine-N-methylacetamide (Figures 2 and  3).In the moderate dispersion case of benzene−water, disp/ elst = 1.1, a slight underbinding appears, but the error is still small: 0.1 kcal/mol at the equilibrium distance (Figure 4).However, for typical high-dispersion nonpolar systems, the binding strength is visibly underestimated.In the worst case of pyridine−pyridine, the error is 1.1 kcal/mol at equilibrium (Figure 5).Note that in this case, all approximate methods, including CCSD, behave similarly.This indicates the importance of connected triples in the pyridine dimer.For the dimer with the largest disp/elst ratio, ethene−pentane, the interaction energy at equilibrium is underestimated by 0.2 kcal/mol (Figure 6).
The rPT2 method exhibits a different pattern.The binding strength of polar dimers tends to be excessive: the errors for methanol-methylamine and methylamine-N-methylacetamide at equilibrium are 0.5 and 0.3 kcal/mol, respectively.The highdispersion cases are described better.The error for pyridine−  pyridine is 0.8 kcal/mol, which is smaller than the errors in both the CCSD method and our approach.The results for the ethene−pentane dimer match the reference to within a few hundredths of a kcal/mol.
The marked difference in the binding strength between our approach and that of rPT2 is related to the use of the semicanonical basis.As seen in the upper panels of Figures 2−6, the E c RPA contribution to the interaction energy is always less attractive in the semicanonical basis.For example, for the ethene−pentane dimer at the equilibrium distance, the direct RPA contribution to the interaction energy is −3.5 kcal/mol in the KS basis and only −2.7 kcal/mol in the semicanonical basis.Also the semicanonical E c SOSEX term is less attractive than its KS-orbital counterpart in polar dimers.At the equilibrium distance of the methanol-methylamine dimer, the KS-orbital E c SOSEX contributes −0.46 kcal/mol, while its semicanonical counterpart contributes only −0.13 kcal/mol.For nonpolar systems, E c SOSEX is small regardless of the orbitals.The contributions of E c rSE and E c 1RDM,lin would be visually indistinguishable on the scale of the presented plots.Therefore, we draw only E c 1RDM,lin in the upper panels of Figures 2−6.It is clear from our data that E c 1RDM,lin is important for all types of dimers.The quadratic mean-field change term is smaller than E c 1RDM,lin but still significant.E c 1RDM,quad is up to

Journal of Chemical Theory and Computation
−0.3 kcal/mol for the polar dimers and up to −0.2 kcal/mol for moderate-and high-dispersion systems.The behavior of the semicanonical RPA correlation energy resembles, to some extent, the well-known underbinding error of RPA on the self-consistent Hartree−Fock reference.However, in the proposed approach, the third-order correction term E c 2g compensates for diminished attraction.In the case of the ethene−pentane dimer at the equilibrium distance, E c 2g contributes −0.55 kcal/mol, which is one-fourth of the total binding effect.E c 2g is comparable in magnitude to E c 1RDM,lin in the moderate-and high-dispersion cases but smaller than E c 1RDM,lin for the polar dimers.
It is clear from our results that only the semicanonicalorbital RPA can accommodate multiple beyond-RPA attractive terms without severe overestimation of the interaction energy.This should be contrasted with the KS-orbital approach, where the interaction energy curve for hydrogen-bonded systems is too deep when both E c rSE and E c SOSEX are present, but the simpler E c RPA + E c rSE approach is remarkably accurate. 233.3.Methane in a Dodecahedral Water Cage.The binding energy of the methane molecule in a dodecahedral water cage is a finite model of the dispersion-dominated interaction in a methane hydrate. 13,52It is an interesting example of a noncovalent system where all local DFT approximations fail due to the inability to account for the interaction nonadditivity. 13,16In this context, developing lowscaling approaches for the n-body interactions, like the method presented in this work, is desirable given their essential role in the incremental energy calculation in molecular crystals. 53,54In what follows, the total CH 4 ...(H 2 O) 20 interaction energy (66) is decomposed into two-, three-, and four-body contributions The nonadditive interactions of n > 4-body clusters are not considered.E int [2] is the sum of 20 interaction energies of CH 4 ...H 2 O dimers where the pairwise interaction energy is E int [3] is the sum of 190 three-body nonadditive energies of where the nonadditive energy of a trimer is the total binding energy minus the pairwise contributions Finally, E int [4] is the sum of 1140 four-body nonadditive energies of CH where the four-body nonadditivity is defined as the tetramer binding energy minus the pairwise energies and three-body nonadditivities (0, , ) (0, , ) The indices i, j, and k in eqs 68−73 correspond to the water molecules; 0 corresponds to the methane molecule.
Let us start by noting that there are large errors in the selfconsistent PBE model, which is the base for the tested RPA methods.It is evident from Table 3 that PBE not only fails to capture the pairwise interaction but also estimates the nonadditive three-and four-body effects to be nearly an order of magnitude larger than the reference values.ndeed, as the cluster size increases from two to three and finally to four molecules, both E c 1RDM,lin and E c 1RDM,quad alternate their signs but do not decay.The absolute value of the E c 1RDM,lin contribution is over 2 times larger than the total three-body nonadditive energy and over 5 times larger than the total fourbody nonadditivity.The E c 1RDM,quad contribution, which we found relatively small for the individual dimer interaction− energy curves, is now important.For instance, the three-body E c 1RDM,quad term amounts to 0.85 kcal/mol out of the total three-body contribution of 1.36 kcal/mol.
In accordance with the results for the dimer curves, the twobody semicanonical direct-ring contribution to the methanebinding energy is less attractive than that of the KS-orbital counterpart (Table 3).The semicanonical two-body contribution of E c RPA amounts to −8.7 kcal/mol while that of the KSorbital counterpart is as large as −11.5 kcal/mol.This is partially corrected by the E c 2g term, which contributes −1.8 kcal/mol to the two-body binding energy.The pattern remains the same in clusters with more than two molecules, where the semicanonical direct-ring contribution is again smaller, but E c 2g corrects the result.The semicanonical SOSEX contribution is small for all types of clusters.The magnitude of the KS-orbital SOSEX term is larger, but the two-body contribution is nearly canceled by the four-body term, which reduces the total SOSEX contribution to only −0.1 kcal/mol.
After summation over all cluster types, the total methanebinding energy in our approach equals −4.7 kcal/mol, which nearly perfectly matches the reference value.This agreement is the result of a fortuitous cancellation between a slightly too attractive two-body contribution and too repulsive nonadditivity.The rPT2 method overshoots the reference total binding energy by about 0.3 kcal/mol, mainly due to the large overestimation of the two-body contribution.Overall, the dominant trend seen in dispersion-driven dimers holds true for CH 4

CONCLUSIONS
We have shown a systematic beyond-RPA approach where the correlation energy approximation is derived from expectationvalue coupled-cluster theory.The base RPA is the ring part of the coupled-cluster doubles expectation value of the normalordered Hamiltonian.The remaining contributions are the beyond-RPA corrections in the form of (i) the coupled-cluster doubles non-ring correction, ΔE c CCD , evaluated with the RPA amplitudes, T 2 , (ii) and the singles correction, expressed with the mean-field one-electron reduced density matrix.The doubles part, ΔE c CCD , takes advantage of Wigner's 2n + 1 rule and is potentially accurate through the third order, but in practice we restrict it to only two least expensive terms, E c SOSEX and E c 2g .The cost of the resulting correlation energy is noniterative N ( ) 4 .The proposed approach includes all beyond-RPA terms of the state-of-the-art rPT2 method, i.e., the linear singles correction and the second-order screened exchange term but supplements those with novel elements: (i) the semicanonical orbital basis, (ii) the non-ring third-order term E c 2g , (iii) and the quadratic 1-RDM (singles) correction E c 1RDM,quad .The semicanonical basis employed in E c RPA and E c SOSEX accounts for the third-and higher-order terms which arise for non-Hartree−Fock Hamiltonians. 18,33,39While those contributions are justified from the MBPT perspective, they greatly diminish the binding strength relative to the conventional KS basis.Fortunately, this effect is counteracted by the proposed beyond-RPA terms E c 2g and, to a lower extent, E c 1RDM,quad .The semicanonical-orbital basis resolves the incompatibility of SOSEX and the singles correction observed by Ren et al. in polar and hydrogen-bonded systems. 23The SOSEX correction computed with semicanonical orbitals is small and no longer causes overbinding in this type of systems.
The third-order doubles term, E c 2g , accounts for a large fraction of the interaction energy in dispersion-dominated dimers, e.g., E c 2g contributes nearly one-third of the binding effect in the CH 4 •••(H 2 O) 20 cluster.From the perspective of noncovalent interactions, E c 2g is a beyond-RPA contribution of primary importance, with a computational cost that is comparable to the well-established SOSEX correction.
The essential role of the singles correction in post-KS beyond-RPA methods is already recognized in the literature, but we have shown that apart from the usual linear singles correction, the quadratic term, E c 1RDM,quad , becomes important in many-body noncovalent systems.This result emphasizes the need to account for the change of the mean-field energy in systems where the underlying DFT exchange−correlation model poorly describes the nonadditive component of intermolecular interactions.
As a general observation, the proposed approach is accurate in polar and mixed-type systems but underbinds highdispersion complexes.This contrasts with the trend seen in the rPT2 results, where the energies of high-dispersion systems are accurate, while polar dimers are bound too strongly.If considered as an approximation within the proposed framework, the performance of rPT2 can be explained by the presence of systematic error cancellation.The rPT2 method neglects both the attractive term E c 2g and the repulsive contributions which arise in MBPT for the non-Hartree− Fock reference.
Finally, we note that the presented scheme allows for the inclusion of further beyond-RPA terms from Table 2, provided that an efficient implementation is developed.Work is currently underway to include additional exchange terms and to survey various third-order contributions to a broader range of molecular properties.

Figure 1
in ref42.See Section 1 of the Supporting Information for the proofs of eq 38, the bound on the eigenvalues, and the description of the algorithm for the N ( ) 4 numerical evaluation of T 2 .2.5.CCD Correction.Consider the CCD expectation value of the full normal-ordered Hamiltonian

Figure 1 .
Figure 1.Expansion of the third-order beyond-RPA contribution E c 2g .Equation 51 accounts for the correct prefactor of the leading diagram, but the higher-order terms are inexact.

3. 1 .
Computational Details.The beyond-RPA code was added on top of the RPA program described in refs 15 and 16.

Figure 2 .
Figure 2. Interaction energy curve of the methanol-methylamine dimer.The upper panel shows the interaction energy contributions corresponding to the single-point energy terms in eqs 64 and 65.The curve denoted as this work is the sum of E HF , E c 1RDM,lin , E c 1RDM,quad , E c RPA , E c SOSEX , and E c 2g , applied with the semicanonical orbitals.The rPT2 energy 23 is the sum of E HF , E c RPA , E c rSE , and E c SOSEX , applied with the KS orbitals.On the scale of the plot, E c rSE of ref 23 and E c 1RDM,lin would be indistinguishable.The reference CCSD(T) and CCSD curves are taken from ref 49.The disp/elst ratio is taken from ref 29

Figure 3 .
Figure 3. Interaction energy curve of the methylamine-N-methylacetamide dimer.See the caption of Figure 2 for details.

Figure 4 .
Figure 4. Interaction energy curve of the benzene-water dimer.See the caption of Figure 2 for details.

Figure 5 .
Figure 5. Interaction energy curve of the pyridine dimer.See the caption of Figure 2 for details.

Figure 6 .
Figure 6.Interaction energy curve of the ethene−pentane dimer.See the caption of 2 for details.
Additional derivations and spreadsheets with raw numerical data (XLSX) Many-body expansion of the binding energy of CH 4 in dodecahedral water cage (XLSX) RPA double-excitation amplitudes, semicanonical basis, MBPT analysis of E c 1RDM , coupled-cluster expectationvalue theory, beyond-RPA terms in the expectation value of H, and adiabatic connection (PDF) ■ AUTHOR INFORMATION Corresponding Author Marcin Modrzejewski − Faculty of Chemistry, University of Warsaw, Warsaw 02-093, Poland; orcid.org/0000-0001-9979-8355;Email: m.m.modrzejewski@gmail.com Journal of Chemical Theory and Computation

Table 1 .
Subset of Second-and Third-Order Correlation Energy Diagrams Which Include the One-Particle Perturbation ac aThe rightmost column indicates the contributing term in the implemented variant of the correlation energy.bEnergycomponents indicated as the origin contribute terms I−VI with the exact prefactors as implied by the diagram translation rules.cOperatordepicted as "X" is the one-particle perturbation δh KS = h − h KS .Note that in the semicanonical basis, the diagonal blocks of the one-particle perturbation are zero.See Figures14a−dof ref 31 for the complete set of third-order diagrams.

Table 2 .
et al. in Linear and Quadratic Contributions to the Cumulant Part of the CCD Expectation Value of the Hamiltonian a

Table 3 .
Many-Body Expansion of the Binding Energy of CH 4 in Dodecahedral Water Cage CH 4 (H 2 O) 20 All energies are denoted in kilocalories per mol.b Sum of 20 interaction energies of CH 4 ...H 2 O. c Sum of 190 three-body nonadditive energies of CH 4 ...(H 2 O) 2 .Renormalized secondorder perturbation theory of Ren et al. 23 E c rSE denotes the renormalized singles energy defined in ref 23.h Self-consistent DFT calculations with the PBE exchange−correlation model.
a d Sum of 1140 four-body nonadditive energies of CH 4 ...(H 2 O) 3 .e Taken from ref 16. f Singlepoint energy contributions defined in eq 64.g i RPA and beyond-RPA energy components are extrapolated (aug-cc-pVTZ → aug-cc-pVQZ) using the method of Halkier et al.