Quasi-variational coupled-cluster theory : Performance of perturbative treatments of connected triple excitations

Quasi-variational coupled-cluster methods are applied to a selection of diatomic molecules. The potential energy curves, spectroscopic constants, and size consistency errors are calculated and compared to those obtained from both single- and multi-reference methods. The effects of connected triple excitations are introduced with either the standard perturbative (T) formulation, or in the renormalised form, and its symmetrised approximation. It is found that the renormalised ansatz is significantly superior to the standard formulation when describing bond breaking and that in most circumstances, the computationally simpler symmetrisation gives nearly identical results.


I. INTRODUCTION
Coupled-cluster (CC) theory [1][2][3] is well established as the ab initio method of choice for the prediction of molecular properties with controlled accuracy when the Hartree-Fock reference state is known to be a reasonable approximation.Geometries, vibrational frequencies, and reaction energy changes can all be computed with errors that are small enough for meaningful predictive chemistry with a cluster operator containing up to double (CCSD), triple (CCSDT), or quadruple (CCSDTQ) excitations.In many circumstances, perturbatively treating the effects of connected triple excitations, CCSD(T), 4 is sufficient.However, it is well understood that standard truncated coupled-cluster methods, like CCSD and CCSD(T), result in significant errors when used to describe chemical systems with strong non-dynamic correlation. 5,6For molecules, this occurs typically when covalent bonds are extended and is the most pronounced for the breaking of double and triple bonds.Normally, such situations have to be investigated with multireference (MR) methods, which means selecting a suitable active space into which the approximate wavefunction can be expanded.The use of a multiconfigurational self-consistent field (MCSCF), with subsequent treatment of dynamical correlation effects through multireference configuration interaction (MRCI), 7 has been very successful in generating global ground and excited state potential energy surfaces for small molecules.Nevertheless, for larger molecules, the difficulty of choosing a suitable active space, the computational scaling with system size, and the lack of a general size-extensive formulation all provide an impetus for seeking alternative approaches.][10][11][12] One such family of methods is quasi-variational coupled-cluster doubles (QVCCD), hereafter collectively denoted as the QV methods, which have been shown to produce accurate potential energy surfaces along dissociation a) Electronic mail: KnowlesPJ@Cardiff.ac.uk paths despite being based on a single-determinant reference function. 11,13,14he first computational implementation of the QVCCD method contained significant inefficiencies, meaning that previous studies were limited in the scope of the problem that could be addressed.We have completely reimplemented the method with near-optimum efficiency within the integrated tensor framework (ITF) 15,16 contained in the Molpro 16 software package, to allow the investigation of larger systems. 17he ITF is a platform that allows developers to write code in a simplified domain-specific language based around tensor manipulations and contractions.This code is then compiled into optimised instructions for a general virtual machine that provides an efficient parallel evaluation of the tensor contractions.
In the first part of this work, the potential energy curves (PECs) and vibrational constants of several diatomic molecules have been investigated and compared with singlereference coupled-cluster methods and multireference configuration interaction (MRCI).Size consistency errors have also been determined by comparing the asymptotic values that the QV methods converge to and the sum of the calculated energies for the separated atoms.
The likely invalidity of the standard (T) triples correction to CCSD and QVCCD in cases where the reference wavefunction is a poor approximation led in our earlier work to the introduction of a renormalised triples correction, 18 based on the formulation presented by Kowalski and Piecuch 19 and by Nooijen and LeRoy. 20The asymmetric formulation of this theory leads to an increased computational cost, and we therefore also proposed a symmetrised approximation, denoted OQVCCDR(T).In this work, we explore, through examples, the effects of the renormalisation and of the symmetric approximation.

