On the distinguishable cluster approximation for triple excitations

TL;DR

This paper extends the distinguishable cluster approximation to triple excitations in coupled cluster methods, significantly enhancing energy accuracy while maintaining computational efficiency and desirable theoretical properties.

Contribution

It introduces a novel application of the distinguishable cluster approximation to triples, improving coupled cluster results with manageable computational scaling.

Findings

Improved absolute and relative energies in coupled cluster calculations.

Method is orbital invariant and size extensive.

Exact for three-electron systems.

Abstract

The distinguishable cluster approximation applied to coupled cluster doubles equations greatly improves absolute and relative energies. We apply the same approximation to the triples equations and demonstrate that it can also improve results of the coupled cluster with singles, doubles and triples. The resulting method has a nominal computational scaling of ${\cal O}(N^7)$ in the real-space representation, is orbital invariant, size extensive, and exact for three electrons.

Tables

Table 1. Table 1: Results for a BeH 2 molecular system with one of the BeH bonds stretched to 100 Å.

	Energy, $E_h$	Error, $\mu E_h$
FCI	-15.662290	0
CCSDT	-15.662155	135
UCCSDTa	-15.662305	-15
DC-CCSDT	-15.662271	19

Table 2. Table 2: Statistical analysis of the absolute deviations of CCSD(T), CCSDT, and DC-CCSDT from CCSDT(Q) results (in kcal/mol).

	MAD	RMSD	MaxD
CCSD(T)	0.178	0.263	0.700
CCSDT	0.160	0.224	0.534
DC-CCSDT	0.114	0.176	0.685
	without reaction 27 from Ref. Huntington et al.,2012
CCSD(T)	0.163	0.239	0.698
CCSDT	0.150	0.208	0.529
DC-CCSDT	0.098	0.135	0.334

Equations

$$R_{abc}^{ijk}$$

$${\hat{\mathcal{P}}}_{qr}^p{X}_{pqr}$$

$$(ld|me)T_{ad}^{il}{T}_{ebc}^{mjk}=$$

$$\int d{\mathbf{r}}_{1}d{\mathbf{r}}_2{\varphi }_l^{\ast }({\mathbf{r}}_1){\varphi }_d({\mathbf{r}}_1)\frac{1}{|{\mathbf{r}}_1-{\mathbf{r}}_2|}{\varphi }_m^{\ast }({\mathbf{r}}_2){\varphi }_e({\mathbf{r}}_2)T_{ad}^{il}{T}_{ebc}^{mjk}$$

$$=\int d{\mathbf{r}}_2{T}_{bc}^{jk}({\mathbf{r}}_2)\int d{\mathbf{r}}_1{T}_a^i({\mathbf{r}}_1)\frac{1}{|{\mathbf{r}}_1-{\mathbf{r}}_2|}$$

$$\mathrm{with}$$

$$T_{bc}^{jk}({\mathbf{r}}_2)={\varphi }_m^{\ast }({\mathbf{r}}_2){\varphi }_e({\mathbf{r}}_2)T_{ebc}^{mjk}$$

$$T_a^i({\mathbf{r}}_1)={\varphi }_l^{\ast }({\mathbf{r}}_1){\varphi }_d({\mathbf{r}}_1)T_{ad}^{il}$$

$$(pq|ai|bj)$$

$$\widetilde{(aq|bj)}\leftarrow -\frac{1}{2}{\delta }_{pi}(pq|ai|bj)$$

$$\widetilde{(ip|bj)}\leftarrow -\frac{1}{2}{\delta }_{qa}(pq|ai|bj)$$

Full text

I Introduction

