How ‘de facto variational’ are fully iterative, approximate iterative, and quasiperturbative coupled cluster methods near equilibrium geometries?

Introduction and statement of the problem

Coupled cluster (CC) theory (for a comprehensive review, see Ref. 1]) has become the tool of choice for accurate wavefunction theory (WFT) electronic structure calculations.

An untruncated coupled cluster ansatz, Ψ = exp ( T ̂ ) ψ 0 , where the cluster operator T ̂ = T ̂ 1 + T ̂ 2 + ⋯ + T ̂ n , is merely a clumsy way of carrying out an FCI (full configuration interaction) calculation, which corresponds to the exact solution within a given finite basis set. However, limited CC, with a truncated cluster operator such as Ψ CCSD = exp ( T ̂ 1 + T ̂ 2 ) ψ 0 or Ψ CCSDT = exp ( T ̂ 1 + T ̂ 2 + T ̂ 3 ) ψ 0 not only converges much more rapidly to the FCI limit, but unlike limited CI is rigorously size extensive. The CCSD(T) method, ^{2 , 3} in particular, is often referred to as “the gold standard of quantum chemistry” (following T. H. Dunning, Jr.), and is widely used as a reference or ‘sanity check’ for low-cost methods like DFT (density functional theory).

On the flip side, limited CC methods are nonvariational and hence (again, unlike limited CI) do not yield guaranteed upper bounds to the system’s total energy. Given that variational energies are upper bounds for the exact energy, and that recently ^{4 , 5} there has been revived interest in establishing lower limits for the exact energy, this offers the tantalizing prospect of a rigorous error bound on approximate energies.

There have been efforts in variational coupled cluster theory (e.g., Refs. [6], [7], [8], [9) as well as in developing so-called ‘quasi-variational coupled cluster’ theory, ¹⁰ but very few researchers have adopted such methods. Prof. Eli Pollak (Weizmann Institute of Science), at the Ph.D. defense of one of us (ES), wondered aloud why. This made the senior author ponder just how serious a practical issue (as distinct from a formal one) the nonvariational character of limited CC truly is. We shall address this question presently for molecular systems near their equilibrium geometries.^[1]

Computational methods

Much of the calculations in this paper were carried out using the MRCC ¹³ general coupled cluster program of Kállay and coworkers. The remainder (particularly the approximate CC methods) were performed using a development version of the CFOUR program system. ¹⁴ The lion’s share of the raw data are already available online in the Supporting Information of Ref. 15]; the CCSDTQ567 energies for a subset of systems can be obtained upon request from the senior author.

The molecules considered are the 140 species from the W4-11 thermochemical benchmark ¹⁶ and its 96-member subset W4-08. ¹⁷ These span a range of inorganic and organic molecules, first-row and second-row, and range from essentially purely dynamical correlation (such as H₂O and SiF₄) to strong nondynamical correlation (such as O₃, S₄, C₂, and BN).

Basis sets considered are the Dunning correlation consistent ^{18 , 19} basis sets. Specifically, we considered cc-pVDZ (correlation consistent polarized double zeta) with the polarization function on hydrogen removed – which we denoted cc-pVDZ(d,s) – and cc-pVDZ with all polarization/angular correlation functions removed, which we denoted cc-pVDZ(p,s). To rule out that some of the conclusions are basis set-specific, additional calculations on a subset of systems were carried out using the cc-pVTZ basis set with the d functions on hydrogen removed, which we denoted cc-pVTZ(f,p).

It was shown a quarter century ago by several groups (e.g., Refs. 12], [20], [21], [22], [23) that coupled cluster energies converge rapidly with substitution rank.

