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 the results of the coupled cluster method with singles, doubles, and triples. The resulting method has a nominal computational scaling of O(N7) in the real-space representation, and is orbital invariant, size extensive, and exact for three electrons.

The coupled cluster hierarchy of methods1 is known to converge very quickly to full configuration interaction (FCI) results with comparably low excitation level for weakly correlated systems. Already triple excitations are usually sufficient to reach the 1 kcal/mol accuracy in relative energies,2 and with inclusion of higher excitations, it is possible to smoothly converge to the exact answer.3–5 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)6–10 and other methods11–23 are usually much better than the results from the coupled cluster singles and doubles (CCSD). Additionally, the DCSD equations can be implemented more efficiently than CCSD by employing density fitting or other integral decomposition techniques.24 Recently, two studies have demonstrated that DCSD can be calculated exactly with the O(N5) computational scaling by changing the representation, whereas for CCSD the O(N6) scaling still remains.25,26 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.9 Here, we investigate the DC approximation for the coupled cluster singles, doubles, and triples (CCSDT) method following 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 singles and doubles amplitude equations of CCSDT untouched. As a 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.

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 amplitudes equations using spin-orbitals can be found, e.g., in Ref. 27. Here, we repeat the most relevant part of the CCSDT equations which contain doubles and triples connected by the two-electron repulsion integrals

RabcijkP^jkiP^bca(ld|me)TadilTebcmjk(ld|me)TaeilTdbcmjkP^jki(ld|me)TdeliTabcmjkP^bca(ld|me)TdalmTebcijkP^ijkP^bca(ld|me)TadijTbeclmkP^jkiP^abc(ld|me)TabilTdecjmk+12P^ijk(ld|me)TdeijTabclmk+12P^abc(ld|me)TablmTdecijk,
(1)

where the permutation operator P^qrp antisymmetrizes a tensor as

P^qrpXpqr=XpqrXqprXrqp.
(2)

Here and in the following indices, i, j, k, l, m denote occupied spin-orbitals, a, b, c, d, e denote virtual, and p, q, r, s are general spin-orbital indices, and we assume the Einstein convention for the repeated indices. Tabij, Tabcijk, 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. (1). The skeleton diagrams corresponding to these terms are shown in Fig. 1.

FIG. 1.

T^2T^3 skeleton diagrams in CCSDT equations.

FIG. 1.

T^2T^3 skeleton diagrams in CCSDT equations.

Close modal