The coupled cluster hierarchy of methods Čížek (1966) is known to converge very quickly to full configuration interaction (FCI) results with comparably low excitation level for weakly correlated systems. Already triple excitation are usually sufficient to reach the 1 kcal/mol accuracy in relative energies, Raghavachari et al. (1989) and with inclusion of higher excitations it is possible to smoothly converge to the exact answer. Tajti et al. (2004); Bomble et al. (2006); Harding et al. (2008) It does not mean, however, that the coupled cluster hierarchy is the most efficient way to approach the FCI limit. On the contrary, results from the distinguishable cluster singles and doubles (DCSD) Kats and Manby (2013a); Kats (2014); Kats et al. (2015); Kats (2016, 2018) and other methods Meyer (1971); Paldus, Čížek, and Takahashi (1984); Piecuch and Paldus (1991); Piecuch, Tobola, and Paldus (1996); Kowalski and Piecuch (2000); Bartlett and Musiał (2006); Nooijen and Le Roy (2006); Neese, Wennmohs, and Hansen (2009); Huntington and Nooijen (2010); Robinson and Knowles (2011); Huntington et al. (2012); Paldus (2016); Black and Knowles (2018) are usually much better than the coupled cluster singles and doubles (CCSD) results. Additionally, the DCSD equations can be implemented more efficiently than CCSD by employing density fitting or other integral decomposition techniques. Kats and Manby (2013b) Recently, two publications have demonstrated that DCSD can be calculated exactly with the ${\cal O}(N^{5})$ computational scaling by changing the representation, whereas for CCSD still the ${\cal O}(N^{6})$ scaling remains. Hummel, Tsatsoulis, and Grüneis (2017); Mardirossian, McClain, and Chan (2018) This raises the question whether there exists a distinguishable cluster hierarchy which would produce more accurate results than the corresponding coupled cluster counterparts also for higher excitations, with a lower nominal scaling.

The distinguishable cluster (DC) doubles equations can be formulated using a screened Coulomb formalism, which can be used to introduce a perturbative triples correction to DCSD that was shown to be able to improve the relative energies from DCSD. Kats (2016) Here we will follow a different approach, more closely related to the original formulation of the DC approximation: we use the coupled cluster amplitude equations, remove (some) exchange terms and restore the exactness for $n$ electrons. Moreover, in this contribution we will investigate the applicability of the DC approximation to the higher order excitation, and therefore we leave the doubles CCSDT amplitude equations untouched. As the result, we do not aim here at reproducing the accuracy of DCSD in the strongly correlated regime, but rather at seeing what gain in the accuracy from the DC approximation can be expected for higher excitations.

II Theory

II.1 DC-CCSDT

As in the original DC formulation, we consider singles to be responsible for orbital relaxation only, and therefore focus our attention only to doubles and triples amplitudes. The CCSDT triple amplitudes equations using spin-orbitals can be found e.g. in Ref. Shavitt and Bartlett, 2009. Here we repeat the most relevant part of the CCSDT triples equations which contain doubles and triples connected by the two-electron repulsion integrals.

[TABLE]

where the permutation operator $\hat{\cal P}^{p}_{qr}$ antisymmetrizes a tensor as

[TABLE]

Here and in the following indices $i,j,k,l,m$ denote occupied spin-orbitals, $a,b,c,d,e$ virtual, and $p,q,r,s$ are general spin-orbital indices, and we assume the Einstein convention for the repeated indices. $T^{ij}_{ab}$, $T^{ijk}_{abc}$, and $(ld|me)$ are the doubles and triples amplitudes, and the two-electron integrals in the chemical (Mulliken) notation, respectively. Note that the integrals here are not antisymmetrized; the amplitudes, however, still possess the proper permutational antisymmetry ensured by the permutation operators in Eq. (II.1). The skeleton diagrams corresponding to these terms are shown in Figure 1.

In the distinguishable cluster approximation we remove diagrams, which correspond to intercluster exchange, in the local particle-hole symmetric fashion, Kats (2016); Kats, Usvyat, and Manby (2018) i.e., the last diagram in Figure 1. By restricting the $i,j,k$ indices in Eq. (II.1) to $i=1,j=2,k=3$ (which is the only possibility if one has only three occupied spin-orbitals and the usual restriction $i\mathchar 12604\relax j\mathchar 12604\relax k$) one can show that the terms originating from this exchange skeleton diagram (the second, seventh, and eighth terms in Eq. (II.1)) cancel with one half of the terms originating from the second and third skeleton diagram in Figure 1 (in particular, the fourth and fifth, third, and sixth terms in Eq. (II.1), respectively) in the case of three electrons. This cancellation is unique (as long as we restrict ourselves to particle-hole symmetric methods). Therefore, removing the exchange diagram and putting a factor one half in front of the two other diagrams yields a local particle-hole symmetric method that is exact for three electrons. This method will be denoted as DC-CCSDT in following. Note that the fourth and the fifth terms in Eq. (II.1) are exactly equal in the case of three electrons, and therefore in principle one can take any linear combination of them to cancel the second term. However, we chose the equally weighted linear combination which corresponds to a local-particle-hole symmetric equations.