In addition, as one of us showed almost two decades ago, ^{24 , 25} high-rank coupled cluster increments, such as E[T ₅] = E[CCSDTQ5] − E[CCSDTQ], converge very rapidly with the basis set, in fact the more rapidly so as the excitation level increases (and increasingly, nondynamical rather than dynamical correlation is sampled). This has been exploited in high-accuracy computational thermochemistry protocols such as W4 theory ^{24 , 26} and HEAT. ^{27 , 28 , 29 , 30} Hence, for the reference correlation energies, we approximated CCSDTQ56/cc-pVDZ(d,s) reference energies as follows:

1. where CCSDTQ5/cc-pVDZ(d,s) is available, E[CCSDTQ5/cc-pVDZ(d,s)]+ E[CCSDTQ56/cc-pVDZ(p,s)]−E[CCSDTQ5/cc-pVDZ(p,s)]

2. otherwise, where CCSDTQ(5)_Λ/cc-pVDZ(d,s) is available, E[CCSDTQ(5)_Λ/cc-pVDZ(d,s)]+E[CCSDTQ5(6)_Λ/cc-pVDZ(p,s)]−E[CCSDTQ(5)_Λ/cc-pVDZ(p,s)]

3. for the largest species, E[CCSDTQ/cc-pVDZ(d,s)]+E[CCSDTQ(5)_Λ/cc-pVDZ(p,s)]-E[CCSDTQ/cc-pVDZ(p,s)]

The approximate CC methods considered include CCSD[T] (a.k.a., CCSD + T(CCSD)); ³¹ CCSD(T); ^{2 , 3} CCSDT-1a and CCSDT-1b; ³² CCSDT-2 ³³ and CCSDT-3; ³³ full CCSDT; ^{31 , 34} CCSDTQ-1 ³⁵ and CCSDTQ-3; ³⁶ full CCSDTQ; ³⁷ CCSDT(Q); ³⁸ CCSD(T)_Λ, ^{39 , 40 , 41 , 42} CCSDT(Q)_Λ ^{36 , 43} ,^[2] and CCSDTQ(5)_Λ; ^{36 , 43} and finally, higher-order fully iterative CC methods ^{21 , 44} such as CCSDTQ5 and CCSDTQ56. It needs to be kept in mind that asymptotic CPU time scaling of fully iterative m-fold CC theory will be O ^m V ^m+2 N _iter (where O and V represent the numbers of occupied and virtual orbitals, respectively, and N _iter is the number of iterations), compared to O ^m−1 V ^m+1 N _iter) followed by a single O ^m V ^m+1 step for a quasiperturbative ‘parentheses’ method, with a Λ approach roughly doubling the first step.^[3]

Some additional approximate coupled cluster methods were newly implemented by one of us (GHJ) in the development version of CFOUR, with total energies verified against published values. Those are:

1. the ‘addition by subtraction’ 3CC method of Bartlett and Musiał, ⁴⁵ recently advocated by Valeev and coworkers ⁴⁶ as an O ³ V ⁵ N _iter scaling approximation to CCSDT(Q) or CCSDTQ;

2. the CCSD(cT) method of Grüneis and coworkers, ⁴⁷ which is in effect a quasiperturbative approximation to CCSDT-3. (See also Ref. 48] for a recent application to noncovalent interactions.)

We also considered the CC3 approach, ⁴⁹ which is normally used more for excited states than for ground states, and its extensions ³⁶ to quadruple and quintuple excitations, CC4 and CC5, respectively.

For open-shell species, UHF references were used throughout. (A few species, namely the two HOOO isomers, FOO, and ClOO, could not reliably be converged to the same UHF solution in both codes, and were hence eliminated from comparisons.) Reference geometries were taken from the supporting information for Ref. 50] and used ‘as is’, without further optimization.

Results and discussion

We first attempted to ensure, for a subset of molecules, that the fully iterative CCSDTQ56 was close enough to full CI. In order to do so, we carried out fully iterative CCSDT567 calculations and assessed E[CCSDTQ567] − E[CCSDTQ56]. The largest correlation energy increments from connected septuple (!) excitations were seen for such troublesome species as C₂ and BN, and even there did not exceed several microhartree; the largest contribution found was 6 μE _h for ozone. For the remaining species, connected septuples contributions were on the order of 1 μE _h or less, meaning that CCSDTQ56 is full CI quality in all but name.