II. THEORY
3][24][25] On dissociation, standard projection coupled-cluster theory can give energies that are unphysically below the exact energy, sometimes even diverging without bound, whereas this cannot happen with VCC, for which the energy is a strict upper bound to the exact Schrödinger eigenvalue.However, VCC, even with only single and double excitation operators (VCCSD), 6 is an impractical theory because of its very high computational cost. 5QVCCD introduces an energy functional that is a closed form based on matrix powers, whose expansion in powers of the cluster operator T agrees with the corresponding expansion of VCCD to third order in T , and which is exact for an assembly of non-interacting electron pairs. 11his is achieved by modifying the linearized coupled electron pair approximation [CEPA(0)] function, 26 which itself is a second-order approximation to VCCD, to use modified cluster amplitudes defined to have the desired properties.These modified amplitudes can be generated via a transformation of the true cluster amplitudes, which remain as the parameters to be varied, where the Einstein summation convention has been used, τ ij permutes labels i, j, and the U tensors are responsible for generating the four unique O(T 3 ) VCCD terms.The meaning of the powers of the four-index U tensors is that they are considered as matrices where the upper and lower index pairs each form a composite single matrix row or column index.The resulting set of amplitudes are denoted transformed amplitudes and pre-labeled by the power of the U tensor. 11n the two-electron limit, the QVCCD energy simply reduces to configuration interaction doubles (CID).In the N-electron case, the series expansion of the U − q 2 tensors generates series of Feynman diagrams, some of which match those arising from the equivalent expansion of VCCD at each power of T ; in the case of non-interacting electron pairs, additional relationships between the various diagrams emerge, and render QVCCD exact. 18Note that in the implementation of QVCCD, this series is never evaluated term by term and therefore need not be truncated.The implicit inclusion of not just O(T 3 ) contributions but also an infinite subset of higher-order terms may suggest why QVCCD can correctly describe bond dissociations.
A simple QVCC formulation involving both single and double excitation operators that is exact for two electrons has not so far been found.However the effects of single excitations can be included via orbital rotations, using the doubles-only energy expression, and auxiliary conditions from either the Brueckner ansatz [27][28][29] or variational minimisation of the correlated energy functional with respect to the reference determinant. 30The latter, denoted OQVCCD, has been used in calculations reported here.
To go beyond OQVCCD and include the effects of connected triple excitations, the standard (T) correction 4 can simply be added to the OQVCCD energy. 13However, it has been argued that when T is large, one should use a formulation involving the transformed amplitudes, as well as additional quadratic terms, in the perturbation expression, [18][19][20] An additional, Asymmetric-Renormalised [AR(T)], approximation omits the term quadratic in T, Equation 5 requires two different matrix elements to be evaluated compared to just one in the standard (T) correction.As a result, the computational time to calculate E AR(T) is nearly double that of E (T) .A computationally simpler approximation can be developed by replacing the different cluster amplitudes by 1 T2 , which we refer to as the symmetric-renormalised triples correction, OQVCCDR(T), and which has the same computational effort as E (T) .These corrections are also related to the non-iterative corrections of quasivariational method-of-moments coupled cluster (QVMMCC) and quadratic method-of-moments coupled cluster (QMMCC), 31,32 as well as the triples corrections designed for full extended singles and doubles coupled cluster (ECCSD) and quadratic extended singles and doubles coupled cluster (QECCSD) methods. 33,34Each try to approximate the VCC norm which occurs in the denominator of the triples corrections.Our corrections attempt this via the use of the transformed amplitudes which have the effect of diving the E (T) correction by the CID norm in the limit of two electrons, whereas QVMMCC uses a truncated VCC norm which includes second-order triples effects in the excitation manifold.

III. RESULTS AND DISCUSSION
All calculations have been carried out using Molpro and a new ITF-based implementation of the QV methods. 17For each diatomic system, a PEC has been calculated, from which spectroscopic constants can be calculated by polynomial fitting.The number of points in the curve and the degree of the polynomial were chosen such that the constants were stable to at least the digits presented below.All energies and spectroscopic quantities have been calculated using CCSD(T), the three QV methods [OQVCCD(T), OQVCCDR(T), and OQVCC-DAR(T)], and MRCI, with and without relaxed quadratic cluster corrections. 35Only valence-shell electrons were included in the correlation treatment.In the multireference calculations, the valence space was taken as the active space, except in the case of Cl 2 , for which an additional set of π u and π g orbitals were included.
A well-established alternative to using QV for simultaneously representing static and dynamic correlation with a closed-shell reference wavefunction is the completely renormalised coupled-cluster (CR-CC) family of methods. 9It is already well documented that CR-CC(2,3) performs well in describing the breaking of at least single covalent bonds. 36revious work has also indicated the need to extend beyond CR-CC(2,3) when connected quadruple excitations assume importance; in some cases, their effects can be introduced through an additive correction. 37,38Calculations have also been performed with the CR-CC(2,3), variant D, and CR-CC(2,3)+Q, variant B, methods from the GAMESS computational package 39,40 and used as a comparison to the QV results.
The cc-pVQZ and cc-pV5Z basis sets 41 were used, with the standard X 3 extrapolation of the correlation energy to the basis-set limit. 42n order to support comparison with empirical values, 43 additional CCSD(T) calculations have been carried out by correlating the core orbitals and including relativistic scalar effects via the second-order Douglas-Kroll-Hess Hamiltonian 44 using the cc-pwCVQZ-DK and cc-pwCV5Z-DK basis sets 45 extrapolated to estimate the limit.
All the calculations carried out used closed-shell Hartree-Fock to define the reference wavefunction.At dissociation, therefore, there will be a large size-consistency error in the reference energy, meaning that the reference energy is significantly higher than the sum of the Hartree-Fock energies of separate ground-state atoms.This error would disappear in a subsequent full configuration interaction treatment, but in any approximate theory there will be only partial cancellation.The size consistency error was estimated by taking the difference of the QV energy calculated at long bond lengths and the sum of the CCSD(T) energies of the separate atoms; in the case of open-shell atoms, the fully spin-restricted RCCSD(T) 46,47 ansatz was used.Open-shell QV methods are currently not available, but they would be expected to give similar energies in the case of the atoms where there is no strong static correlation; this aspect is quantified below.