By applying the same analysis as in Ref. Mardirossian, McClain, and Chan, 2018, it is straightforward to demonstrate that this method has a nominal scaling ${\cal O}(N^{7})$ with the systems size $N$, as opposed to ${\cal O}(N^{8})$ scaling of CCSDT. Let us focus at the first and the second term from Eq. (II.1) to show that the first term can be calculated with ${\cal O}(N^{7})$ scaling by going to the real-space representation, whilst the second term remains ${\cal O}(N^{8})$ scaling. In Eq. (3) a route is sketched how to calculate the first term of Eq. (II.1) with the ${\cal O}(N^{7})$ scaling.

[TABLE]

The same (but more practical) scaling reduction for the $\hat{T}_{2}\hat{T}_{3}$ DC-CCSDT terms can be achieved by employing the density fitting approximation, as has been demonstrated previously for DCSD. Kats and Manby (2013b) However, for some other terms in the triples equations, e.g., the $\hat{T}_{3}$-ladder diagrams, the real-space representation Mardirossian, McClain, and Chan (2018) or the tensor hypercontraction Parrish et al. (2014); Hummel, Tsatsoulis, and Grüneis (2017) have to be used to reduce the scaling to ${\cal O}(N^{7})$ in the same way as for the $\hat{T}_{2}$-ladder diagrams in the doubles amplitude equations.

For the exchange term on the other hand this simplification would not work, since no reduction of the dimensionality with respect to the system size of the amplitudes can be achieved by going to the real space.

II.2 Connection to the screened Coulomb formalism

In this section we briefly discuss a connection to the screened Coulomb formalism. In the screened Coulomb formalism from Ref. Kats, 2016 the screening is done using a direct-random-phase-approximation dressing of the fluctuation densities. The internally dressed Fock matrices have been obtained by switching to a local-particle-hole formulation of the Fock matrix. The same can be done to obtain the $\hat{T}_{2}\hat{T}_{3}$ terms. We have to extend our dressing by introducing a three-body dressed integrals,

[TABLE]

which contribute to effective two-body integrals after an exchange-Fock-like internal contraction

[TABLE]

With this screening all the $\hat{T}_{2}\hat{T}_{3}$ terms in DC-CCSDT can be obtained trough screening of linear terms in the triples equations, see Figure 2. The first diagram in Figure 1 contributes to the first diagram in Figure 2, and the second and third diagrams in Figure 2 contain the second and third diagrams from Figure 1 scaled by a factor one half.

III Test calculations

An initial test implementation of the DC-CCSDT equations from sec. II.1 was obtained by an automated term generation approach, as implemented in the GeCCo program.Köhn et al. The CCSDT, CCSDT(Q), and CCSDTQ results were obtained using either the GeCCo or MRCCKállay et al. implementations, and Molpro Werner et al. (2015) was used for the Hartree-Fock (HF) calculations and for generation of the integrals. We employed the cc-pVDZ basis in all calculations apart from the beryllium dimer calculations, for which the aug-cc-pVTZ basis was used, together with the frozen core approximation.

Our present implementation is not applicable to molecular systems with unpaired electrons, therefore we tested its accuracy for three-electron systems by calculating the BeH2 molecule with one of the hydrogen atoms at a distance of 100 Å from the beryllium atom (the HBeH angle was set to 104 degree and the other BeH bond length to 1.1 Å). The results are presented in Table 1. As expected the new method is very accurate for this quasi-three-electron system. We also see that it is much more accurate than the CCSDT method for breaking the BeH single bond. The error is comparable to the UCCSDT method, where the discrepancy compared to the FCI results comes from the difference in core orbitals between UHF and closed-shell RHF employed in UCCSDT and the rest, respectively.

To verify the higher accuracy of DC-CCSDT in breaking single bonds compared to the normal CCSDT method we calculated the potential energy curve of the F2 molecule. The plot can be found in Figure 3. While the error of CCSDT remains at 1 m$E_{h}$ in the asymptotic region, DC-CCSDT approaches the reference result. At 8 bohr separation the energy difference between the F2 molecular system and two separate fluorine atoms is 0.1 m$E_{h}$ for CCSDTQ and 0.9 m$E_{h}$ for CCSDT, which verifies the CCSDT issues in breaking single bonds.

