Analytic Gradients of Approximate Coupled Cluster Methods with Quadruple Excitations.

The analytic gradient theory for both iterative and non-iterative coupled-cluster approximations that include connected quadruple excitations is presented. These methods include, in particular, CCSDT(Q), which is an analog of the well-known CCSD(T) method which starts from the full CCSDT method rather than CCSD. The resulting methods are implemented in the CFOUR program suite, and pilot applications are presented for the equilibrium geometries and harmonic vibrational frequencies of the simplest Criegee intermediate, CH$_2$OO, as well as to the isomerization pathway between dimethylcarbene and propene. While all methods are seen to approximate the full CCSDTQ results well for"well-behaved"systems, the more difficult case of the Criegee intermediate shows that CCSDT(Q), as well as certain iterative approximations, display problematic behavior.


Introduction
The inclusion of higher-order correlation effects has been recognized as an important contribution to the calculation of accurate molecular energies for more than a decade. [1][2][3][4][5] Such contributions are critical to the accurate evaluation of molecular and reaction enthalpies, [6][7][8] atomization energies, 9,10 barrier heights, 11 and intermolecular interaction energies, 12 and are included in many standard thermochemical model chemistries such as HEAT, [13][14][15][16] Wn, [17][18][19] and ANL-n. 20 However, the effect of higher-order electron correlation on equilibrium molecular geometries, vibrational frequencies, and other properties is less well-understood. Including geometric effects beyond the CCSD(T) 21 or CCSDT 22 level could be desirable in a number of circumstances. For example, Morgan et al. calculated the equilibrium geometry of formaldehyde including CCSDT(Q) 23 and CCSDTQ [24][25][26] higher-order correlation contributions, as well as anharmonic vibrational frequencies including CCSDT(Q) corrections. 27 Puzzarini et al. investigated the effect of higher-order correlation on computed rotational constants. 28 Heckert et al. computed the higher-order correlation contribution to the geometry of several small molecules, and found rather large changes in geometry especially for triply bonded species such as N 2 , HCN, and HCCH, as well as for F 2 . 29, 30 Ruden et al.
investigated the contribution of quadruple excitations to the harmonic frequencies of several diatomics and found corrections on going from CCSDT to CCSDTQ as large as 20 cm −1 . 31 The high computational expense of CCSDTQ naturally limits applicability to small molecules. Approximate coupled cluster methods would ideally extend the range of applicability, especially for non-iterative approximations such as CCSDT(Q). However, for geometric derivatives and especially for harmonic frequencies, the lack of analytic gradients further increases the cost as finite-difference methods must be used. This additional expense furthermore scales with the molecular size rather prohibitively as the finite difference technique depends on the number of degrees of freedom while analytic gradients do not. 32 In order to efficiently study the effect of such higher-order effects on molecular properties, we present here the derivation of analytic gradients for a number of approximate coupled cluster methods that include connected quadruple excitations: CCSDT(Q) 23 (and its A and B variants 33 ), CCSDTQ-1a, -1b, and -3, 34 and CC4. 34 While analytic gradients for general coupled cluster models, including CCSDT 35 and CCSDTQ, 36 have been available for some time, analytic gradients for general approximate coupled cluster methods 34 have not as yet been derived, except for the special cases of CCSDT-n, 37,38 CC3, 38 and of course CCSD(T). 39,40 2 Theory The theory of analytic CCSDT(Q) gradients is developed by first reviewing the basic theory of coupled cluster and its gradients. 32,[41][42][43] Then, the derivation of the CCSDT(Q) energy is reviewed and considerations for non-Hartree Fock references (such as ROHF and QRHF) are discussed. These theories are then combined to derive explicit CCSDT(Q) gradient expressions. Next, iterative approximations to the CCSDTQ energy are reviewed, and finally the corresponding analytic gradient theories are developed.

