Explicitly correlated ring-coupled-cluster-doubles theory.

THE JOURNAL OF CHEMICAL PHYSICS 142, 194106 (2015)

Explicitly correlated ring-coupled-cluster-doubles theory Anna-Sophia Hehn,1 David P. Tew,2 and Wim Klopper1 1

Institute of Physical Chemistry, Karlsruhe Institute of Technology (KIT), Fritz-Haber-Weg 2, D-76131 Karlsruhe, Germany 2 School of Chemistry, University of Bristol, Bristol BSB 1TS, United Kingdom

(Received 13 February 2015; accepted 6 May 2015; published online 18 May 2015) The connection between the random-phase approximation and the ring-coupled-cluster-doubles method bridges the gap between density-functional and wave-function theories and the importance of the random-phase approximation lies in both its broad applicability and this linking role in electronic-structure theory. In this contribution, we present an explicitly correlated approach to the random-phase approximation, based on the direct ring-coupled-cluster-doubles ansatz, which overcomes the problem of slow basis-set convergence, inherent to the random-phase approximation. Benchmark results for a test set of 106 molecules and a selection of 10 organic complexes from the S22 test set demonstrate that convergence to within 99% of the basis-set limit is reached for triple-zeta basis sets for atomisation energies, while quadruple-zeta basis sets are required for interaction energies. Corrections due to single excitations into the complementary auxiliary space reduce the basis-set incompleteness error by one order of magnitude, while contributions due to the coupling of conventional and geminal amplitudes are in general negligible. We find that a non-iterative explicitly correlated correction to first order in perturbation theory exhibits the best ratio of accuracy to computational cost. C 2015 AIP Publishing LLC. [http://dx.doi.org/10.1063/1.4921256] I. INTRODUCTION

The random-phase approximation (RPA)1–3 has proven to be a powerful approach for the description of dispersion interactions, preferable to dispersion-corrected density-functional theory (DFT)4 or spin-component scaled second-order MøllerPlesset perturbation theory (SCS-MP2).5 RPA has the advantage that the method is parameter-free, is suitable for systems with small HOMO-LUMO gaps, and captures non-additive dispersion interactions, which can become dominant in large systems.6 Recent implementations based on resolution-of-theidentity techniques7,8 scale as N 4 log N with the system size N and therefore allow to benchmark density functionals for larger systems where MP2 theory fails due to its perturbative nature. However, a current bottleneck for the application of RPA is its slow basis-set convergence. Eshuis and Furche have shown that the convergence to the basis-set limit is proportional to X −3, where X denotes the cardinal number of the basis set.9 While quadruple-zeta basis sets yield sufficient accuracy for medium- and short-range correlation, basis sets of quintuplezeta size are necessary to obtain reliable estimates for binding energies of dispersion-bound systems. With typical errors of 66% (27%) of the binding energy for triple-zeta (quadruplezeta) basis sets, convergence is as slow as for all other correlated wave-function methods based on orbitals, requiring large basis sets with high angular-momentum functions. The reason for the slow basis-set convergence is the inability of conventional wave-function methods to accurately describe the correlation hole around the electronic cusp.10 Expansions based on atom-centred Gaussian or Slater functions are not well suited to describe features depending on the interelectronic distance. Basis-set convergence can however be improved when parametrising or explicitly including this 0021-9606/2015/142(19)/194106/9/$30.00

dependence in the Hamiltonian or wave-function ansatz. Range-separated RPA approaches11,12 partition the electronelectron interaction into short-range DFT and long-range RPA correlation, resulting in a faster convergence similar to standard Kohn-Sham (KS) calculations. As the reduced basisset dependence is due to a change of the Hamiltonian itself, the improved convergence behaviour comes with a change of the basis-set limit. In contrast, explicitly correlated (F12) methods accelerate convergence to the conventional RPA basis-set limit by introducing geminals into the wave-function expansion which depend explicitly on the interelectronic distance. The so obtained convergence of the energy with basis size is proportional to X −7 13 and in practice near basis-set limit accuracy is obtained with triple-zeta basis sets, as, for example, shown for the explicitly correlated coupled-cluster theory with singles and doubles excitations (CCSD(F12)).14 Starting from CCSD(F12), it is straightforward to set up an explicitly correlated RPA method by exploiting the connection between RPA and coupled-cluster theory. Scuseria showed that RPA is equivalent to direct ring-coupled-clusterdoubles (drCCD) theory,15,16 which is a CCD method, where only those terms that correspond to ring diagrams contribute to the amplitudes equation and the correlation energy. Scuseria et al. presented an implementation based on Cholesky decomposition that scales as N 5; Rekkedal et al. have implemented drCCD gradients with an N 6 scaling.17 First results for a drCCD(F12) method have already been published in Ref. 18, where it was shown for a benzene-imidazole and a benzenepyrrole complex that the explicitly correlated wave-function ansatz improves conventional drCCD results significantly. The results for drCCD(F12) were similar to those obtained with an additive F12 correction, in earlier work denoted RPA+F12.19

142, 194106-1

© 2015 AIP Publishing LLC

This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 202.177.173.189 On: Fri, 22 May 2015 05:25:35

194106-2

Hehn, Tew, and Klopper

J. Chem. Phys. 142, 194106 (2015)

In the following, we present working equations and benchmark results for the explicitly correlated drCCD(F12) method. Working equations for drCCD(F12) can be derived from the already established CCSD(F12) method and are given in Sec. II. Further extensive drCCD(F12) calculations are summarised in Sec. III, examining the contribution of coupling terms as well as the complementary auxiliary basis (CABS)singles correction.

II. METHODS A. Random-phase approximation

Within the random-phase approximation, the correlation energy, defined as the difference between the expectation value of the exact wave-function |Ψ⟩ and the reference waveˆ function |Ψref ⟩ over the Hamilton operator H, Ecorr = ⟨Ψ| Hˆ |Ψ⟩ − ⟨Ψref | Hˆ |Ψref ⟩,

ia

Creation operator aˆ a† creates a particle in the virtual spin orbital a; annihilation operator aˆ i annihilates a particle in the occupied spin orbital i, respectively. With Qˆ † acting onto the RPA ground state, excited determinants are obtained by excitation from the ground state or de-excitation from a doubly excited determinant with the corresponding amplitudes X ai and Yia . In the literature, different derivations for the equations of X ai and Yia can be found, for example, based on response theory,1 canonical transformation,16,20 or the hypervirial theorem,21–23 all leading to the result that X ai and Yia can be obtained by solving the matrix eigenvalue problem,

XTX − YTY = 1.

0 + , −Ω-

ib KS Aia, j b = (ε KS a − ε i )δ i j δ ab + ga j ,

(5)

Bia, j b = giab j ,

(6)

where the two-electron integrals over the spin orbitals are given in Dirac notation, −1 giab j = ⟨ab|r 12 |i j⟩.

(7)

a, b, c, d, . . . denote virtual spin orbitals, i, j, k, l, . . . occupied spin orbitals, and p, q,r, s, . . . the complete basis of virtual and occupied spin orbitals, respectively. The ε KS p denote KS orbital energies. With the calculated excitation energies Ω, the RPA ground-state correlation energy is given as RPA Ecorr =

1 tr[Ω − A]. 2

(8)

(1)

is captured by introducing coupled single excitations in the wave-function expansion via the excitation operator Qˆ †,  Qˆ † = (X ai aˆ a† aˆ i + Yia aˆ i†aˆ a ). (2)

B + *X Y+ *X Y+ *Ω *A = ,−B −A- ,Y X- ,Y X- , 0 subject to the normalisation constraint

The matrices A and B are defined as

(3)

(4)

B. Connection between the random-phase approximation and direct ring-coupled-cluster-doubles theory

Already in early works on the random-phase approximation,24–26 it was shown that the RPA eigenvalue problem (Eq. (3)) can be rewritten as an iterative equation B + TA + AT + TBT = 0,

(9)

with the double-excitation amplitudes T = XY−1. Scuseria15,16 furthermore concluded that this amplitudes equation is equivalent to a CCD approach when neglecting all exchange contributions and retaining only ring diagrams, as illustrated by the first four diagrams in Figure 1. The explicit drCCD equations for closed-shell systems are given, for example, in Ref. 27. By contracting the amplitudes T with the two-electron integrals B, the same correlation energy is obtained as with RPA, 1 RPA tr[BT] = Ecorr . (10) 2 The connection to CC theory was, for example, exploited to incorporate RPA correlation in range-separated DFT drCCD Ecorr =

FIG. 1. Diagrams for the drCCD(F12) residual equations. The first four diagrams correspond to the conventional drCCD amplitudes equation. Diagrams (5)-(7) represent the additional F12 terms for the conventional residual (Eq. (16)). The F12 residual of Eq. (17) is illustrated by diagrams (8)-(11).


194106-3


J. Chem. Phys. 142, 194106 (2015)

functionals.28–30 In the following, we use this connection to set up an explicitly correlated RPA approach based on the certainty that both explicitly correlated RPA and drCCD correlation energies are equivalent in the limit of a complete basis. C. Explicitly correlated direct ring-coupled-cluster-doubles theory

2. Working equations for drCCD(F12)

1. Conventional coupled-cluster and F12 theory

In conventional coupled-cluster-doubles theory, the wavefunction is expanded in terms of excited determinants, |ΨCCD⟩ = exp(Tˆ2)|Ψref ⟩, (11)  † 1 ab with the double-excitation operator Tˆ2 = 4 i j ab t i j aˆ a aˆ i aˆ b† aˆ j , the conventional cluster amplitudes t iab j , and the reference determinant |Ψref ⟩. Within F12 theory, geminals depending on the interelectronic distance r 12 are included in the wavefunction expansion by including double excitations into the infinite basis α, β, γ, δ, . . ., |ΨCCD−F12⟩ = exp(Tˆ2 + Tˆ2′)|Ψref ⟩, (12)   with the analogous excitation operator Tˆ2′ = 14 i jα β x y cixjy w xαyβ aˆα† aˆ i aˆ †β aˆ j and the geminal amplitudes cixjy . x, y, v, w, . . . indicate the orbitals of geminal functions, which are most often taken to be occupied orbitals. The geminal functions, w xαyβ = ⟨α β|Qˆ 12 f 12| xy⟩, include the projection operator Qˆ 12 and the Slater-type correlation factor f 12,31 which depends exponentially on the interelectronic distance r 12,32 f 12 = γ {1 − exp(−γr 12)}. −1

projecting onto the space of doubly excited determinants. The residual equations are solved iteratively yielding the optimised amplitudes and hence, the correlation energy.38 To reduce the computational cost, F12 terms that are either of second and higher order perturbation theory or quadratic in the geminal amplitudes are standardly neglected, an approach denoted as (F12) approximation.39

In contrast to CCD and CCD-F12, the drCCD and analogously the drCCD(F12) equations cannot be derived from an appropriate similarity transformed Hamiltonian due to the fact that the amplitudes do not correspond to an appropriate wavefunction.40 However, starting from the CCD(F12) equations, drCCD(F12) equations are obtained when excluding all terms corresponding to ladder diagrams and neglecting all exchange contributions. Projecting onto the doubly excited determinants, the conventional vector function for the drCCD(F12) method within ansatz 2 is given as   ab KS cb KS ac KS ab Ωab Fac ti j + Fbc t i j − (ε KS i j (A2) = gi j + i + ε j )t i j c  c  ak cb ac k b l k cb + (gic t k j + t ik gc j ) + t ilad gdc tk j kc

+

p ′′k x y

+

3 1 | xy⟩ = ( + Pˆx y )|x y⟩. (14) 8 8 The permutation operator Pˆx y flips the spatial part of orbitals x and y keeping the spin functions σ x and σ y the same, Pˆx y ϕ x σ x ϕ y σ y = ϕ y σ x ϕ x σ y . It should be noted that, in order to satisfy the Coulomb singularity in direct ring-coupledcluster theory, Eq. (14) applies to all spin cases, not just the opposite spin case as for other explicitly correlated approaches.36 Within ansatz 2 (A2), the projection operator Qˆ 12 ensures that the geminals are strongly orthogonal to the basis of occupied orbitals and that the geminal manifold is orthogonal to the standard doubles manifold, Qˆ 12(A2) = 1 − Oˆ 1 Pˆ2′′ − Pˆ1′′Oˆ 2 − Pˆ1 Pˆ2, (15)  where Oˆ 1 = i |i⟩⟨i| is the projector onto the occupied orbitals,  Pˆ1′′ = p′′ |p′′⟩⟨p′′| the projector onto the CABS,37 and Pˆ1  = p |p⟩⟨p| the projector onto the total basis of occupied and virtual orbitals, respectively. p′′, q ′′,r ′′, s ′′, . . . indicate the CABS orbitals. Based on the explicitly correlated wave-function ansatz, xy CCD-F12 equations for the amplitudes t iab j and ci j can be derived starting from a similarity transformed Hamiltonian and

+



′′

(t ilac gcl kp′′ f xpy b ckx jy + cilx y f xayp g pl k′′c t kcbj ) ′′

 p ′′ck l x y

(13)

The exponent γ is optimised depending on the chosen basis set.33–35 In the following, the integral over the Slater geminal y⟩, where is denoted by f xpq y = ⟨pq| f 12 | x

k cl d p b xy x y ap kb (giak p ′′ f x y ck j + cik f x y g p ′′ j ) ′′



′′

Cxaby cixjy .

(16)

(i j)

xy

The vector function for the explicitly correlated amplitudes cixjy can be obtained analogously by projecting onto the determinants with double excitations into the geminal basis,   (i j) ˜ x y v w (i j) x y ab Bv w ci j + Cab t i j Ωixjy (A2) = Vixj y + vw

+



p ′′k c d ( f px′′yd gic tk j

ab cd k p + t ik gdj f cxpy′′). (17) ′′

p ′′c dk KS For canonical orbitals, the KS matrix elements Fpq are given by the corresponding orbital energies, KS Fpq = δ pq ε KS p .

(18)

Within ansatz 2, the two residual equations are coupled xy through the intermediates (i j)Cab and (i j)Cxaby , defined as KS ˆ y⟩. Cxaby = ⟨ab|(Fˆ1KS + Fˆ2KS − ε KS i − ε j )Q 12 f 12 | x

(i j)

(19)

Assuming the extended Brillouin condition, FpKSp′′ ≈ 0, coupling is however avoided. This approach is denoted ansatz 2∗ (A2∗). The intermediates Vixj y and (i j) B˜ vxwy are given by −1 Vixj y = ⟨ xy | f 12Qˆ 12r 12 |i j⟩, (20) KS KS KS KS (i j) ˜ x y Bv w = ⟨ xy | f 12Qˆ 12(Fˆ1 + Fˆ2 − ε i − ε j )Qˆ 12 f 12| v w⟩.

(21)


194106-4


J. Chem. Phys. 142, 194106 (2015)

In our implementation, the construction of the matrices Vixj y and (i j) B˜ vxwy follows Ref. 41 and uses robust density fitting. The KS operator Fˆ KS appears in the matrices B˜ and C, Fˆ KS = Tˆ + Vˆ + Jˆ + Vˆxc ,

(22)

where the exchange-correlation potential Vˆxc replaces the exchange operator Kˆ of Hartree-Fock (HF) theory and commutes with the correlation factor for generalised-gradientapproximation functionals. In this case, approximations A and B13 differ only in the treatment of terms that for the fixedamplitudes ansatz are zero under the generalised Brillouin condition, FiKS ≈ 0. Additionally, the [T + V] and [F + K] p ′′ commutation approximations are related by [F + K] = [T + V + J + V xc ]. xy With the amplitudes t iab j and ci j , the correlation energy can be obtained as 1  ab ab 1  x y x y drCCD(F12) c V t g + Ecorr = 2 aib j i j i j 2 xi y j i j i j 1  xy xy + c¯ Ω . 2 xi y j i j i j

(23)

c¯ixjy is a Lagrangian multiplier for the residual term Ωixjy and is approximated as c¯ixjy = cixjy when using fixed geminal amplitudes. To visualise the correspondence between conventional drCCD and drCCD(F12), Goldstone diagrams for both residual equations, Eqs. (16) and (17), are given in Figure 1. The corresponding nomenclature is explained in Figure 2, where it should be noted that, due to the lack of exchange, the diagrams given here are not antisymmetrised. Based on these diagrams, the comparison of other possible explicitly correlated direct ring-coupled-cluster-doubles approaches is straightforward: direct rCCD[F12] is obtained by neglecting diagram (7) and is equivalent to direct rCCD(F12∗).42 Direct rCCD-F12a and rCCD-F12b approaches43,44 are identical, both neglecting diagrams (6), (7) and (11). The perturbational F12 correction of RPA+F12 within ansatz 2∗ consists solely of diagrams (9) and (10), representing the F12 part of the well-established MP2-F12 method. RPA+F12 within ansatz 2 would require subsequent evaluation of diagram (8) using the RPA doubles amplitudes, which inhibits a simple additive correction. For an extensive diagrammatic analysis of the full CCSD-F12 approach, the reader is furthermore referred to Ref. 45.

3. CABS-singles correction

Within the framework of the CABS-singles correction, single excitations into the CABS basis are included to reduce the error in the expectation value of the KS determinant over the Hamiltonian. Explicit equations for a HF reference are given in Refs. 36 and 46. For a KS reference, the HF matrix is replaced by the KS Fock matrix in the singles residual equation,  ′ KS −1 KS t ia (A2) = − [FaKS (24) ′b ′ − δ a ′b ′ε i ] Fb ′i , b′

with a ′, b′, c ′, d ′, . . . denoting the combined virtual and CABS basis. The CABS-singles correction to the reference energy is ′ obtained by contracting the singles amplitudes t ia with the HF matrix, ′

a E∆CABS = FaHF ′i t i .

(25)

It should be noted that neither the common nor the generalised Brillouin condition is fulfilled because the HF matrix is built from KS orbitals.

III. RESULTS

As pointed out by Eshuis and Furche, basis-set convergence of RPA is quite different for non-covalent and covalent interactions and has to be investigated separately.9 To validate the performance of the drCCD(F12) method, we therefore discuss atomisation energies for a test set of 106 molecules as well as interaction energies for 10 complexes of the S22 test set regarding convergence to the basis-set limit. The 106 molecules of the first test set are compounds containing the elements H, C, N, O, and F representing a wide range of bonding situations used previously to assess F12 methods.47 In contrast, the S22 test set is specifically constructed to contain non-covalent interactions, including complexes with hydrogen bonds, predominant dispersion, and mixed noncovalent contributions.48 All calculations are performed with the drCCD(F12) code as implemented in the 12 module of the TURBOMOLE program package.14,49,50 The total energy is always calculated as the sum of the energy expectation value of the KS determinant and the correlation energy, Etotal = ⟨KS| Hˆ |KS⟩ + Ecorr,

(26)

FIG. 2. Nomenclature for the Goldstone diagrams in Figure 1.


194106-5


J. Chem. Phys. 142, 194106 (2015)

where the KS orbitals are obtained from a self-consistent PBE calculation. In the following, the expectation value of the KS determinant is denoted “PBE-orbital reference” to highlight the difference to energies obtained within density-functional theory. Moreover, error statistics are presented with respect to drCCD(F12) calculations within ansatz 2 for the largest basis set considered, i.e., for the aug-cc-pV5Z51,52 or the jun-ccpV(Q+d)Z basis,53,54 based on the assumption that these calculations are converged to the basis-set limit. Analogous analyses with respect to extrapolated basis-set limits yield comparable errors. To estimate the accuracy of the “converged” drCCD(F12) limits, we compare them to limits obtained by extrapolation. The extrapolation is performed in two steps. First, scaling factors FXY are determined, as proposed by Schwenke,55 by inverting the general extrapolation formula, E(∞) = [E(Y ) − E(X)]FXY + E(X),

(27)

taking the converged drCCD(F12) results as accurate reference E(∞), and corresponding drCCD calculations in two succeeding basis sets of cardinal numbers X and Y as data points for the two-point extrapolation. Reference and correlation contributions are treated separately, so that E(X) and E(Y ) correspond to either PBE-orbital reference or CC correlation energies. In general, drCCD(F12) results within ansatz 2 represent an adequate basis-set limit for both contributions even when compared to conventional drCCD results in the same basis, as the latter neither include F12 terms nor the CABS-singles correction. For the test set of 106 molecules with a drCCD(F12)[A2]+CABS/aug-cc-pV5Z reference, we extrapolated the factor for the correlation part from drCCD results in the aug-cc-pV5Z and aug-cc-pV6Z basis. The calculation of the factor for the reference contribution was however limited to basis sets of quadruple- and quintuplezeta size as the aug-cc-pV5Z and aug-cc-pV6Z results were both converged to the basis-set limit within computational

TABLE I. Scaling factors for the basis-set limit extrapolation according to Eq. (27) and corresponding standard deviations per valence electron (in kJ/mol). Reference energy Test set 106 molecules S22

Correlation energy

F (XY)

σ

F (XY)

σ

0.7628 (Q5) 1.0330 (TQ)

0.0126 0.0010

2.3837 (56) 1.7764 (TQ)

0.0223 0.0050

accuracy and therefore allowed no extrapolation. For the 10 complexes of the S22 test set, both factors were calculated from drCCD(F12) results in the jun-cc-pV(Q+d)Z basis and corresponding drCCD results in the triple- and quadruple-zeta basis sets. The scaling factors as obtained by linear regression over all test molecules are given in Table I. Based on these universal scaling factors, Eq. (27) can, in a second step, be used in reverse to calculate the basis-set limits E(∞), again with drCCD results for the same two basis sets with cardinal numbers X and Y and corresponding energies E(X) and E(Y ). The deviation with respect to the converged drCCD(F12) result reflects the inherent error of the extrapolation scheme and should provide a reasonable estimate for the accuracy of the basis-set limit. The corresponding standard deviations are given in Table I and marked in yellow in Figures 3 and 5. Note that the total error bar of the basis-set limit is, in accordance with the error bars for atomisation and interaction energies, twice the standard deviation σ, from −σ to σ. In the following, errors are reported per valence electron to eliminate the dependence on the system size. It should be noted that the number of valence electrons ranges between 16 and 38 for the complexes of the S22 test set and between 2 and 36 for the test set of 106 molecules, with corresponding mean values of 28 and 18 valence electrons, respectively. We restrict the discussion to mean errors with the corresponding standard deviations as these are sufficient to interpret the

FIG. 3. Test set of 106 molecules: basis-set convergence of the reference and correlation energy contributions to the atomisation energies for the aug-ccpVXZ basis sets. “PBE orb.” denotes the PBE-orbital reference.


194106-6


J. Chem. Phys. 142, 194106 (2015)

results. However, for the sake of completeness, we provide detailed error statistics for both the S22 test set and the test set of 106 molecules including mean absolute deviations, rootmean-square errors, standard deviations, and maximum errors in Tables S1–S4 in the supplementary material.56 A. Atomisation energies

Atomisation energies for the test set of 106 molecules were calculated using the augmented correlationconsistent basis sets aug-cc-pVXZ,51,52 with the aug-ccpwCV(X+1)Z57,58 and aug-cc-pV(X+1)Z59 basis sets as cbas and jkbas, respectively. Furthermore, aug-cc-pVXZ/OPTRI basis sets60 were used as CABS basis for the F12 integrals and the CABS-singles correction. Figure 3 visualises the basis-set convergence as a function of the cardinal number for both the reference and correlation energy contributions to the atomisation energy. On the left, mean errors in the reference energy are plotted with respect to the PBE-orbital reference calculation in the aug-cc-pV5Z basis including the CABS-singles correction for ansatz 2. With 0.01 kJ/mol, the standard deviation of this limit with respect to the extrapolated value obtained from Eq. (27) indicates that PBE-orbital reference as well as PBE-orbital reference +CABS calculations in the largest quintuple- and sextuplezeta basis sets are practically converged to the basis-set limit, which qualifies both the conventional and the CABS-singles corrected result as reference. For the smaller double- and triple-zeta basis sets, the error in the PBE-orbital reference energy is significant, rising up to −2.10 and −0.31 kJ/mol. The CABS-singles correction within ansatz 2∗ accelerates basis-set convergence, with a remaining error of −0.34 kJ/mol for the aug-cc-pVDZ and of −0.08 kJ/mol for the aug-cc-pVTZ basis. Convergence is not smooth due to the perturbative nature of the correction, resulting in a larger error of −0.09 kJ/mol for the quadruple-zeta basis which decreases again to 0.003 kJ/mol for the quintuple-zeta basis. Coupling contributions as introduced within ansatz 2 are negligible; deviations from ansatz 2∗ amount only up to 0.05 kJ/mol. Even though the CABSsingles correction improves PBE-orbital reference errors by one order of magnitude, representing a gain of about one cardinal number for the smaller basis sets, it should be noted that the basis-set error in the KS reference is still of the same order of magnitude and sometimes even larger than the error in the correlation energy. Convergence of the correlation contribution is much faster, as depicted in the graph on the right. The conventional drCCD results converge only slowly to the drCCD(F12)[A2]+CABS/aug-cc-pV5Z limit, with mean errors ranging from −6.11 to −0.33 kJ/mol when going from double- to sextuple-zeta basis sets. In comparison, all explicitly correlated approaches reach the basis-set limit already for the aug-cc-pVDZ basis. Corresponding mean errors of −0.17, 0.26, and −0.05 kJ/mol for ansatz 2∗, ansatz 2, and RPA+F12 are even smaller than those of the conventional aug-cc-pV6Z result. More detailed analysis shows that convergence for drCCD(F12) within ansatz 2∗ is not smooth: double- and quintuple-zeta basis sets underestimate the limit by −0.17 and −0.04 kJ/mol, whereas triple- and quadruple-zeta

FIG. 4. Test set of 106 molecules: normal distributions of the error per valence electron in the atomisation energies.

basis sets overshoot the reference value by 0.15 and 0.05 kJ/mol. A monotonic convergence from above is obtained for drCCD(F12) within ansatz 2 even though the mean errors are about twice as large as those for ansatz 2∗. The additive RPA+F12 approach, which only includes drCCD(F12) terms to first order perturbation theory, yields slightly smaller errors than drCCD(F12) within ansatz 2∗, ranging from −0.05 to 0.10 kJ/mol. Even though the accordance is reasonable as missing higher-order terms are expected to be small, it should be mentioned that the relatively good performance of the perturbative treatment RPA+F12 is due to fortuitous error cancellation, as already pointed out in our previous work.19 For total atomisation energies, it can be summarised that convergence is drastically accelerated: explicitly correlated aug-cc-pVTZ results are closer to the basis-set limit than conventional aug-cc-pV6Z results with error distributions that are both narrower and closer to zero, as depicted in Figure 4. However, it should be noted that the relatively good performance of the aug-cc-pVTZ basis is partly due to error cancellation of the reference and correlation energies. For the drCCD(F12)[A2]+CABS method in the aug-cc-pVDZ basis, errors of 0.26 and −0.28 kJ/mol cancel analogously, resulting in a small deviation of 0.02 kJ/mol, which explains the outstanding performance of ansatz 2 for the smallest basis set. For all other methods, errors of the aug-cc-pVDZ basis add up, leading to comparatively broad normal distributions. The effect on the larger quadruple- and quintuple-zeta basis sets is negligible; both yield nearly equivalent and sufficiently narrow distributions. B. Interaction energies

As RPA is known to perform well for dispersion interactions, we proceed our discussion on drCCD(F12) regarding


194106-7


J. Chem. Phys. 142, 194106 (2015)

FIG. 5. Ten molecules of the S22 test set: basis-set convergence of the interaction energies with and without counterpoise correction (CP) for the jun-ccpV(X+d)Z basis sets.

interaction energies of the S22 test set.48 Due to the N 6 scaling of drCCD(F12), calculations were limited to the smallest 10 complexes of the test set, which nevertheless still include hydrogen-bonded, dispersion-dominated, and mixed complexes in a balanced ratio. Furthermore, we use the comparatively small seasonal basis sets jun-cc-pV(X+d)Z, based on the recommendation of Papajak and Truhlar.53,54 The auxiliary basis sets cbas, jkbas, and CABS were chosen as for the test set of 106 molecules: aug-cc-pwCV(X+1)Z, aug-cc-pV(X+1)Z, and aug-cc-pVXZ/OPTRI, respectively. Basis-set convergence of interaction energies for the selection of 10 complexes is illustrated in Figure 5. The graph is split in two parts: mean errors without the counterpoise (CP) correction 61 are depicted on the left. They are compared to corresponding results including the counterpoise correction on the right. For both treatments, plotted errors are calculated with respect to the counterpoise-corrected drCCD(F12)[A2]+CABS result in the jun-cc-pV(Q+d)Z basis. The extrapolated standard deviation of 0.005 kJ/mol (see Table I) was obtained accordingly from drCCD calculations including the counterpoise correction as this was shown to be crucial for a monotonic convergence behaviour and thus for extrapolation techniques.62 This important aspect of the counterpoise correction is also reflected in the graph on the left showing that the non-corrected results lack a systematic convergence when going from double- via triple- to quadruple-zeta basis-set size. Moreover, most results overshoot the basis-set limit resulting in mean errors with positive sign and convergence to the basis-set limit is rather slow: the errors lie for all methods in the range of 0.017–0.026 kJ/mol for the triplezeta basis while corresponding quadruple-zeta results amount up to 0.020 kJ/mol. Both findings confirm our conclusion that the errors, which are relatively small in comparison to the counterpoise-corrected results, are solely fortuitous. Monotonic convergence is obtained when including the counterpoise correction, as shown in the graph on the right.

The correction is at least about a factor of 2 smaller for the explicitly correlated results compared to the conventional ones, emphasising that F12 methods generally reduce the basis-set superposition error.10,18 However, also for the counterpoise-corrected interaction energies, basis-set convergence is only slightly improved through the explicitly correlated wave-function ansatz and not comparable to the drastically accelerated convergence of atomisation energies. Mean deviations for drCCD amount up to −0.254 kJ/mol for the double-zeta, to −0.097 kJ/mol for the triple-zeta, and to −0.038 kJ/mol for the quadruple-zeta basis. In comparison, drCCD(F12) within ansatz 2⋆ including the CABS-singles

FIG. 6. Ten molecules of the S22 test set: normal distributions of the error per valence electron in the counterpoise-corrected interaction energies.


194106-8


correction shows an error of −0.113 kJ/mol for the doublezeta basis, which decreases to −0.023 and −0.003 kJ/mol for the triple- and quadruple-zeta basis, respectively. All other F12 approaches yield comparable mean errors which deviate by at most 0.006 kJ/mol, indicating that both CABSsingles correction and coupling contributions are negligibly small for the counterpoise-corrected interaction energies. Only the standard deviations for the double-zeta basis of both drCCD(F12)[A2∗] and RPA+F12 are about twice as large as those for drCCD(F12) including the CABS-singles correction. Convergence to within the error estimate of the basis-set limit is reached for the jun-cc-pV(Q+d)X basis with an error of −0.003 kJ/mol for all three methods. The corresponding normal distributions, plotted in Figure 6, are thus similar and highlight the fact that the overall gain through the explicitly correlated ansatz is about one cardinal number. Conventional drCCD results in a basis of cardinal number X are thus comparable to drCCD(F12) calculations performed in a basis of cardinal number X − 1. IV. CONCLUSIONS

As outlined in the beginning, the efficiency of RPA and drCCD approaches can be enhanced by resolution-of-theidentity techniques and Cholesky decomposition. Even though current implementations already scale as N 4 log N or N 5, broad applications on larger systems are hindered due to the slow basis-set convergence requiring, in general, basis sets of quadruple- or even quintuple-zeta size. The central purpose of this work was to outline a possible option, namely, the explicitly correlated drCCD method, to overcome this bottleneck. The presented benchmark results demonstrate the unequal performance of F12 for atomisation and interaction energies: while triple-zeta basis sets yield sextuple-zeta accuracy for atomisation energies, they can only gain one cardinal number for interaction energies, reaching at most quadruple-zeta

J. Chem. Phys. 142, 194106 (2015)

quality. The percentage errors, as visualised in Figure 7, furthermore highlight that triple-zeta results for atomisation energies are converged to within 99% of the basis-set limit, while corresponding interaction energies still lack 7% of the basis-set limit. For the latter, quadruple-zeta basis sets are required to reach the desired 99% limit. These results are in line with earlier findings: for atomisation energies, basis-set convergence of the presented F12 approaches is comparable to the already established MP2F12 method,10 reflecting the close connection to the additive correction of RPA+F12 and the relatively small contribution of drCCD(F12) terms to higher order perturbation theory. This rapid convergence is in contrast to the relatively poor performance of F12 for dispersion interactions, which is inherent to the F12 wave-function ansatz, being designed to capture short-range correlation instead of long-range phenomena.63 Monomers are thus more efficiently described than the complex itself, leading to relatively large errors in the binding energy. In accordance with Eshuis and Furche,9 we therefore conclude that quadruple-zeta basis sets are required for interaction energies, while triple-zeta basis sets yield sufficient accuracy for total energies as well as atomisation energies. Based on the comparison of the presented F12 approaches, we furthermore consider the perturbative treatment RPA+F12 as an efficient alternative to drCCD(F12) due to its more favourable scaling of N 5. We recommend to improve total energies by the CABS-singles correction, as the computational cost is negligible and basis-set incompleteness errors of atomisation energies were shown to be reduced by one order of magnitude. As both drCCD(F12) and RPA+F12 are based on the ring-coupled-cluster formalism, correlation energies are in the limit of a complete basis equivalent to analogous explicitly correlated approaches based on the RPA eigenvalue equations. Basis-set convergence, however, is not comparable due to the differences in the working equations, which hamper

FIG. 7. Mean errors in percent for the atomisation energies and the counterpoise-corrected interaction energies.


194106-9


a direct adoption of the explicitly correlated ansatz. Our results therefore allow no prediction on the presumable convergence behaviour of other possible explicitly correlated RPA approaches, which will be further investigated. ACKNOWLEDGMENTS

A.S.H. gratefully acknowledges support by the Fonds der Chemischen Industrie (FCI) through a Chemiefonds fellowship including a visiting researcher’s fellowship for the University of Bristol. 1H.

Eshuis, J. Bates, and F. Furche, Theor. Chem. Acc. 131, 1084 (2012). Heßelmann and A. Görling, Mol. Phys. 109, 2473 (2011). 3X. Ren, P. Rinke, C. Joas, and M. Scheffler, J. Mater. Sci. 47, 7447 (2012). 4S. Grimme, Wiley Interdiscip. Rev.: Comput. Mol. Sci. 1, 211 (2011). 5S. Grimme, J. Chem. Phys. 118, 9095 (2003). 6R. A. DiStasio, Jr., O. A. von Lilienfeld, and A. Tkatchenko, Proc. Natl. Acad. Sci. U. S. A. 109, 14791 (2012). 7H. Eshuis, J. Yarkony, and F. Furche, J. Chem. Phys. 132, 234114 (2010). 8A. M. Burow, J. E. Bates, F. Furche, and H. Eshuis, J. Chem. Theory Comput. 10, 180 (2014). 9H. Eshuis and F. Furche, J. Chem. Phys. 136, 084105 (2012). 10C. Hättig, W. Klopper, A. Köhn, and D. P. Tew, Chem. Rev. 112, 4 (2012). 11J. Toulouse, I. C. Gerber, G. Jansen, A. Savin, and J. G. Ángyán, Phys. Rev. Lett. 102, 096404 (2009). 12J. Toulouse, W. Zhu, J. G. Ángyán, and A. Savin, Phys. Rev. A 82, 032502 (2010). 13W. Kutzelnigg and W. Klopper, J. Chem. Phys. 94, 1985 (1991). 14D. P. Tew, W. Klopper, C. Neiss, and C. Hättig, Phys. Chem. Chem. Phys. 9, 1921 (2007). 15G. E. Scuseria, T. M. Henderson, and D. C. Sorensen, J. Chem. Phys. 129, 231101 (2008). 16G. E. Scuseria, T. M. Henderson, and I. W. Bulik, J. Chem. Phys. 139, 104113 (2013). 17J. Rekkedal, S. Coriani, M. F. Iozzi, A. M. Teale, T. Helgaker, and T. B. Pedersen, J. Chem. Phys. 139, 081101 (2013). 18S. Ahnen, A.-S. Hehn, K. D. Vogiatzis, M. A. Trachsel, S. Leutwyler, and W. Klopper, Chem. Phys. 441, 17 (2014). 19A.-S. Hehn and W. Klopper, J. Chem. Phys. 138, 181104 (2013). 20A. L. Fetter and J. D. Walecka, Quantum Theory of Many-Particle Systems (Dover Publications, New York, 2003). 21D. J. Rowe, Rev. Mod. Phys. 40, 153 (1968). 22A. E. Hansen and T. D. Bouman, Mol. Phys. 37, 1713 (1979). 23P. Ring and P. Schuck, The Nuclear Many-Body Problem (Springer-Verlag, Berlin, 1980). 24N. Ostlund and M. Karplus, Chem. Phys. Lett. 11, 450 (1971). 25E. A. Sanderson, Phys. Lett. 19, 141 (1965). 26A. Szabo and N. S. Ostlund, J. Chem. Phys. 67, 4351 (1977). 27J. G. Ángyán, R.-F. Liu, J. Toulouse, and G. Jansen, J. Chem. Theory Comput. 7, 3116 (2011). 2A.

J. Chem. Phys. 142, 194106 (2015) 28B.

G. Janesko, T. M. Henderson, and G. E. Scuseria, J. Chem. Phys. 130, 081105 (2009). 29B. G. Janesko, T. M. Henderson, and G. E. Scuseria, J. Chem. Phys. 131, 034110 (2009). 30J. Toulouse, W. Zhu, A. Savin, G. Jansen, and J. G. Ángyán, J. Chem. Phys. 135, 084119 (2011). 31S. Ten-no, Chem. Phys. Lett. 398, 56 (2004). 32D. P. Tew and W. Klopper, J. Chem. Phys. 123, 074101 (2005). 33J. G. Hill, S. Mazumder, and K. A. Peterson, J. Chem. Phys. 132, 054108 (2010). 34J. G. Hill and K. A. Peterson, Phys. Chem. Chem. Phys. 12, 10460 (2010). 35K. E. Yousaf and K. A. Peterson, J. Chem. Phys. 129, 184108 (2008). 36D. P. Tew and W. Klopper, Mol. Phys. 108, 315 (2010). 37E. F. Valeev, Chem. Phys. Lett. 395, 190 (2004). 38T. Helgaker, P. Jørgensen, and J. Olsen, Molecular Electronic-Structure Theory (Wiley, Chichester, 2000). 39H. Fliegl, W. Klopper, and C. Hättig, J. Chem. Phys. 122, 084107 (2005). 40V. Lotrich and R. J. Bartlett, J. Chem. Phys. 134, 184108 (2011). 41R. A. Bachorz, F. A. Bischoff, A. Glöß, C. Hättig, S. Höfener, W. Klopper, and D. P. Tew, J. Comput. Chem. 32, 2492 (2011). 42C. Hättig, D. P. Tew, and A. Köhn, J. Chem. Phys. 132, 231102 (2010). 43T. B. Adler, G. Knizia, and H.-J. Werner, J. Chem. Phys. 127, 221106 (2007). 44G. Knizia, T. B. Adler, and H.-J. Werner, J. Chem. Phys. 130, 054104 (2009). 45J. Noga and W. Kutzelnigg, J. Chem. Phys. 101, 7738 (1994). 46G. Knizia and H.-J. Werner, J. Chem. Phys. 128, 154103 (2008). 47F. A. Bischoff, S. Wolfsegger, D. P. Tew, and W. Klopper, Mol. Phys. 107, 963 (2009). 48P. Jureˇ ˇ cka, J. Šponer, J. Cerný, and P. Hobza, Phys. Chem. Chem. Phys. 8, 1985 (2006). 49TURBOMOLE V6.6; a development of University of Karlsruhe (TH) and Forschungszentrum Karlsruhe GmbH, 1989-2007; TURBOMOLE GmbH, since 2007. Available from http://www.turbomole.com (accessed May 2014). 50R. A. Bachorz, “Implementation and application of the explicitly correlated coupled-cluster method in turbomole,” Ph.D. thesis (Universität Karlsruhe (TH), 2009). 51T. H. Dunning, Jr., J. Chem. Phys. 90, 1007 (1989). 52R. A. Kendall, T. H. Dunning, Jr., and R. J. Harrison, J. Chem. Phys. 96, 6796 (1992). 53E. Papajak and D. G. Truhlar, J. Chem. Theory Comput. 6, 597 (2010). 54E. Papajak and D. G. Truhlar, J. Chem. Theory Comput. 7, 10 (2011). 55D. W. Schwenke, J. Chem. Phys. 122, 014107 (2005). 56See supplementary material at http://dx.doi.org/10.1063/1.4921256 for tables with mean absolute deviations, root-mean-square errors, standard deviations, and maximum errors. 57F. Weigend, A. Köhn, and C. Hättig, J. Chem. Phys. 116, 3175 (2002). 58C. Hättig, Phys. Chem. Chem. Phys. 7, 59 (2005). 59F. Weigend, J. Comput. Chem. 29, 167 (2008). 60K. E. Yousaf and K. A. Peterson, Chem. Phys. Lett. 476, 303 (2009). 61S. Boys and F. Bernardi, Mol. Phys. 19, 553 (1970). 62A. Halkier, W. Klopper, T. Helgaker, P. Jørgensen, and P. R. Taylor, J. Chem. Phys. 111, 9157 (1999). 63D. P. Tew and W. Klopper, J. Chem. Phys. 125, 094302 (2006).


Journal of Chemical Physics is copyrighted by AIP Publishing LLC (AIP). Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. For more information, see http://publishing.aip.org/authors/rights-and-permissions.

Density fitting for three-electron integrals in explicitly correlated electronic structure theory.

Accurate thermochemistry from explicitly correlated distinguishable cluster approximation.

Accurate potential energy curve of the LiH+ molecule calculated with explicitly correlated Gaussian functions.

Prediction of (1)P Rydberg energy levels of beryllium based on calculations with explicitly correlated Gaussians.

Nuclear-electronic orbital reduced explicitly correlated Hartree-Fock approach: Restricted basis sets and open-shell systems.

Linear-scaling explicitly correlated treatment of solids: periodic local MP2-F12 method.

Communication: stochastic evaluation of explicitly correlated second-order many-body perturbation energy.

Embedded correlated wavefunction schemes: theory and applications.

Quantum treatment of protons with the reduced explicitly correlated Hartree-Fock approach.

High angular momentum states of lithium atom, studied with symmetry-projected explicitly correlated Gaussian lobe functions.

Critical Examination of Explicitly Time-Dependent Density Functional Theory for Coherent Control of Dipole Switching.

A comparison of two types of explicitly correlated Gaussian functions for non-Born-Oppenheimer molecular calculations using a model potential.

Characterization of Zn(q+)-imidazole (q = 0, 1, 2) organometallic complexes: DFT methods vs. standard and explicitly correlated post-Hartree-Fock methods.

Direct non-Born-Oppenheimer variational calculations of all bound vibrational states corresponding to the first rotational excitation of D2 performed with explicitly correlated all-particle Gaussian functions.

An algorithm for nonrelativistic quantum-mechanical finite-nuclear-mass variational calculations of nitrogen atom in L = 0, M = 0 states using all-electrons explicitly correlated Gaussian basis functions.

Nonadiabatic time-dependent spin-density functional theory for strongly correlated systems.

Generalized multi-Gaussian correlated Schell-model beam: from theory to experiment.

Accurate and systematically improvable density functional theory embedding for correlated wavefunctions.

Linear response theory for the density matrix renormalization group: efficient algorithms for strongly correlated excited states.

We believe it is better not to explicitly highlight spirituality.

Correlation consistent basis sets for explicitly correlated wavefunctions: pseudopotential-based basis sets for the post-d main group elements Ga-Rn.

An algorithm for non-Born-Oppenheimer quantum mechanical variational calculations of N = 1 rotationally excited states of diatomic molecules using all-particle explicitly correlated Gaussian functions.

Structural Learning in a Visuomotor Adaptation Task Is Explicitly Accessible.

Explicitly targeting explicit sepsis: fair metrics for quality assessment*.