A. Singly bonded molecules
A majority of chemical reactions involve breaking and forming single bonds.Therefore it is important for a quantum chemical method to describe this phenomenon correctly.To start, we investigate the breaking of three singly bonded molecules: HCl, BCl, and Cl 2 .

HCl
For the dissociation of HCl (Fig. 1), we observe the typical behavior of CCSD(T) at long bond lengths: the energy forms a maximum, in this case around 3.0 Å, before falling rapidly.OQVCCD(T) also matches this behavior, illustrating the breakdown of the (T) correction.Eigenvalues of the Fock matrix appear in the denominator of this correction and so can lead to an over-estimation of the effects of the triples when these eigenvalues are close in energy.This problem can be corrected by using a renormalised triples scheme.
OQVCCDR(T) performs well, leading to a flat dissociation limit of around 460.201 hartree, which is bounded by both the MRCI and MRCI+Q energies.CR-CC(2,3) has already been shown to dissociate HCl and reproduce the behavior of CCSDT at bond lengths up to 5 Å. 48This behavior is reproduced by the two CR-CC methods; both predict a curve that is almost identical to MRCI+Q.
The size consistency error for the renormalised QV methods, shown in Table I, is an order of magnitude larger than the MRCI error; however, it is an order of magnitude less than the OQVCCD(T) error.
All the QV methods predict a first vibrational constant in agreement with the CCSD(T) value (Table II); there is a slight overprediction of around 5 cm 1 from the renormalised QV methods.The CR-CC methods underpredict these values compared to CCSD(T); CR-CC(2,3)+Q predicts a constant that is 3 cm 1 lower.

BCl
For BCl (Fig. 2), the error in the CCSD(T) energy comes mainly from the inability to model multi-reference systems as OQVCCD(T) produces an energy that is qualitatively similar to the MR methods.At around 4.0 Å, the CCSD(T) energy sharply drops.The renormalised QV methods do not form this maximum and flatten out to asymptotic limits.The symmetric-and asymmetric-renormalised triples predict a slightly higher energy than the standard triples at long bond lengths; the difference being 0.014 hartree.Surprisingly, the standard triples correction produces a lower size consistency error than the renormalised QV methods by two orders of magnitude.
Both the CR-CC(2,3) fall and CR-CC(2,3)+Q fall appear to follow parallel to the OQVCCDR(T) energy at longer bond lengths.These methods failed to converge past 4.5 Å.The QV and CR-CC methods predict vibrational constants which are close to the CCSD(T) and empirical values; the difference being 1.8 cm 1 between CCSD(T) and OQVCCDAR(T) and 1.3 cm 1 between CCSD(T) and CR-CC(2,3)+Q.

Cl 2
Again, CCSD(T) produces a maximum at long bond lengths for Cl 2 , in part, due to the breakdown of the triples correction (Fig. 3).Both renormalised QV methods follow the MCRI+Q energy, with an overprediction in energy from around 3.0 Å.The energy for both these methods slowly begin to fall at 4.3 Å, the difference between here and 6.5 Å being FIG. 2. Calculated potential energy curves for BCl with extrapolated cc-pVQZ:cc-pV5Z basis set.
1.8 × 10 3 hartree.This can be explained by the non-variational nature of the QV methods; therefore, the energy is not bounded from below by the exact Schrödinger energy.
Both the CR-CC methods are also able to qualitatively dissociate this system and predict an energy curve that is closer to the MRCI+Q energy than OQVCRT(T).However, like the QV results, they also form a maximum, though at the longer bond lengths of 4.8 and 5.3 Å for CR-CC(2,3) and CR-CC(2,3)+Q, respectively.
The renormalised QV methods and MRCI converge to a larger asymptotic limit compared to OQVCCD(T) and MRCI+Q.The size consistency error for all the QV methods is an order of magnitude greater than the MR methods.Again, all the QV and CR-CC methods predict vibrational constants that are close to the CCSD(T) values.Both the renormalised methods slightly overpredict this value by 6 cm 1 .

