Explicitly-correlated ring-coupled-cluster-doubles theory: Including exchange for computations on closed-shell systems

Accepted Manuscript Explicitly-correlated ring-coupled-cluster-doubles theory: Including exchange for computations on closed-shell systems Anna-Sophia Hehn, Christof Holzer, Wim Klopper PII: DOI: Reference:

S0301-0104(16)30620-6 http://dx.doi.org/10.1016/j.chemphys.2016.09.030 CHEMPH 9678

To appear in:

Chemical Physics

Received Date: Accepted Date:

2 August 2016 25 September 2016

Please cite this article as: A-S. Hehn, C. Holzer, W. Klopper, Explicitly-correlated ring-coupled-cluster-doubles theory: Including exchange for computations on closed-shell systems, Chemical Physics (2016), doi: http:// dx.doi.org/10.1016/j.chemphys.2016.09.030

This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

Explicitly-correlated ring-coupled-cluster-doubles theory: Including exchange for computations on closed-shell systems Anna-Sophia Hehn, Christof Holzer, and Wim Klopper Theoretical Chemistry Group, Institute of Physical Chemistry, Karlsruhe Institute of Technology (KIT), P. O. Box 6980, 76049 Karlsruhe, Germany

Abstract Random-phase-approximation (RPA) methods have proven to be powerful tools in electronic-structure theory, being non-empirical, computationally efficient and broadly applicable to a variety of molecular systems including small-gap systems, transition-metal compounds and dispersion-dominated complexes.

Applications are however hindered due to the slow basis-set

convergence of the electron-correlation energy with the one-electron basis. As a remedy, we present approximate explicitly-correlated RPA approaches based on the ring-coupled-cluster-doubles formulation including exchange contributions. Test calculations demonstrate that the basis-set convergence of correlation energies is drastically accelerated through the explicitly-correlated approach, reaching 99% of the basis-set limit with triple-zeta basis sets. When implemented in close analogy to early work by Szabo and Ostlund [J. Chem. Phys. 67 (1977) 4351], the new explicitly-correlated ring-coupled-cluster-doubles approach including exchange has the perspective to become a valuable tool in the framework of symmetry-adapted perturbation theory (SAPT) for the computation of dispersion energies of molecular complexes of weakly interacting closed-shell systems. Keywords: Slater-type geminals, coupled-cluster theory, ring diagrams, exchange diagrams, random-phase approximation, symmetry-adapted perturbation theory, dispersion interaction Preprint submitted to Chemical Physics

September 26, 2016

1. Introduction Wave-function methods based on molecular orbitals suffer from the slow convergence of the correlation energy with respect to the size of the basis set. While the reference energy converges exponentially with increasing cardinal number X [1], much slower asymptotics proportional to X −3 are found for the correlation energy [2], Ec (X) = Ec (∞) + a/X 3 .

(1)

The power law describes the convergence of the correlation energy Ec (X), obtained within a basis of cardinal number X, towards the basis-set limit Ec (∞). It was shown to fit for a variety of perturbation-theory and coupled-cluster methods [3] as well as for the random-phase approximation (RPA), which is a post-Kohn-Sham (KS) approach relying on molecular orbitals and which has recently attracted much attention [4, 5, 6, 7]. RPA has become popular as a computationally efficient and broadly applicable method, capturing long-range dynamic [8, 9] as well as static correlation [10]. Applications on larger molecules are however hindered due to its slow basis-set convergence. In general, basis-sets of quadruple-zeta quality are found to be necessary if at all sufficient to reach satisfying accuracy [11, 12]. However, several strategies to accelerate convergence and thus avoid the use of high angular momentum basis functions exist. Firstly, as convergence is smooth, asymptotic laws like in Eq. (1) allow to set up extrapolation schemes [3, 13, 14, 15]. Fabiano et al. for example applied various extrapolation formula to post-KS RPA calculations [12], investigating test sets of first- and second-row molecules using Dunning’s correlation-consistent basis sets cc-pVXZ [16, 17, 18]. They showed that two-point extrapolation schemes using quintuple- and sextuple-zeta basis sets with optimized parameters achieve an accuracy of 2 mEh for RPA correlation energies, which has to be compared with an error of 10 mEh for standard septuple-zeta results. Extrapolation including quadruple-zeta basis sets is however found to be of too low quality, and 2

Fabiano et al. therefore recommend to use either semiempirical extrapolations or basis-sets with higher cardinal numbers. Secondly, another strategy to circumvent the intrinsically slow basis-set convergence of wave-function methods is given by range separation [19]. Rangeseparated RPA methods [20, 21, 22] partition the Coulomb electron-electron interaction into a short-range density-functional-theory (DFT) and a longe-range RPA contribution. Convergence is accelerated as the short-range part of the correlation hole, whose accurate description requires high-angular momentum basis functions, is now treated within density-functional theory. Franck et al. showed that, for range-separated approaches, both the short- and long-range parts converge exponentially with the maximum angular momentum quantum number [23]. However, modifying the Hamiltonian by range separation always implicates a change of the corresponding basis-set limit. Hence, aiming to accelerate convergence without changing the underlying correlation method requires to go for a third strategy: explicitly-correlated (or F12) wave-function methods [24]. The idea of F12 methods is to improve the description of the Coulomb hole by taking into account geminals which explicitly depend on the inter-electronic distance [25, 26]. By doing so, correlation energies converge proportional to X −7 [27], which allows to investigate basisset limits without falling back onto large basis sets. In our recent paper on explicitly-correlated direct RPA (dRPA) [28], we demonstrate that triple-zeta basis sets are sufficient to reach sextuple-zeta accuracy for atomization energies. For interaction energies, the gain in terms of the possible reduction of the basisset size is about one cardinal number. For interaction energies, for example, quadruple-zeta quality is reached with triple-zeta basis-set size. Our investigations presented in Ref. [28] were restricted to the dRPA approach, in which exchange contributions are neglected. In consequence, the method suffers from the self-interaction error leading to an incorrect description of bond dissociation as well as a strong overestimation of correlation energies. For conventional RPA, different exchange RPA methods have been developed as a remedy, based on perturbation theory [29, 30, 31], the adiabatic connec3