As expected, the DC-CCSDT method shows problems for the simultaneous breaking of two bonds, since we have reintroduced all CCSD quadratic doubles terms in the doubles amplitude equations. Indeed, DC-CCSDT fails for the simultaneous dissociation of two hydrogen atoms from the water molecule, Figure 4. However, the curve in the intermediate region is noticeably closer to the CCSDTQ reference than the CCSDT curve, which demonstrates again the enhanced accuracy of the DC-CCSDT.

The beryllium dimer is another molecular system which is known to require high level theories for a quantitatively accurate description. Lesiuk et al. (2019) We used the aug-cc-pVTZ basis set and compared the accuracy of CCSDT and DC-CCSDT to CCSDTQ (which in the frozen core case, as employed here, is equivalent to FCI). The potential energy curves together with an error curve demonstrating the discrepancy of CCSDT and DC-CCSDT curves to the CCSDTQ results are shown in Figure 5. Clearly the distinguishable cluster approximation improves again upon the CCSDT results. For example, around the equilibrium distance (which in these calculations is around 4.75 bohr) the error of DC-CCSDT is approximately 0.22 m$E_{h}$, while that of CCSDT is approximately 0.35 m$E_{h}$.

Twisting the ethylene molecule is one of the problems where CCSDT is considered to be very accurate. Here we compared the accuracy of CCSDT and DC-CCSDT versus CCSDTQ results. We varied the HCCH angle from 0 to 90 degrees and fixed each of C-CH2 blocks to be planar. The remaining degrees of freedom ($r_{\rm CC}$, $r_{\rm CH}$, and HCH angle) were optimised at the DCSD level, which has been demonstrated to be very accurate for this system. Kats et al. (2015) The plots of the potential energy curve and the errors of CCSDT and DC-CCSDT compared to CCSDTQ can be found in Figure 6. Although the potential energy curves look very similar, the errors of CCSDT reach over 1.5 m$E_{h}$ (nearly 1 kcal/mol) for geometries around 80 degrees twisting. The DC-CCSDT errors in this region are comfortably low, less than 0.5 m$E_{h}$ ( $\mathchar 12604\relax 0.3$ kcal/mol), and even the largest error encountered at 90 degrees is only 0.8 m$E_{h}$ (0.5 kcal/mol).

The accuracy of relative energies is much more important for the applications than the accuracy of absolute energies. In order to investigate it, we have calculated 36 reaction energies involving closed-shell molecules. Huntington et al. (2012); Kats (2018); Kats and Tew (2019) Mean absolute deviations (MAD), root-mean squared deviations (RMSD), and maximal deviations (MaxD) with respect to the CCSDT(Q) results can be found in Table 2. The MAD and RMSD improve only slightly by going from perturbative triples to full CCSDT triples (improvement in the range $10-15\%$) However, DC-CCSDT reaction energies are much closer to CCSDT(Q) results and demonstrate an improvement of around one third upon CCSD(T) results. If one removes the reaction which shows the largest deviation for all three methods considered here (reaction 27 from Ref. Huntington et al., 2012, which is also the most exothermic reaction in the set, with over 100 kcal/mol larger reaction energy than the next one), the results become even more striking: the MAD of DC-CCSDT goes below 0.1 kcal/mol, and the improvement upon CCSD(T) results increases to $40-45\%$. The CCSDT improvement on the other hand remains around $10\%$.

IV Conclusions

We applied the distinguishable cluster approximation to the $\hat{T}_{2}\hat{T}_{3}$ terms of the CCSDT triples amplitude equations, and demonstrated that the DC-CCSDT method is superior to CCSDT in all cases we have investigated so far. Although it lacks the ability of the DCSD method to qualitatively correctly break multiple bonds, the new method is noticeably more accurate than CCSDT in breaking one or two bonds. For example, it seems to dissociate the F2 molecule to the correct limit, whilst CCSDT has some remaining error even at the interatomic distance of 8 bohr. The reaction energies are also much more accurate than from CCSDT, which barely improves upon CCSD(T). All the benchmarks clearly indicate that the improvement upon the CC hierarchy seen on the DCSD level can also be achieved for higher order methods.

Additionally this work demonstrates that a size-extensive orbital-invariant method, exact for three electrons, can be formulated with a nominal computational scaling of ${\cal O}(N^{7})$ with respect to the system size.

Further work is required to make it applicable to strong electron correlation in the similar way as DCSD, and is the focus of our ongoing research.

Acknowledgements.

DK acknowledges financial support from the Max-Planck Society.