AlO −
Figure 4 presents the PECs for AlO , formally a singly bonded anion.The CCSD(T) method fails when the MR nature of the system becomes too large at around 3.4 Å.The QV methods are able to dissociate the molecule and lead to an asymptotic limit.The renormalised triples methods are bounded by the MRCI and MRCI+Q energies, but lead to an asymptotic limit at around 316.985 hartree, which is lower than the MRCI+Q energy by 0.02 hartree.All the QV and MR methods predict a lower asymptotic limit, producing size consistency errors that are comparable in magnitude; however, the cluster correction on MRCI appears to increase this error by an order of magnitude.
All QV methods predict spectroscopic constants in agreement with CCSD(T) though the differences are larger than the previous molecules.OQVCCDR(T) deviates from the CCSD(T) value by 8.8 cm 1 , whereas this difference grows to 12.3 cm 1 when OQVCCDAR(T) is used.CR-CC(2,3)+Q predicts a vibrational frequency that is closer to the CCSD(T) value; this differs by 2.9 cm 1 .

S 2
For the Σ + g excited state of S 2 , the QV methods are able to dissociate the molecule and produce a qualitatively correct PEC (Fig. 5).Renormalised QV methods predict a slightly higher asymptotic value than MCRI, while the standard triples method flattens out to a limit above the MRCI+Q energy.The multireference energies are quicker to reach an asymptotic limit compared to the renormalised QV methods, which is a characteristic of the residual ionic character in the wavefunction at long bond lengths.Interestingly the standard triples QV method does not exhibit this ionic character and follows the MRCI+Q energy well; nevertheless, the energy gradually begins to fall at longer bond lengths.Both CR-CC methods failed to converge for any points along this surface.
The QV and MR methods overpredict the dissociation energy, with OQVCCDAR(T) producing the largest size consistency error of 0.057 hartree.The cluster correction to MRCI reduces this error by an order of magnitude compared to MRCI without the correction.In the case of the sulfur atom, QV calculations can be carried out for the single-determinant mixture of 1 D and 1 S terms, and we find that the OQVCCDR(T) energy lies only 0.003 hartree above RCCSD(T).Thus the majority of the extensivity error cannot be explained by the adoption of RCCSD(T) for the atoms.
All QV methods predict larger vibrational constants compared to the CCSD(T) value; both the renormalised methods deviate from this by around 15 cm 1 .

P 2
To break a triple bond, excitations up to hextuples should be included in the wavefunction expansion.CCSD(T) does not disprove this as the energy shows the unphysical maximum in Fig. 6.The QV results show a marked improvement even for these challenging systems.This is most likely due to the implicit inclusion of higher excitations via the U − q 2 tensors.
The CR-CC methods were not designed to model the breaking of more than double bonds.This is evident in the curves produced for P 2 .The CR-CC(2,3) potential curve is almost identical to CCSD(T), exhibiting unphysical divergence at long bond lengths, and does not display the flattening to an asymptotic energy.The inclusion of the additive quadruples correction does not improve the method for this case.
The QV energy for P 2 , an analog of N 2 , shows it is able to predict a qualitatively correct PEC, even in large nondynamic correlation regimes.The standard triples correction appears to perform better in this case, as the energy closely follows the MRCI energy, whereas the renormalised methods FIG. 6. Calculated potential energy curves for P 2 with extrapolated cc-pVQZ:cc-pV5Z basis set.exhibit a strong ionic character, leading to larger asymptotic limits.
The QV methods all overpredict the first vibrational constant compared to the CCSD(T) value.The smallest deviation occurs for OQVCCD(T) which overpredicts the constant by 10.5 cm 1 .This difference increases to 15.2 cm 1 when OQVCCDAR(T) is applied.The CR-CC methods however produce constants that are closer in value to CCSD(T); CR-CC(2,3)+Q underpredicts this by 3.6 cm 1 .

SiO
Again, at the bond lengths examined, all the QV methods can dissociate SiO without the energy dropping toward negative infinity; however, all the triples schemes predict a higher energy at long bond lengths compared to MRCI (Fig. 7).The energy is no longer bounded by the MRCI and MRCI+Q result.The standard (T) correction produces a very similar answer to the renormalised methods.
The CR-CC methods reproduce the behavior seen in the P 2 system; CR-CC(2,3)+Q clearly forms a maximum at FIG. 7. Calculated potential energy curves for SiO with extrapolated cc-pVQZ:cc-pV5Z basis set.
2.6 Å, while the CR-CC(2,3) curve is still rising in energy before failing to converge at longer bond lengths.
The QV methods produce larger vibrational constants by around 10-12 cm 1 compared with the CCSD(T) result.Surprisingly, given the shape of the PEC, both CR-CC methods drastically underpredict the vibrational constant by around 180-12 cm 1 .This value is stable with respect to the addition and removal of points along the curve and the degree of the polynomial fit.