tion [32, 33], or the equivalence of RPA with ring-coupled-cluster-doubles theory [34, 35, 36, 37, 38, 39]. Among the latter approaches are two approximate rCCD variants developed by Szabo and Ostlund in the 1970s to calculate correlation energies for interacting closed-shell systems [36, 37]. Recently, Toulouse and Mussard et al. adopted Szabo and Ostlund’s ideas to set up analogous rCCD variants with range separation [40, 41]. Their findings, based on the examination of rare-gas dimers, weakly interacting complexes, atomization energies and reaction barrier heights, proved these approximate rCCD approaches as very promising for general chemical applications. Based on the work by Toulouse and Mussard et al., we have developed analogous explicitly-correlated rCCD approaches for closed- and open-shell systems, denoted rCCD(F12). In contrast to the direct counterpart, drCCD(F12) [28], spin-flipped excitations have to be taken into account for open-shell rCCD(F12) theory, and spin adaptation in terms of singlet and triplet amplitudes is required for the corresponding closed-shell formulation based on spatial orbitals (cf. Section 2.1). The thereby entailed triplet and spin-flipped instabilities can however be avoided when introducing the approximations suggested by Szabo and Ostlund (cf. Section 2.2). Furthermore, the geminal amplitudes account for antisymmetry (cf. Section 2.3) in such a way that all rCCD(F12) variants yield the MP2-F12 correlation energy at second-order perturbation theory (cf. Section 2.4). The derived rCCD(F12) approaches were implemented in the Turbomole program package based on the already available drCCD(F12) code [42], allowing for the analysis of the basis-set convergence of correlation and atomization energies for an exemplary test set of small molecules (cf. Section 3). For the future, we plan to employ rCCD(F12)-based approaches for the computation of dispersion-energy contributions in the framework of (KS-based) symmetry-adapted perturbation theory (SAPT) [43, 44, 45, 46]. In some of our recent work [47, 48], we have corrected the basis-set-incompleteness error in SAPT dispersion energies by adding corrections computed at the level of explicitly-correlated second-order perturbation theory (MP2-F12) [49], but rCCD(F12)-based approaches could provide much more accurate corrections 4

for basis-set-incompleteness.

2. Theory 2.1. Symmetry of the ring-coupled-cluster-doubles amplitudes In their 1977 papers [36, 37], Szabo and Ostlund start from a reformulation of the RPA eigenvalue problem that is based on double-excitation amplitudes T and the corresponding residual equation Ω = 0, with Ω = TD + B + AT + TA + TBT .

(2)

All terms in Eq. (2) bear the important characteristic that they solely represent particle-hole interactions which correspond to ring diagrams. The method is therefore, often synonymously to RPA, called ”ring-coupled-cluster-doubles” (rCCD) approach. Depending on the definition of the matrices A and B, it is furthermore common to discriminate between the direct ring-coupled-clusterdoubles (drCCD) and full ring-coupled-cluster-doubles (rCCD) approaches. The direct approach neglects exchange and the matrices A and B therefore only inab clude Coulomb contributions, Aaj ib = haj|ibi and Bij = hab|iji. The rCCD

approach, in contrast, assumes antisymmetrized two-electron integrals, A¯aj ib = haj||ibi = haj|ibi − haj|bii ,

(3)

ab ¯ij B = hab||iji = hab|iji − hab|jii .

(4)

Here and in the following, antisymmetrized matrices are indicated with an overline and, to distinguish the direct-ring CCD amplitudes from the ring CCD ¯ for the ring CCD ones, we use the matrix T for the direct-ring CCD and T amplitudes. All spin orbitals are assumed to be real, and {i, j, . . . } denote occupied spin orbitals while {a, b, . . . } denote virtual spin orbitals, respectively. The matrix D represents for both drCCD and rCCD the zeroth-order term Dijab = εa + εb − εi − εj , with the (Kohn–Sham) molecular orbital energies ε. Note that the above notation should be understood in such a manner that

5

the row (ia) and column (jb) indices of the matrices in Eq. (2) are, for exab ab ample, as follows: (A)ia,jb = Aaj ib , (T)ia,jb = Tij and so on, where Tij is a

direct-ring-coupled-cluster amplitude for the double substitution of the occupied spin-orbital pair ij with the virtual pair ab. The restriction to ring diagrams in the residual equation has the consequence that neither the so-determined drCCD nor the rCCD amplitudes are antisymmetric with respect to the individual exchange of the index i with j or the index ba ab ¯ab ¯ba ¯ab a with b, that is, tab ij 6= −tji = −tij and tij 6= −tji = tij . This “missing anti-

symmetry” imposes important restrictions on the implementation of open- and closed-shell rCCD, as we will see in the following. First of all, the drCCD and rCCD correlation-energy equations, EcdrCCD = 21 tr [BT] ,   ¯T ¯ , EcrCCD = 41 tr B

(5) (6)

are not equivalent and the rCCD correlation energy can — in contrast to standard coupled-cluster schemes like coupled-cluster-doubles, CCD — only be for¯ This furthermore implies mulated in terms of the antisymmetrized matrix B. that open-shell rCCD implementations have to rely on the spin-integrated formulation of Eq. (6), EcrCCD =

1 4

X h

X

i JσIσ 0 Bσ 0 Aσ IσJσ BσAσ IσJσ 0 BσAσ 0 ¯AσBσ , (7) B tJσIσ + BAσBσ 0 tJσIσ 0 − BAσBσ 0 tJσIσ 0

IJAB σ,σ 0 6=σ

requiring to take into account spin-flipped excitations, parameterized by the 0

amplitudes tAσBσ Iσ 0 Jσ [50]. For spin integration, spin orbitals are split into the spin functions σ, σ 0 , . . . and the occupied or virtual spatial orbitals, {I, J, . . . } and {A, B, . . . }. Moreover, the spin-integrated residual equations reveal that sameand opposite-spin amplitudes are no longer mutually related, 0

0

AσBσ AσBσ tAσBσ IσJσ 6= tIσJσ 0 − tJσIσ 0 ,

(8)

a symmetry constraint which otherwise allows spin adaptation in terms of a single closed-shell amplitudes equation [51]. For drCCD, spin-integrated equations can nevertheless be reduced to one closed-shell counterpart, relying on 6

0

AσBσ the equivalence of same- and opposite-spin amplitudes, tAσBσ IσJσ = tIσJσ 0 . Inclu-

sion of exchange, however, lifts this drCCD symmetry relation and closed-shell rCCD schemes therefore require spin adaptation in terms of singlet and triplet amplitudes [32, 40], 0

1 AB tIJ

¯AσBσ = t¯AσBσ IσJσ + tIσJσ 0 ,

3 AB tIJ

¯AσBσ ¯AσBσ = t¯AσBσ IσJσ − tIσJσ 0 = tIσ 0 Jσ .

0

(9) 0

(10)

Closed-shell rCCD implementations are thus based on the corresponding singlet and triplet residuals, s

¯ = s TD + s B ¯ + sA ¯ s T + s Ts A ¯ + s Ts B ¯ sT , Ω