Coupled Cluster Gradients
The coupled cluster energy is conveniently written as a matrix element of the coupled cluster transformed Hamiltonian,H = e −TĤ eT = Ĥ eT c , [44][45][46] whereĤ is the Hamiltonian in the molecular orbital basis, normal-ordered with respect to the reference single-particle wavefunction |0 , and the cluster operatorT is an excitation However, when deriving the theory of coupled cluster gradients and properties, it is more convenient to use the stationary coupled cluster energy functional, 41 Since the coupled cluster equations are satisfied, for all excited determinants P | = S| + D| + T | + Q| + . . . (i.e. single, double, triple, quadruples excitations etc.) up to the number of excitations included in the model, then the energy obtained with this functional is trivially the same as in (1). For CCSDT and CCSDT(Q), P | = S| + D| + T |. However, both derivatives of the energy with respect to some parameter χ and expectation values of an arbitrary propertyÔ can be succinctly written using this same functional, These expressions can be further generalized by constructing one-and two-particle density matrices, assuming that the operatorÔ may have one-electron and/or two-electron components.
The partial derivatives in (9) and (10) are derivatives of the molecular orbital basis quantities, and include contributions from the response of the orbitals to the perturbation. 41,47 These contributions can be separated out and explicit dependence on the perturbation removed through the use of the coupled perturbed Hartree Fock (CPHF) Z-vector equations (similarly to how theΛ operator removes the dependence on the differentiatedT amplitudes). 32 The relationship of the orbital-response-corrected density to the molecular orbital density is entirely independent of the source of the density, and so is not specific to the particular flavor of coupled cluster theory or even to coupled cluster at all. Similarly, as many coupled cluster methods, including CCSDT(Q), require semi-canonical orbitals (i.e. f a b = ǫ a δ ab and f i j = ǫ i δ ij ), the modification of the density matrix to the use of perturbed canonical orbitals is also independent of the source of the density matrix. 39,40 This means that for methods which require (semi-)canonical orbitals, only the diagonal elements of the virtual-virtual and occupied-occupied one-electron density matrices need to be computed.

The CCSDT(Q) Energy
The derivation of the CCSDT(Q) energy is similar in many ways to the derivation of the CCSD(T) energy. In the context of many-body perturbation theory (MBPT), taking the lowest-order correction to the CCSDT energy, and replacing the approximateT 1 CCSDT +H [1] CCSDT +H [2] CCSDT + . . .
is the complement space. Perturbation orders are assigned according to the usual Møller-Plesset partitioning. Taking the lowest (fourth)-order contribution gives the CCSDT(Q) Λ (also called ΛCCSDT(Q)) energy, CCSDT |0 (17) For canonical Hartree Fock references (or other references where f a i = f i a = 0), the convergedΛ 2 andΛ 3 amplitudes are the same asT 2 andT 3 to lowest order. Thus, as in CCSD(T), theΛ amplitudes may be approximated byT † , giving the CCSDT(Q) method, For non-Hartree Fock references such as ROHF and QRHF, though, Kállay notes thatT 3 andΛ 3 are no longer identical at lowest order. 33 By incorporating the additional lowestorder (disconnected) contributions and replacing approximateT 1 andT 2 amplitudes with their converged counterparts at one of two stages, two additional methods, CCSDT(Q)/A and CCSDT(Q)/B may be derived, One difficulty with these methods, though, is that when applied to a canonical reference, they do not reduce to the "normal" CCSDT(Q) energy. While the additional contributions toΛ 3 in the non-HF case are technically at the same order as the portion in common withT 3 , these terms are entirely disconnected, which tends to produce a numerically less significant contribution. On the other hand, the balance between connected and disconnected terms is sometimes necessary to avoid excessive basis set dependency and other problems, as in the case of CCSDT[Q] vs. CCSDT(Q). The gradients of the normal HF CCSDT(Q) method as well as the two non-HF variants will be examined on an equal footing.