A box-and-whiskers plot of errors for various approximate coupled cluster methods is given in Fig. 1. First of all, the fully iterative approaches appear to be de facto variational, albeit not de jure. Also, as can reasonably be expected, CCSDTQ5 is extremely close to the FCI limit: the largest connected sextuples contribution is 85 μE _h for ozone (0.053 kcal/mol). The largest connected quintuples contribution is rather more substantial, skirting the edges of a millihartree for some systems: 0.77 mE _h (0.49 kcal/mol) for S₄, followed by 0.67 mE _h (0.42 kcal/mol) for ozone. That connected quadruples are quite important (e.g., on the order of 6 mE _h , or 3.8 kcal/mol for molecules like S₄, ozone, and FOOF) is well-known by now (e.g., Ref. 24]); the importance of connected triple excitations even more so. (In fact, the percentage of connected triples in the total atomization energy has been proposed almost two decades ago as a diagnostic for nondynamical correlation. ²⁴ )

Fig. 1:

Box plot of errors in W4-11 correlation energies (hartree) relative to the de facto FCI reference data from Ref. 15]. The cc-pVDZ(d,s) basis set was used throughout. The ‘yellow sea’ indicates energies that violate the variational criterion. The box encompasses the middle half of the distribution, i.e., the IQR (interquartile range) between the 25th and 75th percentile (Q1 and Q3, respectively). Whiskers span from Q1-1.5IQR to Q3+1.5IQR. Filled circles are outliers, open circles are ‘extreme outliers’, <Q1-3IQR or > Q3+3IQR.

Let us now turn to approximate triples approaches. The oldest ‘parenthetical’ method, CCSD[T] a.k.a. CCSD + T(CCSD), predictably overshoots the correlation energy in a number of cases, most blatantly for C₂ and BN. This is mostly remedied in the familiar ‘gold standard’ CCSD(T) method, which includes a usually repulsive fifth-order E S T [ 5 ] -like term,^[4] where in our sample only singlet BN diatomic has a CCSD(T) correlation energy that dips below the FCI limit. It is well established ^{23 , 24 , 25 , 27 , 56} that the generally good performance of CCSD(T) for thermochemistry is a consequence of felicitous error compensation between neglect of higher-order triples (typically repulsive) and of connected quadruples (universally attractive). This is perhaps even clearer when the CCSD(T) and CCSDTQ methods are compared against CCSDT rather than FCI (Fig. 2). Also, as seen in Fig. 1, said error compensation can break down severely for molecules with strong static correlation, like O₃ and S₄.

Fig. 2:

Box plot of errors in W4-08 approximate vs. exact CCSDT correlation energies (hartree). The cc-pVTZ(f,p) basis set was used throughout.

In contrast, the ‘lambda quasiperturbative coupled cluster’ approach ^{39 , 40} CCSD(T)_Λ stays above FCI for all species, but is substantially further away from FCI than CCSD(T) for molecules with strong static correlation. Only in an isolated case, the pathologically multireference BN, does the connected triples contribution from CCSD(T)_Λ exceed the CCSDT-CCSD difference.