(11)

with s = 1 for the singlet and s = 3 for the triplet quantities. The singlet matrices are defined as 1 1

A¯AJ IB = 2hAJ|IBi − hAJ|BIi ,

(12)

AB ¯IJ B = 2hAB|IJi − hAB|JIi ,

(13)

while the triplet analogons are given as 3 3

A¯AJ IB = −hAJ|BIi ,

(14)

AB ¯IJ B = −hAB|JIi .

(15)

These singlet and triplet A/B matrices still carry an overline to distinguish 1

them from the direct, Coulomb-only contributions (without overline) AAJ IB = 1

AB = 2hAB|IJi. As summarized in Ref. [40], the symmetry 2hAJ|IBi and BIJ

of the rCCD amplitudes furthermore restricts the spin-adapted representation of the rCCD energy (Eq. (6)) to solely one formulation, EcrCCD = 14 tr

h

i 3 ¯3 B T+3 B T ,

1 ¯1

(16)

annulling the equivalence to another, in coupled-cluster theory commonly used expression, EcrCCD 6= 21 tr

h

1¯

B

7

!1


i T − 3T .

(17)

It should be noted that Eq. (16) is equivalent to the original plasmon formula of McLachlan and Ball [52] and that a neglect of spin-flipped excitations would lead to the energy [50] EcrCCD,noflip = 41 tr

h

i 3 ¯3 B T+ B T .

1 ¯1

(18)

2.2. Approximate variants of ring-coupled-cluster-doubles theory Eqs. (6) and (16) give exact open- and closed-shell formulations for the rCCD correlation energy, requiring on the one hand the calculation of same-spin, opposite-spin and spin-flipped amplitudes, on the other hand the calculation of singlet and triplet amplitudes. The computational effort in comparison to standard coupled-cluster schemes like CCD is increased, because for the latter, the corresponding open- and closed-shell formulations only rely either on same- and opposite-spin amplitudes or just on one type of spin-adapted, closed-shell amplitudes. More cost-efficient rCCD approaches have, however, been introduced by Szabo and Ostlund, suggesting to use Eqs. (5) and (17) as approximate closed-shell rCCD variants,
i T − 3T ,   = 12 tr 1 B1 T .

EcrCCD-SO1 = 12 tr EcrCCD-SO2

h

1¯

B

!1

(19) (20)

Toulouse and co-workers have recently constructed range-separated rCCD schemes based on the above energy expressions, Eqs. (19) and (20). They denoted the methods as rCCD-SO1 and rCCD-SO2 [40]. We adopt this nomenclature. Both approximate schemes, rCCD-SO1 and rCCD-SO2, were shown to be correct to second-order perturbation theory and to describe dispersion at the coupled Hartree-Fock level. Szabo and Ostlund preferred the rCCD-SO1 method, as the related dispersion coefficient is identical to the one given by the Casimir-Polder formula based on RPA polarizabilities. Toulouse et al. however suggested to use the more practical rCCD-SO2 approach as it solely involves singlet excitations and therefore is insensible to triplet instabilities. Furthermore, the spin-adapted formulation of Eq. (20) can be recast into an open-shell 8

analogon, whereas there is no such reformulation possible for the rCCD-SO1 approach. Analogously to the closed-shell variants, open-shell rCCD-SO2 theory reduces the computational effort in comparison to full rCCD as it does not require spin-flipped amplitudes. Based on their investigations for rangeseparated RPA, Toulouse and co-workers finally conclude that the rCCD-SO2 variant is the most accurate approach regarding atomization energies, reaction barrier heights and weak intermolecular interactions [41]. 2.3. Introducing explicitly-correlated double excitations Pursuing the ideas and suggestions of Szabo and Ostlund [36] on the one hand and of Toulouse et al. [40, 41] on the other, explicitly-correlated rCCD approaches can be set up by taking into account double excitations not only into the finite orbital basis but also into the infinite virtual basis {α, β, . . . }, parameterized through 1 xy αβ ˆ ¯ij w ¯xfy = c¯xy xyi , tαβ ij hαβ|Q12 f12 |f ij = 2 c

(21)

with the geminal amplitudes c¯xy ¯xαβ ij and the geminal functions w f y [24]. The orbitals {x, y, . . . } define the geminal space, which is often equal to the space of occupied orbitals. The infinite basis is projected onto the space of geminal functions, which depend on the inter-electronic distance r12 , incorporated through the Slater-type correlation factor f12 (for details, see for example Refs. [53, 54]). Thus, additional contributions have to be considered for the amplitudes equation as well as for the correlation energy, which bare — when schematically summarized as in Eqs. (2) and (6) — close analogy to the conventional counterparts, ¯ ← AT ¯ F12 + TF12 A ¯ + TBT ¯ F12 + TF12 BT ¯ , Ω i h   ¯ †C ¯ F12 = 1 tr V ¯ . Ec ← 1 tr BT 4

4

(22) (23)

The energy contribution is here rewritten in terms of the antisymmetrized in¯ termediate V, x f y ˆ 12 r−1 |ij − jii , V¯ij = hf xy|f12 Q 12

9

(24)

¯ The latter can be kept fixed according to the and the geminal amplitudes C. cusp conditions [55, 56, 57], implying c¯xy ij = δix δjy − δjx δiy .

(25)

Note that the rational generator Sˆxy is included in |f xyi, |f xyi = Sˆxy |xyi = ( 38 + 18 Pˆxy )|xyi .

(26)

The permutation operator Pˆxy acts on the spatial part of the orbital product only, that is, it flips the spatial parts of the geminal orbitals x and y while keeping the spin functions σx and σy unchanged, Pˆxy ϕx σx ϕy σy = ϕy σx ϕx σy . Keeping the geminal amplitudes fixed eliminates the need for an iterative solu¯ F12 . Instead, the constant residual tion of the corresponding residual equation Ω term is added as a constraint to the Lagrangian formulation of the correlation energy, Ec ←

1 h ¯ F12 ¯ LAG i tr Ω C . 4

(27)

The thereby introduced Lagrangian multipliers (¯ cLAG )xy ij are approximated as ¯xy (¯ cLAG )xy ij . ij = c ˆ 12 ensures that the geminals and the Finally, the projection operator Q geminal-excitation manifold are strongly orthogonal to the reference determinant and orthogonal to the conventional double-excitation manifold. Approximating the infinite virtual space as the sum of the conventional virtual basis set and a (finite) complementary auxiliary basis set (CABS) [58], the projection manifold is thus restricted to include double excitations into the CABS as well as mixed excitations into the CABS and virtual basis, ˆ 12 = Q