CCSDT(Q) Gradients
The gradient of the CCSDT(Q) energy can be simply obtained by differentiating the energy expression: and used the identityR χ n = −R nF χR n . We can separate each of these expressions into two parts: one which depends on derivatives of the Hamiltonian multiplied by CCSDT(Q) density matrix elements, and one which depends on derivatives of the coupled cluster amplitudes: given in terms of the density matrices, where ∆ pq = δ pq + δ pa δ qi + δ pi δ qa , and intermediates, where the contraction line indicates that at least one creation/annihilation operator inV must remain uncontracted (i.e.V must have at least one external line in the diagrammatic representation). The factor ∆ pq indicates that the density matrix elements only need be computed for the occupied-virtual and virtual-occupied blocks and for the diagonal of the occupied-occupied and virtual-virtual blocks. Note that the CCSDT(Q) one-particle density matrix does not have an occupied-virtual contribution, and is only valid for canonical (HF) references.
At this point the equations for the CCSDT(Q) and CCSDT(Q)/A,B gradients may be combined in a generic expression, whereŜ X =Ŝ X 1 +Ŝ X 2 +Ŝ X 3 and either X = A, B or is "empty", indicating the CCSDT(Q) gradient and intermediates (withŜ 1 = 0). Now, the definition of the derivative coupled cluster amplitudes (obtained by differentiating the coupled cluster amplitude equations) may be inserted to transform dependence on the derivative amplitudes to dependence on the derivative Hamiltonian, Lastly, the solution of theΛ X equations can be combined with the solution of the CCSDT Λ equations by solving for their linear combination, and the combined density matrices constructed which give the gradient of the total CCSDT(Q)/X energy (total correlation energy; the gradient of the reference energy is computed and added in the standard way), These density matrices may then be processed in the usual way to compute energy gradients and properties as desired.

Iterative Approximations
The CCSDT(Q) energy is a non-iterative correction that is applied after the CCSDT equations have converged. Another approach to approximating the CCSDTQ energy may be obtained by deleting some terms from the CCSDTQ equations, but maintaining the iterative structure of the problem. In particular, we wish to delete at least allT 4 →T 4 terms except those coming fromĤ Depending on the additional terms that are deleted (guided by a mixture of perturbation theory and pragmatism), the following approximations may derived, 34 CCSDTQ-1a (≡ CCSDTQ-1): CCSDTQ-1b: whereP X = |X X| is a projection operator onto the given excitation manifold andV ′ = is almost equal which are excluded by (P TVT3PD )T 3 (since the projection operators and connectivity condition require thatV be connected to the innerT 3 vertex by exactly three indices). These terms are specifically deleted as they would require O(n 10 ) computation.

Gradients in the Iterative Approximation
The derivation of the gradient of the energy in the iterative approximation closely follows the derivation of gradients for canonical CC methods. However, since the cluster equations have been modified to delete specific terms, we express the energy and its derivatives using an effective transformed HamiltonianH (also called the Jacobian), where X is one of the approximate methods above. The effective transformed Hamiltonian may be specified in block form for each of the approximate methods, The stationarity conditions of the energy functional ( (66) and (67)) along with the definition ofH determine the equations forΛ, CCSDTQ-1a: CCSDTQ-1b: CCSDTQ-3: where 0|(1+Λ)H|X CCSDT = 0|(1+Λ 1 +Λ 2 +Λ 3 )H CCSDT |X are the CCSDTΛ equations.
As mentioned by Gauss et al., 38 the equations for theΛ amplitudes may also be derived diagrammatically by "capping" the each cluster amplitude diagram with the appropriatê Λ amplitude vertex to form a closed diagram, and then sequentially deleting each cluster amplitude vertex to give the set of openΛ amplitude diagrams.
As with CCSDT(Q), the gradient and molecular properties are computed through the one-and two-particle density matrices. These are easily derived from (67), are the CCSDT density matrices, but with cluster and Λ amplitudes determined by the modified equations for the specific method in question. As with theΛ amplitudes, the density matrices may be derived diagrammatically starting from the diagrammatic form of the energy functional (66), and then deleting the Hamiltonian vertex. Using the diagrammatic approach, care must be taken with the sign and numerical prefactors such that the definition of the density matrices matches that of (9).