The approximate iterative triples method CCSDT-1b has a median error of about zero for the triples contribution (Fig. 2), but a much wider spread than CCSD(T), and severely overestimates the triples in many cases. Because of an error compensation with neglect of quadruples, its ‘box’ compared to FCI is actually narrower than that of CCSD(T), leaving aside the dismal performance for BN. CCSDT-3, on the other hand, has broader boxes than both CCSDT and CCSD(T), but nowhere dips below FCI: this can be rationalized ⁵⁵ because, while CCSDT-1b only has T ̂ 2 influencing the triples amplitudes, CCSDT-3 has exp ( T ̂ 1 + T ̂ 2 ) , lacking only direct influence of T ̂ 3 on itself (dropping this term reduces CPU time scaling by one power of V. The difference between CCSDT-1b and CCSDT-3 starts at fifth order with E T Q ( 5 ) , which terms were found to be repulsive.

Despite CCSDT-3 being a clear improvement over CCSDT-1b, it still overshoots the triples contribution for some systems such as BN, SiO, and AlF₃, while on the other hand seriously underestimating it for B₂, C₂, P₄, and FOOF. The difference with CCSDT starts out ⁵⁵ at fifth order in perturbation theory with the triples-triples interaction term E T T ( 5 ) . (Incidentally, we compared CCSDT-3 and CCSDT for selected anions from the G2-1 electron affinities dataset ^{57 , 58} and found that CCSDT-3 overestimates the triples for OH⁻ and F⁻. It is possible that this is linked to the behavior for strongly ionic oxides and fluorides.)

Grüneis and coworkers ⁴⁷ proposed CCSD(cT) – in effect, a noniterative approximation to CCSDT-3 – as a remedy for what they termed ‘the infrared catastrophe’, i.e., the tendency of CCSD(T) to overestimate the effect of triples in systems with small HOMO-LUMO gaps. Indeed, CCSD(cT) stays above FCI throughout (Fig. 1); while the box is broader than for CCSD(T), there are basically no outliers. (In other words, the distribution is less leptokurtic.)

The CC3 method, on the other hand, clearly overestimates triples on average. Owing to an error compensation with the missing quadruples, its error distribution relative to FCI (Fig. 1) is much narrower than that of CCSDT, and in fact more akin to the ‘addition by subtraction’ 3CC approach.

We shall now turn to approximate quadruple excitations methods. As expected, the CCSDT(Q) method ³⁸ is much closer to the FCI limit than CCSDT, but it can be seen here to exceed the FCI correlation energy by up to 1.5 millihartree for a number of molecules – not just the usual suspects like BN, C₂, O₃, and S₄, but also N₂O and a number of others. In fact, even for the quintuples-including CCSDTQ(5), there is still a degree of nonvariational character. CCSDT(Q)_Λ, on the other hand, represents an unqualified improvement over CCSDT(Q) – in fact, its error distribution looks more akin to CCSDTQ(5) than to CCSDT(Q). The largest ‘nonvariational outlier’ of CCSDT(Q)_Λ is C₂, at 0.23 millihartree – nearly an order of magnitude less than for CCSDT(Q). The story seen above for CCSDT-1b vs. CCSDT-3 repeats itself for their quadruples siblings: while CCSDTQ-1b is thermochemically much more accurate than CCSDTQ-3, it does overcorrelate for some systems, while CCSDTQ-3 nowhere violates the ‘variational’ criterion.

Consistent with the observations of Kállay and Gauss ³⁶ for a much smaller (and first-row only) sample, namely the HEAT dataset, ²⁷ CCSDT(Q)_Λ, CCSDTQ-1b, and CC4 all perform about equally well.^[5] CC4 ³⁶ slightly overcorrelates many systems though (most pronouncedly BN and C₂), and even CC5 ³⁶ and CCSDTQ(5)_Λ are not entirely immune to this problem, despite otherwise stellar performance. In contrast, fully iterative CCSDTQ5 does behave ‘variationally’. CCSDTQ(5)_Λ having about one-quarter of the distribution below FCI is furthermore less serious in practice than it sounds – considering that the worst case, singlet BN, is overcorrelated by just 61 microhartree (0.038 kcal/mol).

Looking at the quadruples contributions in isolation (Fig. 3), (Q)_Λ evidently has a much narrower distribution than (Q), but does seem to overestimate the quadruples on average. However, it is hard not to notice the similarity between the box-and-whiskers plots of CCSDT(Q)_Λ − CCSDTQ and of CCSDTQ5 − CCSDTQ: note the narrowness of the CCSDT(Q)_Λ − CCSDTQ5 box. It has been argued previously by Stanton and coworkers ⁵⁹ that CCSDT(Q)_Λ is a superior method to CCSDTQ at less cost, and the present results (as well as Ref. 15]) definitely bear this out.