X p0 q 0

|p0 q 0 ihp0 q 0 | +

X

|p0 bihp0 b| +

p0 b

X

|aq 0 ihaq 0 | .

(28)

aq 0

ˆ 12 is usually called ansatz 2. The set {p0 , q 0 , . . . } denotes the This choice for Q CABS orbitals.

10

2.4. Working equations for open- and closed-shell rCCD(F12) Eqs. (22), (23) and (27) sketch the working equations for explicitly-correlated rCCD within ansatz 2∗ using fixed geminal amplitudes. Note that the

∗

in the

ansatz refers to the fact that the extended Brillouin condition (EBC) is assumed to be fulfilled, which avoids the occurence of coupling terms between the conventional amplitude equations and the F12 Lagrangian [49]. More explicit equations for the three approaches rCCD(F12), rCCD(F12)-SO1 and rCCD(F12)SO2 are summarized in Figure 1. Except for rCCD(F12)-SO1, for which only a closed-shell formulation is defined, equations are written out for both closedand open-shell references. The integrals over the correlation factor are hereby x f y κλ defined as fµν = hκλ|f12 |µνi. Furthermore, the F12 intermediate Bvw f, x f y ˆ 12 (Fˆ1 + Fˆ2 − εi − εj )Q ˆ 12 f12 |vwi Bvw xy|f12 Q f , f = hf

(29)

is calculated within approximation A[T + V] (for a definition, see Ref. [49]). A comparison of the full rCCD(F12) approach with the direct rCCD(F12) analogon and the approximate rCCD(F12)-SO2 ansatz can be visualized in terms of Goldstone diagrams, depicted for the open-shell equations in Figure 2. It should be noted that the chosen nomenclature is not the one which is commonly used in coupled-cluster theory, interpreting up- and down-going lines in terms of creation or annihilation operators and allowing to get the algebraic expressions by applying Wick’s theorem. Instead, elements of standard coupled-cluster diagrams are directly associated with the corresponding integrals or amplitudes and special interpretation rules to obtain the algebraic equations are defined for the different rCCD variants, see Figure 3. This is because the rCCD amplitude equations can not be derived starting from an appropriate wave function (expressed in terms of creation and annihilation operators and evaluated by applying Wick’s theorem) [59]. However, when assuming a nomenclature as summarized in Figure 3, all three methods — rCCD(F12), drCCD(F12) and rCCD(F12)-SO2 — can be recast to identically looking diagrams for both residual equations (diagrams (R1) to (R7)) and correlation energy expressions (diagrams (E1) to (E5)): Direct rCCD(F12) always relies on non-antisymmetrized 11

two-electron integrals and geminal amplitudes, while rCCD(F12) requires contrary definitions using solely the antisymmetrized counterparts. rCCD(F12)SO2 is a mixed variant, with antisymmetrized integrals for the residual equations but adopting the drCCD(F12) expression for the correlation energy, being based on non-antisymmetrized two-electron integrals. The geminal amplitudes are antisymmetrized with the exception of diagram (E4), reflecting the special, non-antisymmetrized definition of the rCCD(F12)-SO2 Lagrangian. It should be furthermore noted that all rCCD(F12) variants that include exchange terms reproduce, in accordance with the conventional schemes, the MP2-F12 correlation energy when truncated at second-order perturbation theory. The methods only start to deviate at third-order. A diagrammatic analysis of the third-order expansions of the conventional rCCD and rCCD-SO1 theories is, for example, given in Refs. [60, 30]. An F12 correction as suggested in Ref. [61], which takes into account all F12 contributions up to second-order perturbation theory, is thus equivalent for rCCD, rCCD-SO1 and rCCD-SO2. The corresponding perturbative schemes denoted ”+F12” are represented by energy diagrams (E1) to (E4), depending on the corresponding second-order amplitudes comprising diagrams (R1) to (R3).

3. Results In order to provide reference values for the final implementation of the rCCD(F12)-SO2 method in the Turbomole program package [42], the closedshell (restricted Hartree-Fock/Kohn-Sham, RHF/RKS) and open-shell (unrestricted Hartree-Fock/Kohn-Sham, UHF/UKS) explicitly-correlated rCCD approaches were all implemented in an exclusively for this purpose written program, which is based on the integral and Hartree-Fock codes of the Koala program package [62], and tested on a set of small molecules taken from Ref. [63], + + comprising CH2 , HF, F2 , N2 , CO, C2 H+ 3 , NO , C2 , O3 , CN , BN, C2 H2 , C2 H4 ,

CH4 , CO2 , H2 , H2 O and H2 O2 . Preliminary calculations with this program using the def2 basis sets of Ahlrichs and co-workers (def2-SVP, def2-TZVPP,

12

def2-QZVPP) [64] showed that the rCCD and rCCD-SO1 approaches suffered severely from triplet or spin-flipped instabilities, in contrast to the rCCD-SO2 variant. For some of the test molecules and some of the corresponding atoms, we were not able to converge the calculations, neither when modifying the DIIS algorithm (direct inversion in the iterative subspace) nor adding a shift to the quasi-Newton update. Some calculations may eventually be converged after putting more effort into finding the optimal algorithm, but the problem of instabilities is well-known [50] and provides a strong argument for favoring the rCCD-SO2 approximation over the rCCD-SO1 approximation [40]. Hence, considering the rCCD-SO2 method as the most promising approach, we implemented the rCCD-SO2 and rCCD(F12)-SO2 approaches in the Turbomole program package [42], allowing to investigate the basis-set convergence using the larger augmented correlation-consistent basis sets aug-cc-pVXZ (X=D,T,Q,5,6) of Dunning [16, 65] and the cc-pVXZ-F12 (X=D,T,Q) basis sets of Peterson [66]. The results presented in the following are thus based on two different implementations. Conventional as well as explicitly-correlated rCCD, rCCD-SO1 and rCCD-SO2 calculations were performed using the small standalone program, while only rCCD-SO2 results were obtained with the Turbomole program code. Both codes differ only in terms of efficiency, thereby allowing for a direct comparison of the -SO2 inplementations, as the computed energies are identical. All calculations reported here are based on Hartree–Fock determinants, but note that the calculations could also have been carried out as post-Kohn–Sham approach. For the performance assessment, conventional results are extrapolated using the two-point formula [2], E(XY ) =

X 3 E(X) − Y 3 E(Y ) , X3 − Y 3

(30)

where E(X) and E(Y ) are correlation energies obtained in basis sets of cardinal numbers X and Y .

The extrapolated result E(XY ) is labeled by

writing “(XY )” for the cardinal number.

Results are reported as devia-