Implementation
The total CCSDT(Q), CCSDT(Q)/A, CCSDT(Q)/B, CCSDTQ-1a, CCSDTQ-1b, CCSDTQ-3, and CC4 density matrices have been implemented in the NCC module 50 of the CFOUR program system. 51 As of this time, this implementation only handles closed-shell reference states, including RHF and closed-shell QRHF (or other restricted non-HF orbitals, e.g. Kohn-Sham orbitals). Additionally, the implementation makes use of non-orthogonal spinadaptation techniques 50,52 as well as recent advances in high-performance tensor contraction, 53,54 which allows for high efficiency and relatively compact working equations.
The cost of a CCSDT(Q) gradient calculation is between two and three times that of the corresponding energy calculation. The energy calculation requires the iterative solution of the CCSDT coupled cluster equations with O(n 8 ) cost, and four non-iterative O(n 9 ) steps (multiplication ofT 3 with v ab ci and v ia jk , andT 2 with the three-particle intermediates W abc ije andW abm ijk [see below]). The gradient calculation instead requires two iterative O(n 8 ) procedures (the CCSDT coupled cluster equations and the combinedΛ +Λ equations) and ten non-iterative O(n 9 ) steps. In order to reach this minimal number of O(n 9 ) steps, the contributions to (Γ ′ (Q) ) ab ci and (Γ ′ (Q) ) ci ab must be combined (and similarly for (Γ ′ (Q) ) ia jk and (Γ ′ (Q) ) jk ia ) when computing the symmetrized two-electron density, In the same way, theT where the permutation operator P antisymmetrizes the labels on either side of the slash.
For CCSDT(Q) gradients, the similar intermediate three-particle density matricesΓ aij bcd and Γ ijk abl are also required. These intermediates must either be stored on disk or recalculated as needed. In our implementation we have found that storing them to disk tends to be the most efficient solution as recalculation is often more expensive and can, in extreme circumstances, even increase the formal scaling of the this term. However, as theW abc ije intermediate in particular may be much larger than even theT 3 amplitudes, the disk space requirements of CCSDT(Q) calculations are indeed more strenuous than for CCSDT.
The iterative approximations all scale as O(n 9 ), with the rough order of cost given by CCSDTQ-3 > CC4 > CCSDTQ-1b > CCSDTQ-1a. The main feature that impacts the computational cost is the number of O(n 9 ) steps. For CCSDTQ-1a, the cluster amplitude equations require only four such steps per iteration, while all other methods require six steps due to the directT 4 →T 3 coupling. The difference in theΛ equations is even larger: CCSDTQ-1a again requires four O(n 9 ) steps, but the other methods now require 12 steps per iteration due to the fact the theT 4 amplitudes are required in addition toΛ 4 , and that these amplitudes are generally recomputed rather than stored due to their extremely large size. Additionally, the current implementations of the CCSDTQ-3 and CC4Λ equations are sub-optimal in that they include two additional O(n 9 ) steps per iteration. The construction of the density matrices is not a major bottleneck for the iterative approximations because it occurs only once and not every iteration.
Many terms in the CCSDTQ-3 amplitude equations may be included by defining suitable two-and three-particle intermediates (including the seemingly expensive (P DVT3PS )T 3 term), but in particular the terms t abcd