In the distinguishable cluster approximation, we remove diagrams, which correspond to intercluster exchange, in the local particle-hole symmetric fashion,9,28 i.e., the last diagram in Fig. 1. By restricting the i, j, k indices in Eq. (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 < j < k), one can show that the terms originating from this exchange skeleton diagram [the second, seventh, and eighth terms in Eq. (1)] cancel with one half of the terms originating from the second and third skeleton diagrams in Fig. 1 [in particular, the fourth and fifth, third, and sixth terms in Eq. (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 the following. Note that the fourth and the fifth terms in Eq. (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. 26, it is straightforward to demonstrate that this method has a nominal scaling O(N7) with the systems size N, as opposed to O(N8) scaling of CCSDT. Let us focus at the first and the second terms from Eq. (1) to show that the first term can be calculated with O(N7) scaling by going to the real-space representation, while the second term remains O(N8) scaling. In Eq. (3), a route is sketched on how to calculate the first term of Eq. (1) with O(N7) scaling

(ld|me)TadilTebcmjk=dr1dr2ϕl*(r1)ϕd(r1)1|r1r2|×ϕm*(r2)ϕe(r2)TadilTebcmjk=dr2Tbcjk(r2)dr1Tai(r1)1|r1r2|
(3)

with

Tbcjk(r2)=ϕm*(r2)ϕe(r2)Tebcmjk,Tai(r1)=ϕl*(r1)ϕd(r1)Tadil.

The same (but more practical) scaling reduction for the T^2T^3 DC-CCSDT terms can be achieved by employing the density fitting approximation, as has been demonstrated previously for DCSD.24 However, for some other terms in the triples equations, e.g., the T^3-ladder diagrams, the real-space representation26 or the tensor hypercontraction25,29 has to be used to reduce the scaling to O(N7) in the same way as for the 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.

In this section, we briefly discuss a connection to the screened Coulomb formalism. In the screened Coulomb formalism from Ref. 9, 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 T^2T^3 terms. We have to extend our dressing by introducing three-body dressed integrals

(pq|ai|bj)=(pq|kc)Tcabkij,
(4)

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

(aq|bj)̃12δpi(pq|ai|bj)=12(iq|kc)Tcabkij,(ip|bj)̃12δqa(pq|ai|bj)=12(pa|kc)Tcabkij.
(5)

With this screening all the T^2T^3 terms in DC-CCSDT can be obtained through screening of linear terms in the triples equations (see Fig. 2). The first diagram in Fig. 1 contributes to the first diagram in Fig. 2, and the second and third diagrams in Fig. 2 contain the second and third diagrams from Fig. 1 scaled by a factor one half.

FIG. 2.

Origin of T^2T^3 diagrams in screened DC-CCSDT equations.

FIG. 2.

Origin of T^2T^3 diagrams in screened DC-CCSDT equations.

Close modal

An initial test implementation of the DC-CCSDT equations from Sec. II A was obtained by an automated term generation approach, as implemented in the GeCCo program.30 The CCSDT as well as CCSDT with perturbative quadruples [CCSDT(Q)] and full quadruples (CCSDTQ) results were obtained using either the GeCCo or MRCC31 implementations, and Molpro32 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°, and the other BeH bond length was set to 1.1 Å). The results are presented in Table I. 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 CCSDT method based on unrestricted HF (UHF) orbitals (UCCSDT), 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.

TABLE I.

Results for a BeH2 molecular system with one of the BeH bonds stretched to 100 Å.

Energy, EhError, μEh
FCI −15.662290 
CCSDT −15.662155 135 
UCCSDTa −15.662305 −15 
DC-CCSDT −15.662271 19 
Energy, EhError, μEh
FCI −15.662290 
CCSDT −15.662155 135 
UCCSDTa −15.662305 −15 
DC-CCSDT −15.662271 19 
a

Calculated using the MRCC code.31 

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 Fig. 3. While the error of CCSDT remains at 1 mEh 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 mEh for CCSDTQ and 0.9 mEh for CCSDT, which verifies the CCSDT issues in breaking single bonds.

FIG. 3.

Potential energy curves of the F2 molecule from CCSDT, CCSDTQ, and DC-CCSDT. The difference plot is given with respect to CCSDTQ.

FIG. 3.

Potential energy curves of the F2 molecule from CCSDT, CCSDTQ, and DC-CCSDT. The difference plot is given with respect to CCSDTQ.

Close modal

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 (Fig. 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 method.

FIG. 4.

Potential energy curves of the simultaneous dissociation of two hydrogen atoms from a water molecule from CCSDT, CCSDTQ, and DC-CCSDT. The HOH angle was set to 107.6°. The difference plot is given with respect to CCSDTQ.

FIG. 4.

Potential energy curves of the simultaneous dissociation of two hydrogen atoms from a water molecule from CCSDT, CCSDTQ, and DC-CCSDT. The HOH angle was set to 107.6°. The difference plot is given with respect to CCSDTQ.

Close modal

The beryllium dimer is another molecular system which is known to require high level theories for a quantitatively accurate description.33 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 Fig. 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 bohrs), the error of DC-CCSDT is approximately 0.22 mEh, while that of CCSDT is approximately 0.35 mEh.

FIG. 5.

Potential energy curves of the beryllium dimer from CCSDT, CCSDTQ, and DC-CCSDT. The difference plot is given with respect to CCSDTQ.

FIG. 5.

Potential energy curves of the beryllium dimer from CCSDT, CCSDTQ, and DC-CCSDT. The difference plot is given with respect to CCSDTQ.

Close modal

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 vs CCSDTQ results. We varied the HCCH angle from 0° to 90° and fixed each of C–CH2 blocks to be planar. The remaining degrees of freedom (rCC, rCH, and HCH angle) were optimized at the DCSD level, which has been demonstrated to be very accurate for this system.8 The plots of the potential energy curve and the errors of CCSDT and DC-CCSDT compared to CCSDTQ can be found in Fig. 6. Although the potential energy curves look very similar, the errors of CCSDT reach over 1.5 mEh (nearly 1 kcal/mol) for geometries around 80° twisting. The DC-CCSDT errors in this region are comfortably low, less than 0.5 mEh (<0.3 kcal/mol), and even the largest error encountered at 90° is only 0.8 mEh (0.5 kcal/mol).

FIG. 6.

Potential energy curves of twisting the ethylene molecule from CCSDT, CCSDTQ, and DC-CCSDT. The difference plot is given with respect to CCSDTQ.

FIG. 6.

Potential energy curves of twisting the ethylene molecule from CCSDT, CCSDTQ, and DC-CCSDT. The difference plot is given with respect to CCSDTQ.

Close modal

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.10,21,34 Mean absolute deviations (MADs), root-mean squared deviations (RMSD), and maximal deviations (MaxD) with respect to the CCSDT(Q) results can be found in Table II. The MAD and RMSD improve only slightly by going from perturbative triples to full CCSDT (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. 21, 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%.

TABLE II.

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

MADRMSDMaxD
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. 21  
CCSD(T) 0.163 0.239 0.698 
CCSDT 0.150 0.208 0.529 
DC-CCSDT 0.098 0.135 0.334 
MADRMSDMaxD
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. 21  
CCSD(T) 0.163 0.239 0.698 
CCSDT 0.150 0.208 0.529 
DC-CCSDT 0.098 0.135 0.334 

We applied the distinguishable cluster approximation to the T^2T^3 terms of the CCSDT 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, while CCSDT has some remaining error even at the interatomic distance of 8 bohrs. 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 O(N7) with respect to the system size.

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

D.K. acknowledges financial support from the Max-Planck Society.

1.
J.
Čížek
,
J. Chem. Phys.
45
,
4256
(
1966
).
2.
K.
Raghavachari
,
G. W.
Trucks
,
J. A.
Pople
, and
M.
Head-Gordon
,
Chem. Phys. Lett.
157
,
479
(
1989
).
3.
A.
Tajti
,
P. G.
Szalay
,
A. G.
Császár
,
M.
Kállay
,
J.
Gauss
,
E. F.
Valeev
,
B. A.
Flowers
,
J.
Vázquez
, and
J. F.
Stanton
,
J. Chem. Phys.
121
,
11599
(
2004
).
4.
Y. J.
Bomble
,
J.
Vázquez
,
M.
Kállay
,
C.
Michauk
,
P. G.
Szalay
,
A. G.
Császár
,
J.
Gauss
, and
J. F.
Stanton
,
J. Chem. Phys.
125
,
064108
(
2006
).
5.
M. E.
Harding
,
J.
Vázquez
,
B.
Ruscic
,
A. K.
Wilson
,
J.
Gauss
, and
J. F.
Stanton
,
J. Chem. Phys.
128
,
114111
(
2008
).
6.
D.
Kats
and
F. R.
Manby
,
J. Chem. Phys.
139
,
021102
(
2013
).
7.
D.
Kats
,
J. Chem. Phys.
141
,
061101
(
2014
).
8.
D.
Kats
,
D.
Kreplin
,
H.-J.
Werner
, and
F. R.
Manby
,
J. Chem. Phys.
142
,
064111
(
2015
).
9.
D.
Kats
,
J. Chem. Phys.
144
,
044102
(
2016
).
11.
W.
Meyer
,
Int. J. Quantum Chem. Symp.
5
,
341
(
1971
).
12.
J.
Paldus
,
J.
Čížek
, and
M.
Takahashi
,
Phys. Rev. A
30
,
2193
(
1984
).
13.
P.
Piecuch
and
J.
Paldus
,
Int. J. Quantum Chem.
40
,
9
(
1991
).
14.
P.
Piecuch
,
R.
Tobola
, and
J.
Paldus
,
Phys. Rev. A
54
,
1210
(
1996
).
15.
K.
Kowalski
and
P.
Piecuch
,
J. Chem. Phys.
113
,
18
(
2000
).
16.
R. J.
Bartlett
and
M.
Musiał
,
J. Chem. Phys.
125
,
204105
(
2006
).
17.
M.
Nooijen
and
R. J.
Le Roy
,
J. Mol. Struct.
768
,
25
(
2006
).
18.
F.
Neese
,
F.
Wennmohs
, and
A.
Hansen
,
J. Chem. Phys.
130
,
114108
(
2009
).
19.
L. M. J.
Huntington
and
M.
Nooijen
,
J. Chem. Phys.
133
,
184109
(
2010
).
20.
J. B.
Robinson
and
P. J.
Knowles
,
J. Chem. Phys.
135
,
044113
(
2011
).
21.
L. M. J.
Huntington
,
A.
Hansen
,
F.
Neese
, and
M.
Nooijen
,
J. Chem. Phys.
136
,
064101
(
2012
).
22.
23.
J. A.
Black
and
P. J.
Knowles
,
Mol. Phys.
116
,
1421
(
2018
).
24.
D.
Kats
and
F. R.
Manby
,
J. Chem. Phys.
138
,
144101
(
2013
).
25.
F.
Hummel
,
T.
Tsatsoulis
, and
A.
Grüneis
,
J. Chem. Phys.
146
,
124105
(
2017
).
26.
N.
Mardirossian
,
J. D.
McClain
, and
G. K.-L.
Chan
,
J. Chem. Phys.
148
,
044106
(
2018
).
27.
I.
Shavitt
and
R. J.
Bartlett
,
Many-Body Methods in Chemistry and Physics: MBPT and Coupled-Cluster Theory
(
Cambridge University Press
,
2009
).
28.
D.
Kats
,
D.
Usvyat
, and
F. R.
Manby
,
Mol. Phys.
116
,
1496
(
2018
).
29.
R. M.
Parrish
,
C. D.
Sherrill
,
E. G.
Hohenstein
,
S. I. L.
Kokkila
, and
T. J.
Martínez
,
J. Chem. Phys.
140
,
181102
(
2014
).
30.
A.
Köhn
,
Y. A.
Aoto
,
A.
Bargholz
,
M.
Hanauer
, and
P.
Samanta
, “
The generalized tensor contraction code (GeCCo)
,” for further applications of the code, see, e.g.,
M.
Hanauer
and
A.
Köhn
,
J. Chem. Phys.
131
,
124118
(
2009
);
[PubMed]
M.
Hanauer
and
A.
Köhn
,
J. Chem. Phys.
134
,
204111
(
2011
).
[PubMed]
31.
M.
Kállay
,
Z.
Rolik
,
J.
Csontos
,
I.
Ladjánszki
,
L.
Szegedy
,
B.
Ladóczki
, and
G.
Samu
, “
MRCC, a quantum chemical program suite
,” See also
Z.
Rolik
,
L.
Szegedy
,
I.
Ladjánszki
,
B.
Ladóczki
, and
M.
Kállay
,
J. Chem. Phys.
139
,
094105
(
2013
), as well as: www.mrcc.hu.
32.
H.-J.
Werner
,
P. J.
Knowles
,
G.
Knizia
,
F. R.
Manby
,
M.
Schütz
,
P.
Celani
,
W.
Györffy
,
D.
Kats
,
T.
Korona
,
R.
Lindh
,
A.
Mitrushenkov
,
G.
Rauhut
,
K. R.
Shamasundar
,
T. B.
Adler
,
R. D.
Amos
,
A.
Bernhardsson
,
A.
Berning
,
D. L.
Cooper
,
M. J. O.
Deegan
,
A. J.
Dobbyn
,
F.
Eckert
,
E.
Goll
,
C.
Hampel
,
A.
Hesselmann
,
G.
Hetzer
,
T.
Hrenar
,
G.
Jansen
,
C.
Köppl
,
Y.
Liu
,
A. W.
Lloyd
,
R. A.
Mata
,
A. J.
May
,
S. J.
McNicholas
,
W.
Meyer
,
M. E.
Mura
,
A.
Nicklass
,
D. P.
O’Neill
,
P.
Palmieri
,
D.
Peng
,
K.
Pflüger
,
R.
Pitzer
,
M.
Reiher
,
T.
Shiozaki
,
H.
Stoll
,
A. J.
Stone
,
R.
Tarroni
,
T.
Thorsteinsson
, and
M.
Wang
, molpro, version 2015.1, a package of ab initio programs,
2015
, see http://www.molpro.net.
33.
M.
Lesiuk
,
M.
Przybytek
,
J. G.
Balcerzak
,
M.
Musial
, and
R.
Moszynski
,
J. Chem. Theory Comput.
15
,
2470
(
2019
).
34.
D.
Kats
and
D. P.
Tew
,
J. Chem. Theory Comput.
15
,
13
(
2019
).