tions from the basis-set limit, for which we chose the explicitly-correlated 13

(F12) result obtained in the largest basis set.

For the def2 basis-set se-

ries, for example, rCCD(F12)@def2-QZVPP, rCCD(F12)-SO1@def2-QZVPP and rCCD(F12)-SO2@def2-QZVPP results serve as reference for the remaining (F12) calculations in smaller basis sets, for the conventional rCCD, rCCD-SO1 and rCCD-SO2 calculations as well as for the corresponding perturbative +F12 approaches. To validate the accuracy of the basis-set limit, we performed an extrapolation scheme as suggested in our earlier work on drCCD(F12) [28] based on the extrapolation formula of Schwenke [14]. The standard deviations σ obtained in this manner provide accuracy limits for the reference data. These limits amount to 0.03 / 0.17 / 0.36 kJ/mol per valence electron for the rCCD(F12)SO2@aug-cc-pVQZ / rCCD(F12)-SO2@cc-pVQZ-F12 / rCCD(F12)-SO2@def2QZVPP limits, respectively. They are indicated in the graphs as yellow bars around the basis-set limit, ranging from −σ to σ. Details on the validation scheme are given in the supplementary material. In the following, we restrict the discussion to mean errors, the corresponding standard deviations and mean percentage errors, each reported per valence electron. More detailed statistical measures are given in the supplementary material. 3.1. Comparing explicitly-correlated rCCD, rCCD-SO1 and rCCD-SO2 Basis-set convergence as a function of the cardinal number is visualized for the various rCCD approaches in Figure 4, using the def2-SVP, -TZVPP and -QZVPP basis sets [64]. All calculations are based on Hartree-Fock references, and core orbitals are excluded from the ring-coupled-cluster correlation treatment. The complementary auxiliary basis CABS for the F12 contributions is chosen identically to the auxiliary basis for density fitting (cbas) [67]. The plot shows the correlation contribution to the total energy on the left-hand side and to the atomization energy on the right-hand side, depicting mean errors per valence electron with standard deviations as error bars. For rCCD and rCCDSO1, instabilities hampered the calculations on BN, C2 , C2 H4 , CH2 , CN+ and O3 as well as atomic oxygen and fluorine. For the latter two atoms, this implies that the corresponding atomization energies could neither be computed for the 14

rCCD-SO2 approach, because rCCD-SO2 calculations with UHF or UKS reference require the determination of all doubles amplitudes. The affected molecules were therefore excluded from the statistics and the plots in Figure 4 thus refer to different sets of test molecules depending on the particular approach. We nevertheless consider the comparison between rCCD, rCCD-SO1 and rCCDSO2 to be reliable and representative for first-period main-group compounds, since results for the remaining test molecules all showed analogous convergence behavior with comparable statistical measures. Regarding the errors in the correlation energy on the left-hand side of Figure 4, the slow basis-set convergence of all three conventional approaches, rCCD, rCCD-SO1 and rCCD-SO2, becomes apparent. The mean error for the def2SVP basis amounts up to −35.8 kJ/mol for rCCD. Analogous rCCD-SO1 and rCCD-SO2 results are −32.1 and −29.5 kJ/mol and are of the same magnitude, slightly decreasing from rCCD over rCCD-SO1 to rCCD-SO2. This cascade of decreasing errors when comparing rCCD with the -SO1 and -SO2 approximations is also found for the larger basis sets and the explicitly-correlated results, relating the smallest errors and fastest convergence to rCCD-SO2. However, it should be noted that mean percentage errors are lowest for rCCD, followed by rCCD-SO1 and rCCD-SO2, thus indicating a different ordering. Increasing the basis-set from def2-SVP over def2-TZVPP to def2-QZVPP, mean errors for all three approaches decrease gradually, first by about 20 and then by about 7 kJ/mol. The basis-set limit is finally reached for the extrapolated def2-(TQ)ZVPP result with a remaining deviation of 0.1–0.3 kJ/mol. In comparison, explicitly-correlated (F12) variants converge faster, reducing the error for the def2-SVP basis set to −7.4, −5.2 and −4.3 kJ/mol for rCCD(F12), rCCD(F12)-SO1 and rCCD(F12)-SO2, respectively, and being thus comparable to the conventional def2-QZVPP results. For the triple-zeta basis, mean errors range from −0.8 to −1.3 kJ/mol, reaching 99% of the basis-set limit. The perturbative +F12 schemes show comparable convergence behavior with slightly smaller (rCCD+F12-SO1 and rCCD+F12-SO2) or larger (rCCD+F12) mean errors and can thus be considered as efficient alternatives to the (F12) variants. 15

For rCCD and rCCD-SO2, open-shell (UHF or UKS) equations exist, allowing for a comparison of atomization energies as well, which are depicted on the right-hand side of Figure 4. Conventional rCCD and rCCD-SO2 results differ by at most 1.3 kJ/mol. For the explicitly-correlated +F12 and (F12) variants, the deviations amount to only up to 0.7 kJ/mol. The comparability of rCCD and rCCD-SO2 correlation energies is thus transferred to atomization energies. However, in contrast to correlation energies, the explicitly-correlated approaches do not converge smoothly to the basis-set limit from below, but overshoot the reference by about 0.2–0.9 kJ/mol for the def2-SVP basis. Tripleand quadruple-zeta results show again negative mean errors, reaching 99% of the basis-set limit. 3.2. Basis-set convergence of explicitly-correlated rCCD-SO2 Based on the rCCD-SO2 implementation in the Turbomole program package, the results obtained in the def2 basis sets can be transferred to the larger correlation-consistent cc-pVXZ-F12 and aug-cc-pVXZ basis sets. In both cases, the aug-cc-pwCV(X+1)Z fitting basis sets [68, 69] were used for correlation (cbas) density fitting while the aug-cc-pV(X+1)Z fitting basis sets [70] was used for Coulomb-exchange (jkbas) fitting. The cc-pVXZ-F12/OPTRI [71] and aug-cc-pVXZ-OPTRI [72] basis sets were taken as complementary auxiliary basis set (CABS). Following the cc-pVXZ-F12 series from double-zeta basis-set size to the extrapolated cc-pV(TQ)Z-F12 limit, conventional rCCD-SO2 mean errors decrease from −18.8 to −1.9 kJ/mol per valence electron. The extrapolated result thus comprises 98.1% of the basis-set limit. Both rCCD+F12-SO2 and rCCD(F12)-SO2 results are in the cc-pVDZ-F12 basis even closer to the reference value, with a mean error of −1.2 (−1.0) kJ/mol for rCCD(F12)-SO2 (rCCD+F12-SO2). Both approaches yield equivalent results for the triple- and quadruple-zeta basis sets, showing that triple-zeta results are already converged to 99.8% of the basis-set limit with a corresponding mean error of −0.2 kJ/mol. Comparably fast convergence can be found for the aug-cc-pVXZ series. Here, the explicitly-correlated rCCD(F12)-SO2@aug-cc-pVDZ result is comparable to 16