Dimethylcarbene Isomerization
The most stable conformation of singlet DMC is the C 2 gauche conformation (g-DMC), 60 and the DMC and propene minima are connected by a single chiral transition state (g-TS) along the 1,2-hydrogen shift isomerization pathway. Table 1    The theoretical methods are listed in an abbreviated notation, e.g. Q-1a = CCSDTQ-1a.
(a) g-DMC  indicate that the major geometrical errors in DMC relate to non-bonding intermolecular interactions. As these interactions arise primarily form many-body effects, the fact that connected triple excitations (at least) are necessary to fully describe the structure is not surprising. An alternative hypothesis is that the biradical nature of DMC could account for the greater importance of higher-order correlation. However, a better description of biradical character would be expected to primarily affect the C-C-C angle. An inspection of the detailed data in the Supplemental Information shows that this is not the case, and that  Mean absolute errors for CI are given in Table 2. Right off the bat, an inspection of the errors in the rotational constants, with errors approximately five times larger, show that CI is significantly more sensitive to higher-order correlation effects than the DMC-propene system. Increases in error of a similar magnitude are evident for bond lengths and angles, except that the angle error of DMC still outstrips that of CI at the CCSD(T) level. A major difference between DMC-propene and CI is that a number of higher-order methods fail to improve on CCSD(T) for these geometric quantities, with only CCSDTQ-1b and CC4 (and to some extent CCSDTQ-3) providing a reliable and significant decrease in error. Interestingly, and CCSDTQ-1a further exhibit severe errors (up to 30 cm −1 ) in the out-of-plane motions that are not present at the CCSD(T) level. It is quite difficult to ascribe these errors to a particular feature of the theory, for example using diagrammatic or perturbation arguments.

Criegee Intermediate
Perhaps the simplest, albeit the least satisfying, explanation is a tendency of some of the methods, in particular CCSDT(Q), to "overshoot" the full CCSDTQ results. This can be seen explicitly in the absolute energies, where CCSDT(Q) overshoots CCSDTQ by more than a kcal/mol. On the other hand, CCSDT(Q)/A and /B under-and overestimate the CCSDTQ contribution by roughly the same amount, while CC4 also overshoots CCSDTQ. Thus, the extrapolation of simple energetics to geometric and vibrational parameters is qualitative at best. Perhaps the best example of the disconnect between energetics and geometry is seen for CCSDT(Q), which, as noted above, performs rather poorly for geometries and frequencies, but also halves the error in the absolute energy compared to CCSD(T). CCSDTQ-3 does make some improvement over CCSD(T) for CI, but is again not as accurate as CCSDTQ-1b.

Timings
While the results above show that certain iterative, and in some cases, non-iterative approximations to CCSDTQ can yield substantial accuracy gains compared to CCSD(T) and CCSDT, the computational cost of these methods is a critically important factor in determining when such methods may be used. Full CCSDTQ scales with the tenth power of system size (O(n 10 )), and is applicable only to very small molecular systems, especially when analytic gradients or properties are desired. While the approximate methods described here reduce this scaling to O(n 9 ), they still represent a quite considerable cost increase over CCSDT and especially CCSD(T). Figure 1 shows timings for a single gradient evaluation, as a percentage of the CCSDTQ timing for all methods considered. CCSD(T) only requires 1%. The A and B variants, while increasing the cost again slightly, likewise remain in the 1% range. When these methods can be considered reliable (strongly single-configurational systems with important many-body dynamic correlation), this modest increase in time over CCSDT for systems of the size studied here shows that they can be applied rather routinely.
For CI, however, the non-iterative approximations are not reliable and an iterative approximations is mandatory. Of the iterative methods, CCSDTQ-1a is the cheapest, requiring 3-4% of the CCSDTQ time. As shown above, however, this methods consistently fails to improve upon CCSD(T). CCSDTQ-1b, however, "only" increases the computational time to 4-5% of CCSDTQ, and performs perhaps the best among all of the approximate methods tested. CCSDTQ-3 and CC4 further increase the cost to between 6 and 10% of CCSDTQ.
CCSDTQ-3 is somewhat inconsistent in predicting both geometric and energetic quantities, and due to its increased cost relative to CCSDTQ-1b, the former seems better justified as a standard iterative approximation. CC4 also performs consistently well, but again entails a cost increase over CCSDTQ-1b. Depending on whether an iterative or non-iterative approximation is appropriate for the system in question, a roughly one to two order-of-magnitude reduction in cost relative to full CCSDTQ can be expected, with very nearly the same level of accuracy.