Fig. 3:

Box plot of errors in W4-11 approximate vs. exact CCSDTQ correlation energies (hartree). The cc-pVDZ(d,s) basis set was used throughout.

At the far right of Fig. 1 are the CCSDT(Q)_/A and CCSDT(Q)_/B approximations ⁴³ to CCSDT(Q)_Λ. CCSDT(Q)_/A is clearly a poor substitute, while CCSDT(Q)_/B appears to be preferable to CCSDT(Q) but remains inferior to CCSDT(Q)_Λ.

Besides CCSDT(Q)_Λ, an even more pronouncedly ‘less is more’ method is of course 3CC, ⁴⁵ which omits certain diagrams that are not required for the method to be size-extensive and exact for 3 particles. The omitted diagrams would cancel with quadruples diagrams, and hence 3CC has been argued ^{45 , 46} to be closer to CCSDTQ than is CCSDT (which has the same cost scaling). This is indeed evident from Fig. 1; visually, one could estimate that 3CC bridges about half the gap between CCSD(T) and CCSDT(Q). Intriguingly, the biggest errors are found not just for severely multireference species like S₄ and O₃, but for highly polar molecules such as SiF₄, CF₄, and SO₃.^[6]

What about CC5 vs. CCSDTQ(5)_Λ? The RMS difference between the energies afforded by the two methods is just 11 microhartree, which drops to 7 microhartree if the vexing case of BN is dropped. Hence they are nearly interchangeable, but of course CCSDTQ(5)_Λ is the more economical of the two.

Finally, CCSDT5(6)_Λ (not displayed in Fig. 1) does an exceedingly good job of capturing the connected sextuples: the largest discrepancy with fully iterative CCSDTQ56 is just 5 microhartree for the pathologically multireference BN diatomic.

A reviewer raised the caveat that our conclusions are mostly based on cc-pVDZ calculations, and that the importance of static correlation may be exaggerated in small basis sets owing to coupling between 1-particle and n-particle spaces. However, as comparison of Figs. 1 and 2 reveals, the salient features are the same with the cc-pVTZ basis set. For a detailed analysis of the basis set convergence of high-order coupled cluster contributions, see Karton ⁶⁰ as well as Ref. 15].

Conclusions

While limited coupled cluster theory is formally nonvariational, it is not broadly appreciated whether this is a major issue in practice. Through comparison with de facto full CI energies for a relatively large and diverse set of molecules near equilibrium geometries in a polarized double-zeta basis set, we were able to establish that:

1. Fully iterative limited CC methods such as CCSDT, CCSDTQ, CCSDTQ5 do represent upper bounds to the FCI energy in practice.

2. The approximate triples methods CCSD(cT) and CCSDT-3 appear free from ‘nonvariational’ errors, but especially the latter has a broader error distribution than not just CCSDT but even CCSD(T). The iterative approximate quadruples method CCSDTQ-3 likewise behaves like an upper bound for FCI.

3. In contrast, quasiperturbative approaches such as CCSD(T), and especially CCSDT(Q), may significantly over-correlate molecules if there is significant static correlation.

4. Lambda coupled cluster methods such as CCSDT(Q)_Λ and CCSDTQ(5)_Λ strongly mitigate the issue.

5. The ‘addition by subtraction’ 3CC is about midway in quality between CCSDT and CCSDT(Q)_Λ, but does overcorrelate some systems.

6. Fully iterative CCSDTQ56 as well as CCSDTQ5(6)_Λ are for all intents and purposes of full CI quality, while CCSDTQ(5)_Λ comes very close.