conventional rCCD-SO2@aug-cc-pV6Z with mean errors of −1.4 respectively −1.3 kJ/mol. Errors reduce further to −0.7 kJ/mol for aug-cc-pVTZ and to −0.2 kJ/mol for aug-cc-pVQZ. The perturbative rCCD+F12-SO2 variant deviates for the double-zeta basis by solely −0.3 kJ/mol from the basis-set limit. However, this seemingly better performance is explained by fortuitous error cancellation as the error successively increases to −0.5 kJ/mol for aug-cc-pVTZ. rCCD+F12@aug-cc-pVQZ results are still off by −0.1 kJ/mol.

4. Conclusions In the 1970s, the rCCD-SO1 and rCCD-SO2 approaches were suggested by Szabo and Ostlund as promising rCCD alternatives to accurately describe dispersion in a complex comprised of two interacting closed-shell systems [36, 37]. Based on their ideas, we derived and implemented the rCCD(F12), rCCD(F12)SO1 and rCCD(F12)-SO2 approaches in order to improve the intrinsically slow basis-set convergence of these ring-coupled-cluster correlation methods. For the hierarchy of def2-basis sets, convergence was shown to be equally accelerated for the original rCCD(F12) method as for the approximate rCCD(F12)-SO1 and rCCD(F12)-SO2 approaches. The latter proved to be the most promising approach being by construction insusceptible to triplet and spin-flipped instabilities. Even though convergence of the cluster amplitudes could not be achieved for some atomic calculations based on an open-shell unrestricted UHF or UKS formalism (where all doubles amplitudes are needed), we agree with the conclusion of Toulouse [40] to prefer rCCD-SO2 over rCCD and rCCD-SO1. The overall gain of the explicitly-correlated rCCD-SO2 approach is about 2–4 cardinal numbers for correlation energies obtained within the cc-pVXZ-F12 and aug-cc-pVXZ basis-set series. Convergence to within 99% of the basis-set limit is in both cases achieved already with triple-zeta basis sets. With the proposed rCCD(F12)-SO2 approach, we aim to make a first step towards an explicitly-correlated approach to the theoretical description of weakly interacting systems based on symmetry adapted perturbation theory.

17

The

rCCD(F12)-based approaches for the computation of dispersion-energy contributions have the perspective to be applied successfully in the framework of (Kohn-Sham-based) symmetry-adapted perturbation theory (SAPT). We plan to implement an explicitly-correlated SAPT approach, denoted SAPT(F12), which is based on the computation of intramolecular singlet rCCD(F12)-SO2 amplitudes. These amplitudes enter the equations for the intermolecular amplitudes, from which the dispersion energy is computed [73]. The SAPT(F12) approach shall be applied to weakly interacting closed-shell systems, and the SO2 approximation seems ideally suited for this purpose.

Acknowledgments Helpful discussions with Dr. Bastien Mussard (Universit´e Pierre et Marie Curie, Paris) are gratefully acknowledged. We thank the Deutsche Forschungsgemeinschaft (DFG) for financial support through the Priority Programme SPP 1807 “Control of London dispersion interactions in molecular chemistry” (grant No. KL 721/5-1). We are grateful to Dr. Sebastian H¨ofener (KIT) for making his Koala program code available to us.

References [1] A. Halkier, T. Helgaker, P. Jørgensen, W. Klopper, J. Olsen, Chem. Phys. Lett. 302 (1999) 437. [2] T. Helgaker, W. Klopper, H. Koch, J. Noga, J. Chem. Phys. 106 (1997) 9639. [3] A. Halkier, T. Helgaker, P. Jørgensen, W. Klopper, H. Koch, J. Olsen, A. K. Wilson, Chem. Phys. Lett. 286 (1998) 243. [4] F. Furche, J. Chem. Phys. 129 (2008) 114105. [5] H. Eshuis, J. E. Bates, F. Furche, Theor. Chem. Acc. 131 (2012) 1084. [6] A. Heßelmann, A. G¨orling, Mol. Phys. 108 (2010) 359. 18

[7] X. Ren, P. Rinke, C. Joas, M. Scheffler, J. Mater. Sci. 47 (2012) 7447. [8] J. F. Dobson, J. Wang, Phys. Rev. Lett. 82 (1999) 2123. [9] H. Eshuis, F. Furche, J. Phys. Chem. Lett. 2 (2011) 983. [10] M. Hellgren, F. Caruso, D. R. Rohr, X. Ren, A. Rubio, M. Scheffler, P. Rinke, Phys. Rev. B 91 (2015) 165110. [11] H. Eshuis, F. Furche, J. Chem. Phys. 136 (2012) 084105. [12] E. Fabiano, F. Della Sala, Theor. Chem. Acc. 131 (2012) 1278. [13] D. Feller, K. A. Peterson, J. Grant Hill, J. Chem. Phys. 135 (2011) 044102. [14] D. W. Schwenke, J. Chem. Phys. 122 (2005) 014107. [15] J. M. L. Martin, Chem. Phys. Lett. 259 (1996) 669. [16] T. H. Dunning Jr., J. Chem. Phys. 90 (1989) 1007. [17] A. K. Wilson, T. van Mourik, T. H. Dunning Jr., J. Mol. Struct. (Theochem) 388 (1996) 339. [18] D. E. Woon, T. H. Dunning Jr., J. Chem. Phys. 98 (1993) 1358. [19] J. Toulouse, F. Colonna, A. Savin, Phys. Rev. A 70 (2004) 062505. ´ [20] J. Toulouse, I. C. Gerber, G. Jansen, A. Savin, J. G. Angy´ an, Phys. Rev. Lett. 102 (2009) 096404. [21] B. G. Janesko, T. M. Henderson, G. E. Scuseria, J. Chem. Phys. 130 (2009) 081105. [22] B. G. Janesko, T. M. Henderson, G. E. Scuseria, J. Chem. Phys. 131 (2009) 034110. [23] O. Franck, B. Mussard, E. Luppi, J. Toulouse, J. Chem. Phys. 142 (2015) 074107.

19