Conclusions
The application of approximate quadruples methods to the equilibrium geometries and harmonic vibrational frequencies of the dimethylcarbene to propene isomerization pathway system shows a uniformly good comparison to the full CCSDTQ results, perhaps with the exception of the CCSDTQ-1a method. Since this method alone lacks a direct coupling from the quadruples amplitudes back to the triples, this seems the most likely explanation for this poor behavior. Even for DMC, which might be expected to exhibit mild multi-reference character due to the biradical carbon center, the performance of iterative and non-iterative approximate quadruples method alike remains very good. CCSD(T) does indeed show an increase in error in this case, which can be traced almost exclusively to the interaction between the two gauche hydrogens which, presumably, has important many-body contributions. For these and, we predict, other "well-behaved" molecular systems, the CCSDT(Q) method, or one of the two variants developed to rigorously treat the case of a non-Hartree-Fock reference state, seem good choices that maximize the accuracy-to-computational cost ratio. In particular, based on the effect of non-bonded interactions in DMC, we expect that the inclusion of quadruples effects in the calculation of equilibrium geometries and (harmonic) vibrational frequencies to be important for dispersion-bound complexes and strained molecules.
The application to the simplest Criegee Intermediate is much more mixed, and clear evidence of multi-reference character is seen in e.g. the magnitude of the converged coupled cluster amplitudes. In this case, all methods exhibit a degradation of error with respect to CCSDTQ, but in particular, the CCSDT(Q) method (and to a lesser extent its A and B variants) as well as CCSDTQ-3 show a distinct worsening in accuracy relative to CCSD(T), CCSDT, and the other approximate quadruples methods. CCSDT(Q) performs especially poorly, and in fact worsens the agreement with full CCSDTQ compared to CCSD(T) in every regard except the total energy. Since the vast majority of previous studies so far have focused on the energy alone in evaluating CCSDT(Q) and other approximate quadruples methods, it seems important to consider derived properties in order to gain a full picture of the intrinsic accuracy of such methods. The results presented here suggest caution when applying CCSDT(Q) to moderately multi-reference systems, although further benchmarking is necessary to fully delineate the realm of reasonable applicability. Luckily, the A and B variant show a distinct improvement over "plain" CCSDT(Q) even in the restricted Hartree-Fock case and serve as convenient substitutes. CCSDT(Q)/A can be implemented at essentially no additional cost compared to CCSDT(Q), although the current implementation in CFOUR is not optimal. Likewise CCSDTQ-1b is seen to perform well in all circumstances tested and could in some cases be justified as a higher-cost alternative.
In sum, these results show that the application of post-CCSDT methods to the problem of equilibrium geometries, harmonic frequencies, and potentially to other properties as well is both feasible and potentially worthwhile when very high accuracy is necessary. Because of the computational cost of the CCSDTQ benchmark calculations, only molecules with three first-row atoms were investigated here. However, molecules with as many as six firstrow atoms should be accessible at the double-ζ level using the analytic gradient theory for approximate quadruples methods.

Acknowledgement
This work was supported by a generous start-up grant from SMU, and all calculations were performed on the ManeFrame II system at the SMU Center for Scientific Computation. I would like to specially thank Dr. John Stanton for his unrelenting encouragement, which was instrumental in bringing this work to completion.

Supporting Information Available
The following files are available free of charge. A supporting information file is available which contains individual errors for each geometrical and vibrational parameter as well as the CCSDTQ benchmark values.