C. Asymmetric triples correction
Previous work has not tested the assumptions made in arriving at the renormalised triples corrections, i.e., whether Eq. ( 6) is a good approximation to Eq. ( 5).Table III shows the equilibrium bond lengths and the first and second vibrational constants calculated with both the symmetric and asymmetric renormalised triples corrections.These results show that there is little difference between the methods.The largest differences occur for the first vibrational constant, with OQVCCDAR(T) generally producing larger values.The largest difference occurs for AlO of around 4 cm 1 .
For sake of clarity, the OQVCCDAR(T) PECs have not been included in the graphs as they are indistinguishable from those arising from OQVCCDR(T) on the scale of the graphs.However, in general, the energies are not exactly the same, and OQVCCDAR(T) produces slightly higher energies than OQVCCDR(T).
For singly bonded molecules like Cl 2 (Fig. 8), the differences in energy are around 1 × 10 5 hartree and so are virtually insignificant at all bond lengths.The difference rises to a maximum at around 3.4 Å before falling to an apparent limiting value at longer bond lengths.This behavior is not observed for the triply bonded P 2 molecule (Fig. 9), where the energy difference rises to the order of 1 × 10 3 hartree; the largest difference occurs at 3.8 Å which converts to an energy of 6.8 kJ/mol.After 3.8 Å, it becomes difficult to converge the OQVCCD energy and so it is not possible to establish unambiguously the dissociation limiting value.

IV. CONCLUSIONS
This investigation has provided further evidence of the performance of QV methods when dissociating and breaking apart singly and multiply bonded molecules.The unphysical maximum that is typically encountered when calculating PECs with CCSD(T) is avoided in most cases with the QV methods, apart from OQVCCD(T), which is susceptible to the breakdown of the non-iterative (T) correction.The QV methods also compare well with the CR-CC methods investigated in this paper, which produce similar results for the singly bonded molecules, while avoiding the CCSD(T)-like maximum for diatomics like P 2 .
The QV methods can provide reasonable estimates to the vibrational wavenumber and first anharmonicity constants though, as has previously been recognised, they do not perform as well as CCSD(T).CCSD(T) generally behaves well in the immediate region of PEC minima and is therefore well suited to calculating spectroscopic constants.
Estimates of the QV size consistency error have shown that there are large energy differences between the apparent asymptotic limit and the energy of the separated atoms.From the molecules examined, the largest difference is around 0.057 hartree for OQVCCDAR(T) and S 2 .For the multiply bonded molecules, the renormalised QV methods tend to overpredict the energy at longer bond lengths for molecules like S 2 , P 2 , and SiO, in some cases yielding potential energy curves whose slope goes only very slowly to zero in the limit.This can be interpreted as a characteristic of the residual ionic character in the wavefunction as a consequence of the character of the symmetry-restricted single-determinant reference wavefunction.The near-variational character of OQVCCD means that the energy will for this reason be above the exact energy, even when the proper (renormalised) additive correction for the effects of triple excitations is applied.The unrenormalised (T) correction becomes increasingly negative with bond length and provides a compensation that may sometimes mean that the overall size-consistency error is reduced fortuitously.
The large energy differences for the multiply bonded molecules can possibly be explained by the lack of higherorder connected excitations included in the QV ansatz, which may become more important when dissociating P 2 compared to N 2 , which can be qualitatively dissociated even by OQVCCD(T). 13Other research has shown that a correct description of the connected quadruples may be necessary when describing the breaking of multiply bonded molecules; for example, it was shown that these excitations when added to ECCSD described N 2 more accurately. 33This suggests that the QV ansatz could be improved by including connected quadruple excitations via an additive correction in the same vain as the renormalised quadruples corrected, described by Fan et al. in Ref. 34.
The asymmetric renormalised triples correction has also been investigated.It is apparent that Eq. ( 6) is an excellent approximation to Eq. ( 5) based on the calculation of different diatomic constants and PECs.Differences in energy do become larger for P 2 as it is dissociated and the multireference character of the system becomes more pronounced.In such cases, OQVCCDAR(T) may be preferable, but normally the reduced computational cost of OQVCCDR(T) means that it is the method of choice in its class.
a Douglas-Kroll Hamiltonian and all electrons correlated.