[24] C. H¨attig, W. Klopper, A. K¨ohn, D. P. Tew, Chem. Rev. 112 (2012) 4. [25] W. Kutzelnigg, Theor. Chim. Acta 68 (1985) 445. [26] W. Kutzelnigg, W. Klopper, J. Chem. Phys. 94 (1991) 1985. [27] W. Klopper, C. C. M. Samson, J. Chem. Phys. 116 (2002) 6397. [28] A.-S. Hehn, D. P. Tew, W. Klopper, J. Chem. Phys. 142 (2015) 194106. [29] A. Gr¨ uneis, M. Marsman, J. Harl, L. Schimka, G. Kresse, J. Chem. Phys. 131 (2009) 154115. [30] A. Heßelmann, J. Chem. Phys. 134 (2011) 204107. [31] J. E. Bates, F. Furche, J. Chem. Phys. 139 (2013) 171103. ´ [32] J. G. Angy´ an, R.-F. Liu, J. Toulouse, G. Jansen, J. Chem. Theory Comput. 7 (2011) 3116. ´ [33] B. Mussard, D. Rocca, G. Jansen, J. G. Angy´ an, J. Chem. Theory Comput. 12 (2016) 2191. [34] G. E. Scuseria, T. M. Henderson, D. C. Sorensen, J. Chem. Phys. 129 (2008) 231101. ´ [35] G. Jansen, R.-F. Liu, J. G. Angy´ an, J. Chem. Phys. 133 (2010) 154106. [36] A. Szabo, N. S. Ostlund, J. Chem. Phys. 67 (1977) 4351. [37] A. Szabo, N. S. Ostlund, Int. J. Quantum Chem. S11 (1977) 389. [38] A. Heßelmann, Top Curr. Chem. 97 (2015) 365. [39] A. Heßelmann, Phys. Rev. A 85 (2012) 012517. ´ [40] J. Toulouse, W. Zhu, A. Savin, G. Jansen, J. G. Angy´ an, J. Chem. Phys. 135 (2011) 084119. ´ [41] B. Mussard, P. Reinhardt, J. G. Angy´ an, J. Toulouse, J. Chem. Phys. 142 (2015) 154123. 20

[42] TURBOMOLE V7.0, a development of University of Karlsruhe (TH) and Forschungszentrum Karlsruhe GmbH, 1989-2007; TURBOMOLE GmbH, since 2007. Available from http://www.turbomole.com (accessed April 2016). [43] B. Jeziorski, R. Moszynski, K. Szalewicz, Chem. Rev. 94 (1994) 1887. [44] K. Szalewicz, WIREs Comput. Mol. Sci. 2 (2012) 254. [45] G. Jansen, WIREs Comput. Mol. Sci. 4 (2014) 127. [46] K. U. Lao, J. M. Herbert, J. Phys. Chem. A 119 (2015) 235. [47] S. Ahnen, A.-S. Hehn, K. D. Vogiatzis, M. A. Trachsel, S. Leutwyler, W. Klopper, Chem. Phys. 441 (2014) 17. [48] J. A. Frey, C. Holzer, W. Klopper, S. Leutwyler, Chem. Rev. 116 (2016) 5614. [49] R. A. Bachorz, F. A. Bischoff, A. Gl¨oß, C. H¨attig, S. H¨ofener, W. Klopper, D. P. Tew, J. Comput. Chem. 32 (2011) 2492. [50] W. Klopper, A. M. Teale, S. Coriani, T. B. Pedersen, T. Helgaker, Chem. Phys. Lett. 510 (2011) 147. [51] G. E. Scuseria, A. C. Schreiner, T. J. Lee, J. E. Rice, H. F. Sch¨ afer III, J. Chem. Phys. 86 (1987) 2881. [52] A. D. McLachlan, M. A. Ball, Rev. Mod. Phys. 36 (1964) 844. [53] S. Ten-no, Chem. Phys. Lett. 398 (2004) 56. [54] D. P. Tew, W. Klopper, J. Chem. Phys. 123 (2005) 074101. [55] S. Ten-no, J. Chem. Phys. 121 (2004) 117. [56] S. Ten-no, J. Chem. Phys. 126 (2007) 014108. [57] D. P. Tew, W. Klopper, C. H¨attig, Chem. Phys. Lett. 452 (2008) 326.

21

[58] E. F. Valeev, Chem. Phys. Lett. 395 (2004) 190. [59] V. Lotrich, R. J. Bartlett, J. Chem. Phys. 134 (2011) 184108. [60] D. L. Freeman, M. Karplus, J. Chem. Phys. 64 (1976) 2641. [61] A.-S. Hehn, W. Klopper, J. Chem. Phys. 138 (2013) 181104. [62] S. H¨ofener, J. Comput. Chem. 35 (2014) 1716. [63] S. H¨ofener, D. P. Tew, W. Klopper, T. Helgaker, Chem. Phys. 356 (2009) 25. [64] F. Weigend, R. Ahlrichs, Phys. Chem. Chem. Phys. 7 (2005) 3297. [65] R. A. Kendall, T. H. Dunning Jr., R. J. Harrison, J. Chem. Phys. 96 (1992) 6796. [66] K. A. Peterson, T. B. Adler, H.-J. Werner, J. Chem. Phys. 128 (2008) 084102. [67] F. Weigend, M. H¨aser, H. Patzelt, R. Ahlrichs, Chem. Phys. Lett. 294 (1998) 143. [68] F. Weigend, A. K¨ohn, C. H¨attig, J. Chem. Phys. 116 (2002) 3175. [69] C. H¨attig, Phys. Chem. Chem. Phys. 7 (2005) 59. [70] F. Weigend, J. Comput. Chem. 29 (2008) 167. [71] K. E. Yousaf, K. A. Peterson, J. Chem. Phys. 129 (2008) 184108. [72] K. E. Yousaf, K. A. Peterson, Chem. Phys. Lett. 476 (2009) 303. [73] R. Moszynski, B. Jeziorski, K. Szalewicz, Int. J. Quantum Chem. 45 (1993) 409. [74] I. Shavitt and R. J. Bartlett, Many-Body Methods in Chemistry and Physics:

MBPT and Coupled-Cluster Theory, (Cambridge University

Press, New York, 2009), p. 296. 22

Open-shell equations implying c¯xy ij = δix δjy − δiy δjx :   1  †  1  F12  ¯ + tr Ω ¯ ¯ C ¯ EcrCCD(F12) = 41 tr BT + 4 tr V¯ C 4  †  1  F12  rCCD(F12)-SO2 1 1 ¯ + tr Ω C ¯ Ec = 2 tr [BT] + 2 tr V C 2 ab ac ¯kb ac ¯ kl db ¯ ab ¯ ab ¯ak cb Ω ij = Bij + Dijab tij + Aic tkj + tik Acj + tik Bcd tlj 0

