Method of Moments of Coupled-cluster Equations: Externally Corrected Approaches Employing Configuration Interaction Wave Functions

A new approach to the many-electron correlation problem, termed the method of moments of coupled-cluster equations (MMCC), is further developed and tested. The main idea of the MMCC theory is that of the noniterative energy corrections which, when added to the energies obtained in the standard coupled-cluster calculations, recover the exact (full configuration interaction) energy. The MMCC approximations require that a guess is provided for the electronic wave function of interest. The idea of using simple estimates of the wave function, provided by the inexpensive configuration interaction (CI) methods employing small sets of active orbitals to define higher–than–double excitations, is tested in this work. The CI-corrected MMCC methods are used to study the single bond breaking in HF and the simultaneous breaking of both O–H bonds in H 2 O.


Introduction
One of the most important problems in coupled-cluster (CC) theory [1][2][3][4][5][6][7] is extension of the standard single-reference CC (SRCC) method to quasi-degenerate electronic states. Potential energy surfaces (PESs) involving bond breaking represent a particularly challenging problem in this area.
The main idea of the ground-state MMCC formalism [7,[16][17][18]22] is that of the noniterative energy correction which, when added to the energy obtained in the standard SRCC calculations, such as CCSD or CCSDT (the CC approach with singles, doubles, and triples), recovers the exact (full configuration interaction or full CI) energy. It has been demonstrated that the MMCC formalism allows us to renormalize the existing noniterative SRCC approximations, such as CCSD [T], CCSD(T), CCSD(TQ f ), and CCSDT(Q f ) [12], so that they can correctly describe entire groundstate PESs in situations where the standard arguments based on MBPT, on which the CCSD(T), CCSD(TQ f ), and similar approximations are based, completely fail [7,[16][17][18][19][20][21][22]. In our view, the MMCC formalism represents an interesting idea in the area of new CC methods for molecular PESs that needs to be developed and tested further.
The MMCC energy correction is expressed in terms of the generalized moments of CC equations, i.e., the SRCC equations projected on the excited configurations whose excitation level exceeds that defining a given SRCC approximation. For example, if we would like to improve the results of the CCSD calculations by adding the noniterative MMCC correction to the CCSD energy, we must calculate the generalized moments of the CCSD equations, i.e., the CCSD equations projected on higher-than-doubly excited configurations. In addition, in order to calculate the MMCC energy correction, we have to suggest an approximate form of the electronic wave function |Ψ . In our earlier studies of the ground-state MMCC approach, we focused on very simple guesses for |Ψ , obtained in low-order MBPT calculations. This allowed us to formulate the renormalized and completely renormalized CCSD(T), CCSD(TQ f ), and CCSDT(Q f ) approaches [7,[16][17][18][19][20][21][22]. In this paper, we consider a possibility of approximating the wave function |Ψ that enters the MMCC energy correction by wave functions obtained in inexpensive CI calculations. In analogy to the externally corrected CC methods [6,10,[39][40][41][42][43][44][45][46][47][48][49][50], in which non-CC wave functions are used to provide information about higher-than-doubly excited clusters, we call the resulting MMCC schemes the externally corrected MMCC methods. The CI-corrected MMCC approaches discussed in this work are tested in pilot calculations for the potential energy curve of HF and simultaneous breaking of both O-H bonds in H 2 O.

The Method of Moments of Coupled-Cluster Equations: An Overview of the Ground-State Formalism
In the SRCC theory, we represent the ground-state wave function of an N -electron system, described by the Hamiltonian H, in the following way: where T is the cluster operator and |Φ is an independent-particle-model reference configuration (usually, the Hartree-Fock determinant). In the exact theory, T is a sum of all many-body cluster components that can be written for a given N -electron system. In the standard SRCC approximations, the many-body expansion of cluster operator T is truncated at some, conveniently chosen, excitation level. Thus, if A represents the standard SRCC approximation and if m A (m A < N) is the corresponding excitation level that defines method A, the formula for the cluster operator T (A) used in calculations with method A is as follows: where T n , n = 1, . . . , m A , are the many-body components of T (A) . The CCSD method is obtained by setting m A = 2; in the CCSDT method, m A = 3, etc. The cluster operator T (A) of method A is obtained by solving the standard SRCC equations, where subscript C designates the connected part of the corresponding operator expression and Q (A) is the projection operator onto the subspace of all excited configurations described by T (A) , i.e., with Q n representing the projection operator onto the subspace of n-tuply excited configurations relative to reference |Φ . Once the system of equations, Eq. (3), is solved for T (A) , the energy is calculated using the well-known formula The main idea of the ground-state MMCC theory is that of the noniterative energy correction which, when added to energy E (A) , obtained in the standard CC calculations with method A, recovers the corresponding exact (full CI) energy E. The main purpose of the approximate MMCC calculations is to estimate correction δ (A) , so that the resulting MMCC energy, is close to the corresponding exact energy E.
We have recently demonstrated that correction δ (A) can be expressed in terms of the generalized moments of the SRCC equations, i.e., the SRCC equations projected on the excited configurations whose excitation level exceeds that defining a given SRCC approximation [7,[16][17][18]22]. The precise mathematical definition of the generalized moments of the SRCC equations is as follows: where |Φ (j) J designate the j-tuply excited configurations relative to |Φ . The generalized moments M . We obtain [7,[16][17][18]22], where represents the (n − j)-body component of the wave operator e T (A) , defining the SRCC approximation A, and All this means that if we, for example, want to correct the results of the CCSD calculations (the m A = 2 case) and recover the full CI energy by adding the noniterative correction δ (A) to the CCSD energy, we have to calculate the generalized moments of the CCSD equations, i.e., the CCSD equations projected on triply, quadruply, pentuply, and hextuply excited configurations, or, symbolically, the quantities where T 1 and T 2 are the singly and doubly excited clusters resulting from the CCSD calculations and j = 3 − 6 (M j (2)|Φ = 0 for j > 6). The original proof of Eq. (9) has been based on the Fundamental Theorem of the Formalism of β-Nested Equations [7], which describes mathematical relationships between multiple solutions of the SRCC equations representing different levels of CC theory (CCSD, CCSDT, etc.). An elementary derivation of Eq. (9), based on a simple asymmetric energy expression, termed the MMCC functional, i.e., introduced in Ref. 16, has been given in Appendix A of Ref. 16. The possibility of applying the MMCC functional in direct calculations of noniterative corrections to standard SRCC energies has recently been examined in Refs. 73, 74.
Equation (9) is the basic equation of the ground-state MMCC formalism. There are two main elements of Eq. (9): the aforementioned generalized moments of CC equations and the wave function |Ψ , which in the exact MMCC theory represents the exact ground-state wave function. The generalized moments of CC equations can easily be determined, once we know the corresponding cluster operator T (A) . The remaining issue is what do we do with |Ψ . This issue is addressed in the next subsection.

The MMCC(m A , m B ) Approximations. The CI-Corrected MMCC(2,3) and MMCC(2,4) Methods
The MMCC formula for the energy correction δ (A) , Eq. (9), is expressed in terms of the exact ground-state wave function |Ψ , which we usually do not know (if we knew the exact |Ψ , we would not have to perform any calculations !). In approximate MMCC methods, designated as the MMCC(m A , m B ) schemes, wave functions |Ψ are evaluated by performing, for example, limited CI calculations.
In the MMCC(m A , m B ) approximations, we use wave functions |Ψ , whose CI expansions do not contain higher-than-m B -tuply excited components (m B < N). This assumption reduces the summation over n in Eq. (9) to are calculated according to the following formula [7,[16][17][18]22]: where When m B = N and when |Ψ is exact, we obtain the exact MMCC theory described in the previous subsection. Clearly, we must assume that m B > m A to obtain a nonzero value of δ(m A , m B ). The lowest-order approximations belonging to the MMCC(m A , m B ) hierarchy are the MMCC(2,3) and MMCC(2,4) approaches [7,[16][17][18]22]. In the MMCC(2,3) and MMCC(2,4) methods, we use corrections δ(2, 3) and δ (2,4) to improve the results of the CCSD calculations.
The MMCC(2,3) energy expression has the following form [7,[16][17][18]22]: where E CCSD is the CCSD energy and with T 1 and T 2 representing the singly and doubly excited clusters obtained in the CCSD calculations. According to Eq. (11), the M 3 (2)|Φ quantities entering Eq. (17) can be expressed in terms of the projections of the CCSD equations on triply excited configurations. We obtain, where with |Φ abc ijk representing the triexcited configurations relative to |Φ . Here and elsewhere in the present paper, we use the standard notation in which i, j, k, l, . . . represent the spin-orbitals occupied in |Φ and a, b, c, d, . . . are the unoccupied spin-orbitals. In the MMCC(2,3) approach, we assume that the CI expansion of the wave function |Ψ that enters Eq. (17) does not contain higher-than-triexcited components.
The MMCC(2,4) method represents a natural extension of the MMCC(2,3) scheme in which, in addition to the triexcited moments of the CCSD equations, M abc ijk (2), Eq. (19), we also consider the quadruply excited moments, where |Φ abcd ijkl are the quadruply excited configurations. In the MMCC(2,4) approximation, we assume that the CI expansion of the wave function |Ψ does not contain higher-than-quadruply excited components. The MMCC(2,4) energy is calculated as follows [7,16,17,22]: where The M 3 (2)|Φ quantities are the same as in the MMCC(2,3) approximation [cf. Eq. (18)]. The M 4 (2)|Φ quantities can be expressed in terms of the quadruply excited moments of the CCSD equations defined by Eq. (20). We have, Different types of the MMCC(2,3) and MMCC(2,4) approximations are obtained by making different choices for |Ψ in Eqs. (17) and (22) [7,[16][17][18]22]. We can, for example, choose the loworder MBPT expressions to represent |Ψ , obtaining the so-called renormalized and completely renormalized CCSD[T], CCSD(T), and CCSD(TQ) methods mentioned in the Introduction. As shown in Refs. 7, 16-22, these new approaches represent promising computational tools that remove the pervasive failing of the standard CCSD[T], CCSD(T), and CCSD(TQ f ) approximations at larger internuclear separations, while preserving the "black-box" character of the noniterative SRCC methods. Slightly modified versions of the original renormalized and completely renormalized CC methods, introduced by us in Refs. 7, 16, 17, have recently been examined in Ref. 75. The analytic gradients for the renormalized CCSD(T) approach have been developed by Stanton et al. [76].
We can also think of using the limited CI methods to define wave functions |Ψ in Eqs. (17) and (22). An interesting possibility (suggested, for the first time, in Ref. 7 and examined, in detail, in this paper) is offered by the active-space CISDt and CISDtq approximations, in which higher-than-doubly excited components of |Ψ are defined through active orbitals. We use the CISDt approach to calculate wave function |Ψ for the MMCC(2,3) calculations, while the CISDtq method is used to determine |Ψ for the MMCC(2,4) calculations.
In order to calculate the CISDt and CISDtq wave functions, we divide the available spin-orbitals into core spin-orbitals (i, j, k, l, . . . ), active spin-orbitals occupied in |Φ (I, J, K, L, . . . ), active spin-orbitals unoccupied in |Φ (A, B, C, D, . . . ), and virtual spin-orbitals (a, b, c, d, . . . ). Once active orbitals are selected, we define the CISDt and CISDtq wave functions, |Ψ CISDt and |Ψ CISDtq , respectively, as follows [7,71,72]: where C 0 |Φ , C 1 |Φ , and C 2 |Φ are the reference, singly excited, and doubly excited components of |Ψ CISDt and |Ψ CISDtq and Thus, in the CISDt method, we construct wave function |Ψ by including all singles and doubles from |Φ and a relatively small set of internal and semi-internal triples containing at least one active occupied and at least one active unoccupied spin-orbital indices. In the CISDtq approach, we also include a relatively small set of quadruples containing at least two active occupied and at least two active unoccupied spin-orbital indices [see Eqs. (26) and (27)]. As in all CI approaches, the CI expansion coefficients defining |Ψ CISDt and |Ψ CISDtq are determined variationally. Choosing the CISDt and CISDtq wave functions |Ψ in Eqs. (17) and (22) should be particularly useful in the context of studies of bond breaking. The CISDt approach should provide us with a qualitatively correct description of the single bond breaking (as, e.g., in HF), whereas the CISDtq method seems appropriate for a double bond breaking or for the simultaneous breaking of two single bonds (e.g., two O-H bonds in H 2 O).
It should be noted that the CISDt and CISDtq approaches, which are used to generate wave functions |Ψ for the MMCC(2,3) and MMCC (2,4) calculations, can be regarded as the CI analogs of the active-space CC approaches, such as CCSDt [19,61,62] or SSCCSD(T) [14,[51][52][53][54][55][56][57][58][59][60] and CCSDtq [61] or SSCCSD(TQ) [54][55][56][57], which are known to provide excellent description of bond breaking [14,19,57,61,62]. Although the numbers of triples and quadruples used in the CISDt and CISDtq calculations are identical to the numbers of triples and quadruples used in the CCSDt and CCSDtq calculations (provided that we use the same active space in all these calculations) and although the most expensive N o N u n 2 o n 4 u and N 2 o N 2 u n 2 o n 4 u steps of the CISDt and CISDtq calculations are of the same type as the most expensive steps characterizing the CCSDt and CCSDtq methods, the CISDt and CISDtq approaches are less expensive than the CCSDt and CCSDtq methods. As is always the case, when we are comparing the CI and CC approaches, the prefactors defining the N o N u n 2 o n 4 u and N 2 o N 2 u n 2 o n 4 u scalings of the CISDt and CISDtq methods are smaller than the prefactors defining the analogous scalings of the CCSDt and CCSDtq approaches. The absence of nonlinear terms in the CISDt and CISDtq equations and their presence in the CCSDt and CCSDtq equations also contribute to the lower cost of the CISDt and CISDtq calculations, when compared with the corresponding CCSDt and CCSDtq calculations. In consequence, the CISDtbased MMCC(2,3) method and the CISDtq-based MMCC(2,4) approach are less expensive than their CCSDt and CCSDtq counterparts. It is also very easy to adapt the CISDt and CISDtq equations to the spin and spatial symmetries (a difficult thing to do in the CCSDt, CCSDtq, and, as a matter of fact, many CC approximations) and our considerable experience with the CISDt/CISDtq and CCSDt/CCSDtq calculations indicates that it is easier to converge the CISDt/CISDtq equations compared to the CCSDt/CCSDtq equations at larger internuclear separations (although a considerable progress has been made in improving the convergence of the CCSDt/CCSDtq equations [61,62]). It is true that the CISDt and CISDtq methods are much less accurate than the CCSDt and CCSDtq approaches. However, the CISDt and CISDtq methods are sufficiently good to provide us with relatively inexpensive and yet quite reasonable choices of spin-and symmetryadapted wave functions |Ψ for the MMCC(2,3) and MMCC(2,4) calculations. Another advantage of using the CISDt and CISDtq methods (not examined in this work, but also important) is an easy access to excited electronic states, which can be subsequently used to construct the highly successful MMCC(2,3) and MMCC (2,4) corrections to excited-state energies [22,71,72]. Thus, the CISDt-based MMCC(2,3) approach and the CISDtq-based MMCC(2,4) method, considered in this work, represent an interesting new way of incorporating the triples and quadruples effects in the SRCC calculations, which might be viewed as a very useful alternative to the existing active-space CCSDt and CCSDtq methods. The CI-corrected MMCC(2,3) and MMCC (2,4) approaches are also the alternatives to the recently proposed renormalized and completely renormalized CCSD(T) and CCSD(TQ) approximations [7,[16][17][18][19][20][21][22]. Finally, we can view the CI-corrected MMCC(2,3) and MMCC(2,4) methods as noniterative approximations to the reduced MRCCSD (RMRCCSD) method of Refs. 43-49. In the RMRCCSD approach, the ground-state multi-reference CI (MRCI) wave function is used to extract information about triply and quadruply excited clusters of the SRCC theory, which are, in turn, employed to correct the CCSD equations. In the CI-corrected MMCC(2,3) and MMCC(2,4) approaches, we solve the standard CCSD equations in the absence of terms containing T 3 and T 4 clusters, and use the triply and quadruply excited moments of the CCSD equations and simple, MRCI-like, CISDt and CISDtq wave functions to construct noniterative a posteriori corrections δ(2, 3) and δ (2.4), which are subsequently added to the CCSD energy.

The Remaining Computational Details
The CI-corrected MMCC(2,3) and MMCC(2,4) methods have been implemented within our system of CC programs. The required CCSD calculations were performed with the orthogonally spin-adapted [10,[77][78][79][80] CCSD code described elsewhere [35]. The triply and quadruply excited moments of the CCSD equations, used in the CISDt-based MMCC(2,3) and CISDtq-based MMCC(2,4) calculations, M abC Ijk (2) and M abCD IJkl (2), respectively, were calculated with the new code obtained by modifying the computer program described in Ref. 16. The CISDt and CIS-Dtq calculations were performed with the computer programs described elsewhere [7,71,72]. All our programs are interfaced with the gamess code [81], which is used to generate the restricted Hartree-Fock (RHF) molecular orbitals (MOs) and to transform atomic integrals to MO basis.
The correctness of the CISDt-based MMCC(2,3) and CISDtq-based MMCC(2,4) codes was tested in various ways. The most essential series of tests were the MMCC(2,4) calculations for a few four-electron systems, using the full CI wave functions |Ψ in Eq. (22)

Numerical Examples
In order to illustrate the performance of the CI-corrected MMCC approaches described in Sec. 2 in studies of bond breaking, we applied the MMCC(2,3) and MMCC(2,4) methods to potential energy curves of HF and H 2 O. As in our original studies of the renormalized and completely renormalized CCSD(T) and CCSD(TQ) methods [16] (see, also, Ref. 7), we used a double zeta (DZ) basis set [82], for which the exact, full CI energies [14,83,84] and many other useful results, including the full CCSDT and CCSDTQ energies [14,85,86] and their standard and renormalized CCSD(T) and CCSD(TQ) analogs [16], are available. In all calculations reported in this work, the ground-state RHF determinant was used as a reference.
In our discussion, we compare the MMCC(2,3) and MMCC(2,4) results with the exact, full CI results and a variety of CC results. We also compare the MMCC(2,3) and MMCC(2,4) results with the results of the CISDt and CISDtq calculations, which provide wave functions |Ψ for computing corrections δ(2, 3) and δ (2,4). The purpose of the latter comparison is to demonstrate how much the results improve when the relatively poor (at best, qualitative) CISDt and CISDtq wave functions are inserted into the MMCC(2,3) and MMCC(2,4) energy expressions.

The HF Molecule
The results of our CI-based MMCC calculations for the potential energy curve of HF are shown in Table 1. In this case, there is a 1.634 millihartree difference between the CCSD and full CI energies at the equilibrium geometry, R = R e (R is the H-F internuclear separation), which increases to 12.291 millihartree at R = 5R e (for all practical purposes, R = 5R e can be regarded as a dissociation limit). The large differences between the CCSD and full CI energies at larger values of R are primarily caused by the absence of the connected T 3 clusters in the CCSD wave function. Indeed, the full CCSDT method, which includes these clusters, reduces large errors in the CCSD results, relative to full CI, to as little as 0.173 millihartree at R = R e and 0.431 millihartree at R = 5R e .
The The completely renormalized CCSD(T) and CCSD(TQ) approaches of Refs. 7, 16-18 (cf., also, Refs. [19][20][21][22], designated as the CR-CCSD(T) and CR-CCSD(TQ) methods, considerably improve the results, eliminating, in particular, the unphysical humps on the CCSD(T) and CCSD(TQ f ) potential energy curves. The CR-CCSD(T) and CR-CCSD(TQ) methods reduce the 53.183 and 35.078 millihartree errors in the CCSD(T) and CCSD(TQ f ) results at R = 5R e to 1.650 and 0.454 millihartree, respectively [16] (see Table 1). The errors in the CR-CCSD(T) energies, relative to full CI, do not exceed 2 millihartree over the entire range of R values. The errors in the CR-CCSD(TQ) results are even smaller (they do not exceed 0.5 millihartree; cf. Table 1). This clearly demonstrates that the MMCC theory, on which the CR-CCSD(T) and CR-CCSD(TQ) approaches are based, represents a powerful new formalism, which is capable of removing the pervasive failing of the standard CCSD(T) and CCSD(TQ f ) approximations at larger internuclear separations.
The CI-based MMCC(2,3) and MMCC(2,4) results shown in Table 1 correspond a choice of the highest occupied and lowest unoccupied σ orbitals, 3σ and 4σ, respectively, and the valence π orbitals as active orbitals in the related CISDt and CISDtq calculations. This is a natural choice of active space for the description of bond breaking in HF, since for larger internuclear separations R the full CI wave function of HF is dominated by the ground-state RHF configuration, the doubly excited configuration, corresponding to the (3σ) 2 → (4σ) 2 excitation, and the 3σ → 4σ singly excited configuration. In principle, we could consider a smaller active space, consisting of only two valence σ orbitals [7], but we decided to use a somewhat more complete description here, in which all MOs that correlate with the 1s and 2p shells of the H and F atoms, respectively, are included in active space.  [17]). g The active space consisted of the 3σ, 1π, 2π, and 4σ orbitals.
The above choice of active space guarantees that the description of the potential energy curve of HF by the CISDt and CISDtq methods is qualitatively correct. However, the CISDt and CISDtq results (particularly, the former ones) are quantitatively rather poor. This remark applies to all internuclear separations. Indeed, at the equilibrium geometry, both methods give the > 5 millihartree errors relative to full CI and the situation only worsens as we approach the dissociation region. For example, the 5.783 millihartree error in the CISDt result at R = R e increases to 33.627 millihartree at R = 5R e . In the case of the CISDtq method, the error increase is less dramatic, but the 7.468 millihartree error in the CISDtq result at R = 5R e is still relatively large.
In spite of the relatively poor performance of the CISDt and CISDtq methods and in spite of the large errors in the CCSD results at larger R values, the CISDt-based MMCC(2,3) results and their CISDtq-based MMCC (2,4) analogs are very good. The errors in the MMCC(2,3) results vary between 1.195 millihartree at R = R e and 3.669 millihartree at R = 3R e . Although these results are somewhat worse than the results of the CR-CCSD(T) calculations, they are much better than the results of the CISDt and CCSD calculations, which are used to construct corrections δ (2,3). The fact that we can use an inexpensive CISDt method to construct correction δ(2, 3) and reduce in this way the 33.627 and 12.291 millihartree errors in the CISDt and CCSD results at R = 5R e to 3.255 millihartree is very encouraging. We are not gaining a lot in this case by performing the more expensive MMCC(2,4) calculations (which is a consequence of the fact that the H-F bond is a single bond), but it is good to see that the MMCC theory behaves in a systematic manner. The use of the better CISDtq wave function in the MMCC [MMCC(2,4)] calculations gives a slightly better description of the potential energy curve of HF than that provided by the CISDt-based MMCC(2,3) approach. The MMCC(2,4) curve is virtually parallel to the exact (full CI) curve. The small (∼ 1 millihartree) errors in the MMCC(2,4) results in the equilibrium region slowly (and monotonically) increase with R, reaching a maximum value of ∼ 3 millihartree in the dissociation region.
The fact that the CI-based MMCC methods and their MBPT-based CR-CCSD(T) and CR-CCSD(TQ) analogs provide comparable improvements in the CCSD results for the potential energy curve of HF indicates that the MMCC theory is a robust formalism, in which very good results can be obtained independent of the quality of the wave function |Ψ used to construct correction δ (A) . As in the case of the CR-CCSD(T) and CR-CCSD(TQ) methods [16], one of the main reasons of the excellent performance of the CI-based MMCC approaches at large internuclear separations, in spite of the relatively poor description of the potential energy curve of HF by the CISDt and CISDtq methods, is the presence of the Ψ|e T 1 +T 2 |Φ denominators in the MMCC energy expressions. These denominators increase their values from ∼ 1.0 at R = R e to 2.3-2.4 at R = 5R e , damping the corrections due to triples [in the MMCC(2,3) case] or triples and quadruples [in the MMCC(2,4) case], which are considerably overestimated by the traditional CCSD (T) and CCSD(TQ f ) approaches and would be overestimated by the MMCC(2,3) and MMCC (2,4) methods if the Ψ|e T 1 +T 2 |Φ denominators were not present in the MMCC expressions.

The H 2 O Molecule
The CI-based MMCC(2,3) and MMCC(2,4) results for the water molecule are shown in Table 2 Table 2; notice the change of sign of the T 4 energy contribution, compared to the R = 1.5R e case).
It is difficult to describe the large T 3 and T 4 effects with the noniterative SRCC approximations. Indeed, the quasi-degenerate nature of the ground-state wave function of H 2 O for larger values of R leads to a complete failure of the CCSD(T) and CCSD(TQ f ) methods. At R = 2R e , the unsigned errors in the CCSD(T) and CCSD(TQ f ) results, relative to full CI, are 7.699 and 5.914 millihartree, respectively. Even the full CCSDT approach gives a negative, −2.211 millihartree, error at R = 2R e , which might be an indication of the beginning of the breakdown of the CCSDT approximation (which lacks important T 4 clusters) at larger R distances. The complete incorporation of T 4 contributions via the full CCSDTQ approach considerably improves this situation, reducing the 2.211 millihartree error in the CCSDT result at R = 2R e to 0.108 millihartree (see Table 2).
The performance of the MMCC-based CR-CCSD(T) and CR-CCSD(TQ) methods, which use the second-order MBPT expressions to represent wave functions |Ψ in the MMCC(2,3) and MMCC(2,4) energy formulas [7,[16][17][18]22], is excellent. The CR-CCSD(T) and CR-CCSD(TQ) methods reduce the 7.699 and 5.914 millihartree errors in the CCSD(T) and CCSD(TQ f ) results at R = 2R e to 1.830 and 1.461 millihartree, respectively. At the same time, the CR-CCSD(T) and CR-CCSD(TQ) results at R = R e are virtually identical to the very good results obtained with the standard CCSD(T) and CCSD(TQ f ) approaches. As pointed out in Ref. 16, the CR-CCSD(TQ) method provides very good estimates of T 4 effects for R ≤ 1.5R e . It is interesting to examine if the estimates of T 4 effects provided by the CISDtq-based MMCC (2,4) approach are as good as those obtained with the CR-CCSD(TQ) method.
Two different choices of active orbitals for the CISDt-based MMCC(2,3) and CISDtq-based MMCC (2,4) calculations are investigated in this paper (see Table 2). The natural choice of and singly and doubly excited configurations relative to |Φ involving the 3a 1 , 1b 2 , 4a 1 , and 2b 2 orbitals. The somewhat better choice of active orbitals, which guarantees a more uniform description of the equilibrium and bond breaking regions is provided by the 1b 1 , 3a 1 , 1b 2 , 4a 1 , 2b 1 , and 2b 2 orbitals [57]. Both choices of active orbitals for the MMCC(2,3) and MMCC (2,4) calculations and the related CISDt and CISDtq calculations are considered in this work. As shown in Table 2, the results of the CISDt-based MMCC(2,3) calculations and the CISDtqbased MMCC(2,4) calculations, employing the above two active spaces, are excellent. In view of the small computer effort associated with constructing the MMCC(2,3) correction δ (2,3), the CISDt-based MMCC(2,3) results are particularly impressive. As can be seen in Table 2, the MMCC(2,3) method is capable of reducing the huge, ∼ 50 millihartree, errors in the CISDt energies and the sizable, 9.333 millihartree, error in the result of the CCSD calculations at R = 2R e to less than 2 millihartree. Similar error reductions are observed at R = 1.5R e , where the ∼ 19 millihartree errors in the CISDt results are reduced to 2-3 millihartree, when the CISDt wave function is inserted into the MMCC(2,3) energy expression. The CISDtq-based MMCC(2,4) results at R = 2R e seem to be slightly worse than the corresponding CISDt-based MMCC(2,3) results, but the MMCC(2,4) approach offers a more balanced description of the simultaneous breaking of both O-H bonds in H 2 O. Indeed, the small, 0.5-1 millihartree, errors in the MMCC(2,4) results at R = R e slowly and monotonically increase with R, whereas the errors in the MMCC(2,3) results initially increase, as we go from the R = R e region to R = 1.5R e , and then decrease, as we approach the R > 1.5R e region, which might be the first sign of the breakdown of the MMCC(2,3) approximation at very large distances R (notice that a very similar pattern is observed, when we compare the CR-CCSD(T) and CR-CCSD(TQ) energies; let us recall that the CR-CCSD(T) method is an example of the MMCC(2,3) approximation and that the CR-CCSD(TQ) approach is a special case of the MMCC(2,4) approximation [7,16,17]). The small, 0.501, 0.942, and 2.416 millihartree, errors in the CISDtq-based MMCC(2,4) results at R = R e , 1.5R e , and 2R e , respectively, obtained with the active space consisting of only 6 valence orbitals, and the fact that the MMCC(2,4) and full CI potential energy curves are nearly parallel are the very encouraging findings from the point of view of future applications of the CI-based MMCC methods to PESs involving bond breaking. Although it is possible that the CISDt-based MMCC(2,3) approximation eventually breaks down for R 2R e , the fact that this simple and inexpensive approach provides the 1-3 millihartree errors for all R values ranging between R e and 2R e is very encouraging, too.
The fact that the CISDtq-based MMCC(2,4) approach provides excellent description of the simultaneous breaking of both O-H bonds in water is, at least to some extent, a consequence of the relatively good description of this bond breaking by the CISDtq method (see Table 2). The incorporation of the quadruply excited moments of the CCSD equations in the MMCC(2,4) calculations is responsible for some improvements, too, particularly in the region of the intermediate R values (the R ≈ 1.5R e region). However, bulk of the improvement in the CCSD results by the MMCC theory is already achieved at the lowest MMCC(2,3) level, which uses a relatively poor CISDt wave function to construct the relevant correction δ (2,3). This is mostly related to the dominant role of T 3 clusters in describing the O-H bond breaking in H 2 O [these clusters are already included in approximate manner in the MMCC(2,3) calculations] and, as in the case of calculations for HF, to the presence of the Ψ|e T 1 +T 2 |Φ denominators in the MMCC energy expressions. These denominators increase their values from ∼ 1.0 at R = R e to 1.5-1.7 at R = 2R e , damping the large negative corrections due to triples that are grossly overestimated by the traditional noniterative approaches, such as CCSD(T). Without the presence of the Ψ|e T 1 +T 2 |Φ denominator in the formula for correction δ (A) , the MMCC(2,3) and MMCC(2,4) results would be much worse.
We have already mentioned that the CI-based MMCC approaches describe bulk of the T 3 effects. Interestingly enough, the MMCC(2,4) approximation employing the CISDtq wave function |Ψ seems to provide a very reasonable description of the T 4 cluster components. The "true" values of the T 4 contributions to the energy, as measured by the differences between the CCSDTQ and CCSDT energies, are −0.419 millihartree at R = R e , −1.332 millihartree at R = 1.5R e , and 2.319 millihartree at R = 2R e . The corresponding differences between the MMCC(2,4) and MMCC(2,3) energies for the active space consisting of the 1b 1 , 3a 1 , 1b 2 , 4a 1 , 2b 1 , and 2b 2 orbitals are −0.310, −1.465, and 0.785 millihartree, respectively, in reasonable agreement with the values obtained by forming the differences between the CCSDTQ and CCSDT energies. Remarkably enough, the change of sign of the T 4 energy contribution at R = 2R e is correctly reproduced by the MMCC(2,4) theory. The CISDtq-based MMCC(2,4) approach is better in this regard than the MBPT-based CR-CCSD(TQ) approximation, which gives the negative value of the T 4 effect at R = 2R e , if we subtract the CR-CCSD(T) energy from the CR-CCSD(TQ) energy (see Table 2). On the other hand, the CR-CCSD(TQ) results are very good for all values of R (including R = 2R e ), so that one should not view this particular feature of the renormalized CC methods as detrimental to the high quality of the results that these methods offer. We must not forget that the second-order MBPT wave functions |Ψ , used to construct the noniterative triples and quadruples corrections in the renormalized CC theories, such as CR-CCSD(TQ), are much worse in the region of larger R values than the CISDtq wave functions used in the CI-based MMCC(2,4) calculations.
Finally, it is very encouraging to observe the systematic improvements in the MMCC(2,3) and MMCC (2,4) results, when the smaller active space consisting of the 3a 1 , 1b 2 , 4a 1 , and 2b 2 orbitals is replaced by the active space consisting of the 1b 1 , 3a 1 , 1b 2 , 4a 1 , 2b 1 , and 2b 2 orbitals. For each of the three geometries of water considered in this work, the MMCC(2,3) results improve by ∼ 0.3 millihartree. The MMCC(2,4) results improve by 0.6-0.7 millihartree, when the active space consisting of the 3a 1 , 1b 2 , 4a 1 , and 2b 2 orbitals is replaced by the larger active space consisting of the 1b 1 , 3a 1 , 1b 2 , 4a 1 , 2b 1 , and 2b 2 orbitals. As in the case of the HF molecule, all MMCC energies are above the corresponding full CI energies. Thus, the MMCC approaches are capable of eliminating the non-variational behavior of the conventional CCSD(T) and CCSD(TQ f ) methods at larger internuclear separations.

Summary and Concluding Remarks
We tested the CI-based variants of the MMCC theory, in which "trial" wave functions |Ψ that enter the MMCC energy expressions are obtained with the relatively simple CISDt and CISDtq approaches. We demonstrated that the MMCC(2,3) and MMCC (2,4) (2,3) and MMCC(2,4) approximations, referred to as the CR-CCSD(T) and CR-CCSD(TQ) approaches [7,[16][17][18]22]. This implies that the MMCC theory is a robust formalism, in which very good results can be obtained independent of the quality of the wave function |Ψ used to construct the noniterative energy corrections δ (A) . Another advantage of the MMCC methods is the fact that the MMCC results systematically improve, when we switch from the MMCC(2,3) approach to the MMCC(2,4) approximation and when we increase active space used in the CISDt and CISDtq calculations. As in the case of the CR-CCSD(T) and CR-CCSD(TQ) methods, the energies obtained with the CI-based MMCC(2,3) and MMCC(2,4) approaches are above the corresponding full CI energies. Thus, the CI-based MMCC(2,3) and MMCC(2,4) methods eliminate the pervasive failing of the standard CCSD(T) and CCSD(TQ f ) approximations, which suffer from the highly non-variational behavior in the region of larger internuclear separations.
Clearly, the preliminary findings reported in this work need to be investigated further by studying other molecular examples and by performing the CI-based MMCC(2,3) and MMCC(2,4) calculations with larger basis sets. We also plan to incorporate higher-than-quadruply excited moments of the CCSD equations, corresponding to the projections of the CCSD equations on pentuply and hextuply excited configurations [the M 5 (2)|Φ and M 6 (2)|Φ quantities in Eq. (15), where m A is fixed at 2], in our computer codes. The resulting MMCC (2,5) and MMCC(2,6) methods are expected to improve the very good results reported here even further. More importantly, the MMCC(2,6) method, in which noniterative corrections due to triples, quadruples, pentuples, and hextuples will be added to the CCSD energies, should accurately describe the most difficult cases of bond breaking, including the triple bond breaking in N 2 [87].
Another important area that needs to be explored is the extension of the CI-based MMCC methods to excited electronic states. Our preliminary calculations indicate that the excited-state analogs of the CI-corrected MMCC approaches discussed in this work provide highly accurate results for excited-state PESs [22,71,72].