Efficient and accurate treatment of weak pairs in local CCSD(T) calculations. II. Beyond the ring approximation Martin Schütz, Oliver Masur, and Denis Usvyat Citation: The Journal of Chemical Physics 140, 244107 (2014); doi: 10.1063/1.4884156 View online: http://dx.doi.org/10.1063/1.4884156 View Table of Contents: http://scitation.aip.org/content/aip/journal/jcp/140/24?ver=pdfcov Published by the AIP Publishing Articles you may be interested in Orbital-optimized MP2.5 and its analytic gradients: Approaching CCSD(T) quality for noncovalent interactions J. Chem. Phys. 141, 204105 (2014); 10.1063/1.4902226 Efficient and accurate treatment of weak pairs in local CCSD(T) calculations J. Chem. Phys. 139, 164116 (2013); 10.1063/1.4826534 An efficient linear-scaling CCSD(T) method based on local natural orbitals J. Chem. Phys. 139, 094105 (2013); 10.1063/1.4819401 Efficient and accurate approximations to the local coupled cluster singles doubles method using a truncated pair natural orbital basis J. Chem. Phys. 131, 064103 (2009); 10.1063/1.3173827 Efficient and accurate local approximations to coupled-electron pair approaches: An attempt to revive the pair natural orbital method J. Chem. Phys. 130, 114108 (2009); 10.1063/1.3086717
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: 137.149.200.5 On: Sat, 29 Nov 2014 12:58:47
THE JOURNAL OF CHEMICAL PHYSICS 140, 244107 (2014)
Efficient and accurate treatment of weak pairs in local CCSD(T) calculations. II. Beyond the ring approximation Martin Schütz,a) Oliver Masur, and Denis Usvyat Institut für Physikalische und Theoretische Chemie, Universität Regensburg, Regensburg D-93040, Germany
(Received 24 March 2014; accepted 6 June 2014; published online 25 June 2014) In order to arrive at linear scaling of the computational cost with molecular size, local coupled cluster methods discriminate pairs of local molecular orbitals according to the spatial separation R of the latter. Only strong pairs are treated at the full coupled cluster level, whereas for weak pairs a lower level of theory (usually Møller-Plesset perturbation theory of second order, MP2) is used. Yet an MP2 treatment of weak pairs is inadequate in certain situations (for example, for describing π -stacking), which calls for an improved but still inexpensive method for dealing with the weak pairs. In a previous contribution, we proposed as a substituent for MP2 the LrCCD3 method, which is based on ring coupled cluster doubles (ring-CCD) and includes all third-order diagrams with energy contributions decaying not quicker than R−6 . In the present work, we explore a still more accurate method, which is based on the same principles. It turned out to be essential to abandon the restriction to ring-CCD, i.e., to include further CCD diagrams beyond the ring approximation. The occurring intermediates turn out to be formally very similar to LMP2 density matrices, such that an efficient evaluation of these non-ring CCD diagrams is possible. Furthermore, a computationally cheap a posteriori estimate for the fourth-order singles contribution to the weak pair energy, which also exhibits a decay behavior of R−6 , is introduced. The resulting method, denoted as LCCD[S]-R−6 , indeed provides a substantial improvement in accuracy over the previous LrCCD3 method at a relatively modest additional computational cost. © 2014 AIP Publishing LLC. [http://dx.doi.org/10.1063/1.4884156] I. INTRODUCTION
Coupled cluster singles and doubles [CCSD] theory with a perturbative evaluation of connected triples excitations [CCSD(T)] reaches chemical accuracy (provided that large enough basis sets are employed) for chemical systems that are dominated by the (Hartree-Fock) reference determinant, i.e., systems that do not exhibit multireference character. Hence CCSD(T) is frequently referred to as the “gold standard” of numerical quantum chemistry. The computational complexity of canonical CCSD(T) is still rather high; the scaling with respect to molecular size N is O(N 6 ) for the CCSD iterations, and O(N 7 ) for the a posteriori perturbative triples estimate. In order to expand the application range of CCSD(T) theory to big molecular systems, local CCSD(T) methods have been developed during the past 15 years.1–10 The term “local” here implies the formulation of CCSD(T) in terms of spatially localized occupied and virtual molecular orbitals (MOs), rather than the commonly used delocalized canonical orbitals. Occupied (in the reference determinant) local orbitals (LMOs) are usually obtained from a localization procedure such as Pipek-Mezey11 or Boys12 applied to the occupied canonical orbitals. As virtual orbitals, projected atomic orbitals (PAOs),13–18 orbital specific virtuals (OSVs),5, 19–21 or pair natural orbitals (PNOs)9, 10, 22–29 have been employed, either exclusively, or in combination. In local coupled cluster theory the number of relevant determinants and related amplitudes is drastically reduced such a) [email protected]
0021-9606/2014/140(24)/244107/9/$30.00
that eventually O(N ) scaling of the computational complexity is reached.1, 2 This reduction is based on (i) truncations of the virtual space, and (ii) pair approximations.1, 2, 4 While the truncation of the virtual space is motivated by the exponential decay of the related amplitudes with respect to the distance R between the local orbitals, the pair approximations are only motivated by an R−6 decay and thus a more critical approximation. Pair approximations assign the individual LMO pairs to the distinct pair classes strong, weak, and very distant based on the distance R between these two LMOs. The full local coupled cluster (LCC) treatment is reserved exclusively for strong pairs, while weak pairs are treated at a lower level, typically with local second-order Møller-Plesset perturbation theory (LMP2), and very distant pairs are entirely neglected. The weak pair class is subdivided into the class of the close pairs, and the rest. Close pairs were originally introduced in the context of local triples;1, 3 they contribute together with the strong pairs to the triples residual, whereas the remaining weak pairs are omitted. Later, close pairs were also optionally included in the LCC residuals (cf. Eq. (B27) in Ref. 4), which turned out to be essential for explicitly correlated LCC methods,30 and also for a proper treatment of intermolecular interactions.31, 32 In Paper I,31 we investigated the accuracy of the pair approximation for a set of intermolecular complexes and clusters and proposed an inexpensive alternative to the local MP2 (LMP2) treatment of weak pairs, denoted as local rCCD3 (LrCCD3). LrCCD3 is based on ring-CCD33, 34 plus additional exchange terms (such that the antisymmetry of
140, 244107-1
© 2014 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: 137.149.200.5 On: Sat, 29 Nov 2014 12:58:47
244107-2
Schütz, Masur, and Usvyat
the amplitudes in the spin-orbit formalism remains intact). Furthermore, only diagrams up to third-order with energy contributions not decaying quicker than R−6 are maintained. Note that here and in the following we mean by the “order of a diagram” its leading contribution to the correlation energy (within the Møller-Plesset partitioning) [due to the iterative procedure of the individual methods discussed here each diagram generates apart from its leading energy contribution further contributions of higher (essentially infinite) order]. Similarly, the “decay rate of a diagram” corresponds to the slowest decaying (with respect to the inter-orbital separation) energy component, originating from that diagram. The LrCCD3 treatment of weak pairs turned out to be clearly superior to the LMP2 treatment at an only moderate increase in the computational cost.31 On the other hand, a full LrCCD treatment of weak pairs (including additional fourthorder diagrams) was performing worse than LrCCD3. In Ref. 31, we speculated that this disappointing behavior of LrCCD is due to the neglect of further fourth-order CCD diagrams absent in rCCD, which also decay only as R−6 (diagrams D13–D16 in Fig. 5 of Ref. 31). In this contribution, we demonstrate that this is indeed the case. The resulting method which includes these diagrams on top of LrCCD, denoted for further reference as LCCD-R−6 , is superior to LrCCD3 in the sense that it provides similar accuracy as the previous LrCCD3 method, but in contrast to the latter, not due to fortuitous cancellation of errors. Moreover, the diagrams D13–D16 of Ref. 31 can efficiently be evaluated by (i) pre-contraction to a matrix reminiscent of a LMP2 density matrix, and (ii) accumulation of the latter to the Fock matrix in the (orbital invariant) LMP2 residual calculation. Yet in order to reach a substantially higher accuracy than that of the previous LrCCD3 method it is furthermore necessary to include the fourth-order singles contribution to the close pair energy. This is done via a computationally inexpensive a posteriori perturbative estimate, which is reusing the already calculated intermediates available from above mentioned close pair feedback to the LCC singles residual. For the final method, we introduce the acronym LCCD[S]-R−6 . The corresponding LCCSD|LCCD[S]-R−6 hybrid method indeed provides interaction energies in excellent agreement with a full (and much more expensive) LCCSD calculation. Finally, we note that the classification of the LCCSD(T)|LCCD[S]-R−6 diagrams according to their decay rate and the leading order of their contributions to the energy conforms to the following pair hierarchy: (i) strong pairs are treated at the full LCCSD level, (ii) close pairs involve all the diagrams of 4th order decaying as R−6 ; they contribute to the (4th order) triples1 and via 3rd and 4th order diagrams to the LCC residuals,4 (iii) weak/distant pairs involve all the diagrams of 3rd or 2nd order decaying as R−6 , and (v) very distant pairs are disregarded. In the local triples treatment, the leading contribution is of the 4th order and decays as R−6 , which dictates that the strong and close pairs have to be included. This is in accord with the usual L(T) treatment.1, 3
J. Chem. Phys. 140, 244107 (2014)
In Sec. II A below, we just provide the algebraic expressions of the fourth-order R−6 decaying LCCD diagrams absent in LrCCD (diagrams D13–D16 in Fig. 5 of Ref. 31). Furthermore, it turns out to be necessary to include also an estimate for the fourth-order R−6 decaying singles estimate absent in LCCD. We calculate such a fourth-order singles correction a posteriori, as discussed in Sec. II B. As usually, indices i, j, k, l denote in the following localized occupied molecular orbitals (LMOs), while indices r, s, t, u denote localized virtual orbitals, i.e., PAOs. For other choices of localized virtuals such as OSVs or PNOs, the formulae are very similar and straightforwardly obtained from those presented here for PAOs. Furthermore, the integrals are written in Mulliken’s notation and implicit summation of repeated indices (Einstein convention) is employed. Boldfaced letters refer to the matrix form of the respective quantities. A. Additional CCD diagrams
The contributions of the diagrams D13 and D15 (cf. Fig. 5 of Ref. 31) to the doubles residual yield together jl
ij Rrs = −PSrr Trik s Ss s (kt|lu)T˜tu
= −PSrr Trik s Dkj Ss s = −S Dki Tkj + Tik Dkj S,
(1)
where P = 1 + (ij )(rs)(r s ) is a permutation operator, S is the PAO overlap matrix, T˜ ij = 2Tij − Tj i , and Dki = (kt|lu)T˜tuil
(2)
is an object formally reminiscent of the internal LMP2 density matrix,35, 36 but with the amplitudes replaced by the respective electron repulsion integrals. Analogously, the expression for the diagram D14 and D16 contributions to the doubles residual reads as ij ij = −PSrr Tr t (kt|lu)T˜uslk Ss s Rrs ij
ij
= −Dtr Tts Ss s − Srr Tr t Dts = −D† Tij S − STij D,
(3)
Dtr = (kt|lu)T˜urlk Sr r ,
(4)
where
which in turn formally resembles the external LMP2 density matrix. Equations (1) and (3) suggest to accumulate the density alike objects after construction to the internal and external Fock matrix when computing the diagrams corresponding to the MP2 residual contribution [cf. Eq. (5) in Ref. 16]. This way, the diagrams D14–D16 can be included in the method efficiently. A similar approach for evaluating such terms was also used in Refs. 21 and 37. The LrCCD method presented in Ref. 31, augmented with the additional diagrams D13– D16 beyond the ring approximation, defines the LCCD-R−6 method. B. Perturbative fourth-order singles correction
II. THEORY
Theory and residual equations of the LrCCD3 and LrCCD methods are given in Ref. 31 and not reiterated here.
As shown below, the LCCSD|LCCD-R−6 method by itself does not provide a substantial improvement in accuracy over LCCSD|LrCCD3 (though the good performance of the
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: 137.149.200.5 On: Sat, 29 Nov 2014 12:58:47
244107-3
Schütz, Masur, and Usvyat
A
B
A
J. Chem. Phys. 140, 244107 (2014)
B
FIG. 1. The singles contributions to the close-pair doubles residual, which contribute in the fourth-order to the correlation energy (not included in the LCCD-R−6 method). The sum of these two diagrams decays as R−6 .32
latter turns out to be due to a fortuitous cancellation of errors). The reason is the neglect of the fourth-order singles contributions to the close-pair doubles residual, which are displayed in Fig. 1. The sum of these two diagrams have a decay rate of R−6 ,32 and consequently they should be included in our treatment of weak pairs. Unlike all the other diagrams, one of the diagrams of Fig. 1 involves integrals with three external indices, which are computationally quite expensive. To keep the computational cost as low as possible, we calculate an a posteriori perturbational estimate of the energy contribution originating from these diagrams, rather than including them in the evaluation of the close-pair doubles residual. This perturbational estimate is defined similarly as the [T] correction,38 in the sense that an MP4-like energy correction is computed from the converged LCCSD|LCCD-R−6 singles and close-pair doubles. Consequently, we denote this correction as “[S].” Furthermore, when writing the corresponding diagrams upside down, as shown in Fig. 2, it becomes evident that the required intermediates, namely, the contraction of the close pair doubles with the corresponding three external integrals, is already available (vide infra). The algebraic expression for the diagrams D18–D21 is E 4s = 2tri T˜stik (rs|kt) − Srr T˜rkl s (sl|ki) = 2tri 2K(T ik )rk − K(T ki )rk − Srr T˜rkl s (sl|ki) , (5) ij
with K(T ij )rk = Tst (rs|kt). According to our pair hierarchy mentioned at the end of Sec. I, we restrict the doubles amplitudes in Eq. (5) to close pairs. Since the K(Tij )rk operators [cf. Eq. (31) in Ref. 4], constructed with close pair amplitudes, are already available from the calculation of the close pair coupling, i.e., from the feedback of the close pair doubles to the LCCSD singles residual [cf. Eq. (B27) in Ref. 4], the additional effort for calculating E 4s is very small. The LCCD-R−6 method together with the addition of E 4s defines the LCCD[S]-R−6 method. We point out that the non-iterative [S]-treatment only refers to the energy component that originates from the sin-
A
B
(D18)
A
B
(D19)
A (D20)
B
A
B
(D21)
FIG. 2. The diagrams of the a posteriori perturbational estimate of the fourth-order singles contribution to the weak pair correlation energy of the LCCD[S]-R−6 method.
gles’ contribution to the close-pair doubles, while the singles themselves are still fully iterated in the framework of the LCCSD|LCCD-R−6 scheme. This energy component reflects mainly the effect of induction on dispersion (for that reason it is most important for H-bonded systems, cf. Sec. III). For exotic cases where this effect becomes dominant the non-iterative [S]-treatment may fail, but for all closedshell systems we could think of we expect it to work properly. III. TEST CALCULATIONS
In the following, five hybrid methods with different approximate weak-pair treatment are compared: (i) the standard LCCSD(T)|LMP2 method;2, 4 (ii) LCCSD(T)|LrCCD3, with the diagrams D1, D3, and D6 of Ref. 31 included in the weak pair residual, in addition to the LMP2 diagrams; (iii) LCCSD(T)|LrCCD, with the diagrams D2, D4, D5, D7–D9 added to the weak pair residual on top of those of LrCCD3; (iv) LCCSD(T)|LCCD-R−6 , with the additional diagrams D13–D16 on top of LrCCD, and (v) LCCSD(T)|LCCD[S]-R−6 , with an additional a posteriori perturbative estimate of the fourth-order singles contribution to the close pair energy (D18–D21 of Fig. 2 on top of LCCSD(T)|LCCD-R−6 ). All five methods include the feed-back of the closepair amplitudes into the LCCSD singles and doubles residual equations (close pair coupling, cf. Ref. 4), as the uncoupled treatment leads to very large errors.31 Most of the diagrams of the weak pair residual (D2–D5, D7–D9, D15, and D16) are automatically obtained without any computational effort just by taking the contravariant amij plitudes T˜rs (defined above) instead of the covariant ampliij tudes Trs . Furthermore, according to Fig. 3 of Ref. 31 and the related discussion, additional savings without significant loss in accuracy is achieved by restricting the pairs of the occupied indices in the two-electron integrals of diagrams D6–D9 to strong pairs only. In order to assess the individual accuracies of the five hybrid methods a test set of 25 intermolecular complexes was compiled from the S6639 and S2240 benchmark data sets, augmented by some further systems of the JSCH-2005 set40 neither included in the S66, nor the S22 sets. The selected systems represent quite a balanced mixture between large and small hydrogen-bonded (9), dispersion dominated (11), and mixed (5) intermolecular complexes; a detailed compilation of the complexes is given in the supplementary material.41 In all these calculations, the intra- and intermolecular pairs were assigned to the strong and close pair classes, respectively (the close and weak pair lists thus are identical). Assigning all intramolecular pairs to the strong pair class is usually not done in routine calculations (due to higher computational cost), yet for the present purpose it is more convenient, since it trivially allows for the partitioning of the interaction energy into intra- and intermolecular components. Boughton-Pulay domains with a completeness criterion of 0.985, determined at large intermolecular separation and kept fixed in the actual dimer calculations42 were employed throughout, and all intermolecular interaction energies are counterpoise corrected. For the triples correction, the L(T0)
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: 137.149.200.5 On: Sat, 29 Nov 2014 12:58:47
244107-4
Schütz, Masur, and Usvyat
J. Chem. Phys. 140, 244107 (2014)
FIG. 3. Deviations (in kcal/mol) of the correlation contributions to the interaction energies Ec from the respective full LCCSD values, for LCCSD|LMP2 (blue), LCCSD|LrCCD3 (red), LCCSD|LrCCD (yellow), LCCSD|LCCD-R−6 (green), and LCCSD|LCCD[S]-R−6 (brown).
approximation1 was utilized. Furthermore, all calculations were performed in the aVTZ basis set (cc-pVTZ43 on hydrogen atoms, aug-cc-pVTZ44 on all other atoms). As auxiliary basis for density fitting, the MP2FIT basis set of Weigend et al. corresponding to the aVTZ AO basis set45 was employed. In the following, we use the “full” LCCSD or LCCSD(T) methods, i.e., with all pairs treated at the LCCSD level and
entering the perturbative triples correction as the corresponding references. As shown in Tables VI and VII of the supplementary material,41 the deviations of full LCCSD and LCCSD(T) interaction energies from the respective canonical values do not exceed 0.05 kcal/mol, if the MP2 domain error correction4, 5, 20 is applied. Figure 3 displays, for the five above mentioned hybrid methods without L(T) correction, the respective
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: 137.149.200.5 On: Sat, 29 Nov 2014 12:58:47
244107-5
Schütz, Masur, and Usvyat
J. Chem. Phys. 140, 244107 (2014)
FIG. 4. Deviations (in kcal/mol) of the correlation contributions to the interaction energies Ec from the respective full LCCSD(T) values, for LCCSD(T)|LMP2 (blue), LCCSD(T)|LrCCD3 (red), LCCSD(T)|LrCCD (yellow), LCCSD(T)|LCCD-R−6 (green), and LCCSD(T)|LCCD[S]-R−6 (brown).
deviations of the correlation contribution to the interaction energy Ec from the full LCCSD result. Furthermore, Fig. 4 shows an analogous plot for the hybrid methods with L(T) correction1 relative to the full LCCSD(T) result. As is evident from Fig. 3, the inclusion of the diagrams D13–D16 has a considerable effect on the LCCSD|LrCCD interaction energies, generally increasing them (making them less neg-
ative). For example, for the π -stacked benzene dimer Ec increases by 0.70 kcal/mol, for guanine-cytosine dimer even by 1.02 and 1.50 kcal/mol, for the Watson Crick and the stacked configuration, respectively. This means generally a substantial improvement compared to LCCSD|LrCCD. Furthermore, compared to LCCSD|LMP2 and LCCSD|LrCCD3, the LCCSD|LCCD(R−6 ) model is indeed more accurate
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: 137.149.200.5 On: Sat, 29 Nov 2014 12:58:47
244107-6
Schütz, Masur, and Usvyat
J. Chem. Phys. 140, 244107 (2014)
TABLE I. RMS and maximum deviations in Ec (total correlation contribution to interaction energy) and ET (triples component thereof) of the five different hybrid methods from the full (all pairs strong) LCCSD and LCCSD(T) references, respectively. All values in kcal/mol.
X|LMP2 X|LrCCD3 X|LrCCD
TABLE II. Comparison of pure pair theories LCCD|LrCCD3 and LCCD|LCCD(R−6 ) to full LCCD. Given are the correlation contribution to the interaction energy Ec , and its intra- and intermolecular components δEintra and δEinter 31, 42, 46 (all values in kcal/mol).
X|LCCD- X|LCCD[S]R−6 R−6
δ RMS (Ec ) δ MAX (Ec )
0.637 1.626
0.293 0.673
X = LCCSD 0.550 0.240 1.212 0.556
0.033 0.082
δ RMS (Ec ) δ MAX (Ec ) δ RMS (ET ) δ MAX (ET )
0.816 2.130 0.181 0.504
0.384 0.900 0.092 0.227
X = LCCSD(T) 0.712 0.299 1.599 0.693 0.162 0.059 0.387 0.138
0.067 0.160 0.059 0.138
for van der Waals dominated and mixed cases. Yet for H-bonded cases, where LCCSD|LMP2 and LCCSD|LrCCD3 are both already very accurate, LCCSD|LCCD(R−6 ) is generally somewhat worse, systematically underestimating the Ec of full LCCSD. Yet even in this case, LCCSD|LCCD(R−6 ) appears to provide better relative energies of different conformers, as the example of the two Guanine-Cytosine dimers shows [deviation from full LCCSD result: 1.57 kcal/mol, 0.22 kcal/mol, and 0.07 kcal/mol for LCCSD|LMP2, LCCSD|LrCCD3, and LCCSD|LCCD(R−6 )], respectively]. The systematic underestimation of the interaction energies of H-bonded systems is remarkably remedied by our proposed perturbative fourth-order singles correction, as is evident from Figs. 3 and 4: the LCCSD|LCCD[S]-R−6 interaction energies are in excellent agreement with those of full LCCSD, and the same also applies for LCCSD(T)|LCCD[S]R−6 vs. full LCCSD(T), though the deviations are slightly larger for the latter case (we note in passing that the doubles amplitudes entering the triples correction are identical to those obtained with LCCSD|LCCD-R−6 ). Table I compiles RMS and maximum deviations of Ec for the individual hybrid methods. Notably, the maximum deviation decreases from 1.63 over 0.67 to 0.08 kcal/mol on going from LCCSD|LMP2 over LCCSD|LrCCD3 to LCCSD|LCCD[S]-R−6 , and from 2.13 over 0.90 to 0.16 kcal/mol when the perturbative triples are added. The RMS deviation amounts to 0.82, 0.38, and 0.07 kcal/mol, respectively, for the latter three methods. It is clear from these numbers that LrCCD3 constitutes a substantial improvement over LMP2, and LCCD[S]-R−6 , in turn, a substantial improvement over LrCCD3. From Table I it is also clear that inclusion of the fourthorder singles correction to the close pair energy is essential to achieve a considerably higher accuracy than LCCD|LrCCD3; as already mentioned above, the LCCSD|LCCD-R−6 alone performs even worse than LCCD|LrCCD3 for some of the complexes, in particular hydrogen bonded systems. This, in turn, also implies that the LCCD|LrCCD3 method must benefit from some fortuitous error cancellation between the missing singles correction and the fourth-order pair diagrams included at the level of LCCSD|LCCD-R−6 . In order to investigate this further, we performed for a critical subset of the test molecules, additional LCCD|LrCCD3 and LCCD|LCCD-
LCCD
LCCD |LrCCD3
LCCD |LCCD-R−6
δEinter δEintra Ec
− 2.228 1.519 − 0.709
Water-water − 2.322 1.529 − 0.793
− 2.204 1.515 − 0.688
δEinter δEintra Ec
− 8.762 6.927 − 1.835
Uracil-AcNH2 − 9.267 6.996 − 2.271
− 8.741 6.924 − 1.817
δEinter δEintra Ec
Guanine-cytosine (Watson Crick) − 13.319 − 14.118 10.086 10.208 − 3.232 − 3.910
− 13.285 10.088 − 3.197
δEinter δEintra Ec
Pyridine-pyridine (π -π ) − 7.772 − 8.350 2.263 2.373 − 5.509 − 5.976
− 7.800 2.269 − 5.530
δEinter δEintra Ec
Guanine-cytosine (stacked) − 17.821 − 18.802 8.954 9.129 − 8.867 − 9.673
− 17.765 8.967 − 8.798
R−6 calculations for comparison with full LCCD, i.e., comparing pure pair theories. Table II collects the results of these calculations. It turns that LCCSD|LCCD-R−6 indeed is in excellent agreement with full LCCD, with a maximum deviation of 0.07 kcal/mol (for the stacked Guanine-Cytosine dimer). LCCD|LrCCD3, on the other hand, has considerably larger deviations, 0.7 and 0.8 kcal/mol for guanine-cytosine in Watson Crick and stacked configuration, respectively. These larger deviations originate primarily from the intermolecular component δEinter of Ec . The new LCCSD|LCCD[S]-R−6 hybrid method was also tested for four geometries of the (CS2 )3 trimer introduced in Paper I,31 i.e., in parallel and tilted geometry (tilting angle α = 28.2◦ ), each with two different intermolecular C– C distances of R = 3.811 and R = 3.65 Å, respectively. This system turned out to be particularly sensitive with respect to the weak pair approximation.31 Table III compiles interaction energies Eint , the correlation contributions Ec thereof, as well as the three-body correlation contributions E3c . Evidently, the LCCSD|LCCD[S]-R−6 Ec values are in good agreement with those of a full LCCSD calculation, deviating at most by 0.17 kcal/mol, which compares to deviations of 0.72, 1.62, and even 2.61 kcal/mol for LCCSD|LrCCD3, LCCSD|LrCCD, and LCCSD|LMP2, respectively. When contemplating the individual three-body correlation contributions E3c to Ec a striking improvement relative to LCCSD|LrCCD3 manifests, reducing the relative error from 46%–58% to less than 3%. The interaction energies of the four (CS2 )3 geometries were also calculated at the LCCSD(T) level; Table IV compiles the resulting Eint values, along with their triples
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: 137.149.200.5 On: Sat, 29 Nov 2014 12:58:47
244107-7
Schütz, Masur, and Usvyat
J. Chem. Phys. 140, 244107 (2014)
TABLE III. Intermolecular interaction energies Eint , correlation contributions Ec to Eint , and three-body correlation contributions E3c to Ec (all in [kcal/mol]) of equilateral (CS2 )3 for different intermolecular distances R (in Å). Two different tilting angles α were applied: α = 28.2◦ , R = 3.811 Å is the experimental geometry reported in Ref. 47, and α = 0◦ specifies a parallel arrangement of the three CS2 molecules. R = 3.65 Å corresponds to the LMP2 minimum.
TABLE IV. Intermolecular interaction energies Eint and correlation contributions of the local triples correction ET (all in [kcal/mol]) of (CS2 )3 for two different intermolecular distances R (in [Å]) and two different tilting angles α. R = 3.811 Å, α = 28.2◦ is the experimental geometry reported in Ref. 47, R = 3.65 Å corresponds to the LMP2 minimum, and α = 0◦ to a parallel arrangement of the three CS2 molecules. In parentheses the relative stabilities of the four geometries, i.e., the differences of the individual Eint to that of the experimental geometry are given.
LCCSD R
Eint Eint Ec Ec E3c E3c Eint Eint Ec Ec E3c E3c
Full
|LMP2
|LrCCD3
|LrCCD
3.81 3.65 3.81 3.65 3.81 3.65
Equilateral (CS2 )3 , α = 28.2◦ − 3.018 − 4.913 − 3.604 − 4.247 − 1.899 − 4.266 − 2.579 − 3.383 − 8.297 − 10.192 − 8.883 − 9.526 − 10.371 − 12.739 − 11.051 − 11.855 0.172 0.052 0.272 0.258 0.254 0.097 0.386 0.369
3.81 3.65 3.81 3.65 3.81 3.65
Equilateral (CS2 )3 , α = 0◦ − 2.404 − 0.946 − 1.662 − 0.269 1.619 0.722 − 11.338 − 9.880 − 10.596 − 14.332 − 12.443 − 13.341 0.110 0.460 0.442 0.203 0.660 0.638
− 0.328 2.340 − 9.261 − 11.723 0.303 0.452
LCCSD(T)
|LCCD[S]-R−6 R − 3.107 − 2.011 − 8.386 − 10.484 0.168 0.248 − 0.468 2.169 − 9.402 − 11.893 0.309 0.461
|LMP2
Full
|LrCCD3
|LCCD[S]-R−6
Eint Eint ET ET
3.81 3.65 3.81 3.65
Equilateral (CS2 )3 , α = 28.2◦ −4.948 (0.00) −7.424 (0.00) −5.696 (0.00) −5.031 (0.00) −4.267 (0.68) −7.362 (0.06) −5.138 (0.56) −4.369 (0.66) −1.930 −2.511 −2.093 −1.925 −2.368 −3.096 −2.559 −2.358
Eint Eint ET ET
3.81 3.65 3.81 3.65
Equilateral (CS2 )3 , α = 0◦ −2.534 (2.41) −5.317 (2.11) −3.330 (2.37) −2.656 (2.38) −0.394 (4.55) −3.889 (3.54) −1.326 (4.37) −0.541 (4.49) −2.206 −2.913 −2.384 −2.188 −2.734 −3.620 −2.946 −2.710
agreement with the reference values. The same applies to the interaction energies Eint themselves, and, most notably, to the relative stabilities of the four geometries, given in parenthesis in Table IV. The max errors in the latter are 1.01, 0.18, and 0.06 kcal/mol for the three methods, respectively. Note that the LCCSD|LrCCD values compiled in Table III differ somewhat from those of Table 4 in Ref. 31, which is due to the fixing of a bug in a fourth-order diagram, which we became aware of in the course of the present work. Furthermore, there are also slight differences in the LCCSD|LrCCD3 and LCCSD(T)|LrCCD3 results compared
correction contributions ET of full LCCSD(T) and of the hybrid methods LCCSD(T)|LMP2, LCCSD(T)|LrCCD3, and LCCSD(T)|LCCD[S]-R−6 . ET clearly is important for this system, it amounts to 2–3 kcal/mol and thus is of similar size as Eint itself. It is evident from Table IV that there is a clear improvement in the ET values on going from LCCSD(T)|LMP2 over LCCSD(T)|LrCCD3 to LCCSD(T)|LCCD[S]-R−6 , the latter again being in excellent
TABLE V. Computational details for the calculations on the T-shaped indole-benzene complex (total number of pairs: 703, number of intermolecular pairs: 330). In the hybrid calculations, the intramolecular pairs may be strong, close, or weak (number of strong/close/weak pairs: s = 181, c = 413, w = 109). Timings for the transformation of the 3-external and 4-external integrals, for evaluating the individual groups of diagrams, and the time for converging the (hybrid) coupled cluster calculations are given in seconds. The file sizes of the 3-external and 4-external integral distributions are given in Gbytes. The interaction energies Eint and their inter- and intramolecular components δEinter , δEintra are also provided (in kcal/mol). The calculations were performed on seven AMD Opteron 6180 SE cores. LCCSD| Part LMP2 LrCCD3 LrCCD LCCD-R−6 [S] Sum parts Average iteration time Number of iterations Integral transformation (3-ext) Integral transformation (4-ext) 3-ext file size 4-ext file size δEinter δEintra Eint a
Diagramsa
Full
LMP2 92
LrCCD3
75 86 347
LCCD-R−6 76 83 348 4
LCCD[S]-R−6
3404 10 5265 9281 90.646 263.293
92 1382 9 3830 3501 64.271 61.860
181 1603 9 3809 3382 64.271 61.860
508 1908 9 3968 3419 64.271 61.860
511 1875 10 3772 3384 64.271 61.860
76 83 348 4 0.3 511 1875 10 3772 3384 64.271 61.860
− 7.440 3.236 − 3.909
− 8.223 3.504 − 4.431
− 7.837 3.314 − 4.234
− 8.160 3.393 − 4.478
− 7.372 3.204 − 3.871
− 7.403 3.204 − 3.902
D1,D3,D6 D2,D4,D5,D7-9 D13-16 D18-21
92 89
LrCCD
Additional diagrams beyond LMP2 according to Figs. 1, 2, and 5 of Ref. 31, and Fig. 2.
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: 137.149.200.5 On: Sat, 29 Nov 2014 12:58:47
244107-8
Schütz, Masur, and Usvyat
to the earlier publication, which is due to the fact that the feedback of the close pair amplitudes to the strong pair residual via three-external integrals was now calculated with the actual (iterated) close pair amplitudes, rather than the LMP2 close pair amplitudes, as done previously. In Table V, we provide computational details for the calculations on the T-shaped indole-benzene complex, one of the bigger dimers of our test set. Here, unlike the previous calculations, we allow for close and weak intramolecular pairs (as is usually done in routine calculations, since this reduces substantially the number of strong pairs; for the present case 181 instead of 373), yet still all intermolecular pairs are specified as close. An analogous table for the calculations with all intramolecular pairs being treated as strong, is provided in the supplementary material.41 Evidently, the interaction energies, as well as their intra- and intermolecular components,31, 42, 46 which are also provided in Table V, are hardly affected. Notably, the LCCSD|LCCD[S]-R−6 values are in excellent agreement with those of the full LCCSD calculation. The timings for the iterations increase by 16% on going from LCCSD|LMP2 to LCCSD|LrCCD3, and by further 17% on going from LCCSD|LrCCD3 to LCCSD|LCCD[S]R−6 . As can be seen, the evaluation of the non-ring CCD doubles diagrams and the a posteriori perturbative estimate of the fourth-order singles cost virtually nothing, implying that an LCCSD|LCCD[S]-R−6 calculation comes about at the same price as an LCCSD|LrCCD calculation. Moreover, as is also evident from Table V, no additional cost in integral transformation time or integral file size is required when replacing LCCSD|LMP2 by LCCSD|LCCD[S]-R−6 . IV. CONCLUSIONS
In local coupled cluster methods, the individual pairs of occupied localized MOs are assigned to distinct pair classes (strong, weak, and very distant), e.g., according to their spatial separation, or based on an inexpensive estimate of the corresponding pair energies. In such a hierarchical pair approximation only strong pairs are treated at the full coupled cluster level, whereas weak pairs are dealt with at a lower level of theory (typically MP2), and very distant pairs are entirely neglected. Yet an MP2 treatment of weak pairs may not provide sufficient accuracy. For example, if the system comprises π stacked aromatic rings, or in other situations where van der Waals contributions play a significant role, MP2 provides just qualitative rather than quantitative accuracy. Hence the LMP2 part of the LCCSD(T)|LMP2 hybrid method may compromise the high accuracy of the LCCSD(T) part. Recently,31 we proposed the more accurate, and at the same time still inexpensive LrCCD3 method for the treatment of the weak pairs. It includes all the ring-CCD diagrams plus additional exchange terms up to third-order with energy contributions that do not decay quicker than R−6 . In particular, the ladder diagrams, which altogether lead to an energy contribution decaying as quickly as R−9 ,31, 32 are discarded. In the present paper, we presented an even more accurate fourth-order method along the same lines. It turned out to be essential to abandon the ring-CCD approximation and to include further (non-ring) CCD diagrams, and simultane-
J. Chem. Phys. 140, 244107 (2014)
ously also to add a correction for the fourth-order singles contribution to the correlation energy, which all decay also as R−6 . These additional diagrams can still be evaluated efficiently. The intermediates of the non-ring CCD doubles diagrams formally resemble LMP2 density matrices, which in turn appear as “Fock matrices” in a LMP2 residual alike expression. The fourth-order singles contribution, on the other hand, can be evaluated via an a posteriori perturbational estimate, which reuses quantities already available for the close pair feedback to the strong pair residual. In particular, no further 3-external and 4-external electron repulsion integrals are required in the proposed LCCSD(T)|LCCD[S]-R−6 method beyond those necessary for LCCSD(T)|LMP2. The timings of the calculations performed for the T-shaped indole-benzene complex show an additional cost of 17% in the (hybrid) coupled cluster iteration times relative to LCCSD(T)|LrCCD3, and of 36% relative to LCCSD(T)|LMP2. On the other hand, the LCCD[S]-R−6 approach establishes a major improvement in accuracy over LrCCD3 and LMP2; the RMS deviations of LCCSD(T)|LMP2, LCCSD(T)|LrCCD3, and LCCSD(T)|LCCD[S]-R−6 from the full LCCSD(T) reference calculation (all pairs strong) declines from 0.82 over 0.38 to 0.07 kcal/mol for our test set of dimers. It thus appears that replacing LMP2 (or LrCCD3) by LCCD[S]-R−6 in local coupled cluster calculations with pair approximations is highly recommendable. ACKNOWLEDGMENTS
The authors would like to thank Dr. Daniel Kats for useful discussions. Financial support of the Deutsche Forschungsgemeinschaft (DFG), Grant Nos. US-103/1-1 and SCHU 1456/8-1, is gratefully acknowledged. 1 M.
Schütz, J. Chem. Phys. 113, 9986 (2000). Schütz and H.-J. Werner, J. Chem. Phys. 114, 661 (2001). 3 M. Schütz, J. Chem. Phys. 116, 8772 (2002). 4 H.-J. Werner and M. Schütz, J. Chem. Phys. 135, 144116 (2011). 5 M. Schütz, J. Yang, G. K. L. Chan, F. R. Manby, and H.-J. Werner, J. Chem. Phys. 138, 054109 (2013). 6 P. E. Maslen, A. Dutoi, M. S. Lee, Y. H. Shao, and M. Head-Gordon, Mol. Phys. 103, 425 (2005). 7 W. Li, P. Piecuch, J. R. Gour, and S. Li, J. Chem. Phys. 131, 114109 (2009). 8 Z. Rolik and M. Kállay, J. Chem. Phys. 135, 104111 (2011). 9 Z. Rolik, L. Szegedy, I. Ladjanszki, B. Ladoczki, and M. Kallay, J. Chem. Phys. 139, 094105 (2013). 10 C. Riplinger, B. Sandhoefer, A. Hansen, and F. Neese, J. Chem. Phys. 139, 134101 (2013). 11 J. Pipek and P. G. Mezey, J. Chem. Phys. 90, 4916 (1989). 12 S. F. Boys, “Localized orbitals and localized adjustment functions,” in Quantum Theory of Atoms, Molecules, and the Solid State, edited by P. O. Löwdin (Academic Press, New York, 1966), pp. 253–262. 13 P. Pulay, Chem. Phys. Lett. 100, 151 (1983). 14 S. Saebø and P. Pulay, Chem. Phys. Lett. 113, 13 (1985). 15 P. Pulay and S. Saebø, Theor. Chim. Acta 69, 357 (1986). 16 M. Schütz, G. Hetzer, and H.-J. Werner, J. Chem. Phys. 111, 5691 (1999). 17 D. Kats, D. Usvyat, and M. Schütz, Phys. Chem. Chem. Phys. 10, 3430 (2008). 18 C. Pisani, M. Schütz, S. Casassa, D. Usvyat, L. Maschio, M. Lorenz, and A. Erba, Phys. Chem. Chem. Phys. 14, 7615 (2012). 19 J. Yang, Y. Kurashige, F. R. Manby, and G. K. L. Chan, J. Chem. Phys. 134, 044123 (2011). 20 J. Yang, G. K.-L. Chan, F. R. Manby, M. Schütz, and H.-J. Werner, J. Chem. Phys. 136, 144105 (2012). 21 D. Kats and F. R. Manby, J. Chem. Phys. 138, 144101 (2013). 2 M.
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: 137.149.200.5 On: Sat, 29 Nov 2014 12:58:47
244107-9
Schütz, Masur, and Usvyat
J. Chem. Phys. 140, 244107 (2014)
22 W.
36 D.
23 W.
37 J.
Meyer, Int. J. Quantum Chem. 5(S5), 341 (1971). Meyer, J. Chem. Phys. 58, 1017 (1973). 24 R. Ahlrichs, F. Driessler, H. Lischka, V. Staemmler, and W. Kutzelnigg, J. Chem. Phys. 62, 1235 (1975). 25 F. Neese, F. Wennmohs, and A. Hansen, J. Chem. Phys. 130, 114108 (2009). 26 F. Neese, A. Hansen, and D. G. Liakos, J. Chem. Phys. 131, 064103 (2009). 27 A. Hansen, D. G. Liakos, and F. Neese, J. Chem. Phys. 135, 214102 (2011). 28 C. Hättig, D. P. Tew, and B. Helmich, J. Chem. Phys. 136, 204105 (2012). 29 C. Krause and H.-J. Werner, Phys. Chem. Chem. Phys. 14, 7591 (2012). 30 T. B. Adler and H.-J. Werner, J. Chem. Phys. 135, 144117 (2011). 31 O. Masur, D. Usvyat, and M. Schütz, J. Chem. Phys. 139, 164116 (2013). 32 D. Usvyat, “On the physics of intermolecular interactions in local correlation approaches,” (unpublished). 33 G. E. Scuseria, T. M. Henderson, and D. C. Sorensen, J. Chem. Phys. 129, 231101 (2008). 34 J. Toulouse, W. Zhu, A. Savin, G. Jansen, and J. Angyan, J. Chem. Phys. 135, 084119 (2011). 35 M. Schütz, H.-J. Werner, R. Lindh, and F. R. Manby, J. Chem. Phys. 121, 737 (2004).
Usvyat and M. Schütz, J. Phys.: Conf. Ser. 117, 012027 (2008). J. Shepherd, T. M. Henderson, and G. E. Scuseria, J. Chem. Phys. 140, 124102 (2014). 38 M. Urban, J. Noga, S. J. Cole, and R. J. Bartlett, J. Chem. Phys. 83, 4041 (1985). 39 J. Rezᡠˇ c, K. E. Riley, and P. Hobza, J. Chem. Theory Comput. 7, 2427 (2011). 40 P. Jureˇ ˇ cka, J. Šponer, J. Cerný, and P. Hobza, Phys. Chem. Chem. Phys. 8, 1985 (2006). 41 See supplementary material at http://dx.doi.org/10.1063/1.4884156 for a detailed compilation of our test set of complexes, and for detailed numerical results of the test calculations. 42 M. Schütz, G. Rauhut, and H.-J. Werner, J. Phys. Chem. A 102, 5997 (1998). 43 T. H. Dunning, Jr., J. Chem. Phys. 90, 1007 (1989). 44 R. A. Kendall, T. H. Dunning, and R. J. Harrison, J. Chem. Phys. 96, 6796 (1992). 45 F. Weigend, M. Häser, H. Patzelt, and R. Ahlrichs, Chem. Phys. Lett. 294, 143 (1998). 46 D. Usvyat, K. Sadeghian, L. Maschio, and M. Schütz, Phys. Rev. B 86, 045412 (2012). 47 M. Rezaei, J. Norooz Oliaee, N. Moazzen-Ahmadi, and A. R. W. McKellar, Phys. Chem. Chem. Phys. 13, 12635 (2011).
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: 137.149.200.5 On: Sat, 29 Nov 2014 12:58:47
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: 137.149.200.5 On: Sat, 29 Nov 2014 12:58:47
THE JOURNAL OF CHEMICAL PHYSICS 140, 244107 (2014)
Efficient and accurate treatment of weak pairs in local CCSD(T) calculations. II. Beyond the ring approximation Martin Schütz,a) Oliver Masur, and Denis Usvyat Institut für Physikalische und Theoretische Chemie, Universität Regensburg, Regensburg D-93040, Germany
(Received 24 March 2014; accepted 6 June 2014; published online 25 June 2014) In order to arrive at linear scaling of the computational cost with molecular size, local coupled cluster methods discriminate pairs of local molecular orbitals according to the spatial separation R of the latter. Only strong pairs are treated at the full coupled cluster level, whereas for weak pairs a lower level of theory (usually Møller-Plesset perturbation theory of second order, MP2) is used. Yet an MP2 treatment of weak pairs is inadequate in certain situations (for example, for describing π -stacking), which calls for an improved but still inexpensive method for dealing with the weak pairs. In a previous contribution, we proposed as a substituent for MP2 the LrCCD3 method, which is based on ring coupled cluster doubles (ring-CCD) and includes all third-order diagrams with energy contributions decaying not quicker than R−6 . In the present work, we explore a still more accurate method, which is based on the same principles. It turned out to be essential to abandon the restriction to ring-CCD, i.e., to include further CCD diagrams beyond the ring approximation. The occurring intermediates turn out to be formally very similar to LMP2 density matrices, such that an efficient evaluation of these non-ring CCD diagrams is possible. Furthermore, a computationally cheap a posteriori estimate for the fourth-order singles contribution to the weak pair energy, which also exhibits a decay behavior of R−6 , is introduced. The resulting method, denoted as LCCD[S]-R−6 , indeed provides a substantial improvement in accuracy over the previous LrCCD3 method at a relatively modest additional computational cost. © 2014 AIP Publishing LLC. [http://dx.doi.org/10.1063/1.4884156] I. INTRODUCTION
Coupled cluster singles and doubles [CCSD] theory with a perturbative evaluation of connected triples excitations [CCSD(T)] reaches chemical accuracy (provided that large enough basis sets are employed) for chemical systems that are dominated by the (Hartree-Fock) reference determinant, i.e., systems that do not exhibit multireference character. Hence CCSD(T) is frequently referred to as the “gold standard” of numerical quantum chemistry. The computational complexity of canonical CCSD(T) is still rather high; the scaling with respect to molecular size N is O(N 6 ) for the CCSD iterations, and O(N 7 ) for the a posteriori perturbative triples estimate. In order to expand the application range of CCSD(T) theory to big molecular systems, local CCSD(T) methods have been developed during the past 15 years.1–10 The term “local” here implies the formulation of CCSD(T) in terms of spatially localized occupied and virtual molecular orbitals (MOs), rather than the commonly used delocalized canonical orbitals. Occupied (in the reference determinant) local orbitals (LMOs) are usually obtained from a localization procedure such as Pipek-Mezey11 or Boys12 applied to the occupied canonical orbitals. As virtual orbitals, projected atomic orbitals (PAOs),13–18 orbital specific virtuals (OSVs),5, 19–21 or pair natural orbitals (PNOs)9, 10, 22–29 have been employed, either exclusively, or in combination. In local coupled cluster theory the number of relevant determinants and related amplitudes is drastically reduced such a) [email protected]
0021-9606/2014/140(24)/244107/9/$30.00
that eventually O(N ) scaling of the computational complexity is reached.1, 2 This reduction is based on (i) truncations of the virtual space, and (ii) pair approximations.1, 2, 4 While the truncation of the virtual space is motivated by the exponential decay of the related amplitudes with respect to the distance R between the local orbitals, the pair approximations are only motivated by an R−6 decay and thus a more critical approximation. Pair approximations assign the individual LMO pairs to the distinct pair classes strong, weak, and very distant based on the distance R between these two LMOs. The full local coupled cluster (LCC) treatment is reserved exclusively for strong pairs, while weak pairs are treated at a lower level, typically with local second-order Møller-Plesset perturbation theory (LMP2), and very distant pairs are entirely neglected. The weak pair class is subdivided into the class of the close pairs, and the rest. Close pairs were originally introduced in the context of local triples;1, 3 they contribute together with the strong pairs to the triples residual, whereas the remaining weak pairs are omitted. Later, close pairs were also optionally included in the LCC residuals (cf. Eq. (B27) in Ref. 4), which turned out to be essential for explicitly correlated LCC methods,30 and also for a proper treatment of intermolecular interactions.31, 32 In Paper I,31 we investigated the accuracy of the pair approximation for a set of intermolecular complexes and clusters and proposed an inexpensive alternative to the local MP2 (LMP2) treatment of weak pairs, denoted as local rCCD3 (LrCCD3). LrCCD3 is based on ring-CCD33, 34 plus additional exchange terms (such that the antisymmetry of
140, 244107-1
© 2014 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: 137.149.200.5 On: Sat, 29 Nov 2014 12:58:47
244107-2
Schütz, Masur, and Usvyat
the amplitudes in the spin-orbit formalism remains intact). Furthermore, only diagrams up to third-order with energy contributions not decaying quicker than R−6 are maintained. Note that here and in the following we mean by the “order of a diagram” its leading contribution to the correlation energy (within the Møller-Plesset partitioning) [due to the iterative procedure of the individual methods discussed here each diagram generates apart from its leading energy contribution further contributions of higher (essentially infinite) order]. Similarly, the “decay rate of a diagram” corresponds to the slowest decaying (with respect to the inter-orbital separation) energy component, originating from that diagram. The LrCCD3 treatment of weak pairs turned out to be clearly superior to the LMP2 treatment at an only moderate increase in the computational cost.31 On the other hand, a full LrCCD treatment of weak pairs (including additional fourthorder diagrams) was performing worse than LrCCD3. In Ref. 31, we speculated that this disappointing behavior of LrCCD is due to the neglect of further fourth-order CCD diagrams absent in rCCD, which also decay only as R−6 (diagrams D13–D16 in Fig. 5 of Ref. 31). In this contribution, we demonstrate that this is indeed the case. The resulting method which includes these diagrams on top of LrCCD, denoted for further reference as LCCD-R−6 , is superior to LrCCD3 in the sense that it provides similar accuracy as the previous LrCCD3 method, but in contrast to the latter, not due to fortuitous cancellation of errors. Moreover, the diagrams D13–D16 of Ref. 31 can efficiently be evaluated by (i) pre-contraction to a matrix reminiscent of a LMP2 density matrix, and (ii) accumulation of the latter to the Fock matrix in the (orbital invariant) LMP2 residual calculation. Yet in order to reach a substantially higher accuracy than that of the previous LrCCD3 method it is furthermore necessary to include the fourth-order singles contribution to the close pair energy. This is done via a computationally inexpensive a posteriori perturbative estimate, which is reusing the already calculated intermediates available from above mentioned close pair feedback to the LCC singles residual. For the final method, we introduce the acronym LCCD[S]-R−6 . The corresponding LCCSD|LCCD[S]-R−6 hybrid method indeed provides interaction energies in excellent agreement with a full (and much more expensive) LCCSD calculation. Finally, we note that the classification of the LCCSD(T)|LCCD[S]-R−6 diagrams according to their decay rate and the leading order of their contributions to the energy conforms to the following pair hierarchy: (i) strong pairs are treated at the full LCCSD level, (ii) close pairs involve all the diagrams of 4th order decaying as R−6 ; they contribute to the (4th order) triples1 and via 3rd and 4th order diagrams to the LCC residuals,4 (iii) weak/distant pairs involve all the diagrams of 3rd or 2nd order decaying as R−6 , and (v) very distant pairs are disregarded. In the local triples treatment, the leading contribution is of the 4th order and decays as R−6 , which dictates that the strong and close pairs have to be included. This is in accord with the usual L(T) treatment.1, 3
J. Chem. Phys. 140, 244107 (2014)
In Sec. II A below, we just provide the algebraic expressions of the fourth-order R−6 decaying LCCD diagrams absent in LrCCD (diagrams D13–D16 in Fig. 5 of Ref. 31). Furthermore, it turns out to be necessary to include also an estimate for the fourth-order R−6 decaying singles estimate absent in LCCD. We calculate such a fourth-order singles correction a posteriori, as discussed in Sec. II B. As usually, indices i, j, k, l denote in the following localized occupied molecular orbitals (LMOs), while indices r, s, t, u denote localized virtual orbitals, i.e., PAOs. For other choices of localized virtuals such as OSVs or PNOs, the formulae are very similar and straightforwardly obtained from those presented here for PAOs. Furthermore, the integrals are written in Mulliken’s notation and implicit summation of repeated indices (Einstein convention) is employed. Boldfaced letters refer to the matrix form of the respective quantities. A. Additional CCD diagrams
The contributions of the diagrams D13 and D15 (cf. Fig. 5 of Ref. 31) to the doubles residual yield together jl
ij Rrs = −PSrr Trik s Ss s (kt|lu)T˜tu
= −PSrr Trik s Dkj Ss s = −S Dki Tkj + Tik Dkj S,
(1)
where P = 1 + (ij )(rs)(r s ) is a permutation operator, S is the PAO overlap matrix, T˜ ij = 2Tij − Tj i , and Dki = (kt|lu)T˜tuil
(2)
is an object formally reminiscent of the internal LMP2 density matrix,35, 36 but with the amplitudes replaced by the respective electron repulsion integrals. Analogously, the expression for the diagram D14 and D16 contributions to the doubles residual reads as ij ij = −PSrr Tr t (kt|lu)T˜uslk Ss s Rrs ij
ij
= −Dtr Tts Ss s − Srr Tr t Dts = −D† Tij S − STij D,
(3)
Dtr = (kt|lu)T˜urlk Sr r ,
(4)
where
which in turn formally resembles the external LMP2 density matrix. Equations (1) and (3) suggest to accumulate the density alike objects after construction to the internal and external Fock matrix when computing the diagrams corresponding to the MP2 residual contribution [cf. Eq. (5) in Ref. 16]. This way, the diagrams D14–D16 can be included in the method efficiently. A similar approach for evaluating such terms was also used in Refs. 21 and 37. The LrCCD method presented in Ref. 31, augmented with the additional diagrams D13– D16 beyond the ring approximation, defines the LCCD-R−6 method. B. Perturbative fourth-order singles correction
II. THEORY
Theory and residual equations of the LrCCD3 and LrCCD methods are given in Ref. 31 and not reiterated here.
As shown below, the LCCSD|LCCD-R−6 method by itself does not provide a substantial improvement in accuracy over LCCSD|LrCCD3 (though the good performance of the
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: 137.149.200.5 On: Sat, 29 Nov 2014 12:58:47
244107-3
Schütz, Masur, and Usvyat
A
B
A
J. Chem. Phys. 140, 244107 (2014)
B
FIG. 1. The singles contributions to the close-pair doubles residual, which contribute in the fourth-order to the correlation energy (not included in the LCCD-R−6 method). The sum of these two diagrams decays as R−6 .32
latter turns out to be due to a fortuitous cancellation of errors). The reason is the neglect of the fourth-order singles contributions to the close-pair doubles residual, which are displayed in Fig. 1. The sum of these two diagrams have a decay rate of R−6 ,32 and consequently they should be included in our treatment of weak pairs. Unlike all the other diagrams, one of the diagrams of Fig. 1 involves integrals with three external indices, which are computationally quite expensive. To keep the computational cost as low as possible, we calculate an a posteriori perturbational estimate of the energy contribution originating from these diagrams, rather than including them in the evaluation of the close-pair doubles residual. This perturbational estimate is defined similarly as the [T] correction,38 in the sense that an MP4-like energy correction is computed from the converged LCCSD|LCCD-R−6 singles and close-pair doubles. Consequently, we denote this correction as “[S].” Furthermore, when writing the corresponding diagrams upside down, as shown in Fig. 2, it becomes evident that the required intermediates, namely, the contraction of the close pair doubles with the corresponding three external integrals, is already available (vide infra). The algebraic expression for the diagrams D18–D21 is E 4s = 2tri T˜stik (rs|kt) − Srr T˜rkl s (sl|ki) = 2tri 2K(T ik )rk − K(T ki )rk − Srr T˜rkl s (sl|ki) , (5) ij
with K(T ij )rk = Tst (rs|kt). According to our pair hierarchy mentioned at the end of Sec. I, we restrict the doubles amplitudes in Eq. (5) to close pairs. Since the K(Tij )rk operators [cf. Eq. (31) in Ref. 4], constructed with close pair amplitudes, are already available from the calculation of the close pair coupling, i.e., from the feedback of the close pair doubles to the LCCSD singles residual [cf. Eq. (B27) in Ref. 4], the additional effort for calculating E 4s is very small. The LCCD-R−6 method together with the addition of E 4s defines the LCCD[S]-R−6 method. We point out that the non-iterative [S]-treatment only refers to the energy component that originates from the sin-
A
B
(D18)
A
B
(D19)
A (D20)
B
A
B
(D21)
FIG. 2. The diagrams of the a posteriori perturbational estimate of the fourth-order singles contribution to the weak pair correlation energy of the LCCD[S]-R−6 method.
gles’ contribution to the close-pair doubles, while the singles themselves are still fully iterated in the framework of the LCCSD|LCCD-R−6 scheme. This energy component reflects mainly the effect of induction on dispersion (for that reason it is most important for H-bonded systems, cf. Sec. III). For exotic cases where this effect becomes dominant the non-iterative [S]-treatment may fail, but for all closedshell systems we could think of we expect it to work properly. III. TEST CALCULATIONS
In the following, five hybrid methods with different approximate weak-pair treatment are compared: (i) the standard LCCSD(T)|LMP2 method;2, 4 (ii) LCCSD(T)|LrCCD3, with the diagrams D1, D3, and D6 of Ref. 31 included in the weak pair residual, in addition to the LMP2 diagrams; (iii) LCCSD(T)|LrCCD, with the diagrams D2, D4, D5, D7–D9 added to the weak pair residual on top of those of LrCCD3; (iv) LCCSD(T)|LCCD-R−6 , with the additional diagrams D13–D16 on top of LrCCD, and (v) LCCSD(T)|LCCD[S]-R−6 , with an additional a posteriori perturbative estimate of the fourth-order singles contribution to the close pair energy (D18–D21 of Fig. 2 on top of LCCSD(T)|LCCD-R−6 ). All five methods include the feed-back of the closepair amplitudes into the LCCSD singles and doubles residual equations (close pair coupling, cf. Ref. 4), as the uncoupled treatment leads to very large errors.31 Most of the diagrams of the weak pair residual (D2–D5, D7–D9, D15, and D16) are automatically obtained without any computational effort just by taking the contravariant amij plitudes T˜rs (defined above) instead of the covariant ampliij tudes Trs . Furthermore, according to Fig. 3 of Ref. 31 and the related discussion, additional savings without significant loss in accuracy is achieved by restricting the pairs of the occupied indices in the two-electron integrals of diagrams D6–D9 to strong pairs only. In order to assess the individual accuracies of the five hybrid methods a test set of 25 intermolecular complexes was compiled from the S6639 and S2240 benchmark data sets, augmented by some further systems of the JSCH-2005 set40 neither included in the S66, nor the S22 sets. The selected systems represent quite a balanced mixture between large and small hydrogen-bonded (9), dispersion dominated (11), and mixed (5) intermolecular complexes; a detailed compilation of the complexes is given in the supplementary material.41 In all these calculations, the intra- and intermolecular pairs were assigned to the strong and close pair classes, respectively (the close and weak pair lists thus are identical). Assigning all intramolecular pairs to the strong pair class is usually not done in routine calculations (due to higher computational cost), yet for the present purpose it is more convenient, since it trivially allows for the partitioning of the interaction energy into intra- and intermolecular components. Boughton-Pulay domains with a completeness criterion of 0.985, determined at large intermolecular separation and kept fixed in the actual dimer calculations42 were employed throughout, and all intermolecular interaction energies are counterpoise corrected. For the triples correction, the L(T0)
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: 137.149.200.5 On: Sat, 29 Nov 2014 12:58:47
244107-4
Schütz, Masur, and Usvyat
J. Chem. Phys. 140, 244107 (2014)
FIG. 3. Deviations (in kcal/mol) of the correlation contributions to the interaction energies Ec from the respective full LCCSD values, for LCCSD|LMP2 (blue), LCCSD|LrCCD3 (red), LCCSD|LrCCD (yellow), LCCSD|LCCD-R−6 (green), and LCCSD|LCCD[S]-R−6 (brown).
approximation1 was utilized. Furthermore, all calculations were performed in the aVTZ basis set (cc-pVTZ43 on hydrogen atoms, aug-cc-pVTZ44 on all other atoms). As auxiliary basis for density fitting, the MP2FIT basis set of Weigend et al. corresponding to the aVTZ AO basis set45 was employed. In the following, we use the “full” LCCSD or LCCSD(T) methods, i.e., with all pairs treated at the LCCSD level and
entering the perturbative triples correction as the corresponding references. As shown in Tables VI and VII of the supplementary material,41 the deviations of full LCCSD and LCCSD(T) interaction energies from the respective canonical values do not exceed 0.05 kcal/mol, if the MP2 domain error correction4, 5, 20 is applied. Figure 3 displays, for the five above mentioned hybrid methods without L(T) correction, the respective
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: 137.149.200.5 On: Sat, 29 Nov 2014 12:58:47
244107-5
Schütz, Masur, and Usvyat
J. Chem. Phys. 140, 244107 (2014)
FIG. 4. Deviations (in kcal/mol) of the correlation contributions to the interaction energies Ec from the respective full LCCSD(T) values, for LCCSD(T)|LMP2 (blue), LCCSD(T)|LrCCD3 (red), LCCSD(T)|LrCCD (yellow), LCCSD(T)|LCCD-R−6 (green), and LCCSD(T)|LCCD[S]-R−6 (brown).
deviations of the correlation contribution to the interaction energy Ec from the full LCCSD result. Furthermore, Fig. 4 shows an analogous plot for the hybrid methods with L(T) correction1 relative to the full LCCSD(T) result. As is evident from Fig. 3, the inclusion of the diagrams D13–D16 has a considerable effect on the LCCSD|LrCCD interaction energies, generally increasing them (making them less neg-
ative). For example, for the π -stacked benzene dimer Ec increases by 0.70 kcal/mol, for guanine-cytosine dimer even by 1.02 and 1.50 kcal/mol, for the Watson Crick and the stacked configuration, respectively. This means generally a substantial improvement compared to LCCSD|LrCCD. Furthermore, compared to LCCSD|LMP2 and LCCSD|LrCCD3, the LCCSD|LCCD(R−6 ) model is indeed more accurate
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: 137.149.200.5 On: Sat, 29 Nov 2014 12:58:47
244107-6
Schütz, Masur, and Usvyat
J. Chem. Phys. 140, 244107 (2014)
TABLE I. RMS and maximum deviations in Ec (total correlation contribution to interaction energy) and ET (triples component thereof) of the five different hybrid methods from the full (all pairs strong) LCCSD and LCCSD(T) references, respectively. All values in kcal/mol.
X|LMP2 X|LrCCD3 X|LrCCD
TABLE II. Comparison of pure pair theories LCCD|LrCCD3 and LCCD|LCCD(R−6 ) to full LCCD. Given are the correlation contribution to the interaction energy Ec , and its intra- and intermolecular components δEintra and δEinter 31, 42, 46 (all values in kcal/mol).
X|LCCD- X|LCCD[S]R−6 R−6
δ RMS (Ec ) δ MAX (Ec )
0.637 1.626
0.293 0.673
X = LCCSD 0.550 0.240 1.212 0.556
0.033 0.082
δ RMS (Ec ) δ MAX (Ec ) δ RMS (ET ) δ MAX (ET )
0.816 2.130 0.181 0.504
0.384 0.900 0.092 0.227
X = LCCSD(T) 0.712 0.299 1.599 0.693 0.162 0.059 0.387 0.138
0.067 0.160 0.059 0.138
for van der Waals dominated and mixed cases. Yet for H-bonded cases, where LCCSD|LMP2 and LCCSD|LrCCD3 are both already very accurate, LCCSD|LCCD(R−6 ) is generally somewhat worse, systematically underestimating the Ec of full LCCSD. Yet even in this case, LCCSD|LCCD(R−6 ) appears to provide better relative energies of different conformers, as the example of the two Guanine-Cytosine dimers shows [deviation from full LCCSD result: 1.57 kcal/mol, 0.22 kcal/mol, and 0.07 kcal/mol for LCCSD|LMP2, LCCSD|LrCCD3, and LCCSD|LCCD(R−6 )], respectively]. The systematic underestimation of the interaction energies of H-bonded systems is remarkably remedied by our proposed perturbative fourth-order singles correction, as is evident from Figs. 3 and 4: the LCCSD|LCCD[S]-R−6 interaction energies are in excellent agreement with those of full LCCSD, and the same also applies for LCCSD(T)|LCCD[S]R−6 vs. full LCCSD(T), though the deviations are slightly larger for the latter case (we note in passing that the doubles amplitudes entering the triples correction are identical to those obtained with LCCSD|LCCD-R−6 ). Table I compiles RMS and maximum deviations of Ec for the individual hybrid methods. Notably, the maximum deviation decreases from 1.63 over 0.67 to 0.08 kcal/mol on going from LCCSD|LMP2 over LCCSD|LrCCD3 to LCCSD|LCCD[S]-R−6 , and from 2.13 over 0.90 to 0.16 kcal/mol when the perturbative triples are added. The RMS deviation amounts to 0.82, 0.38, and 0.07 kcal/mol, respectively, for the latter three methods. It is clear from these numbers that LrCCD3 constitutes a substantial improvement over LMP2, and LCCD[S]-R−6 , in turn, a substantial improvement over LrCCD3. From Table I it is also clear that inclusion of the fourthorder singles correction to the close pair energy is essential to achieve a considerably higher accuracy than LCCD|LrCCD3; as already mentioned above, the LCCSD|LCCD-R−6 alone performs even worse than LCCD|LrCCD3 for some of the complexes, in particular hydrogen bonded systems. This, in turn, also implies that the LCCD|LrCCD3 method must benefit from some fortuitous error cancellation between the missing singles correction and the fourth-order pair diagrams included at the level of LCCSD|LCCD-R−6 . In order to investigate this further, we performed for a critical subset of the test molecules, additional LCCD|LrCCD3 and LCCD|LCCD-
LCCD
LCCD |LrCCD3
LCCD |LCCD-R−6
δEinter δEintra Ec
− 2.228 1.519 − 0.709
Water-water − 2.322 1.529 − 0.793
− 2.204 1.515 − 0.688
δEinter δEintra Ec
− 8.762 6.927 − 1.835
Uracil-AcNH2 − 9.267 6.996 − 2.271
− 8.741 6.924 − 1.817
δEinter δEintra Ec
Guanine-cytosine (Watson Crick) − 13.319 − 14.118 10.086 10.208 − 3.232 − 3.910
− 13.285 10.088 − 3.197
δEinter δEintra Ec
Pyridine-pyridine (π -π ) − 7.772 − 8.350 2.263 2.373 − 5.509 − 5.976
− 7.800 2.269 − 5.530
δEinter δEintra Ec
Guanine-cytosine (stacked) − 17.821 − 18.802 8.954 9.129 − 8.867 − 9.673
− 17.765 8.967 − 8.798
R−6 calculations for comparison with full LCCD, i.e., comparing pure pair theories. Table II collects the results of these calculations. It turns that LCCSD|LCCD-R−6 indeed is in excellent agreement with full LCCD, with a maximum deviation of 0.07 kcal/mol (for the stacked Guanine-Cytosine dimer). LCCD|LrCCD3, on the other hand, has considerably larger deviations, 0.7 and 0.8 kcal/mol for guanine-cytosine in Watson Crick and stacked configuration, respectively. These larger deviations originate primarily from the intermolecular component δEinter of Ec . The new LCCSD|LCCD[S]-R−6 hybrid method was also tested for four geometries of the (CS2 )3 trimer introduced in Paper I,31 i.e., in parallel and tilted geometry (tilting angle α = 28.2◦ ), each with two different intermolecular C– C distances of R = 3.811 and R = 3.65 Å, respectively. This system turned out to be particularly sensitive with respect to the weak pair approximation.31 Table III compiles interaction energies Eint , the correlation contributions Ec thereof, as well as the three-body correlation contributions E3c . Evidently, the LCCSD|LCCD[S]-R−6 Ec values are in good agreement with those of a full LCCSD calculation, deviating at most by 0.17 kcal/mol, which compares to deviations of 0.72, 1.62, and even 2.61 kcal/mol for LCCSD|LrCCD3, LCCSD|LrCCD, and LCCSD|LMP2, respectively. When contemplating the individual three-body correlation contributions E3c to Ec a striking improvement relative to LCCSD|LrCCD3 manifests, reducing the relative error from 46%–58% to less than 3%. The interaction energies of the four (CS2 )3 geometries were also calculated at the LCCSD(T) level; Table IV compiles the resulting Eint values, along with their triples
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: 137.149.200.5 On: Sat, 29 Nov 2014 12:58:47
244107-7
Schütz, Masur, and Usvyat
J. Chem. Phys. 140, 244107 (2014)
TABLE III. Intermolecular interaction energies Eint , correlation contributions Ec to Eint , and three-body correlation contributions E3c to Ec (all in [kcal/mol]) of equilateral (CS2 )3 for different intermolecular distances R (in Å). Two different tilting angles α were applied: α = 28.2◦ , R = 3.811 Å is the experimental geometry reported in Ref. 47, and α = 0◦ specifies a parallel arrangement of the three CS2 molecules. R = 3.65 Å corresponds to the LMP2 minimum.
TABLE IV. Intermolecular interaction energies Eint and correlation contributions of the local triples correction ET (all in [kcal/mol]) of (CS2 )3 for two different intermolecular distances R (in [Å]) and two different tilting angles α. R = 3.811 Å, α = 28.2◦ is the experimental geometry reported in Ref. 47, R = 3.65 Å corresponds to the LMP2 minimum, and α = 0◦ to a parallel arrangement of the three CS2 molecules. In parentheses the relative stabilities of the four geometries, i.e., the differences of the individual Eint to that of the experimental geometry are given.
LCCSD R
Eint Eint Ec Ec E3c E3c Eint Eint Ec Ec E3c E3c
Full
|LMP2
|LrCCD3
|LrCCD
3.81 3.65 3.81 3.65 3.81 3.65
Equilateral (CS2 )3 , α = 28.2◦ − 3.018 − 4.913 − 3.604 − 4.247 − 1.899 − 4.266 − 2.579 − 3.383 − 8.297 − 10.192 − 8.883 − 9.526 − 10.371 − 12.739 − 11.051 − 11.855 0.172 0.052 0.272 0.258 0.254 0.097 0.386 0.369
3.81 3.65 3.81 3.65 3.81 3.65
Equilateral (CS2 )3 , α = 0◦ − 2.404 − 0.946 − 1.662 − 0.269 1.619 0.722 − 11.338 − 9.880 − 10.596 − 14.332 − 12.443 − 13.341 0.110 0.460 0.442 0.203 0.660 0.638
− 0.328 2.340 − 9.261 − 11.723 0.303 0.452
LCCSD(T)
|LCCD[S]-R−6 R − 3.107 − 2.011 − 8.386 − 10.484 0.168 0.248 − 0.468 2.169 − 9.402 − 11.893 0.309 0.461
|LMP2
Full
|LrCCD3
|LCCD[S]-R−6
Eint Eint ET ET
3.81 3.65 3.81 3.65
Equilateral (CS2 )3 , α = 28.2◦ −4.948 (0.00) −7.424 (0.00) −5.696 (0.00) −5.031 (0.00) −4.267 (0.68) −7.362 (0.06) −5.138 (0.56) −4.369 (0.66) −1.930 −2.511 −2.093 −1.925 −2.368 −3.096 −2.559 −2.358
Eint Eint ET ET
3.81 3.65 3.81 3.65
Equilateral (CS2 )3 , α = 0◦ −2.534 (2.41) −5.317 (2.11) −3.330 (2.37) −2.656 (2.38) −0.394 (4.55) −3.889 (3.54) −1.326 (4.37) −0.541 (4.49) −2.206 −2.913 −2.384 −2.188 −2.734 −3.620 −2.946 −2.710
agreement with the reference values. The same applies to the interaction energies Eint themselves, and, most notably, to the relative stabilities of the four geometries, given in parenthesis in Table IV. The max errors in the latter are 1.01, 0.18, and 0.06 kcal/mol for the three methods, respectively. Note that the LCCSD|LrCCD values compiled in Table III differ somewhat from those of Table 4 in Ref. 31, which is due to the fixing of a bug in a fourth-order diagram, which we became aware of in the course of the present work. Furthermore, there are also slight differences in the LCCSD|LrCCD3 and LCCSD(T)|LrCCD3 results compared
correction contributions ET of full LCCSD(T) and of the hybrid methods LCCSD(T)|LMP2, LCCSD(T)|LrCCD3, and LCCSD(T)|LCCD[S]-R−6 . ET clearly is important for this system, it amounts to 2–3 kcal/mol and thus is of similar size as Eint itself. It is evident from Table IV that there is a clear improvement in the ET values on going from LCCSD(T)|LMP2 over LCCSD(T)|LrCCD3 to LCCSD(T)|LCCD[S]-R−6 , the latter again being in excellent
TABLE V. Computational details for the calculations on the T-shaped indole-benzene complex (total number of pairs: 703, number of intermolecular pairs: 330). In the hybrid calculations, the intramolecular pairs may be strong, close, or weak (number of strong/close/weak pairs: s = 181, c = 413, w = 109). Timings for the transformation of the 3-external and 4-external integrals, for evaluating the individual groups of diagrams, and the time for converging the (hybrid) coupled cluster calculations are given in seconds. The file sizes of the 3-external and 4-external integral distributions are given in Gbytes. The interaction energies Eint and their inter- and intramolecular components δEinter , δEintra are also provided (in kcal/mol). The calculations were performed on seven AMD Opteron 6180 SE cores. LCCSD| Part LMP2 LrCCD3 LrCCD LCCD-R−6 [S] Sum parts Average iteration time Number of iterations Integral transformation (3-ext) Integral transformation (4-ext) 3-ext file size 4-ext file size δEinter δEintra Eint a
Diagramsa
Full
LMP2 92
LrCCD3
75 86 347
LCCD-R−6 76 83 348 4
LCCD[S]-R−6
3404 10 5265 9281 90.646 263.293
92 1382 9 3830 3501 64.271 61.860
181 1603 9 3809 3382 64.271 61.860
508 1908 9 3968 3419 64.271 61.860
511 1875 10 3772 3384 64.271 61.860
76 83 348 4 0.3 511 1875 10 3772 3384 64.271 61.860
− 7.440 3.236 − 3.909
− 8.223 3.504 − 4.431
− 7.837 3.314 − 4.234
− 8.160 3.393 − 4.478
− 7.372 3.204 − 3.871
− 7.403 3.204 − 3.902
D1,D3,D6 D2,D4,D5,D7-9 D13-16 D18-21
92 89
LrCCD
Additional diagrams beyond LMP2 according to Figs. 1, 2, and 5 of Ref. 31, and Fig. 2.
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: 137.149.200.5 On: Sat, 29 Nov 2014 12:58:47
244107-8
Schütz, Masur, and Usvyat
to the earlier publication, which is due to the fact that the feedback of the close pair amplitudes to the strong pair residual via three-external integrals was now calculated with the actual (iterated) close pair amplitudes, rather than the LMP2 close pair amplitudes, as done previously. In Table V, we provide computational details for the calculations on the T-shaped indole-benzene complex, one of the bigger dimers of our test set. Here, unlike the previous calculations, we allow for close and weak intramolecular pairs (as is usually done in routine calculations, since this reduces substantially the number of strong pairs; for the present case 181 instead of 373), yet still all intermolecular pairs are specified as close. An analogous table for the calculations with all intramolecular pairs being treated as strong, is provided in the supplementary material.41 Evidently, the interaction energies, as well as their intra- and intermolecular components,31, 42, 46 which are also provided in Table V, are hardly affected. Notably, the LCCSD|LCCD[S]-R−6 values are in excellent agreement with those of the full LCCSD calculation. The timings for the iterations increase by 16% on going from LCCSD|LMP2 to LCCSD|LrCCD3, and by further 17% on going from LCCSD|LrCCD3 to LCCSD|LCCD[S]R−6 . As can be seen, the evaluation of the non-ring CCD doubles diagrams and the a posteriori perturbative estimate of the fourth-order singles cost virtually nothing, implying that an LCCSD|LCCD[S]-R−6 calculation comes about at the same price as an LCCSD|LrCCD calculation. Moreover, as is also evident from Table V, no additional cost in integral transformation time or integral file size is required when replacing LCCSD|LMP2 by LCCSD|LCCD[S]-R−6 . IV. CONCLUSIONS
In local coupled cluster methods, the individual pairs of occupied localized MOs are assigned to distinct pair classes (strong, weak, and very distant), e.g., according to their spatial separation, or based on an inexpensive estimate of the corresponding pair energies. In such a hierarchical pair approximation only strong pairs are treated at the full coupled cluster level, whereas weak pairs are dealt with at a lower level of theory (typically MP2), and very distant pairs are entirely neglected. Yet an MP2 treatment of weak pairs may not provide sufficient accuracy. For example, if the system comprises π stacked aromatic rings, or in other situations where van der Waals contributions play a significant role, MP2 provides just qualitative rather than quantitative accuracy. Hence the LMP2 part of the LCCSD(T)|LMP2 hybrid method may compromise the high accuracy of the LCCSD(T) part. Recently,31 we proposed the more accurate, and at the same time still inexpensive LrCCD3 method for the treatment of the weak pairs. It includes all the ring-CCD diagrams plus additional exchange terms up to third-order with energy contributions that do not decay quicker than R−6 . In particular, the ladder diagrams, which altogether lead to an energy contribution decaying as quickly as R−9 ,31, 32 are discarded. In the present paper, we presented an even more accurate fourth-order method along the same lines. It turned out to be essential to abandon the ring-CCD approximation and to include further (non-ring) CCD diagrams, and simultane-
J. Chem. Phys. 140, 244107 (2014)
ously also to add a correction for the fourth-order singles contribution to the correlation energy, which all decay also as R−6 . These additional diagrams can still be evaluated efficiently. The intermediates of the non-ring CCD doubles diagrams formally resemble LMP2 density matrices, which in turn appear as “Fock matrices” in a LMP2 residual alike expression. The fourth-order singles contribution, on the other hand, can be evaluated via an a posteriori perturbational estimate, which reuses quantities already available for the close pair feedback to the strong pair residual. In particular, no further 3-external and 4-external electron repulsion integrals are required in the proposed LCCSD(T)|LCCD[S]-R−6 method beyond those necessary for LCCSD(T)|LMP2. The timings of the calculations performed for the T-shaped indole-benzene complex show an additional cost of 17% in the (hybrid) coupled cluster iteration times relative to LCCSD(T)|LrCCD3, and of 36% relative to LCCSD(T)|LMP2. On the other hand, the LCCD[S]-R−6 approach establishes a major improvement in accuracy over LrCCD3 and LMP2; the RMS deviations of LCCSD(T)|LMP2, LCCSD(T)|LrCCD3, and LCCSD(T)|LCCD[S]-R−6 from the full LCCSD(T) reference calculation (all pairs strong) declines from 0.82 over 0.38 to 0.07 kcal/mol for our test set of dimers. It thus appears that replacing LMP2 (or LrCCD3) by LCCD[S]-R−6 in local coupled cluster calculations with pair approximations is highly recommendable. ACKNOWLEDGMENTS
The authors would like to thank Dr. Daniel Kats for useful discussions. Financial support of the Deutsche Forschungsgemeinschaft (DFG), Grant Nos. US-103/1-1 and SCHU 1456/8-1, is gratefully acknowledged. 1 M.
Schütz, J. Chem. Phys. 113, 9986 (2000). Schütz and H.-J. Werner, J. Chem. Phys. 114, 661 (2001). 3 M. Schütz, J. Chem. Phys. 116, 8772 (2002). 4 H.-J. Werner and M. Schütz, J. Chem. Phys. 135, 144116 (2011). 5 M. Schütz, J. Yang, G. K. L. Chan, F. R. Manby, and H.-J. Werner, J. Chem. Phys. 138, 054109 (2013). 6 P. E. Maslen, A. Dutoi, M. S. Lee, Y. H. Shao, and M. Head-Gordon, Mol. Phys. 103, 425 (2005). 7 W. Li, P. Piecuch, J. R. Gour, and S. Li, J. Chem. Phys. 131, 114109 (2009). 8 Z. Rolik and M. Kállay, J. Chem. Phys. 135, 104111 (2011). 9 Z. Rolik, L. Szegedy, I. Ladjanszki, B. Ladoczki, and M. Kallay, J. Chem. Phys. 139, 094105 (2013). 10 C. Riplinger, B. Sandhoefer, A. Hansen, and F. Neese, J. Chem. Phys. 139, 134101 (2013). 11 J. Pipek and P. G. Mezey, J. Chem. Phys. 90, 4916 (1989). 12 S. F. Boys, “Localized orbitals and localized adjustment functions,” in Quantum Theory of Atoms, Molecules, and the Solid State, edited by P. O. Löwdin (Academic Press, New York, 1966), pp. 253–262. 13 P. Pulay, Chem. Phys. Lett. 100, 151 (1983). 14 S. Saebø and P. Pulay, Chem. Phys. Lett. 113, 13 (1985). 15 P. Pulay and S. Saebø, Theor. Chim. Acta 69, 357 (1986). 16 M. Schütz, G. Hetzer, and H.-J. Werner, J. Chem. Phys. 111, 5691 (1999). 17 D. Kats, D. Usvyat, and M. Schütz, Phys. Chem. Chem. Phys. 10, 3430 (2008). 18 C. Pisani, M. Schütz, S. Casassa, D. Usvyat, L. Maschio, M. Lorenz, and A. Erba, Phys. Chem. Chem. Phys. 14, 7615 (2012). 19 J. Yang, Y. Kurashige, F. R. Manby, and G. K. L. Chan, J. Chem. Phys. 134, 044123 (2011). 20 J. Yang, G. K.-L. Chan, F. R. Manby, M. Schütz, and H.-J. Werner, J. Chem. Phys. 136, 144105 (2012). 21 D. Kats and F. R. Manby, J. Chem. Phys. 138, 144101 (2013). 2 M.
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: 137.149.200.5 On: Sat, 29 Nov 2014 12:58:47
244107-9
Schütz, Masur, and Usvyat
J. Chem. Phys. 140, 244107 (2014)
22 W.
36 D.
23 W.
37 J.
Meyer, Int. J. Quantum Chem. 5(S5), 341 (1971). Meyer, J. Chem. Phys. 58, 1017 (1973). 24 R. Ahlrichs, F. Driessler, H. Lischka, V. Staemmler, and W. Kutzelnigg, J. Chem. Phys. 62, 1235 (1975). 25 F. Neese, F. Wennmohs, and A. Hansen, J. Chem. Phys. 130, 114108 (2009). 26 F. Neese, A. Hansen, and D. G. Liakos, J. Chem. Phys. 131, 064103 (2009). 27 A. Hansen, D. G. Liakos, and F. Neese, J. Chem. Phys. 135, 214102 (2011). 28 C. Hättig, D. P. Tew, and B. Helmich, J. Chem. Phys. 136, 204105 (2012). 29 C. Krause and H.-J. Werner, Phys. Chem. Chem. Phys. 14, 7591 (2012). 30 T. B. Adler and H.-J. Werner, J. Chem. Phys. 135, 144117 (2011). 31 O. Masur, D. Usvyat, and M. Schütz, J. Chem. Phys. 139, 164116 (2013). 32 D. Usvyat, “On the physics of intermolecular interactions in local correlation approaches,” (unpublished). 33 G. E. Scuseria, T. M. Henderson, and D. C. Sorensen, J. Chem. Phys. 129, 231101 (2008). 34 J. Toulouse, W. Zhu, A. Savin, G. Jansen, and J. Angyan, J. Chem. Phys. 135, 084119 (2011). 35 M. Schütz, H.-J. Werner, R. Lindh, and F. R. Manby, J. Chem. Phys. 121, 737 (2004).
Usvyat and M. Schütz, J. Phys.: Conf. Ser. 117, 012027 (2008). J. Shepherd, T. M. Henderson, and G. E. Scuseria, J. Chem. Phys. 140, 124102 (2014). 38 M. Urban, J. Noga, S. J. Cole, and R. J. Bartlett, J. Chem. Phys. 83, 4041 (1985). 39 J. Rezᡠˇ c, K. E. Riley, and P. Hobza, J. Chem. Theory Comput. 7, 2427 (2011). 40 P. Jureˇ ˇ cka, J. Šponer, J. Cerný, and P. Hobza, Phys. Chem. Chem. Phys. 8, 1985 (2006). 41 See supplementary material at http://dx.doi.org/10.1063/1.4884156 for a detailed compilation of our test set of complexes, and for detailed numerical results of the test calculations. 42 M. Schütz, G. Rauhut, and H.-J. Werner, J. Phys. Chem. A 102, 5997 (1998). 43 T. H. Dunning, Jr., J. Chem. Phys. 90, 1007 (1989). 44 R. A. Kendall, T. H. Dunning, and R. J. Harrison, J. Chem. Phys. 96, 6796 (1992). 45 F. Weigend, M. Häser, H. Patzelt, and R. Ahlrichs, Chem. Phys. Lett. 294, 143 (1998). 46 D. Usvyat, K. Sadeghian, L. Maschio, and M. Schütz, Phys. Rev. B 86, 045412 (2012). 47 M. Rezaei, J. Norooz Oliaee, N. Moazzen-Ahmadi, and A. R. W. McKellar, Phys. Chem. Chem. Phys. 13, 12635 (2011).
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: 137.149.200.5 On: Sat, 29 Nov 2014 12:58:47