0

0

0

ap ¯ lk cb pb ap ¯kb ac ¯ kl xy p b ¯lj fxfy + c¯xy ¯xy ¯xy + A¯ak ip0 c il fx kj fx ik fx f y Bp0 c tkj f y Ap0 j + tik Bcp0 c f y +c 0

0

0

0

x f y dc ¯kp x f y ¯p k cd ¯ xfy ¯ F12 )xy = B xfy c¯vw (Ω ij vw f ij + Vij + fp0 d Aic tkj + fdp0 tik Acj x f y dc kp x f y x f y p k cd x f y vw (ΩF12 )xy ij = Bvw f cij + Vij + fp0 d Aic tkj + fdp0 tik Acj

Closed-shell equations implying 1 XY c¯IJ

1 XY = 2δIX δJY − δIY δJX , 3 c¯XY IJ = −δIY δJX , cIJ = 2δIX δJY :

i i h 3¯ 1¯ 3 ¯3 3 1 + 3 V¯ † C B T+3 B T + 41 tr V¯ † C h i 1 ¯ F12 1 ¯ 3 ¯ F12 3 ¯ + 14 tr Ω C+3 Ω C h h  i
i 1¯ 1 1¯ 3¯ 1 T − 3 T + 21 tr V¯ † C − C EcrCCD(F12)-SO1 = 12 tr B h  i 1 ¯ F12 1 ¯ 3¯ C− C + 12 tr Ω h i h i   1¯ 1¯ + 1 tr 1 ΩF12 C EcrCCD(F12)-SO2 = 1 tr 1 B1 T + 1 tr 1 V † C EcrCCD(F12) = 14 tr

h

1 ¯1

2

s ¯ AB ΩIJ

2

2

s ¯AK s CB s ¯ AB s AC s ¯KB s AB = B IJ + DIJAB tIJ + AIC tKJ + tIK ACJ s ¯AK s XY P 0 B s ¯ KL s DB s XY AP 0 s ¯KB ¯IK fXY ¯KJ fXY + s tAC IK BCD tLJ + AIP 0 c g AP 0 J g + c s ¯ KL s XY P 0 B s XY AP 0 s ¯ LK s CB ¯IL fXY ¯LJ fXY + s tAC IK BCP 0 c g BP 0 C tKJ g + c

s

g s VW g g s ¯P 0 K s CD g s DC s ¯KP 0 s XY XY ¯ F12 )XY = B XY (Ω c¯IJ + V¯IJ + fPXY tKJ + fDP 0 D AIC 0 tIK ACJ IJ Vg W

1

1 P K 1 CD XY 1 V W XY XY 1 DC 1 KP + fPXY (ΩF12 )XY cIJ + 1 VIJ tKJ + fDP 0 D AIC 0 tIK ACJ IJ = BVg W

g

3 ¯ XY g VIJ

g

g

0

g

1 XY g ˆ 12 r−1 |(2IJ − JI)i, g |f12 Q with V¯IJ = hXY 12

g g |f12 Q ˆ 12 r−1 |JIi and 1 V XY g ˆ −1 = −hXY 12 IJ = 2hXY |f12 Q12 r12 |IJi.

Figure 1: Equations for rCCD(F12), rCCD(F12)-SO1 and rCCD(F12)-SO2.

23

0

(R1)

(R2)

(R4)

(R5)

(R6)

(R7)

(E1)

(E2)

(E4)

(E5)

(R3)

(E3)

Figure 2: Goldstone diagrams for rCCD(F12) variants within ansatz 2∗ using fixed geminal amplitudes. The residual equations comprise diagrams (R1) to (R7), the energy is represented by diagrams (E1) to (E5).

24

hκ|Fˆ |λi hκλ|µνi for drCCD hκλ||µνi for rCCD and rCCD-SO2 hκλ|µνi for drCCD and rCCD-SO2 hκλ||µνi for rCCD

tab ij

αβ cxy ij wx f y for drCCD xy αβ c¯ij wxfy for rCCD and rCCD-SO2

αβ cxy ij wx f y for drCCD and rCCD-SO2 xy αβ c¯ij wxfy for rCCD

Figure 3: Nomenclature for the Goldstone diagrams of Figure 2, deviating for the different rCCD(F12) methods in the definition of the two-electron interaction and the geminal amplitudes. {κ, λ, . . . } indicate the complete infinite basis, given as the sum of the occupied and infinite virtual basis, {i, j, . . . } and {α, β, . . . }. Standard interpretation rules for antisymmetrized Goldstone diagrams representing spin orbitals (see e.g. Ref. [74]) need to be modified for the given rCCD diagrams in such a way that a) an additional factor of 2 needs to be associated with non-antisymmetrized two-electron integrals or geminal amplitudes and b) distinct permutations of inequivalent external lines are accounted for by including the ab with P ˆ ab Aab = Aab + Aba . permutation operator Pˆij ij ij ij ji

25

Correlation contribution to total energy

atomization energy 5

Mean error per valence electron [kJ/mol]

0

-10

0

-20 -5 -30

rCCD-SO2 rCCD(F12)-SO2 rCCD+F12-SO2 rCCD-SO1 rCCD(F12)-SO1 rCCD+F12-SO1 rCCD rCCD(F12) rCCD+F12

-40

-10

rCCD-SO2 rCCD(F12)-SO2 rCCD+F12-SO2 rCCD rCCD(F12) rCCD+F12

-15

-50 SVP

TZVPP QZVPP

(TQ)

SVP

TZVPP QZVPP

(TQ)

Basis set def2-X

Figure 4: Mean error per valence electron in the correlation contribution to the total and the atomization energy for rCCD, rCCD-SO1 and rCCD-SO2 as well as the corresponding explicitly-correlated (F12)- and +F12-variants [in kJ/mol].

26

cc-pVXZ-F12

Mean error per valence electron [kJ/mol]

0

aug-cc-pVXZ

def2-X

0

0

-5

-5

-5

-10

-10

-15

-10

-15 -20 -20

-15

-25 -25

-30

-20 -30 -25

-35

-35

-40 -45

-40

-30

5

6

(56)

)

P

PP

P

Cardinal number X

Q

Q

(T

T

ZV

D

VP

(TQ)

Q

Q

TZ

T

SV

D

rCCD-SO2 rCCD(F12)-SO2 rCCD+F12-SO2

Figure 5: Mean error per valence electron in rCCD-SO2, rCCD(F12)-SO2 and rCCD+F12SO2 correlation energies for the series of cc-pVXZ-F12, aug-cc-pVXZ and def2-basis sets [in kJ/mol].

27