Scalable Electron Correlation Methods. 4. Parallel ... - ACS Publications

Sep 12, 2017 - In accordance with previous findings, it is demonstrated that the F12 ... Accelerating the coupled-cluster singles and doubles method u...
0 downloads 0 Views 970KB Size
Subscriber access provided by UNIV OF CAMBRIDGE

Article

Scalable Electron Correlation Methods. 4. Parallel Explicitly Correlated Local Coupled Cluster with Pair Natural Orbitals (PNO-LCCSD-F12) Qianli Ma, Max Schwilk, Christoph Köppl, and Hans-Joachim Werner J. Chem. Theory Comput., Just Accepted Manuscript • DOI: 10.1021/acs.jctc.7b00799 • Publication Date (Web): 12 Sep 2017 Downloaded from http://pubs.acs.org on September 16, 2017

Just Accepted “Just Accepted” manuscripts have been peer-reviewed and accepted for publication. They are posted online prior to technical editing, formatting for publication and author proofing. The American Chemical Society provides “Just Accepted” as a free service to the research community to expedite the dissemination of scientific material as soon as possible after acceptance. “Just Accepted” manuscripts appear in full in PDF format accompanied by an HTML abstract. “Just Accepted” manuscripts have been fully peer reviewed, but should not be considered the official version of record. They are accessible to all readers and citable by the Digital Object Identifier (DOI®). “Just Accepted” is an optional service offered to authors. Therefore, the “Just Accepted” Web site may not include all articles that will be published in the journal. After a manuscript is technically edited and formatted, it will be removed from the “Just Accepted” Web site and published as an ASAP article. Note that technical editing may introduce minor changes to the manuscript text and/or graphics which could affect content, and all legal disclaimers and ethical guidelines that apply to the journal pertain. ACS cannot be held responsible for errors or consequences arising from the use of information contained in these “Just Accepted” manuscripts.

Journal of Chemical Theory and Computation is published by the American Chemical Society. 1155 Sixteenth Street N.W., Washington, DC 20036 Published by American Chemical Society. Copyright © American Chemical Society. However, no copyright claim is made to original U.S. Government works, or works produced by employees of any Commonwealth realm Crown government in the course of their duties.

Page 1 of 74

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Theory and Computation

Scalable Electron Correlation Methods. 4. Parallel Explicitly Correlated Local Coupled Cluster with Pair Natural Orbitals (PNO-LCCSD-F12) Qianli Ma, Max Schwilk, Christoph Köppl, and Hans-Joachim Werner∗ Institut für Theoretische Chemie, Universität Stuttgart, Pfaffenwaldring 55, D-70569 Stuttgart, Germany E-mail: [email protected]

Abstract We present an efficient explicitly correlated pair natural orbital local coupled cluster (PNOLCCSD-F12) method. The method is an extension of our previously reported PNO-LCCSD approach (Schwilk et al., J. Chem. Theory Comput. 2017, 13, 3650–3675). Near linear scaling with the size of molecule is achieved by using pair, domain, and projection approximations, local density fitting, local resolution of the identity, and by exploiting the sparsity of the local molecular orbitals as well as of the projected atomic orbitals. The effect of the various domain approximations is tested for a wide range of chemical reactions and intermolecular interactions. In accordance with previous findings it is demonstrated that the F12 terms significantly reduce the domain errors. The convergence of the reaction and interaction energies with respect to the parameters that determine the domain sizes and pair approximations is extensively tested. The results obtained with our default thresholds agree within a few tenths of a kcal mol−1 with the ones computed with very tight options. For cases where canonical calculations are still feasible, the relative energies of local and canonical calculations agree within similar error bounds. The PNO-LCCSD-F12 method needs only 25-40% more computer time than a

1 ACS Paragon Plus Environment

Journal of Chemical Theory and Computation

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

corresponding PNO-LCCSD calculation, while greatly improving the accuracy. Our program is well parallelized and capable of computing accurate correlation energies for molecules with more than 150 atoms, using augmented triple-ζ basis sets and more than 5000 basis functions. Using several nodes of a small computer cluster such calculations can be carried out within a few hours.

1 Introduction The steep scaling of computational resources with the molecule size limits the application of conventional wave function based quantum chemistry methods to relatively small molecules. The gold standard of quantum chemistry methods for single-reference systems, the coupled cluster method with single-, double-, and perturbative triple-excitations [CCSD(T)], scales as O(N 7 ) in terms of the CPU time, and O(N 4 ) in terms of memory and storage requirements, where N is a measure of the molecule size. The steep scaling creates a “scaling wall” that cannot be overcome even with the use of massive parallelization techniques on supercomputers. This problem can be significantly alleviated by exploiting the short-range character of electron correlation using local molecular orbitals (LMOs) and domain approximations. 1–6 A large variety of local correlation methods have been published in the past, 7–60 but a scalable method with well-controlled errors is still to be developed. Ideally, the cost of a “scalable” method should scale linearly with the number of correlated electrons and inverse-linearly with the number of CPU cores. Toward this direction, we reported efficient pair natural orbital local second-order Møller-Plesset perturbation theory 61 (PNO-LMP2) and coupled cluster methods 62 (PNO-LCCSD) in earlier papers of this series. Although the PNO-LCCSD method is not completely scalable due to some slowly-decaying contributions in the coupled-cluster equations and significant communication requirements between computer nodes, it scales nearly linearly with the molecular size and is well parallelized up to 100-200 compute cores. The local approximations include truncations in the external spaces (domain approximations) and of the list of correlated LMO pairs (pair approximations), as well as a technique of reducing the number of required integral blocks 2 ACS Paragon Plus Environment

Page 2 of 74

Page 3 of 74

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Theory and Computation

(projection approximations). 49,50,55 In our previous work 62 all these approximations have been described in detail, and this is not repeated here. We have shown that the errors caused by the PNO and projection approximations are likely below 1 kJ mol−1 for relative energies. Nevertheless, the differences to canonical CCSD may still be larger, which is mainly caused by using domains of projected atomic orbitals (PAOs) for determining the PNOs. The correlation energies are sometimes slowly convergent with respect to the sizes of these domains. We believe that this is mostly due to intramolecular basis set superposition effects (BSSE), which are largely removed when small domains are used, 15,63 but re-introduced when the domains are increased in order to reduce other domain errors. These combined effects of BSSE and domain errors, which may either partly cancel or augment each other, make it difficult to predict the quality of a local correlation calculation. This problem is largely removed in the present work through the inclusion of explicitly correlated terms. The very slow basis set convergence of the correlation energy in wave function based methods, and related to this the BSSE problem, are due to the fact that expansions in products of one-electron functions do not describe the shape of the wave function well at small to intermediate values of the interelectronic distance, which should satisfy the electronic wave function cusp conditions. 64,65 Including terms that are explicitly dependent on the interelectronic distances greatly accelerates the basis-set convergence. 61,66–122 In particular, the CCSD(T)-F12x (x=a,b) methods 90,103,108 have been demonstrated as highly cost-effective ways of obtaining not only accurate intramolecular correlation energies, but also those of van der Waals interactions. 123–125 The F12 methods are known to reduce basis set incompleteness errors and BSSE, as well as domain errors in local correlation methods, 112,113 making the combination of local and explicitly correlated methods particularly interesting and useful. A number of explicitly correlated local correlation methods have been developed in the past, and these methods give excellent accuracy. Previously, our group has developed PAO-based explicitly correlated local methods denoted PAO-LMP2-F12 and PAO-LCCSD(T)-F12, using efficient local density fitting (DF) techniques. 126,127 Schmitz et al. 54,128–130 have developed reduced scal-

3 ACS Paragon Plus Environment

Journal of Chemical Theory and Computation

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

ing PNO-based MP2-F12, MP3-F12, and several CCSD(T)-F12 methods, with both iterative and perturbative treatments of the explicit correlation. Recently, Pavoševi´c et al. have developed nearlinear scaling DLPNO-LMP2-F12 122 and DLPNO-CCSD(T)-F12 131 methods. All these methods are only moderately parallel and have a rather long time-to-solution for calculations on large molecules using quality basis sets. It should be noted that in all aforementioned F12 methods the perturbative (T) triples correction does not include F12 terms. The development of well-parallelized coupled-cluster methods is a great challenge due to the large amount of data to be stored and communicated over the network. Many well parallelized canonical CCSD programs have been developed in the last decades, 132–142 and some methods have shown impressive parallel efficiency using thousands of CPU cores. Many of these programs used an integral-direct strategy to reduce the communication cost, in which the 2-electron integrals with 3 or 4 virtual indices are computed on-the-fly in each iteration. Some of the programs also employed optimized tensor contraction frameworks to alleviate the communication problem. These approaches are less helpful in PNO-based local coupled-cluster methods, in which computing and transforming the AO integrals in each iteration would lead to very high computational redundancy, and the use of pair-specific virtuals leads to numerous small tensor contractions that are not suitable for parallelized tensor frameworks. In this paper we report an efficient parallel PNO-LCCSD-F12x (x=a,b) method. In the following, we simply refer to the method as PNO-LCCSD-F12. The theoretical framework used is similar to that of the PAO-LCCSD-F12 method 126 developed earlier in our group, but here we employ much more accurate domain approximations using PNOs, more accurate local approximations for weak pairs, and more efficient integral transformation routines with advanced parallelization techniques to improve the accuracy and efficiency of the method. As in our preceding works, we aim at developing efficient implementations that can handle three-dimensional molecules with 100–200 atoms using augmented triple-ζ basis sets, and are parallelized with reasonable efficiency on up to a few hundred CPU cores. The paper is organized as follows: In Section 2 we briefly outline the PNO-LCCSD-F12 the-

4 ACS Paragon Plus Environment

Page 4 of 74

Page 5 of 74

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Theory and Computation

ory. In section 3 we discuss the local approximations in the PNO-LCCSD-F12 method and present benchmark results on these approximations. Some technical details on the parallel implementation of the method are given in Section 4. We explain in Section 5 the integral transformation algorithms developed for the method, which are critical to the performance and scalability of the PNO-LCCSD-F12 method. In Section 6 we present comprehensive results on several test sets. A discussion in Section 7 concludes the paper.

2 PNO-LCCSD-F12 Theory For the sake of clarity, we show in Table 1 the abbreviations and index notations for different orbital spaces. In the present work, the occupied molecular orbitals are localized using the intrinsic bond orbital (IBO) method. 143 Core and valence orbitals are localized separately. While the core orbital localization has no effect in LCCSD calculations, it slightly affects the F12 energies if domain approximations are applied to the projector, as described below. The PNO-LCCSD-F12 theory is adapted from the PAO-LCCSD-F12 theory developed in our group, 126 and we only provide a brief outline here. Details of the local approximations applied are given in Section 3. Table 1: Index Notation Used for Different Orbital Spaces orbital space indices localized valence molecular orbitals i, j localized occupied molecular orbitals m, n pair natual orbitals (PNOs) ai j , bi j , ... projected atomic orbitals (PAOs) r, s density fitting (DF) basis functions A, B orthonormal RI functions α, ¯ β¯ formally complete orthonormal orbital basis α, β atomic orbital (AO) basis functions µ, ν The PNO-LCCSD-F12 wave function is defined as   ΨLCCSD-F12 = exp Tˆ 1 + Tˆ 2 |0i

5 ACS Paragon Plus Environment

(1)

Journal of Chemical Theory and Computation

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 6 of 74

where |0i is the Hartree–Fock reference wave function, and the excitation operators are defined as Tˆ 1 =

X X i

tai Eˆ ia

(2)

a∈[ii]PNO

and 1 Tˆ 2 = 2 +

X

X

i j ˆ ab Ei j T ab

i j a,b∈[i j]PNO

1 X X X ij T kl hkl|F12 Qˆ i12j |αβi Eˆ iαβj 2 i j αβ kl

(3)

ij where tai and T ab are conventional singles and doubles amplitudes, T kli j are the amplitudes of the

explicit correlation, and Eˆ ia (and similarly Eˆ iα ) are the spin-summed excitation operators. The notation [i j]PNO represents a domain of PNOs for the LMO pair i j (cf. Section 3). Throughout this work, the domains of PNOs are indicated by the range of summation or shown explicitly after the equations, and the labels a, b, ... are used as the PNO indices. We use the fixed amplitudes determined from the first-order wave function cusp conditions 79,80 3 1 T kli j = δik δ jl + δil δ jk 8 8

(4)

F12 is a Slater-type geminal that is approximated by 6 Gaussians 6

1 X −αi r2 1 ci e 12 F12 = − e−γr12 ≈ − γ γ i=1

(5)

Details about this fit are described in ref 88. This correlation factor depends on a length-scale parameter γ. For valence-shell correlation γ = 1 a−1 0 is usually a good choice and is used throughout this work. Qˆ i12j is a projector that keeps the explicitly correlated configurations orthogonal to the reference Hartree-Fock function as well as to the conventional singly and doubly excited configuˆ ab rations |Φai i = Eˆ ia |0i and |Φab i j i = E i j |0i. In canonical methods the projector is independent of the

6 ACS Paragon Plus Environment

Page 7 of 74

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Theory and Computation

pair label i j, and the so-called Ansatz 3 77,88 Qˆ 12 = (1 − oˆ 1 )(1 − oˆ 2 ) − vˆ 1 vˆ 2 is used, where oˆ =

P

m

|mi hm| and vˆ =

P

a

(6)

|ai ha| are one-electron projection operators onto the

occupied and virtual orbitals, respectively. In local correlation methods the projector becomes pair-specific, since the terms vˆ 1 vˆ 2 are restricted to the domain of the corresponding pairs. This can be considered as an “ansatz”, which is discussed in more detail in Section 3.3. Furthermore, local resolution of the identity (RI) approximations 115 are applied in order to achieve linear scaling. The projector then takes the form Qˆ i12j =1 −

X

|abi hab| +

X

|mni hmn|

m,n∈[i j]LMO

a,b∈[i j]PNO



X

(|mαi ¯ hmα| ¯ + |αmi ¯ hαm|) ¯

(7)

m∈[i j]LMO ,α∈[i ¯ j]RI

where [i j]LMO represents a list of LMOs that are spatially close to i or j, and [i j]RI represents a domain of RI basis functions associated to the pair i j, as described in ref 117. The F12 operator, projector, and amplitudes are the same as defined previously in our PNO-LMP2-F12 method. We assume that the RI basis includes the occupied space, but do not use the complementary auxiliary basis set (CABS) formulation. 82 Exactly the same results as with the CABS approach can then be obtained if the union of an AO basis set and an RI basis set is used to approximate the complete space, provided the union of both sets is not linearly dependent. This can be achieved with the so-called OPTRI basis sets of Yousaf and Peterson. 144,145 However, the resulting RI basis sets are rather large and have therefore not been used in the current work. Many tests have shown that equally accurate results are obtained by using the somewhat smaller JKFIT auxiliary basis sets 146 for the RI. In fact, these often yield better CABS singles corrections 90,103 of the Hartree-Fock (HF) energies, in particular if diffuse functions are included.

7 ACS Paragon Plus Environment

Journal of Chemical Theory and Computation

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 8 of 74

The amplitudes in eq 1 are determined by solving the usual coupled LCCSD equations ˆ

ˆ

∀a ∈ [ii]PNO

(8)

ˆ

ˆ

∀a, b ∈ [i j]PNO

(9)

˜ ai |e−T He ˆ T |0i = 0, rai = hΦ −T ˆ T ˜ ab Riabj = hΦ i j |e He |0i = 0,

where the contravariant configurations are defined as 1 ˆa E |0i 2 i

(10)

1  ˆ a ˆ b ˆ b ˆ a 2Ei E j + Ei E j |0i 6

(11)

˜ ai i = |Φ ˜ ab |Φ ij i =

These amplitude equations are solved iteratively, as described in ref 62. A formally complete treatment of the CCSD-F12 method introduces terms in which the unity operator in the projector must be approximated with a double RI. This is expensive and limits the application of the method to very small molecules. In addition, it may also introduce additional errors since the double RI approximation is slowly convergent with the RI basis set size. In the PNO-LCCSD-F12 method we use the so-called CCSD-F12x (x=a,b) approximations, in which no double RI approximations are needed. Furthermore, only contributions that are linear in the F12 amplitudes are considered beyond LMP2-F12. The CCSD-F12x method has been extensively benchmarked, 103,108,147,148 and has been found to give similar or even better results than more complete treatments. 85,99 This is due to certain systematic error cancellations and the avoidance of multiple RI approximations. 103 In addition, we also use the approximation 3*A 88 in PNO-LCCSDF12. The PNO-LCCSD-F12 and the PNO-LCCSD residual equations 62 then differ only slightly: ∆Riabj

ij =V˜ ab −

X   ij ij k t¯b + t¯ak V˜ kb V˜ ak

∀a, b ∈ [i j]PNO

(12)

k∈[i j]CC LMO

∆rai =

X X j

[ii,i j] S ab (2V˜ bi jj − V˜ bjij ) ∀a ∈ [ii]PNO

b∈[i j]PNO

8 ACS Paragon Plus Environment

(13)

Page 9 of 74

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Theory and Computation

where [i j,kl] S ab = hai j |bkl i ij −1 ˆ i j V˜ ab = hab|r12 Q12 F12 |i˜ji

(14) (15)

[i j,i j] ij Since the PNOs for a given pair are orthonormal, S ab = δab . The quantity V˜ ak is defined similarly.

The modified kets |i˜ji =

3 1 |i ji + | jii 8 8

(16)

account for the fixed amplitudes defined in eq 4. For convenience we also define projected amplitude vectors and matrices t¯aki j =

X

[i j,kk] k S ac tc

(17)

c∈[kk]PNO

T¯ ai klj bmn =

[kl,i j] i j [i j,mn] S ac T cd S cb

X

(18)

c,d∈[i j]PNO

This avoids the explicit occurrence of overlap matrices in the LCCSD-F12 amplitude equations. 62 The PNO-LCCSD-F12a energy expression differs from the PNO-LCCSD one only by adding MP2 the MP2-F12 energy correction ∆EF12 :

LCCSD-F12a LCCSD MP2 Ecorr = Ecorr + ∆EF12

(19)

Note, however, that due to the additional terms in the residuals, the amplitudes are different, and LCCSD therefore the Ecorr contribution is different than for LCCSD. The F12b approximation adds fur-

9 ACS Paragon Plus Environment

Journal of Chemical Theory and Computation

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 10 of 74

ther terms which introduce additional couplings with the conventional amplitudes: LCCSD-F12b Ecorr

 X  X LCCSD-F12a  =Ecorr +  ij

+

X 

2V˜ ai jj

a,b∈[i j]PNO



V˜ ajij

a∈[i j]PNO



  ij ji ˜ i j Vab 2Dab − Dab

 X   j  ji i j t¯ai + 2V˜ bi − V˜ bi t¯b 

(20)

b∈[i j]PNO

where ij Diabj = T ab + t¯ai t¯bj ,

a, b ∈ [i j]PNO

(21)

The F12b approximation is more rigorous than F12a. Nevertheless, it has been found that in canonical calculations the F12a approximation sometimes gives more accurate energy differences if small AO basis sets are used. 103 This is due to a systematic overestimation of the F12a energy correction, which favorably compensates for basis set incompleteness errors. In addition to F12a and F12b, we also tested the cruder approximation of simply adding the F12 correction obtained from PNO-LMP2-F12 calculations to the PNO-LCCSD energies. This approach is denoted ∆F12. As will be shown later, it overestimates the F12 correction even more than the F12a approximation. ij The quantity V˜ ab is computed as 113,115

ij −1 V˜ ab = hab|r12 F12 |i˜ji −

X

c,d∈[i j]PNO



X

,c∈[i j]PNO m∈[i j]CC LMO



cd ˜ i j Fcd − Kab

X

mn ˜ i j Kab Fmn

m,n∈[i j]CC LMO

mc ˜ i j cm ˜ i j Kab Fmc + Kab Fcm



(22)

where ij F˜ xy = hxy|F12 |i˜ji

(x, y representing PNOs or LMOs)

(23)

xy ij −1 and the other two-electron integrals are defined as Kab = hab|r12 |xyi. We note that F˜ xy should

be distinguished from the dressed Fock matrices defined in the preceding paper. 62 In eq 22 the 10 ACS Paragon Plus Environment

Page 11 of 74

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Theory and Computation

summation over the occupied MOs is restricted to slightly different domains as used in eq 7 to reduce the cost, and the summation over the virtual orbitals in the last term is restricted to the PNO pair domain. These two approximations are discussed in detail in Section 3.3, and the errors ij introduced were found to be small. The V˜ ak operators can be computed in a similar fashion:

ij −1 V˜ ak = hak|r12 F12 |i˜ji −

X

c,d∈[i j]PNO



X

m∈[i j]CC ,c∈[i j]PNO LMO



cd ˜ i j Fcd − Kak

X

mn ˜ i j Kak Fmn

m,n∈[i j]CC LMO

mc ˜ i j cm ˜ i j Fmc + Kak Fcm Kak



(24)

As discussed in the following sections, in the PNO-LCCSD-F12 method, the number of pairs treated with F12 correction increases linearly with the molecule size. All summations are restricted to domains that asymptotically have a constant size. Using the local density fitting technique for the integral evaluations and transformations, the extra cost of the F12 terms scales linearly with molecular size.

3 Local Approximations In this section we discuss the local approximations used in the PNO-LCCSD-F12 method. These include pair approximations, domain approximations, and projection approximations. In order to test these approximations we used the “Auamin” and “Isomer4” reactions (Figure 1) as the benchmark systems. Auamin is a gold(I)-aminonitrene complex (AuC41 H45 N4 P) taken from ref 149. The Isomer4 reaction is the fourth reaction from the ISOL24 benchmark of Huenerbein et al., 150 which is known to be strongly affected by dispersion interactions. Both reactions are difficult cases for local correlation methods 62 and yet not exceedingly large so that benchmark calculations can be performed with very tight thresholds within a reasonable amount of time. All results presented in this section were computed with the VTZ-F12 basis set. 151 The computational details can be found in Section 6.1. Benchmark calculations were also carried out for the so-called FH test set developed by Friedrich 11 ACS Paragon Plus Environment

Journal of Chemical Theory and Computation

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 12 of 74

and Hänchen. 152 This set contains small to medium sized molecules so that a direct comparison with canonical results is possible. The last 4 reactions involving rather large molecules in the set are not included here since the canonical CCSD-F12 calculations for them are very difficult to carry out. Isomer4 reaction H O HO

H Auamin reaction

P N N

N

Au N N

Au N N

+

N

P

Auamin

Figure 1: Benchmark systems for testing the PNO-LCCSD-F12 method.

3.1 Pair Approximations As in the PNO-LCCSD method, nondistant LMO pairs are classified into strong, close, and weak pairs based on their LMP2 pair energies, using by default the thresholds of T close = 10−4 Eh and T weak = T close /10. We only treat the strong pairs with the full LCCSD-F12 theory, i.e., we apply the F12 correction to the doubles residuals (eq 12) for strong pairs, and restrict the summation over j in eq 13 to strong pairs of i j. For close and weak pairs, the LMP2-F12 correction is simply added to the LCCSD pair energies. The F12 correction of distant and very distant pairs is completely neglected, and the same PNO-LMP2 treatment with multipole approximations as described in ref 62 is used. Figure 2 shows the PNO-LMP2 pair correlation energies and the F12 corrections as a function of the orbital distances Ri j = |ri − r j |, where ri and r j are the charge centers of LMOs i and j, respectively, for the Auamin molecule. Each point in the figure corresponds to one pair. For pairs 12 ACS Paragon Plus Environment

Page 13 of 74

Journal of Chemical Theory and Computation

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 ACS Paragon Plus Environment

Journal of Chemical Theory and Computation

a function of T close for PNO-LCCSD and the three F12 treatments. The figure nicely demonstrates 46

(a) Isomer4 reaction reaction energy / kcal mol−1

64.5 reaction energy / kcal mol−1

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 14 of 74

64.0 LCCSD-F12a LCCSD-F12b 63.5

63.0

(b) Auamin reaction

45 LCCSD-F12a LCCSD-F12b 44

43

LCCSD LCCSD

42 10−3

10−4 10−5 Tclose / Eh

10−6

10−3

10−4 10−5 Tclose / Eh

10−6

Figure 3: PNO-LCCSD and PNO-LCCSD-F12 reaction energies of the Isomer4 (left) and Auamin (right) reactions as a function of the pair selection thresholds T close . Default parameters and T weak = T close /10 are used. The LCCSD values include the default domain corrections. Basis set: VTZ-F12. that the pair approximations are justified and that the reaction energies smoothly converge to the pure strong pair limits, in which these approximations are absent. Both the F12a and the F12b approaches show a very similar dependence on T close . With our default threshold T close = 10−4 Eh , the error is around 0.1 kcal mol−1 , which is less than the intrinsic error of the CCSD(T) method, and also less than the differences of the F12a and F12b approximations. In both systems, the reaction energy decreases in the order of LCCSD-F12a, LCCSD-F12b, LCCSD. The differences between the F12a and F12b results typically amount to a few tenth of a kcal mol−1 . We consider the F12b results as most reliable. They are not only theoretically the most accurate, but also show the smoothest convergence with respect to domain and basis set sizes, as is demonstrated later.

3.2 Domain Approximations The virtual space in the PNO-LCCSD-F12 method is represented by PNOs as in the PNO-LCCSD method. The details for the PNO-generation procedures are reported in previous publications. 61,62 In short, after localizing the occupied molecular orbitals, we compute the PAOs and create orbital 14 ACS Paragon Plus Environment

Page 15 of 74

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Theory and Computation

domains [i]PAO . Each orbital domain is determined by a list of centers (atoms), and all PAOs that originate from basis functions at the selected centers are included. The list is obtained by first selecting a primary domain, which contains centers whose partial IBO charges of i are greater than 0.2 e, and then extending it by adding all centers which are within IEXT layers of atoms or within a radius REXT=(2 × IEXT+1) a0 from any atom in the primary domain (default IEXT=2). Following that, orbital specific virtual orbitals (OSVs) are created for each LMO, and the OSV domains [i]OSV are selected by an occupation number threshold T OSV (default 10−9 ). Finally, the PNOs are generated from the OSVs for each LMO pair, with [i j]OSV = [i]OSV ∪ [ j]OSV . The PNOs are truncated with an occupation number threshold T PNO (default 10−8 ). However, if the semicanonical MP2 PNO pair energy is still less than the corresponding OSV pair energy multiplied by a threshold T en (default 0.997), extra PNOs are added. This additional criterion is important to ensure the PNO domains of weak pairs are sufficiently large. 62 The cost and accuracy of the local correlation methods depend heavily on the domain sizes. In fact, the finite PAO domain sizes are probably causing the largest errors in PNO-LCCSD calculations. Figure 4 shows the dependence of the reaction energies of the two test reactions on IEXT for PNO-LCCSD and the three F12 treatments. Without the F12 terms, the convergence is painfully slow. For both reactions, the difference in the PNO-LCCSD reaction energies from IEXT=2 to IEXT=3 is still larger than 1 kcal mol−1 . With IEXT=3, the maximum PAO pair domain sizes may exceed 1600 orbitals and calculations for large molecules become prohibitively expensive. This problem is much alleviated by using either of the F12 approximations. The convergence is smoothest with the F12b treatment, while the F12a approximation overshoots the correlation energy and overcompensate the inadequate description of the external space with small PAO domains. With F12b, our default IEXT=2 leads to well converged results for both cases, and even IEXT=1 leads to only 0.03 and 0.3 kcal mol−1 error for the Isomer4 and Auamin reactions, respectively. As already mentioned in the introduction, the slow convergence with the PAO domain sizes is at least partly caused by BSSE, which are largely removed in local correlation methods with small PAO domains, but re-introduced with increasing PAO domain sizes. The F12 terms are known

15 ACS Paragon Plus Environment

Journal of Chemical Theory and Computation

66

(a) Isomer4 reaction

(b) Auamin reaction 46

LCCSD-F12a

64

reaction energy / kcal mol−1

reaction energy / kcal mol−1

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 16 of 74

LCCSD-F12b 62 LCCSD 60 58 56

LCCSD-F12a 44 LCCSD-F12b 42 LCCSD 40

38 0

1

2

3

0

1

2

3

IEXT

IEXT

Figure 4: Reaction energies of the Isomer4 (left) and Auamin (right) reactions as a function of the PAO domain parameter IEXT obtained from PNO-LCCSD and different F12 treatments. Default PNO domains and REXT = (2 × IEXT + 1) a0 is used. to reduce the BSSE and domain errors significantly. 112,113 In our previous work on PNO-LMP2F12, 117 we found that this effect reduced the errors caused by small PAO domains by an order of magnitude for a cluster of 100 H2 O molecules, using the VTZ-F12 basis set. From Figure 4 we see that the reduction of PAO domain errors is also very effective in PNO-LCCSD-F12. More results on this topic are given in the benchmark of intermolecular interactions in Section 6.6. We now consider the PNO domain truncation. Figure 5 shows the reaction energies as a function of T PNO using different values of T en . The convergence is rather slow for the PNO-LCCSD method without any domain corrections (dotted lines). This is significantly improved by using the semi-canonical energy corrections or using the F12 treatments. Consistent with previous studies, 131 the F12 terms are less effective in correcting the PNO domain errors for very small PNO domains (T en = 0, T PNO < 10−7 Eh ). This may have to do with the rapid dropping of the PNO domain sizes with increasing orbital distances. 62 In the worst case scenario, the PNO domain sizes can become zero for weak pairs, and the PNO-LMP2-F12 energy for these pairs is then zero. The PNO energy threshold prevents this problem, and therefore improves the convergence significantly. Similar to the composite approach used in the PNO-LCCSD method, 62 it is possible to use smaller PNO domains in PNO-LCCSD-F12 calculations than in the preceding PNO-LMP2-F12 calculations without significant loss of accuracy. This saves a lot of computational resources (time, 16 ACS Paragon Plus Environment

Page 17 of 74

64 62 60 reaction energy / kcal mol−1

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Theory and Computation

LCCSD-F12a LCCSD-F12b LCCSD

58 56 54 10−6

10−8

10−9

Isomer4, Ten = 0.99

Isomer4, Ten = 0.9

Isomer4, Ten = 0 10−7

10−6

10−7

10−8

10−6

10−7

10−9 10−6

10−7

10−9

10−8

10−9

46 44 42 40 38 36 10−6

Auamin, Ten = 0 10−7

10−8

10−9 10−6

Auamin, Ten = 0.99

Auamin, Ten = 0.9 10−7 10−8 TPNO

10−8

10−9

Figure 5: PNO-LCCSD(-F12) reaction energies of the Isomer4 (top panels) and Auamin (lower CC CC panels) as a function of T PNO for various values of T en . T PNO and T en are set to T PNO and T en , respectively. The full and dotted black lines represent the PNO-LCCSD energies with and without the semi-canonical domain corrections, respectively. memory, disk) and is therefore done by default in our program. The small domains are determined CC CC by an occupation number threshold T PNO (default 10−7 ) and an energy threshold T en (default 0.99).

First an LMP2-F12 calculation is performed with the large domains (determined by the thresholds ij ij ij T PNO and T en ), keeping the Kab , Fab , and Uab integrals (see ref 117) in GAs. Next, an LMP2 calcu-

lation with small domains is carried out and the transformation matrices from the large to the small pseudo-canonical domains are obtained. The integrals are then transformed to the small domains and the 2-external contributions to the LMP2-F12 energies are recomputed. Other contributions to the F12 correction are invariant to the PNO domain sizes and can be taken from the large-domain calculation. The LCCSD-F12 calculations are carried out with the small domains. Finally, the correction ∆ELMP2-F12(small) = ELMP2-F12(large) − ELMP2-F12(small) is added to the final LCCSD-F12 energy computed with the small PNO domains.

17 ACS Paragon Plus Environment

(25)

Journal of Chemical Theory and Computation

66 64 62 60 reaction energy / kcal mol−1

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 18 of 74

LCCSD-F12a LCCSD-F12b LCCSD

58 56 54 10−6

CC = 0 Isomer4, Ten

10−7

10−8

10−9 10−6

CC = 0.9 Isomer4, Ten

10−7

10−8

10−9 10−6

CC = 0.99 Isomer4, Ten

10−7

10−8

10−9

46 44 42 40 38 −6 10

CC = 0 Auamin, Ten

10−7

10−8

10−9 10−6

CC = 0.9 Auamin, Ten

10−7 CC 10−8 TPNO

10−9 10−6

CC = 0.99 Auamin, Ten

10−7

10−8

10−9

Figure 6: PNO-LCCSD(-F12) reaction energies of the Isomer4 (top panels) and Auamin (lower CC CC panels) as a function of T PNO for various values of T en . LMP2(-F12) domain corrections from −9 large PNO domains (T PNO = 10 , T en = 0.997) are applied (see text). The full and dotted black lines represent the PNO-LCCSD energies with and without the semi-canonical domain correction of the large domain PNO-LMP2 calculations, respectively. The results with these corrections are shown in Figure 6. The convergence with respect to CC T PNO is again improved by the F12 terms, in particular for PNO-LCCSD-F12b. The PNO domain

correction (eq 25) tends to overshoot with very coarse domain thresholds, but the results obtained with our default thresholds T PNO = 10−7 , T en = 0.99 are very close to the converged results. Further reducing the PNO domain sizes would force us to use less aggressive projections if the same level of accuracy is desired. This would then increase the computational cost despite smaller PNO domains. 62 It is also possible to use smaller PAO domains (i.e. IEXT=1) in the LCCSD-F12 calculation and apply the PAO domain corrections at the LMP2-F12 level with eq 25. For the majority of applications that we tested the errors were negligible (< 0.1 kcal mol−1 ). However, since this approach requires two independent PNO-LMP2-F12 calculations, it is not used by default.

18 ACS Paragon Plus Environment

Page 19 of 74

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Theory and Computation

3.3 Local Approximations in the strong-orthogonality projector The fact that the F12 terms correct for domain errors is due to the restriction a, b ∈ [i j]PNO in the pair-specific projector (eq 7), which projects out just the domain [i j]PNO that is included in the conventional calculation. Implicitly, the F12 terms can then describe excitations into the virtual P ij orbitals outside the domain with amplitudes T αβ = k,l T kli j hkl|F12 Qˆ i12j |αβi (cf. eq 3). 112 This is an approximation, since these amplitudes are neither optimized nor coupled to the conventional ones. Nevertheless, as demonstrated in Figures 5 and 6, the approximation normally works very well. Also the domain correction in eq 25 is usually much smaller than the corresponding correction without F12. The reduction of domain errors by the F12 terms is particularly effective for the PAO domains. 112–116 However, we already noticed in our previous work 61 that the F12 terms are somewhat less effective for correcting the PNO domain errors than for the PAO domain errors, and this has also been mentioned by other authors. 153 In particular for small molecules, which contain atoms with lone electron pairs, the PNO-LMP2-F12 energies sometimes converge very slowly and not monotonically with decreasing threshold T PNO toward the full domain (OSV-LMP2-F12 or PAOLMP2-F12) limit. We have therefore re-investigated this problem using the FH test set. 152 Some results are presented in Figure 7. The red lines in the top panels of this figure show the root-means square (RMS) errors of the PNO-LMP2-F12 correlation energies for the 104 molecules of the FH set relative to OSV-LMP2-F12 results. Results obtained both without (left panel) and with (right panel) the energy threshold T en are presented. While the PNO-LMP2-F12 correlation energies converge rather smoothly with decreasing T PNO , for very small thresholds the convergence becomes slow, and noticeable errors still remain even with very small values of T PNO . These errors are more clearly seen when the reaction energies are considered (lower panels of the figure). In this case the errors are not reduced beyond T PNO = 10−8 , and a residual RMS error of about 0.05 kcal mol−1 remains even for T PNO = 10−11 . This small error only disappears if all PNOs are included. A closer inspection showed that a large fraction of the error results from contributions of PNOs with near-zero occupation numbers (niaj < 10−10 ), which do not affect the PNO-LMP2 energy at all, but 19 ACS Paragon Plus Environment

Journal of Chemical Theory and Computation

still contribute to the F12 correction if they are not included in the projector. The slow convergence when using PNOs in the projector have also been discussed by Pavoševi´c et al., 153 and they proposed the use of so-called geminal-spanning orbitals (GSOs) in the projector.

correlation energy errors / mEh

4

PNO projector OSV projector OSV projector w/ sc-PNO

3

(a) Ten = 0, correlation energy 2 1 0 10−7

10−8

0.6

10−9 TPNO

10−10

10−11

(c) Ten = 0, reaction energy

0.4

0.2

0.0 10−7

10−8

10−9 TPNO

10−10

10−11

reaction energy errors / kcal mol−1

correlation energy errors / mEh

However, in their later work they used again the standard PNO projector. 122,131

reaction energy errors / kcal mol−1

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 20 of 74

2.0

(b) Ten = 0.997, correlation energy

1.5

1.0

0.5

0.0 −7 10 0.15

10−8

10−9 TPNO

10−10

10−11

(d) Ten = 0.997, reaction energy

0.10

0.05

0.00 −7 10

10−8

10−9 TPNO

10−10

10−11

Figure 7: RMS errors on correlation energies and reaction energies of PNO-LMP2-F12 for the FH test set due to PNO domain truncations. The reference values are obtained from calculations including all PNOs (which is equivalent to OSV-LMP2-F12). Basis set: VTZ-F12. As an alternative treatment, one can use the OSVs in the F12 strong orthogonality projector, i.e., to replace the summation a, b ∈ [i j]PNO by a, b ∈ [i j]OSV . In this case the F12 terms do not correct for PNO domain errors at all, but they still correct for the PAO and OSV domain errors. In order to reduce the remaining PNO domain errors, the semi-canonical PNO (sc-PNO) domain correction (i.e., the difference of the semi-canonical OSV-LMP2 and PNO-LMP2 energies) can be added to the resulting PNO-LMP2-F12 correlation energy. In Figure 7 the black lines represent the energies obtained with the OSV projector without the sc-PNO correction, and the blue ones include the sc-PNO correction. In both cases the errors converge smoothly to zero. As expected, the sc20 ACS Paragon Plus Environment

Page 21 of 74

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Theory and Computation

PNO correction significantly reduces the errors. However, for our default thresholds T PNO = 10−8 , T en = 0.997 the overall errors of this method are quite similar to those obtained with the standard PNO projector. Only for smaller thresholds the OSV projector leads to more accurate results, while for larger thresholds, in particularly when the energy threshold is zero, the PNO projector tends to produce more accurate results. The OSV projector cannot easily be used in the additional terms of PNO-LCCSD-F12. This would require much larger sets of 3-external and 4-external integrals, making the LCCSD-F12 calculations prohibitively expensive. Probably for the same reasons, Pavoševi´c et al. did also not use their GSOs in the CCSD-F12 terms. 153 Therefore, the only possibility of applying the OSV projector in LCCSD-F12 calculations is to keep the standard PNO projector in the corresponding terms of eqs 22 and 24, but to use the OSV projector in the PNO-LMP2-F12 contribution. Since the terms in eqs 22 and 24 give only a relatively small contribution to the overall correlation energy this is a good approximation. Additional approximations are needed in the summations over LMOs m, n in eqs 22 and 24 in order to achieve linear scaling. Similar to the PNO-LMP2-F12 method, 117 we restrict the summations to a domain [i j]LMO , which includes LMOs m that are close to i or j. For LMP2-F12, this list of LMOs is comprised of a list of valence LMOs and a list of core LMOs. The former list is selected using a pair energy threshold T val , i.e., an orbital of m is included in the list if the PNOLMP2 pair energies Eim or E jm are above T val . The list of core LMOs is selected by distance and connectivity criteria, as described previously. 117 For LCCSD-F12, we neglect the contributions of core orbitals, which has been found to introduce very small errors. 113 Furthermore, by default the LMO domains are linked to the pair classes with [i j]CC LMO = ({i} s ∪ { j} s ) ∩ ({i}c ∩ { j}c )

(26)

in eqs 22 and 24, where {i} s represents a domain of LMOs m where im is strong, and {i}c for a domain of m where im is strong or close. This choice is made mostly for cost considerations, given

21 ACS Paragon Plus Environment

Journal of Chemical Theory and Computation

Isomer4 reaction reaction energy / kcal mol−1

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

67

66

Tval = 10−5 Eh Tval = 10−4 Eh CC Tval = Tval

65

64

63

10−3 CC Tval

10−4 / Eh

10−5

Figure 8: PNO-LCCSD-F12a (dotted lines) and PNO-LCCSD-F12b (solid lines) reaction energies CC of the Isomer4 reaction as a function of T val . T close is set to 10−6 Eh in these calculations, and all other parameters are program defaults. The red circle represents the LCCSD-F12b reaction energy obtained using the default settings. Basis set: VTZ-F12. that we would like to make the most accurate approximations for the F12 terms possible without computing 3- and 4-external integrals other than those needed in the PNO-LCCSD calculations anyway. Furthermore, the LMO domains are determined by the thresholds T close and T weak , and thus there is no need to introduce extra parameters. Solely for the purpose of testing the effects of the LMO domains [i j]CC LMO in the projector, we CC implemented an additional parameter T val that controls the summations in eqs 22 and 24 similarly

to T val for the LMP2-F12 calculations, independent of T close and T weak . The dependence of the CC Isomer4 reaction energy with respect to T val and T val is shown in Figure 8. Clearly, the reaction

energies are much more sensitive to the LMO domains at the LMP2-F12 level than at the LCCSDF12 level. Our default threshold T val = 10−4 Eh leads to reasonably converged results. Note that the LMO domains at the LCCSD-F12 level are defined by eq 26 with T close = 10−4 Eh and T weak = CC 10−5 Eh by default. These domains are slightly smaller than those obtained with T val = 10−4 Eh .

The F12b reaction energy obtained with the default settings is marked with a red circle in Figure 8. In a strict F12 treatment, the last summations in eqs 22 and 24 should run over the complete virtual space, including the CABS part. Previous studies have shown that restricting the summation to domains of PAOs leads to very small errors. 113,115 This is confirmed by the excellent convergence 22 ACS Paragon Plus Environment

Page 22 of 74

Page 23 of 74

of the PNO-LCCSD-F12 method with respect to the PAO domain sizes, as shown in Figure 4. In the present method, the summation is further restricted to PNO domains, in order to avoid the very expensive computation of 3-external PAO integrals. The effect of this restriction is shown in Figure 9. Clearly, this approximation leads to only negligible errors with our default PNO domains CC CC (T en = 0.99, T PNO = 10−7 ). Even with smaller PNO domains the error is still in the 0.1 kcal mol−1

range, and is well below the difference between various F12 treatments.

reaction energy / kcal mol−1

64.5 (a) Isomer4 reaction

64.0

LCCSD-F12a LCCSD-F12b

63.5

63.0

0.90

0.92

0.94

0.96

0.98

1.00

CC Ten

45.5 reaction energy / kcal mol−1

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Theory and Computation

(b) AuAmin reaction

45.0

LCCSD-F12a

44.5

LCCSD-F12b 0.90

0.92

0.94

0.96

0.98

1.00

CC Ten

Figure 9: Dependence of the Isomer4 and Auamin reaction energies on the PNO energy threshold CC CC T en (T PNO = 10−7 ). All other parameters are kept at program defaults. The dotted lines represent the results obtained using c ∈ [i j]PC-PAO instead of c ∈ [i j]PNO in the last two term of eq 22. Finally, we note that some additional projection approximations have been introduced for the F12 terms. For example, the summation over b in eq 13 should be in [ii]PNO in the exact form. ij Here we use the projected singles amplitudes (eq 17) so that the V˜ ak quantities in eq 13 are a subset

of those in eq 12. This approximation reduces the computational cost and the complexity of the implementation. This term represents the coupling of the F12 terms with the singles amplitudes and 23 ACS Paragon Plus Environment

Journal of Chemical Theory and Computation

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

has a rather small contribution to the F12 energy. Similar projection approximations are applied to the last two terms in eq 20 for the computation of the F12b energies. Otherwise, the same projection approximations as in the PNO-LCCSD method are applied to PNO-LCCSD-F12. We found that the effect of these projection approximations in the LCCSD equations is nearly independent of the F12 terms. In order to demonstrate the overall accuracy of the various approximations in the F12 terms, we present in Figure 10 the errors of the PNO-LCCSD-F12 reaction energies for the FH benchmark set. Default settings were employed for these calculations. Panel (a) shows the deviation of the PNO-LCCSD results from the CCSD ones without F12 terms, and panel (b) shows the deviation of the PNO-LCCSD-F12b results from the CCSD-F12b ones. Both local and canonical F12 calculations used the same ansatz, i.e., 3*A with fixed amplitudes. In the canonical calculations density fitting was only used to compute the F12 terms, as described in refs 88 and 103. In the PNOLCCSD calculations as well as in the PNO-LCCSD-F12 calculations with the OSV projector, the sc-PNO domain corrections were included. From the figure we can see that, with both the standard PNO projector and the OSV projector, the errors are significantly smaller than for PNO-LCCSD. This is attributed to the reduction of the PAO domain error by the F12 terms. This improvement would be much more pronounced if the scPNO domain correction had not been applied in the PNO-LCCSD calculations. A corresponding figure without this correction for PNO-LMP2-F12 can be found in ref 117. Furthermore, it can be seen that with the default domains the errors are of comparable size for both F12 projectors. We have also computed the PNO-LCCSD-F12 reaction energies of a few larger systems with both projectors with default domains, and did not observe significant difference either. A statistical analysis of the errors of the PNO-LCCSD-F12 reaction energies of the FH benchmark set yields an RMS error is below 0.1 kcal mol−1 (using the standard PNO projector and the VTZ-F12 basis set, see Section 6.4). However, there are still a few cases with larger errors, up to 0.22 kcal mol−1 . These errors can only be reduced by using the OSV projector and very small domain thresholds (T PNO ≤ 10−9 ). Such small thresholds strongly increase the computation time

24 ACS Paragon Plus Environment

Page 24 of 74

Page 25 of 74

Journal of Chemical Theory and Computation

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 ACS Paragon Plus Environment

Journal of Chemical Theory and Computation

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

4 Parallelization Our program is parallelized with the help of MPI and the GlobalArrays (GA) toolkit. 154 The GA toolkit provides interfaces for one-sided remote memory access which is helpful in many of our communication tasks. For most quantities, we use standard GAs comprised of equal-sized sharedmemory chunks from each node. However, for the PNO overlap matrices we use irregular GAs and make sure that all communication of these quantities during the iterations is intra-node. These GAs within a node are implemented using shared memory, and the intra-node access of the overlap matrices is then almost as efficient as keeping them redundantly in local memory of each CPU core. In addition to GAs, the standard MPI point-to-point communication is used in the redistribution of non-duplicated data. We use both static and dynamic parallelization techniques. In the static parallelization algorithms, we estimate the cost of the tasks and assign them to various processors. We use this approach in the iterative solution of the PNO-LCCSD-F12 equations. This allows to keep most data (integrals, residuals, and amplitudes) in memory, in node-specific GAs, or on local disks, which strongly reduces the communication through the network. In order to minimize data redundancy, we try to treat residuals that share a lot of data on the same processor. An optimal assignment of the tasks in a way that minimizes communication and memory usage, and at the same time leads to good load balance of all processors, is highly non-trivial. We use the METIS graph partitioning program 155 for this purpose, but our current implementation is still far from being perfect. The disadvantage of the static parallelization is the unavoidable waste of computational power due to uneven task distribution. In the dynamic parallelization approach the tasks are allocated to the processors on a first-come first-serve basis to minimize the waiting time. We use this approach in the integral transformation step, since the 3-index integrals are mostly shared among many tasks. Communication in this step is significant and unavoidable, otherwise a lot of duplicated work would have to be done. The communication requirements make it very difficult to predetermine the cost of each task, and a dynamic parallelization scheme is therefore highly desired.

26 ACS Paragon Plus Environment

Page 26 of 74

Page 27 of 74

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Theory and Computation

5 Integral Transformations When PNOs are used as virtual orbitals, the computational cost of an LCCSD-F12 calculation is largely migrated from solving the amplitude equations to the transformation of 2-electron integrals. In this section we describe our general approach to obtain the PNO-based integrals in an efficient and parallel linear-scaling algorithm, followed by details on specific classes of integrals.

5.1 General Strategy Similar to our previous PNO-LMP2-F12 program, we use the local density fitting (LDF) technique 23 to evaluate the two-electron integrals. In general, the 2-electron Coulomb integrals are obtained from (uv|wx) ≈

fit X

A duv (A|wx)

(27)

A

where the indices u, v, w, x represent any orbitals, including LMOs, PAOs, and PNOs. The fitting A coefficients duv are obtained by solving linear equation systems fit X B (A|uv) = (A|B)duv

(28)

B

The 2- and 3-index integrals are defined as

(A|B) = (A|wx) =

Z Z

R3

R3

dr1

Z

−1 dr2 χA (r1 )r12 χB (r2 )

(29)

dr1

Z

−1 dr2 χA (r1 )r12 χw (r2 )χ x (r2 )

(30)

R3

R3

Table 2 shows all Coulomb integrals used in the PNO-LCCSD-F12 program and our route of integral transformations. The LCCSD equations have been presented in the previous publication, 62 and the only additional integrals used in the F12 terms are those that appear in eqs 22 and 24. The fitting domains we used for eqs 27 and 28 are also shown in Table 2, where a pair domain [i j]fit represents the union of the fitting domain for the LMOs i and j, determined similarly to the PAO

27 ACS Paragon Plus Environment

Journal of Chemical Theory and Computation

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 28 of 74

Table 2: Route of Transformation of the 2-electron Coulomb Integrals in the PNO-LCCSD-F12 Method and the Fitting Domains Used integral route of transformation Fitting domains default storage 4-external (A|ai j bi j ) → (ai j bi j |ci j di j ) [i j]fit local diskc ij ij ij ij ij ij ij i j i j ik symmetric 3-external (A|rk), (A|a b ) → (a b |rk) → (a b |c k), (a b |c k) [i j]fit local diskc ij ij ij ij ij jj ij 3-external for singles residuals (A|a r), (A|c i) → (a r|c i) → (a b |c i) [i j]fit local diskc ij ij i j ik j j a asymmetric 3-external (A|sk), (A|a r) → (a r|sk) → (a b |c k) [i j]fit or [ik]fit GA symmetric 2-external Coulomb (A|kl), (A|ai j bi j ) → (ai j bi j |kl) [i j]fit GA asymmetric 2-external Coulomb (A| jk), (A|ai j r) → (ai j r| jk) → (ai j bik | jk) [i j]fit or [ik]fit a GA 2-external exchange (A|ri) → (ri|s j) → (akl i|bmn j) [i j]fit GA 1-external (A|rk), (A|lm) → (rk|lm) → (ai j k|lm) [i j]fit GA 0-external (A|ik), (A| jl) → (ik| jl)i j b [i j]fit GA a b c

The fitting domain of the weaker pair is used by default. The integrals are pair-specific and do not have the usual permutation symmetry due to the choice of fitting domains. Communicated over the network and sorted once before the PNO-LCCSD-F12 iterations.

domains. 61 With the exception of the 0-external integrals, we have chosen fitting domains that keep the usual permutation symmetry of the Coulomb integrals. Therefore we do not distinguish symmetry-equivalent integrals used in different terms of the amplitude equations. As a general strategy, we first compute the 3-index integrals of the types (A|i j), (A|ri), and (A|rs), and store these integrals in GAs. This requires the evaluation of the 3-index integrals in the AO basis and then transforming them twice to the desired PAO or LMO bases. This step is parallelized over blocks of DF functions and the details are described in our previous publications. 61,117 The LMO and PAO indices i, j, r, s are limited to domains which depend on the center at which the fitting function A is located. Proper determination and counting of the required blocks of these integrals is complicated but is essential to obtain a linear scaling method. In the next step we loop over pairs of LMOs and perform further integral assemblies and transformations. Depending on the integral type, some indices are transformed from PAO to PNO only after the 4-index integrals are assembled, while in other cases the 3-index integrals are partially or fully transformed before the assembly. This order of the integral transformation, as presented in Table 2, is determined from cost considerations. For example, the communication of the (A|rs)type integrals is very expensive. To minimize the communication cost while keeping memory usage moderate, the required integrals are fetched from the GA once for each pair i j, transformed partially to PNO [i.e., (A|ai j r)], and kept in memory for the assembly of various types of 4-index 28 ACS Paragon Plus Environment

Page 29 of 74

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Theory and Computation

integrals.

5.2 4-external Integrals The only term in the canonical CCSD amplitude equations involving the 4-external integrals is Riabj =

X

ij (ac|bd)T ab + ...

(31)

cd

which scales as O(N 6 ) and is normally the most demanding step in canonical CCSD calculations. In PNO-LCCSD this term reads Riabj =

X

ij (ai j ci j |bi j di j )T cd + ..., a, b ∈ [i j]PNO

(32)

cd∈[i j]PNO

Thus, for each pair i j a distinct set of integrals is needed, but nevertheless the total number of integrals is significantly reduced when using PNOs. For an average pair domain of 50 PNOs, the 4-external integrals take about 6 MB of storage per pair, and in-memory assembly and contraction is feasible. Still, in very accurate calculations the size of some domains may exceed 200 PNOs, and the fourth power scaling may then lead to memory bottlenecks. We therefore implemented fallback algorithms to block-wise assemble the integrals and store them on disks or in GAs. A corresponding block-wise contraction algorithm is also implemented. In our static pair distribution scheme of the PNO-LCCSD-F12 iterations, the 4-external integrals needed by different processors do not overlap. This allows us to store the 4-external on local disks of each node. For performance reasons, we assemble and store the 4-external integrals using a dynamical pair distribution, and the integrals are sorted before the iterations using standard MPI point to point communications. We found that the sorting step takes only a small fraction of the total elapsed time. The large disk space usage in accurate PNO-LCCSD-F12 calculations is inevitable. The 4-external integrals require about 150 GB of disk space in total for a calculation of the Auamin molecule (shown in Figure 1) using the VTZ-F12 basis and the default settings.

29 ACS Paragon Plus Environment

Journal of Chemical Theory and Computation

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 30 of 74

However, this requirement is distributed over multiple nodes, and the usage on each node is quite moderate.

5.3 3-external Integrals The 3-external integrals occur in several places in the PNO-LCCSD-F12 amplitude equations. From the perspective of integral transformations, the most demanding term is Riabj = −

1 X 2 c∈[ik]

PNO

X

ik ik i j j j ¯j T¯ ac (c b |d k)td + ...,

d∈[ j j]PNO

a, b ∈ [i j]PNO

(33)

For the most accurate calculations, we assemble the integrals as written for each pair i j using the shared in-memory-intermediate (A|rbi j ) and fetching the required (A|sk) integrals from GAs on the fly. The computed integrals are transformed fully to PNO basis before being stored in GAs. However, even if we only include this term for strong pairs i j and restrict k ∈ [i j]strong in the summation over k, the number of 3-external integrals is normally much larger than that of 4external integrals. This creates significant CPU-time and memory bottlenecks. This problem can be largely alleviated by a projection approximation (denoted as project_je in ref 62): Riabj = −

1 X 2 c,d∈[i j]

ik i j i j i j T¯ ac (c b |d k)tdj + ...,

a, b ∈ [i j]PNO

(34)

PNO

The error introduced by this approximation is negligible and the number of 3-external integrals is significantly reduced. In addition, they can be computed from the shared-intermediate (A|ai j bi j ), which is a lot easier to compute and handle than the (A|ai j r) integrals without this projection approximations. The integral blocks (ci j bi j |di j k) are symmetric to the permutation of ci j and bi j , allowing them to be stored in packed format. Moreover, this type of integrals is also required by

30 ACS Paragon Plus Environment

Page 31 of 74

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Theory and Computation

other terms in the PNO-LCCSD-F12 amplitude equations or intermediates, for example ij V˜ ab =

X

−(mai j |bi j ci j )(im|Fˆ 12 | jci j ) + ...,

a, b ∈ [i j]PNO

(35)

c∈[i j]PNO

With the default settings, we approximate most terms involving 3-external integrals by projection approximations, in order to use the (ai j bi j |ci j k) type of integrals only. Detailed discussion on these approximations can be found in ref 62. In this scenario, the the 3-external integrals needed by different processors do not overlap and we can use disk-based algorithms as for the 4-external integrals, significantly reducing the memory requirement of the program.

5.4 2-external Coulomb Integrals The 2-external Coulomb integrals [of the type (ab|i j), to be referred to as the “J-integrals” hereafter] occur in the linear terms of the PNO-LCCSD-F12 amplitude equations, for example Riabj = −

1 X ¯ ik ik i j T ac (c b |k j) + ..., 2 c∈[ik]

a, b ∈ [i j]PNO

(36)

PNO

These terms are needed in our LCCSD residuals for all nondistant pair classes. 33,62 Without using projection approximations, we assemble these integrals from the shared-intermediate (A|ai j r) where r is restricted to a united PAO domain containing all [ik]PAO . The assembled integrals (ai j r|k j) are then transformed to the target integrals in the PNO space. The J-integrals have the permutation symmetry

(aik bi j | jk)i j = (ai j bik | jk)ik ,

(37)

where the subscript indicates the pair residual in which the integral is used. To keep the cost low, one can either use [ik]DF or [i j]DF as the fitting domain, which lifts the permutation symmetry shown above. In practice, we choose the fitting domain of the weaker pair among ik and i j to keep the permutation symmetry. The LMOs in the weaker pair tend to be more distant and therefore the 31 ACS Paragon Plus Environment

Journal of Chemical Theory and Computation

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 32 of 74

fitting domain is expected to cover more centers. This saves a factor of 2 in computing and storing the J-integrals. We found the difference of the two approaches in terms of correlation energy negligible: For VTZ-F12 calculations of the FH test set 152 the largest difference in the correlation energies was found to be 2 µEh , and in most cases the difference was well below 0.1 µEh . If the 3-external integrals are projected as shown in eq 34, all other integrals only require the intermediates (A|ai j r) with r ∈ [i j]PAO . Therefore, the added cost for the unprojected J-integrals are quite high. To alleviate this problem, we introduced a projection approximation by default (denoted project_j in ref 62): Riabj = −

1 X ¯ ik i j i j T ac (c b |k j) + ..., 2 c∈[i j]

a, b ∈ [i j]PNO

(38)

PNO

Despite projecting a linear term, this approximation works surprisingly well, provided that reasonably large PNO domains are used.

5.5 Other PNO-LCCSD Integrals The 2-external exchange integrals as well as the 1-external integrals are assembled in the PAO space from the pre-computed 3-index integrals as shown in Table 2. These integrals are relatively cheap to compute, except for the 2-external exchange-type integrals (ai j k|bmn l), which are needed in the non-linear terms of the PNO-LCCSD equations. The number of these integral blocks can be drastically reduced using projection approximations as proposed by Neese et al. 49,50,55 By default, we use these projections in all non-linear terms of the PNO-LCCSD equations, even though they can be optionally avoided. All these transformed integrals are stored in GAs. We used the fitting domain of [i j]DF for integrals of the type (iakl |bmn j), (kai j |lm), and (ik|l j). The permutation symmetries of the 1- and 2-external integrals are automatically kept with this choice of fitting domains and we exploited this symmetry to reduce the number of integrals. The 0-external (ik|l j) integrals are cheap to evaluate and store, so that we did not exploit the permutation

32 ACS Paragon Plus Environment

Page 33 of 74

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Theory and Computation

symmetry of i, k or l, j. This ensures consistent fitting domains for terms such as kl Riabj = (ik|l j)T¯ ab + ...,

a, b ∈ [i j]PNO

(39)

5.6 F12 Terms The density-fitting assembly of the F12-specific integrals takes a more complicated form since robust density fitting algorithms are needed. 78 The integrals required in the PNO-LCCSD-F12 method involve no more than two virtual indices and are evaluated in a similar way as in our PNO-LMP2-F12 program. 117 ij ij We compute the intermediates V˜ ab and V˜ ak (denoted V-terms in the following) using eqs 22

and 24 before the PNO-LCCSD-F12 iterations and store them in GAs. The 3-index F12 integrals are computed first and stored in GAs. This step is performed before the assembly of all 4-index integrals and parallelized over blocks of fitting functions. Subsequently, in a loop over pairs, the F12-specific integrals are assembled and stored in memory. After the 3- and 4-external integrals for the LCCSD calculations are computed, they are immediately contracted with the F12 integrals and the contributions are added to the V-terms, provided that sufficient memory is available. Otherwise the 3- or 4-external integrals are stored on local disks first, and subsequently loaded and block-wise contracted with the F12 integrals. The 2-external contributions to the V-terms are then evaluated. Integrals not needed for the LCCSD iterations are immediately discarded once the V-terms for a pair are computed. Using this strategy, the only additional requirements for the F12 terms in the PNO-LCCSD iterations are the trivial contributions to the residuals shown in eqs 13 and 12, and the storage of the F12 V-terms. Both the added CPU time and the GA usage are insignificant.

33 ACS Paragon Plus Environment

Journal of Chemical Theory and Computation

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

6 Benchmark Calculations 6.1 Computational Details In this section we show results of benchmark calculations with the PNO-LCCSD-F12 method. We primarily used the VDZ-F12 and VTZ-F12 basis sets of Peterson et al. 151 that are specially designed for F12 calculations. It has been found previously that diffuse functions in the basis sets are more important for the explicit correlation methods than for the conventional ones. 83,103,106,148 F12 calculations using the standard Dunning basis sets like cc-pVnZ provide rather small gains in the basis-set convergence in comparison with the standard calculations. The VnZ-F12 basis sets are specially optimized for F12 methods 151 and perform much better. They are more diffuse [their s, p bases correspond to the aug-cc-pV(n+1)Z bases] and yield significantly better Hartree-Fock energies as well. We recommend using the VTZ-F12 basis set in combination with our PNOLCCSD-F12 method when kJ mol−1 accuracy is desired. An exception are intermolecular interaction energies, for which diffuse polarization functions are important for an accurate description of the dispersion interaction. For such systems the aug-cc-pVnZ basis sets seem more appropriate. 125,156 To avoid basis-set linear-dependence in some systems, we also used a basis set denoted as aVTZ, which consists of the cc-pVTZ basis sets 157 for H atoms, and the aug-cc-pVTZ basis sets 158 for C, N, and O atoms. Due to the presence of diffuse polarization functions, this basis set is also well suited for calculations of intermolecular interactions. 125 For all transition metal atoms, the Stuttgart effective core potentials ECPnMDF along with the the aug-cc-pVDZ-PP and aug-ccpVTZ-PP basis sets 159–162 were used. The number of core electrons n in the ECP is 10, 28, and 60, respectively, for the 1st (Cu), 2nd (Ru, Pd, Ag), and 3rd (Pt, Au) block transition metals. For computing the Fock matrix the Molpro aug-cc-pVTZ/JKFIT basis was used. This basis was derived from the cc-pVTZ/JKFIT basis set 146 by adding for each angular momentum another shell of diffuse functions in an even tempered manner. Furthermore, we found that the fitting errors of the Hartree-Fock energies obtained with this basis were still large in cases where sulfur atoms are present. For example, the error amounts to 0.28 kcal mol−1 for the reaction energy of reaction 34 ACS Paragon Plus Environment

Page 34 of 74

Page 35 of 74

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Theory and Computation

4 of the FH test set (Section 6.4). We found that this problem is avoided if one tight f -function −2 (exponent 11.0 a−2 0 ) and one g-function (exponent 1.6 a0 ) are added to the JKFIT basis for S. The

error is then reduced by one-order of magnitude to 0.03 kcal mol−1 . Probably, this is related to the fact that the cc-pVTZ/JKFIT was optimized for the cc-pVTZ orbital basis, but it is well known that an additional tight d function in the orbital basis is required to obtain accurate HF energies [yielding the cc-pV(n+d)Z basis sets]. 163 Such tight functions are also present in the VTZ-F12 basis set. We also tested the effect of additional tight f and g functions in the JKFIT basis for phosphorous. This increased the HF reaction energy of the by Auamin reaction (Figure 1) by about 0.09 kcal/mol. In view of this rather small effect we did not repeat the PNO-LCCSD-F12 calculations with the larger fitting basis (the problem was detected only after all calculations had been completed). Similarly, additional f and g functions on chlorine atoms had a negligible effect on chlorine-containing reactions of the FH benchmark set. For all other integral types the aug-cc-pVTZ/MP2FIT basis set 164 was used. The aug-ccpVTZ/JKFIT basis set was employed as RI basis. The CABS singles correction 90,103 of the Hartree-Fock energies was included in all F12 calculations. Unless otherwise noted, the tripleζ auxiliary basis sets were also used for double-ζ calculations. From our preliminary tests, the aug-cc-pVDZ/JKFIT basis set is inadequate for RI when sub-kcal mol−1 accuracy is desired. For some test systems we performed canonical MP2-F12 and CCSD-F12 calculations. Unless we add the prefix “DF-” to the name of the methods, DF approximations are not used for these canonical calculations, except for the evaluation of F12-specific integrals. Details about the canonical MP2-F12 and CCSD-F12x methods can be found in previous publications. 88,90,103,108,165 In all F12 calculations, fixed amplitudes 79,80 and a geminal exponent of 1 a−1 0 have been used. All benchmark calculations were performed using the development version 166 of Molpro 167 compiled with the GNU compiler collections 6.2.0. The Intel Math Kernel Library 11.2, OpenMPI 2.0.2, and GlobalArrays 5.5 were used. Each node on our computer cluster has 2 Xeon E5-2680 v2 @ 2.8 GHz processors with a total of 20 physical cores and 256 GB of memory. The nodes are connected by a QDR-Infiniband network. All elapsed times reported in this work were obtained

35 ACS Paragon Plus Environment

Journal of Chemical Theory and Computation

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 36 of 74

using 20 cores per node, and one MPI process per core was created.

6.2 General Performance Table 3: Reaction Energies in kcal mol−1 Obtained with PNO-LCCSD-F12 Calculations Reaction Isomer4 Coronene dimera Auamin Androstendione Organocatalysis ∆ETS1 Organocatalysis ∆ETS2 Organocatalysis ∆Er a b c

Basis HF+CABS LMP2 LMP2-F12 LCCSD LCCSD+∆F12 LCCSD-F12a LCCSD-F12b VDZ-F12 19.07 79.63 80.69 64.79 65.47 65.05 63.90 VTZ-F12 18.69 78.55 79.83 62.91 64.19 63.81 63.51 VDZ-F12 −15.98 33.38 35.95 10.06 13.03 12.66 11.99 VDZ-F12b −15.97 34.48 36.52 11.22 13.13 12.94 12.30 aVTZc −15.97 35.38 36.92 11.80 13.22 12.97 12.59 VDZ-F12 22.01 53.82 60.47 38.88 45.42 45.20 44.38 VTZ-F12 22.03 57.25 60.27 42.26 45.20 45.01 44.57 VDZ-F12 −1.30 4.26 4.58 4.59 4.82 4.77 4.47 VTZ-F12 −1.38 4.35 4.45 4.56 4.59 4.55 4.46 VDZ-F12 16.79 13.59 13.21 14.12 13.74 13.80 13.93 VTZ-F12 16.82 13.20 13.30 13.88 13.98 13.98 13.99 VDZ-F12 14.28 10.06 9.51 11.11 10.56 10.68 10.91 VTZ-F12 14.32 9.65 9.42 10.82 10.58 10.65 10.74 VDZ-F12 −14.33 −21.91 −23.74 −19.40 −21.22 −21.08 −20.95 VTZ-F12 −14.32 −23.27 −23.73 −20.70 −21.15 −21.08 −21.05

Dissociation energies with counterpoise correction. Even tempered d-type diffuse functions are added to carbon atoms. We were unable to perform VTZ-F12 calculations due to basis-set linear dependence.

Table 3 shows the reaction energies computed with various methods and basis sets for a few large molecular systems described in the preceding paper. 62 In all cases the reaction energy differences from double-ζ to triple-ζ basis sets are reduced when F12 terms are included, both at the LMP2 and LCCSD levels. The basis set effects are somewhat smaller for LCCSD-F12 as compared with LMP2-F12. Among the F12 treatments, the most rigorous one F12b tends to behave best in correcting the basis set errors. In general, already the VDZ-F12 basis set yields quite accurate results. The largest error compared with triple-ζ basis sets is found for the coronene dimer and amounts to 0.6 kcal mol−1 . However, this is a special case, since the VDZ-F12 basis appears not to be diffuse enough to well describe the dispersion interaction. If a diffuse d-function is added to each carbon atom, the agreement becomes much better and the difference to the aVTZ results is reduced to 0.3 kcal mol−1 . Without F12 terms, the Auamin reaction energies differ by 3.4 kcal mol−1 for the two basis sets, while with the F12 methods the difference is only 0.2 kcal mol−1 . 36 ACS Paragon Plus Environment

Page 37 of 74

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Theory and Computation

We also performed the Auamin calculation with the VTZ-F12 basis using larger domains and CC CC tighter thresholds (project_j = close, project_ke = 1, T en = 0.997, T PNO = 10−8 , T close =

10−5 , T weak = 10−6 ). The resulting PNO-LCCSD-F12b reaction energy of 44.52 kcal mol−1 is in excellent agreement with the result of 44.57 kcal mol−1 obtained with default settings. This value can be compared to the zero-point energy (ZPE) corrected experimental value 149 of 47.0 ± 2.7 kcal mol−1 obtained from collision-induced dissociation (L-CID) measurements. 168 Of course, triple excitations are still missing in the calculations, and including them might further improve the agreement. Earlier studies for canonical methods have shown that the F12a approximation yields some favorable error cancellation when small basis sets (e.g. VDZ-F12) are used. The more rigorous F12b treatment was preferred when larger basis sets are used. 103,148 With augmented triple-ζ basis sets, the two methods were found to be approximately of similar quality. This assessment was based on the accurate Ansatz 3C in the F12 theory. 77,88 In the present work, we used the simpler Ansatz 3*A, which is known to slightly overestimate the correlation energy. Therefore, we recommend to use F12b even for smaller basis, given that F12a also systematically overshoots the correlation energy. 90,103,147 The simple ∆F12 approximation gives surprisingly good results in some cases, but with our current program it yields only relative small (∼ 10%) overall savings of the computational time. We therefore recommend using the more rigorous PNO-LCCSD-F12b treatment for all calculations. Unless otherwise noted, the PNO-LCCSD-F12 results reported in the following sections were obtained with F12b. There have been publications on local MP2 and CCSD methods using the formally more accurate Ansätze 3B or 3C. 118,119,121,122,129,131 These treatments are usually preferred in canonical calculations where domain approximations do not apply, although Ansatz 3*A may bring favorable error cancellations that lead to better basis-set convergence in some cases. However, in the local cases such methods have several disadvantages: First, the Fock matrix in the RI basis is needed, which cannot easily be truncated to yield linear scaling, without introducing noticeable errors. Second, higher levels of theory require double RI approximations, where the truncation of the RI space

37 ACS Paragon Plus Environment

Journal of Chemical Theory and Computation

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

introduces additional errors, in particular when local RI approximations are used. In contrast, approximation 3*A does not require any RIs between the Fock matrix and the F12 correlation factor, and it is linear in the RI basis. In sections 6.5 and 6.6 we will compare results obtained with the canonical MP2-F12/3C and MP2-F12/3*A methods, and it is shown that the differences are quite small, in particular when triple-ζ basis sets are used. Table 4: Elapsed Times and GA Usages of PNO-LCCSD(-F12) Calculations Using the VTZ-F12 Basis Set and Default Options number of molecule atoms AO functions nodes Isomer4 reactant 81 2543 3 Isomer4 product 81 2543 3 Coronene dimer 72 2544 4 AuAmin 92 3345 4 (gly)32 227 7306 4 a b c d

LMP2 1.8 1.0 3.4 3.7 5.8

elapsed time / min GA usage / gigaworda b c d LMP2-F12 LCCSD LCCSD-F12 F12 extra LCCSD LCCSD-F12 F12 extra 13.9 72.0 91.6 27% 17.9 18.8 5% 5.5 19.1 26.5 39% 8.6 9.0 5% 19.2 135.5 175.5 30% 29.0 29.2 1% 23.2 129.1 177.3 37% 36.1 37.5 4% 15.5 65.5 81.7 25% 17.7 18.8 6%

1 gigaword = 8 × 109 bytes. Including the preceding LMP2 calculation. Including the preceding LMP2-F12 calculation. Extra cost of LCCSD-F12 calculations compared to the corresponding LCCSD calculations.

Table 4 shows a comparison of timings and GA usages for the PNO-LCCSD and PNO-LCCSDF12 methods (excluding the preceding Hartree-Fock calculations). The PNO-LCCSD-F12 calculations typically take only 25–40% more elapsed time than the corresponding PNO-LCCSD calculations. The majority of the additional cost comes from PNO-LMP2-F12 relative to PNO-LMP2. The additional effort in the LCCSD-F12 part is small, since the F12 terms do not require any extra 3- and 4-external integrals. The increase in computation time tends to be more noticeable when the molecule is relatively small or less spatially compact. We note that, similar to PNO-LCCSD, 62 the PNO-LCCSD-F12 program can run on a single node, exclusively using shared disk files to avoid GA usage entirely. Using this option it is possible to perform most of the calculations reported in this work on one node using 4–6 GB of memory per CPU core.

6.3 Illustration of Scalability In this section we demonstrate the scaling of the PNO-LCCSD-F12 method with respect to the number of correlated electrons for linear glycine polypeptides (gly)n . The results are shown in 38 ACS Paragon Plus Environment

Page 38 of 74

Page 39 of 74

4000

(a) comparison of LCCSD and LCCSD-F12

3000 PNO-LCCSD(-F12) total

integrals

2000

1000 iterations elapsed time / s

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Theory and Computation

0

4

8

12

16

20

24

28

32

1200 (b) integral transformations 1000

3-, 4-external and J-integrals

800 3-index PAO integrals

600

F12 specific integrals

400 200 other integrals 0

4

8

12 24 16 20 number of glycine units n

28

32

Figure 11: Elapsed time of various steps in the PNO-LCCSD-F12 (solid lines) and PNO-LCCSD (dashed lines) calculations for the glycine polypeptides (gly)n . The calculations were performed with 80 cores on 4 nodes using the aVTZ basis set and default parameters as described in the text. 10 iterations were performed for each calculation. The actual number of iterations to converge the energy to within 10−6 Eh varies between 9 and 13. Figure 11, which is directly comparable to the corresponding study for PNO-LCCSD published in the preceding paper. 62 The dotted and full lines show the elapsed times of calculations without and with the F12 terms, respectively. From the top panel, we can see that the scaling behavior is very similar with or without F12 terms, and the increased cost of these terms comes almost exclusively from the integral transformations. During the iterations, the F12 terms only add trivial updates to the residuals as shown in eqs 12 and 13, and the additional cost is invisible on the scale of Figure 11. The bottom panel of Figure 11 shows step-wise scaling curves of the PNO-LCCSD-F12 integrals. The red line in the panel shows all extra cost for the F12 terms. This includes the 3-index F12 integrals, the assembly of the extra 0–2 external integrals, and the computation of the V terms.

39 ACS Paragon Plus Environment

Journal of Chemical Theory and Computation

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

These steps show nearly perfect linear scaling with the molecular size and cost significantly less than the assembly of the 3-external, 4-external, and the J-integrals. It is interesting to see that the latter steps are faster in a PNO-LCCSD-F12 calculation than in the corresponding calculation without the F12 terms (dotted lines), while exactly the same quantities are computed. This is likely due to better load-balancing, since the F12 terms add a noticeable amount of CPU time, while the communication requirements are rather low. Since different processors can be working on different steps, the added CPU time leads to a less congested network and hence more efficient integral transformations. The scaling with the number of processors of the PNO-LCCSD-F12 method is very similar to that reported for the PNO-LCCSD method. 62 The parallel efficiency is excellent on a single node and starts to drop slightly on multiple nodes due to the unavoidable communication between nodes. The parallelization is also limited by the number of pairs (in the assembly step) and the number of fitting blocks in the calculation and transformation of the 3-index integrals, which is proportional to the number of centers. For example, in a PNO-LCCSD-F12 calculation of the coronene dimer with the VTZ-F12 basis set, we observed a speed-up of 2.5 when using 120 cores (83% efficiency) relative to the timings obtained with 40 cores. When using 240 cores, the speed up was 4.0 (67% efficiency). The speed-up is identical to those reported for the PNO-LCCSD method (cf. Figure 12 of ref 62). These timings do not include the preceding PNO-LMP2-F12 calculation, which has a slightly better parallel efficiency due to its lower communication demands. 117

6.4 FH Test Set In Table 5 we present some statistical data for the FH test set. 152 Some of the underlying calculations and results for the individual reactions have already been described in section 3.3. The first block of the table shows the local errors of the PNO methods with respect to the corresponding canonical calculations. Comparing the results obtained from different basis sets, we observe that the local errors are slightly larger for the smaller VDZ-F12 basis, likely due to the smaller PAO domains of the external space. The difference is significantly reduced when F12 terms are included. 40 ACS Paragon Plus Environment

Page 40 of 74

Page 41 of 74

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Theory and Computation

The local errors are slightly reduced from LMP2 to LCCSD, which is probably due to some fortuitous error cancellations. On the other hand, the errors of the LCCSD-F12 method are somewhat larger than those of the LMP2-F12 method, since we have used more drastic approximations at the LCCSD-F12 level to minimize the cost (cf. Section 3.3). Table 5: Error Statistics in kcal mol−1 for the Reaction Energies of the FH Test Set Obtained with Various Methods method PNO-LMP2 PNO-LMP2-F12 PNO-LCCSD PNO-LCCSDd PNO-LCCSDe PNO-LCCSD-F12 PNO-LCCSD-F12d PNO-LCCSD-F12e

VDZ-F12 reference MAXa RMSb MADc DF-MP2 1.21 0.28 0.19 DF-MP2-F12 0.19 0.06 0.04 DF-CCSD 1.13 0.27 0.17 DF-CCSD 1.13 0.27 0.18 DF-CCSD 1.12 0.27 0.18 f DF-CCSD-F12 0.25 0.10 0.07 f DF-CCSD-F12 0.25 0.08 0.05 DF-CCSD-F12f 0.20 0.07 0.05

DF-HF DF-MP2 DF-CCSD

HF MP2 CCSD

0.03 0.01 0.03 0.01 0.03 0.01

0.01 0.01 0.01

0.04 0.02 0.05 0.01 0.04 0.01

0.01 0.01 0.01

PNO-LMP2 PNO-LMP2-F12 PNO-LCCSD PNO-LCCSD-F12

MP2 MP2-F12 CCSD CCSD-F12

1.20 0.18 1.13 0.27

0.18 0.04 0.16 0.08

0.78 0.17 0.54 0.22

0.15 0.03 0.12 0.06

a b c d

e f

0.27 0.06 0.26 0.10

VTZ-F12 MAXa RMSb MADc 0.79 0.21 0.15 0.16 0.05 0.03 0.54 0.16 0.12 0.58 0.17 0.13 0.57 0.17 0.13 0.23 0.08 0.06 0.19 0.07 0.05 0.18 0.07 0.05

0.21 0.05 0.16 0.08

Maximum deviation. Root-mean-squares deviation. Mean average deviation. CC CC With large domains (T PNO = 10−8 , T en = 0.997, and T val = 10−5 Eh ) and reduced projections (project_j=weak and project_ke=0). Same as d, but without any pair approximations. Estimated from ECCSD-F12b − ECCSD + EDF-CCSD .

Another possible source of error is the absence of density fitting in the CCSD-F12 reference calculations. Presently we do not have an efficient DF-CCSD-F12 program using exactly the same ansatz as in the local calculations. The reference values were therefore obtained by adding the F12 correction (ECCSD-F12 − ECCSD ) to the DF-CCSD results. Since the F12 terms in the CCSD-F12 method are computed using density fitting, this should be a very good approximation. The RMS 41 ACS Paragon Plus Environment

Journal of Chemical Theory and Computation

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

error of the PNO-LCCSD-F12 method is below 0.1 kcal mol−1 . The results obtained without any pair approximations are also shown in Table 5, but the error statistics are almost identical to those with default pair approximations. This is partly due to our tight default pair thresholds that give rather converged results (cf. Figure 3), and partly due to the relatively small size of the molecules in the test set. We now consider the error from density fitting itself. The second block of Table 5 shows the fitting errors of the canonical DF-HF, DF-MP2, and DF-CCSD methods. The calculations for both AO basis sets used the same triple-ζ fitting bases. It was found previously 32 that the fitting errors of correlation energies in PAO-LCCSD calculations are significantly larger than those in PAO-LMP2 calculations, unless a fitting basis with one cardinal number higher than the orbital basis was used. However, these errors canceled largely out when relative energies were computed. Consistent with these findings we do not see this effect on the relative energies of the FH test set, and the fitting errors are much smaller than the other sources of errors. In the last block of Table 5 we present the deviation of the PNO methods from the standard canonical methods, in which density fitting is used for F12-specific integrals only. The MAX and RMS errors are found to be very similar to the corresponding ones relative to the canonical methods with density fitting. The individual results can be found in the Supporting Information.

6.5 Isomerization Energies In order to test our method for systems that are known to be significantly affected by correlation effects and long-range dispersion forces we used the ISOL24 benchmark set of Huenerbein et al. 150 This set contains 24 real-life isomerization reactions of molecules with 24–81 atoms, which have been used in the past to test the performance of DFT methods with numerous different functionals. 150,169,170 The original reference values were computed by SCS-LMP2 and SCS-LMP3 using basis set extrapolations based on TZ2P(2d,2p) and QZ3P(3d2f, 3p2d) basis sets. 150 One year later, the benchmark set has been refined by Luo et al. 169 For the 6 smallest systems, CCSD(T)F12a/aug-cc-pVDZ calculations were carried out, while for the others a hybrid method denoted 42 ACS Paragon Plus Environment

Page 42 of 74

Page 43 of 74

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Theory and Computation

MCQCISD-MPW 171 was used. The subset of 6 reactions, which is sometimes denoted ISO6/11, was also used by Liakos and Neese to benchmark their DLPNO-CCSD(T) method. 172 Table 6: CCSD-F12x and PNO-LCCSD-F12x Isomerization Energies (in kcal mol−1 ) for the ISOL6/11 Benchmark Set 150,169 Using the aug-cc-pVDZ Basis Seta No. 3 9 10 13 14 20 a

HF MP2 LMP2 +CABS -F12/3*A -F12/3*A 12.19 10.80 10.74 19.42 22.24 22.34 0.34 9.05 9.10 29.86 36.91 36.96 −0.25 5.62 5.57 5.02 5.96 5.89

SCS-MP2 SCS-LMP2 CCSD- LCCSDCCSD- LCCSD- CCSD-F12/3*A -F12a/3*A F12a/3*A F12a/3*A F12b/3*A F12b/3*A F12b/3C 9.99 9.91 11.36 11.33 11.23 11.18 10.94 20.95 21.08 20.66 20.77 20.73 20.80 21.23 7.03 7.09 5.50 5.61 5.52 5.60 5.87 35.18 35.25 32.50 32.59 32.55 32.59 32.83 4.24 4.19 4.15 4.20 4.03 4.09 4.11 5.45 5.36 5.35 5.37 5.49 5.51 5.23

The aug-cc-pVTZ/JKFIT basis set was used for RI and computing the Fock matrix, and the aug-cc-pVTZ/MP2FIT basis for all other integrals.

In Table 6 we compare our PNO-LMP2 and PNO-LCCSD-F12 calculations for the 6 smallest systems with the corresponding canonical counterparts, using the aug-cc-pVDZ basis set. We have recomputed all canonical MP2-F12, SCS-MP2-F12, and CCSD-F12 values, using aug-cc-pVTZ fitting and RI basis sets. Luo et al. 169 apparently used the cc-pVDZ/JKFIT basis for the RI and the Fock matrix, which is the current default in Molpro for the aug-cc-pVDZ orbital basis. However, we found that in some cases the use of larger auxiliary basis sets has a significant effect on the results (up to 0.3 kcal/mol for reaction 20, which contains fluorine atoms). This is mainly due to missing diffuse functions in the cc-pVDZ/JKFIT basis. The canonical and local results are in outstanding agreement, the differences between the canonical and local CCSD-F12b values are always smaller than 0.1 kcal mol−1 . It should be noted that slight differences can arise from the fact that the canonical values were computed without density fitting (except for the F12 integrals and the CABS singles correction), while in the local methods the Fock matrix and all integrals were density fitted. Since the local LMP2 and LCCSD calculations in our program are always based on the MP2F12/3*A approximation (using fixed amplitudes and a geminal exponent of 1.0 a−1 0 ), we computed the canonical CCSD-F12b/3*A as well as the CCSD-F12b/3C values for comparison. The differences are quite large, up to 0.5 kcal mol−1 for reaction 9. For the VTZ-F12 or aug-cc-pVTZ orbital 43 ACS Paragon Plus Environment

Journal of Chemical Theory and Computation

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

basis sets, which have been used for most other calculations in this work, the differences are significantly smaller (see below). Somewhat unexpectedly, the difference between the aug-cc-pVDZ results in Table 6 and the corresponding VTZ-F12 results in Table 7 are mostly smaller for the 3*A approximation than for the 3C approximation, and F12b yields most consistent results. Overall, the various F12 approximations (3C vs 3*A and F12a vs F12b) have in most cases a larger effect than the local approximations. Unfortunately, it is not possible for the larger molecules to compare with reliable extrapolated complete basis set (CBS) results, and therefore the overall accuracy relative to CCSD/CBS results is difficult to estimate. However, the fact that the differences between PNO-LCCSD-F12b calculations with double-ζ and triple-ζ basis sets are small indicates that the triple-ζ results should be well converged. Some comparisons of reaction energies obtained by PNO-LMP2-F12/VTZ-F12 and [45]-extrapolated DF-MP2 values can be found in ref 33, and the agreement was outstanding. In Table 7 we present results for the whole ISOL24 benchmark set using the VTZ-F12 basis set. The reference values of Huenerbein et al. 150 and Luo et al. 169 are included for comparison, but since these have been computed with much smaller basis sets and at lower computational levels without F12 corrections, they cannot be used as a reference for testing the local methods. Our PNOSCS-LMP2 values can on the one hand be compared with corresponding canonical SCS-MP2-F12 calculations that we performed, and on the other hand with the extrapolated “CBS” values of Huenerbein et al. The agreement of our local and canonical SCS-MP2-F12/3*A results is excellent. In most cases the differences of the isomerization energies is within 0.1 kcal mol−1 ; only in one case (system 4) it is 0.25 kcal mol−1 and in one other case (system 19) it amounts to 0.18 kcal mol−1 . This demonstrates that the domain errors in our local calculations are negligible. They are smaller than the differences between the MP2-F12/3C and MP2-F12/3*A approximations (which are also quite small, the largest difference being 0.6 kcal mol−1 for system 2). Normally, the domain errors are very similar at the LMP2-F12 and LCCSD-F12 levels, and we can therefore safely assume that the PNO-LCCSD-F12 values are also converged to within a few tenth of a kcal mol−1 . The pair and projection approximations can only be tested by checking the convergence with tighter

44 ACS Paragon Plus Environment

Page 44 of 74

Page 45 of 74

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Theory and Computation

Table 7: PNO-LMP2-F12 and PNO-LCCSD-F12 Isomerization Energies in kcal mol−1 for the ISOL24 Benchmark Set 150 using the VTZ-F12 Basis No. 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 b c d e f

HF +CABS 74.41 19.61 12.28 18.69 22.73 23.51 7.22 32.44 19.29 0.15 51.51 0.08 29.76 −0.50 −3.83 19.25 4.47 21.71 19.24 5.19 5.53 −13.19 17.39 12.87

MP2 SCS-MP2 SCS-MP2 -F12/3C -F12/3C -F12/3*A —f —f 71.43 43.63 37.51 36.87 10.68 9.89 9.79 —f —f 64.72 37.56 33.41 33.11 28.01 24.11 24.04 20.52 17.14 16.77 22.93 23.75 24.24 22.38 21.18 21.04 9.14 7.13 7.03 33.58 37.01 37.14 0.78 0.86 0.88 36.96 35.21 35.17 5.50 4.12 4.02 −0.66 2.11 1.91 27.74 22.63 22.70 12.60 10.07 9.84 26.64 26.67 26.52 16.32 17.56 17.65 5.51 5.06 5.05 11.75 10.94 10.50 −0.92 −2.87 −3.24 29.17 26.41 26.24 15.87 15.34 15.25

SCS-LMP2 -F12/3*A 71.48 36.94 9.82 64.97 33.10 23.99 16.85 24.23 21.05 7.07 37.14 0.87 35.21 4.07 1.81 22.70 9.81 26.56 17.83 5.09 10.45 −3.15 26.15 15.31

SCS-MP2b 69.3 41.8 10.9 77.4 35.9 26.2 19.3 18.8 22.8 8.5 36.8 0.9 35.9 5.3 3.1 23.0 11.4 27.2 16.5 4.8 13.3 0.0 28.2 16.1

PNO-LCCSD-F12b/3*A MCQCISDdefaults accuratec accurated MPWe 70.64 70.61 70.81 69.17 36.85 36.84 36.76 37.54 11.31 11.32 11.32 9.77 63.51 63.50 63.32 66.43 33.09 33.11 33.04 32.84 25.06 25.02 25.00 25.51 16.93 16.93 16.87 17.37 20.15 20.19 20.28 22.34 20.88 20.89 20.84 21.76 5.52 5.53 5.48 6.82 41.98 41.92 42.14 37.87 0.31 0.31 0.32 0.20 32.63 32.65 32.65 33.52 3.97 3.97 3.90 5.30 2.87 2.88 2.77 3.06 21.79 21.80 21.90 22.78 11.00 11.00 10.89 10.33 24.69 24.69 24.57 22.57 17.72 17.73 17.77 18.25 5.17 5.14 5.16 4.66 12.05 12.06 11.90 11.21 −1.19 −1.21 −1.20 0.77 24.93 24.91 25.01 23.43 15.25 15.27 15.20 14.94

From ref 150. With T close = 10−5 Eh , project_j=weak, and project_ke=2. With T PNO = 10−8 , T en = 0.997, project_j=weak, and project_ke=2. From ref 169. We were unable to perform these calculations due to disk space requirements.

pair thresholds and less aggressive projections. For this purpose we performed PNO-LCCSD-F12 calculations with tighter thresholds, and found that the results agree well within 0.1 kcal mol−1 (RMS error 0.02 kcal mol−1 ). We also performed calculations with larger PNO domains and found the difference is within 0.1–0.2 kcal mol−1 (RMS error 0.10 kcal mol−1 ) of the results obtained from default settings. We therefore expect the PNO-LCCSD-F12b results with the default parameters to be very close to the canonical counterparts. Moreover, the PNO-LCCSD-F12b results for the 6 smallest systems agree well with the aug-cc-pVDZ values in Table 6. This indicates that basis set incompleteness errors are small as well. In fact, benchmarks for small molecules have shown that with the CCSD-F12b method better than quintuple-ζ accuracy is achieved with the VTZ-F12 basis 45 ACS Paragon Plus Environment

Journal of Chemical Theory and Computation

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

set. The agreement of our SCS-MP2-F12 and SCS-LMP2-F12 results with the extrapolated CBS values of Huenerbein et al. 150 is at best fair, in some cases poor. The worst case is system 4, where the results differ by more than 4 kcal mol−1 . This demonstrates that the basis sets used in ref 150 are too small, and their error estimates are somewhat too optimistic. Also the MCQCISD-MPW reference values reported in ref 169 show quite large differences to our PNO-LCCSD-F12 results, with an RMS deviation of 1.50 kcal mol−1 and the largest deviation being 4.11 kcal mol−1 . We believe that our results should be much closer to canonical CCSD-F12 (or CCSD/CBS) results than the MCQCISD-MPW values. But the still missing triple excitations would be needed to obtain reference values with < 1 kcal mol−1 accuracy.

6.6 Intermolecular Interactions Table 8: Comparison of Canonical CCSD-F12 and Local PNO-LCCSD-F12 Calculations for Intermolecular Interaction Energiesa HF molecule Benzene dimer (stacked) Benzene dimer (T-shaped) Uracil dimer (H-bonded) Uracil dimer (stacked) Pyrazine dimer a b c d e f

5.36 1.51 −16.29 0.06 4.16

3C −4.93 −3.62 −20.39 −11.09 −6.88

MP2-F12 3*A loc. −4.95 −4.89 −3.64 −3.61 −20.47 −20.45 −11.14 −11.06 −6.91 −6.86

CCSD-F12a can.b loc.c −1.21 −1.18 −2.01 −2.00 −19.73 −19.72 -7.93 −7.86 −2.73 −2.68

CCSD-F12b can.b loc.c loc.d loc.e loc.f −1.12 −1.06 −1.08 −1.06 −1.11 −1.96 −1.94 −1.95 −1.93 −1.96 −19.65 −19.60 −19.61 −19.60 −19.65 -7.81 −7.71 −7.70 −7.69 −7.79 −2.62 −2.54 −2.56 −2.55 −2.61

In kcal mol−1 , using the aug-cc-pVTZ basis set. All values are counterpoise corrected and include the CABS singles correction. All CCSD-F12 calculations employ the MP2-F12/3*A approximation. Using default settings. As c, but with reduced projections: project_j=weak and project_ke=0. As d, but T close = 10−5 Eh . CC = 10−9 , and T = T CC = 0.997. As e, but tight domain thresholds: IEXT=3, T OSV = 10−10 , T PNO = T PNO en en

In this section we apply our method to some dimers from the S22 benchmark set 173 for which canonical CCSD-F12x values are available from previous work in our group. 124 Table 8 shows a comparison of interaction energies computed with the canonical and local CCSD-F12x methods for four of the larger dimers, which could still be treated with the canonical CCSD-F12 method. Default thresholds have been used in the PNO-LCCSD-F12 calculations. The basis sets and methods are essentially the same as those used in ref 124. The only difference is that in the previous calculations the rigorous MP2-F12/3C approximation has been applied, while here we used the 46 ACS Paragon Plus Environment

Page 46 of 74

Page 47 of 74

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Theory and Computation

more approximate and simpler MP2-F12/3*A approximation. To demonstrate the difference of these approximations, we present in Table 8 canonical MP2-F12/3*A as well as the MP2-F12/3C results. For the hydrogen-bonded uracil dimer, the difference in the interaction energy amounts to 0.08 kcal mol−1 out of a total correlation contribution of ∼ 4 kcal mol−1 . In the other four dispersion-dominated cases, it amounts to 0.05 kcal mol−1 or less. This difference is nearly the same for MP2-F12 and CCSD-F12. The agreement of the local and the corresponding canonical results is excellent. With the default settings, the largest difference amounts to 0.10 kcal mol−1 for the stacked uracil dimer. We further show the results with reduced projections, tighter pair and domain thresholds, and then the interaction energies match even better with the canonical values. It is anticipated that the pair approximations affect the interaction energies. We see that these are reduced by 0.01–0.02 kcal mol−1 when using a tighter pair threshold, which slightly deteriorates the agreement with the canonical results. However, near-perfect agreement with canonical results can be obtained if the domain thresholds are also tightened. This shows that with the default parameters there is a slight error compensation between pair and domain approximations. However, all these effects are below 0.1 kcal mol−1 and normally negligible. Other sources of error, in particular the omission of triple excitations, are certainly much larger. We now consider the coronene dimer from the L7 benchmark of Sedlak et al. 174 Obviously, due to its delocalized electronic structure this is a rather difficult system for local correlation methods. The authors presented an interaction energy of 24.36 kcal mol−1 . This value results from an MP2/CBS estimate obtained by extrapolating MP2/aug-cc-pVDZ and MP2/aug-cc-pVTZ energies, plus a QCISD(T) correction EQCISD(T) − EMP2 , computed with a mixed aug-cc-pVDZ and cc-pVDZ basis. In Table 9 we present results obtained with canonical MP2 and SCS-MP2 and their local counterparts, as well as PNO-LCCSD and PNO-LCCSD-F12b results (using default values) as a function of the PAO domain size parameter IEXT. We obtained 38.1 kcal/mol using DF-MP2F12/aVTZ [counterpoise (CP) corrected], which is probably closer to the MP2/CBS value than the extrapolated value of 39.0 kcal/mol of Sedlak et al. The corresponding result with the cc-pVDZ-

47 ACS Paragon Plus Environment

Journal of Chemical Theory and Computation

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

F12 basis is 37.4 kcal/mol. Very recently, Pavoševi´c et al. 131 carried out DLPNO-CCSD(T)-F12 calculations, and their best result (not CP corrected) is 19.14 kcal/mol, using the VDZ-F12 basis and “verytight” program options. Their triples correction amounts to 5.25 kcal mol−1 , yielding a DLPNO-CCSD-F12 value of 13.9 kcal mol−1 . This can be compared to our PNO-LCCSD-F12 results of 12.8 kcal mol−1 (VDZ-F12), or 12.9 kcal mol−1 (aVTZ) obtained with IEXT=3 (Table 9). We believe that the remaining difference of about 1 kcal mol−1 is likely outside our error bounds for the local approximations. A difference could still arise from the different F12 approximations (MP2-F12/3*A vs MP2-F12/3C and CCSD-F12b vs CCSDF12 ). Unfortunately, we were unable to compare canonical MP2-F12/3C and MP2-F12/3*A calculations for the coronene dimer, due to excessive disk space requirements of the MP2-F12/3C calculation. For the benzene dimer, however, the difference of the MP2-F12/3C and MP2-F12/3*A interaction energies is only 0.02 kcal/mol using the aug-cc-pVTZ basis set (Table 8). The interaction energy of the coronene dimer obtained with PNO-LCCSD-F12a amounts to 13.41 (12.92) kcal/mol (CP corrected result in parenthesis), and these values are larger than the corresponding PNO-LCCSD-F12b values of 12.76 (12.25) kcal/mol by 0.65 (0.67) kcal/mol. For the larger aVTZ basis the differences between F12a and F12b results are reduced to 0.05 (0.35) kcal/mol. The CP corrections are also significantly smaller for the aVTZ basis set. Therefore an error of 1 kcal/mol caused by the CCSD-F12b approximation seems unlikely. In any case it appears that the CBS estimate of Sedlak et al. is too high. This might be due to the strong overshooting of MP2, which overestimates the correlation contribution to the binding energy by about 50%. In Table 9 it can be seen that the CP corrections in the local calculations are much smaller than in the canonical ones. As expected, they increase with increasing PAO domain size. The two coronene molecules are separated by about 6.5 a0 . Since our primary domains are extended by REXT = (2 × IEXT + 1) a0 , the PAO domains are located only on each individual monomer for IEXT=1 and 2, but extend over both molecules for IEXT=3 (REXT=7 a0 ). This explains the steep increase of the BSSE from IEXT=2 to IEXT=3 by nearly a factor of 2. Accordingly, the CPuncorrected local results (without F12 corrections) converge very slowly to the canonical values

48 ACS Paragon Plus Environment

Page 48 of 74

Page 49 of 74

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Theory and Computation

Table 9: Computed Dissociation Energies (in kcal mol−1 ) of the Coronene Dimer for Various PAO Domain Sizesa Method HF DF-MP2 DF-SCS-MP2

VDZ-F12 IEXT ∆E ∆ECP b ∆CPc −14.6 −16.1 −1.5 45.6 35.4 −10.2 34.4 24.0 −10.4

DF-MP2-F12 DF-SCS-MP2-F12

∆E −14.9 46.5 35.4

aVTZ ∆ECP b −16.0 37.4 25.7

∆CPc −1.1 −9.1 −9.7

38.3 26.5

37.4 25.9

−0.9 −0.6

37.7 25.8

38.1 26.4

0.4 0.6

PNO-LMP2 PNO-LMP2 PNO-LMP2

1 2 3

33.4 36.3 40.0

31.3 33.4 34.1

−2.1 −2.9 −5.6

35.9 37.6 40.4

34.3 35.4 36.2

−1.5 −2.2 −4.2

PNO-LMP2-F12 PNO-LMP2-F12 PNO-LMP2-F12

1 2 3

34.6 36.3 37.4

34.4 36.0 36.7

−0.2 −0.4 −0.7

35.9 37.0 37.4

35.9 36.9 37.5

−0.0 −0.1 0.1

PNO-SCS-LMP2 PNO-SCS-LMP2 PNO-SCS-LMP2

1 2 3

23.1 25.5 29.0

20.9 22.6 23.1

−2.2 −2.9 −6.0

24.9 26.4 29.2

23.3 24.1 24.8

−1.6 −2.3 −4.5

PNO-SCS-LMP2-F12 PNO-SCS-LMP2-F12 PNO-SCS-LMP2-F12

1 2 3

23.9 25.2 26.0

23.7 24.9 25.5

−0.2 −0.3 −0.5

24.8 25.6 25.9

24.7 25.6 26.0

−0.1 −0.0 0.1

PNO-LCCSD PNO-LCCSD PNO-LCCSD

1 2 3

11.8 13.5 16.4

9.6 10.6 10.8

−2.3 −2.9 −5.6

13.0 14.0 16.1

11.4 11.8 12.2

−1.6 −2.1 −3.9

PNO-LCCSD-F12 PNO-LCCSD-F12 PNO-LCCSD-F12

1 2 3

11.8 12.4 12.8

11.5 12.0 12.3

−0.3 −0.4 −0.4

12.4 12.8 12.9

12.2 12.6 12.8

−0.2 −0.2 −0.1

a

b c

All calculations are carried out with default options, except IEXT and REXT = (2 × IEXT + 1) a0 . Counterpoise corrected values. Counterpoise corrections.

with increasing PAO domain size. The CP-corrected values converge much better, but even for IEXT=3 the error still amounts to 1.2 kcal mol−1 for LMP2, and 0.9 kcal mol−1 for SCS-LMP2. It remains unclear whether these differences are due to inaccuracies of the CP corrections or other 49 ACS Paragon Plus Environment

Journal of Chemical Theory and Computation

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

effects. Part of them may be due to missing ionic and exchange contributions in the local methods. 15 Fortunately, these problems are to a large extent eliminated in the local explicitly correlated methods. Table 9 demonstrates that the PNO-LCCSD-F12 results are much less dependent on IEXT, and the values obtained with our default value IEXT=2 are very close to those with IEXT=3. Furthermore, the CP corrections are reduced by an order of magnitude. As recently pointed out by Brauer et al., 175 there may be a favorable error compensation if the CP correction is not applied for small basis sets such as cc-pVDZ-F12. This is consistent with the data in Table 9. For larger basis sets (aVTZ) the CP corrections are very small can be neglected anyway in local F12 calculations. Further evidence for the accuracy is provided by the fact that the MP2-F12 and PNO-LMP2-F12 values are in close agreement. As already mentioned, similar effects can also be caused by intramolecular BSSE effects in other molecules. This makes the convergence of local results toward the canonical ones very slow (see for example Figure 4). Depending on whether the BSSE speeds up the convergence to the CBS limit (typical for molecular dissociation energies) or slows it down (typical for intermolecular interactions) the local results may be closer or further away from the CBS limit than the canonical ones, and the choice of the PAO domain size has a large effect. We therefore strongly recommend only to use the F12 methods.

6.7 Transition-metal coordination reactions Organometallic compounds play a crucial role in many fields of chemistry such as homogenous catalysis, photochemistry, and redox reactions in general. The metal center introduces significant correlation effects whose local approximations (pair, domain, and projection approximations) may behave differently from pure main group element compounds. It is therefore necessary to also consider organometallic compounds when testing the robustness of local methods. We present in this section preliminary results on the WCCR10 benchmark set of transitionmetal coordination reactions, which have been studied previously by Weymuth et al., 176 both ex50 ACS Paragon Plus Environment

Page 50 of 74

Page 51 of 74

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Theory and Computation

perimentally and computationally. The molecules in this set have comparable sizes to the Auamin molecule (see Figure 1) discussed in previous sections. An exception is the reaction 4 of the set, which involves a Ru-complex of 174 atoms. In all the reactions, a positively charged organometallic complex dissociates into two fragments. All ten reactions in this set are assumed to have predominantly closed-shell single reference character, and they are therefore challenging test cases for our PNO-LCCSD-F12 method. Weymuth et al. 176 tested a number of widely used DFT functionals by comparison with their experimental dissociation energies (obtained by L-CID experiments 168 ) with and without several versions of Grimme’s dispersion corrections. 177,178 They reported mean absolute deviations (MAD) between 4.9 and 10.7 kcal mol−1 as well as maximum absolute deviations (MAX) between 9.4 and 21.9 kcal mol−1 . They observed in general very large deviations between the results obtained with different functionals. Since a systematic improvement of the results by dispersion corrections could not be established, they recommended to use the structures optimized using the BP86 179,180 functional and the ZPE corrections obtained by the corresponding frequency calculations. We followed this recommendation in the current calculations. Our results with default options are presented in Table 10 for double-ζ and for triple-ζ basis sets along with the experimental results. It can be seen that the basis set convergence for LCCSDF12 is generally good (the difference of double-ζ and triple-ζ results amounts at most to 0.70 kcal mol−1 ). Also the differences between LCCSD-F12a and LCCSD-F12b are mostly also far below 1 kcal mol−1 . In Table 11 we compare the LCCSD-F12b/VTZ-F12 results obtained with default options to the ones with much tighter pair, domain, and projection approximations. The deviations only amount to 0.2 kcal mol−1 or less. Thus, we believe that our results should be close to the corresponding CCSD-F12 values. However, the agreement with experimental values is fair for some reactions, but rather poor in most cases. The largest differences between the experimental and computed results are found for reactions 4 and 6, and these amount to 12.7 and 16.3 kcal mol−1 , respectively. It is still not known where these disappointingly large differences come from. From our tests we

51 ACS Paragon Plus Environment

Journal of Chemical Theory and Computation

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Table 10: PNO-LMP2-F12 and PNO-LCCSD-F12 Reaction Energies (in kcal mol−1 ) of the WCCR10 Set Using the VDZ-F12 and the VTZ-F12 Basis with the Default Options No.a exptl.b basis 1 25.89 ± 0.65 VDZ-F12 VTZ-F12 2 47.60 ± 1.20 VDZ-F12 VTZ-F12 3 48.24 ± 0.70 VDZ-F12 VTZ-F12 4 34.41 ± 0.91 VDZ-F12c VTZ-F12d 5 44.53 ± 1.20 VDZ-F12 VTZ-F12 6 48.23 ± 1.61 VDZ-F12 VTZ-F12 7 50.13 ± 1.30 VDZ-F12 VTZ-F12 8 44.56 ± 2.01 VDZ-F12 VTZ-F12 9 38.73 ± 1.01 VDZ-F12 VTZ-F12 10 22.77 ± 0.60 VDZ-F12 VTZ-F12 a b c d

HF HF+CABS LMP2 LMP2-F12 LCCSD LCCSD-F12a LCCSD-F12b 25.72 24.55 23.99 22.42 26.13 24.53 24.42 24.59 24.43 22.25 21.73 24.91 24.41 24.43 43.20 42.75 60.91 59.78 57.47 56.27 55.64 42.89 42.63 58.79 59.17 55.52 55.80 55.34 43.60 43.15 61.44 60.25 57.76 56.51 55.92 43.28 43.02 59.31 59.57 55.90 56.07 55.59 33.67 31.48 51.48 52.36 46.06 47.17 46.41 31.40 31.12 51.21 52.55 46.25 47.62 47.06 24.43 22.52 57.20 55.03 43.85 41.52 40.45 22.50 22.34 56.12 55.70 42.32 41.62 40.87 49.25 49.02 70.92 72.54 63.06 64.46 64.50 48.95 48.98 71.98 72.74 64.00 64.54 64.55 40.45 40.16 68.23 70.06 57.70 59.35 58.66 40.12 40.11 69.30 70.00 58.25 58.86 58.28 31.19 31.05 52.34 54.17 45.80 47.06 46.80 31.01 31.01 53.59 54.31 46.52 46.76 46.56 10.94 10.44 47.49 48.78 32.35 33.58 32.38 10.39 10.34 48.88 49.01 33.03 32.98 32.11 9.553 9.624 25.86 26.69 20.25 21.07 21.19 9.70 9.62 26.63 26.90 20.95 21.17 21.33

The reaction numbering is as given in ref 176. The zero point energy is substracted from the given experimental values. The cc-pVDZ basis set was used for hydrogen atoms in order to avoid basis set linear dependencies. The cc-pVTZ basis set was used for hydrogen atoms in order to avoid basis set linear dependencies.

can almost certainly exclude that basis set effects or local approximations cause the discrepancies. Triple (or higher) excitations, which are still missing in our calculations, could have a significant effect, but it is rather unlikely that this amounts to 10 kcal mol−1 or more. Accounting for the relativistic effects using ECPs is also very unlikely to cause an error of this magnitude. Other sources of errors could be a breakdown of the CCSD approximation due to strong multi-reference character, wrong geometries, or a problem in the experimental analysis. We were able to reproduce the geometries and various DFT results in ref 176 using Molpro. Since the geometries obtained with different DFT functionals differ, we also carried out PNOLCCSD-F12 calculations for selected reactions at various different DFT geometries, but the results remained very similar. We also computed distinguishable cluster singles and doubles 35,181,182

52 ACS Paragon Plus Environment

Page 52 of 74

Page 53 of 74

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Theory and Computation

(PNO-LDCSD-F12b) energies, which are known to be often somewhat more accurate than the CCSD ones, in particular when the multi-reference character becomes strong. Results of the two methods are compared in Table 11. The LDCSD-F12b values do not lie systematically closer to the experimental results and do not deviate more than 2.1 kcal mol−1 from the LCCSD-F12b values. Furthermore, we probed the potential multireference character of the molecules of the test set with the D1 183 and T 1 184 diagnostics, which are also presented in Table 11. In their original publication, Janssen and Nielsen 183 concluded in their studies on small and medium size molecules that the D1 diagnostic performs better than the T 1 diagnostic in terms of size-intensitivity. They conclude that values of D1 (CCSD) < 0.05 would indicate a well performing single-reference coupled cluster ansatz. This criterion is fulfilled for all molecules in the WCCR10 test set (cf. Table 11). Lee and Taylor 184 conclude in their original study on small molecules that values of T 1 (CCSD) > 0.02 indicate a potential multireference case. All T 1 values for the molecules in this test set lie at least slightly below 0.02. A more thorough study of the potential multireference character of the molecules in the test set would require multireference calculations, but this is beyond the scope of the current paper. We can conclude that our local approximations in the PNO-LCCSD-F12 method with default options are very well converged even for these challenging transition metal cases. However, the inclusion of perturbative triples, leading to a PNO-LCCSD(T)-F12 method, will be important to achieve the high intrinsic accuracy of the single reference gold standard method CCSD(T), and only then a more reliable comparison with experimental values will be possible. A method to compute the perturbative triples correction is under development and will be the subject of a future publication.

7 Summary and Conclusions In this work we presented a parallel implementation of the PNO-LCCSD-F12 method, as well as extensive benchmarks on reaction energies and intermolecular interaction energies to illustrate

53 ACS Paragon Plus Environment

Journal of Chemical Theory and Computation

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 54 of 74

Table 11: D1 and T 1 Diagnostics, PNO-LCCSD-F12b, and PNO-LDCSD-F12b Reaction Energies (in kcal mol−1 ) of the WCCR10 Set Using the VTZ-F12 Basis D1 diagnosticd T 1 diagnosticd No. exptl. LCCSD-F12b LCCSD-F12b LDCSD-F12b compl. fragm. 1 fragm. 2 compl. fragm. 1 fragm. 2 1 25.89 ± 0.65 24.43 24.48 24.20 0.028 0.028 0.017 0.016 0.016 0.009 2 47.60 ± 1.20 55.34 55.41 55.79 0.028 0.029 0.028 0.013 0.013 0.012 3 48.24 ± 0.70 55.59 55.74 56.02 0.028 0.029 0.028 0.013 0.013 0.012 4e 34.41 ± 0.91 47.06 47.16 46.85 0.044 0.032 0.021 0.014 0.016 0.010 5 44.53 ± 1.20 40.87 40.89 42.49 0.023 0.022 0.027 0.013 0.017 0.011 6 48.23 ± 1.61 64.55 64.45 65.46 0.030 0.030 0.041 0.013 0.013 0.018 7 50.13 ± 1.30 58.28 58.30 59.87 0.031 0.037 0.041 0.013 0.013 0.018 8 44.56 ± 2.01 46.56 46.53 47.50 0.031 0.037 0.041 0.013 0.013 0.018 9 38.73 ± 1.01 32.11 31.91 33.99 0.034 0.027 0.035 0.013 0.012 0.015 10 22.77 ± 0.60 21.33 21.37 22.11 0.028 0.029 0.030 0.015 0.015 0.013 a

a b c d e

b

c

d

c

The reaction numbering is as given in ref 176. The zero point energy is substracted from the given experimental values. Computed using default options. CC CC = 0.997, T PNO = 10−8 , project_j = close, Computed using tight options (T close = 10−5 Eh , T weak = 10−6 Eh , T en project_ke = 1). The cc-pVTZ basis set was used for hydrogen atoms in order to avoid basis set linear dependencies.

its accuracy and performance. Aside from improving the basis-set convergence, the F12 terms strongly reduce the local domain errors, which are mainly caused by the PAO domains. Probably, a large fraction of these errors is due to intra- or intermolecular BSSE effects. The errors due to the PNO domain truncations are also quite effectively corrected by the F12 terms. The exclusive use of explicitly correlated methods is therefore strongly recommended. As has been demonstrated and discussed in section 3.3, it is difficult to reduce the remaining local errors below ≈1 kJ mol−1 (0.25 kcal mol−1 ) without using very large PNO domains. This would make the calculations exceedingly expensive. Nevertheless, with our default options, the local errors are probably in most cases smaller than the intrinsic error of the canonical CCSD(T)F12 method. Apart from the remaining domain errors, the accuracy obtained with triple-ζ basis sets is mainly limited by basis set incompleteness errors, which should be similar to those of the corresponding canonical CCSD-F12 method. This is indicated by the fact that CCSD-F12a and CCSD-F12b relative energies with the VTZ-F12 basis set may differ by up to ∼0.5 kcal mol−1 , and approximations in the F12 ansatz (3*A vs 3C) can also have effects of similar size. The benchmarks in this pa-

54 ACS Paragon Plus Environment

Page 55 of 74

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Theory and Computation

per indicate that the PNO-LCCSD-F12b results are more consistent and reliable, as also expected theoretically. The PNO-LCCSD-F12 method is designed to minimize computational redundancy by computing nearly all types of integrals only once for the whole molecule. As a trade-off, the communication requirements are rather significant, and this limits its efficient parallelization to a few hundred CPU cores. The F12 treatment that is highly cost-effective. Typically, a PNO-LCCSDF12 calculation takes only 25–40% longer than a corresponding PNO-LCCSD calculation, and most of the additional time is spent for computing the F12 correction of the LMP2 energy. Other additional resources required for the F12 terms, such as memory and disk space, are negligible. Our method shows good parallel efficiency using up to ∼200 CPU cores, and can perform PNO-LCCSD-F12/VTZ-F12 calculations for molecules with 100–200 atoms within a few hours of elapsed time. Our current program does not yet include the perturbative triples energy correction, which is crucial for the application of the PNO-LCCSD-F12 method in real-life chemistry. Furthermore, many applications require the extension of the current method to open-shell cases. Work in these directions is in progress. The application of local methods for geometry optimizations and molecular property calculations requires the implementation of analytic gradients, which is challenging even at the MP2 level due to the additional side conditions in the Lagrangian. Some recent developments on this topic include the (approximate) analytical calculation of first-order properties for DLPNO-CCSD by Datta et al., 185 and the PNO-MP2 analytic gradients by Frank et al. 186 In our group, analytical energy gradients have been implemented quite long time ago for closed-shell PAO-LMP2, 8,24 and recently also for open-shell PAO-LRMP2. 187

Acknowledgement This work has been supported by the European Research Council (ERC) Advanced Grant 320723

55 ACS Paragon Plus Environment

Journal of Chemical Theory and Computation

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

(ASES).

Supporting Information Available • fh_summary.txt: Reaction energies for the test set of Friedrich and Hänchen computed with various methods.

References (1) Pulay, P. Localizability of dynamic electron correlation. Chem. Phys. Lett. 1983, 100, 151–154. (2) Saebø, S.; Pulay, P. Local Configuration-Interaction: An efficient approach for larger molecules. Chem. Phys. Lett. 1985, 113, 13–18. (3) Pulay, P.; Saebø, S. Orbital-invariant formulation and second-order gradient evaluation in MøllerPlesset perturbation theory. Theor. Chim. Acta 1986, 69, 357–368. (4) Saebø, S.; Pulay, P. Fourth-order Møller-Plesset perturbation theory in the local correlation treatment. I. Method. J. Chem. Phys. 1987, 86, 914–922. (5) Saebø, S.; Pulay, P. The local correlation treatment. II. Implementation and tests. J. Chem. Phys. 1988, 88, 1884–1890. (6) Saebø, S.; Pulay, P. Local treatment of electron correlation. Annu. Rev. Phys. Chem. 1993, 44, 213– 236. (7) Hampel, C.; Werner, H.-J. Local treatment of electron correlation in coupled cluster theory. J. Chem. Phys. 1996, 104, 6286–6297. (8) ElAzhary, A.; Rauhut, G.; Pulay, P.; Werner, H.-J. Analytical energy gradients for local second-order Møller-Plesset perturbation theory. J. Chem. Phys. 1998, 108, 5185–5193. (9) Maslen, P. E.; Head-Gordon, M. Non-iterative local second order Moller-Plesset theory. Chem. Phys. Lett. 1998, 283, 102–180.

56 ACS Paragon Plus Environment

Page 56 of 74

Page 57 of 74

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Theory and Computation

(10) Maslen, P. E.; Head-Gordon, M. Noniterative local second order Moller-Plesset theory: Convergence with local correlation space. J. Chem. Phys. 1998, 109, 7093–7099. (11) Ayala, P. Y.; Scuseria, G. E. Linear scaling second-order Møller-Plesset theory in the atomic orbital basis for large molecular systems. J. Chem. Phys. 1999, 110, 3660–3671. (12) Scuseria, G. E.; Ayala, P. Y. Linear scaling coupled cluster and perturbation theories in the atomic orbital basis. J. Chem. Phys. 1999, 111, 8330–8343. (13) Hetzer, G.; Pulay, P.; Werner, H.-J. Multipole approximation of distant pair energies in local MP2 calculations. Chem. Phys. Lett. 1998, 290, 143–149. (14) Hetzer, G.; Schütz, M.; Stoll, H.; Werner, H.-J. Low-order scaling local correlation methods II: Splitting the Coulomb operator in linear scaling local second-order Møller-Plesset perturbation theory. J. Chem. Phys. 2000, 113, 9443–9455. (15) Schütz, M.; Rauhut, G.; Werner, H.-J. Local treatment of electron correlation in molecular clusters: Structures and stabilities of (H2 O)n , n = 2-4. J. Phys. Chem. A 1998, 102, 5997–6003. (16) Schütz, M.; Hetzer, G.; Werner, H.-J. Low-order scaling local electron correlation methods. I. Linear scaling local MP2. J. Chem. Phys. 1999, 111, 5691–5705. (17) Schütz, M.; Werner, H.-J. Local perturbative triples correction (T) with linear cost scaling. Chem. Phys. Lett. 2000, 318, 370–378. (18) Schütz, M. Low-order scaling local electron correlation methods. III. Linear scaling local perturbative triples correction (T). J. Chem. Phys. 2000, 113, 9986–10001. (19) Schütz, M.; Werner, H.-J. Low-order scaling local electron correlation methods. IV. Linear scaling local coupled-cluster (LCCSD). J. Chem. Phys. 2001, 114, 661–681. (20) Schütz, M. Linear scaling local connected triples beyond local (T): Local CCSDT-1b with O(N) scaling. J. Chem. Phys. 2002, 116, 8772–8785. (21) Schütz, M. A new, fast, semi-direct implementation of linear scaling local coupled cluster theory. Phys. Chem. Chem. Phys. 2002, 4, 3941–3947.

57 ACS Paragon Plus Environment

Journal of Chemical Theory and Computation

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

(22) Schütz, M.; Manby, F. R. Linear scaling local coupled cluster theory with density fitting. I : 4-external integrals. Phys. Chem. Chem. Phys. 2003, 5, 3349–3358. (23) Werner, H.-J.; Manby, F. R.; Knowles, P. Fast linear scaling second-order Møller-Plesset perturbation theory (MP2) using local and density fitting approximations. J. Chem. Phys. 2003, 118, 8149–8160. (24) Schütz, M.; Werner, H.-J.; Lindh, R.; Manby, F. R. Analytical energy gradients for local second-order Møller-Plesset perturbation theory using density fitting approximations. J. Chem. Phys. 2004, 121, 737–750. (25) Auer, A.; Nooijen, M. Dynamically screened local correlation method using enveloping localized orbitals. J. Chem. Phys. 2006, 125, 024104. (26) Hrenar, T.; Rauhut, G.; Werner, H.-J. Impact of local and density fitting approximations on harmonic vibrational frequencies. J. Phys. Chem. A 2006, 110, 2060–2064. (27) Werner, H.-J.; Pflüger, K. On the selection of domains and orbital pairs in local correlation treatments. Ann. Reports in Comput. Chem. 2006, 2, 53–80. (28) Mata, R.; Werner, H.-J. Calculation of smooth potential energy surfaces using local electron correlation methods. J. Chem. Phys. 2006, 125, 184110. (29) Mata, R.; Werner, H.-J. Local correlation methods with a natural localized molecular orbital basis. Mol. Phys. 2007, 105, 2753–2761. (30) Mata, R. A.; Werner, H.-J.; Thiel, S.; Thiel, W. Towards accurate barriers for enzymatic reactions: QM/MM case study on p-hydroxybenzoate hydroxylase. J. Chem. Phys. 2008, 128, 025104. (31) Mata, R.; Werner, H.-J.; Schütz, M. Correlation regions within a localized molecular orbital approach. J. Chem. Phys. 2008, 128, 144106. (32) Werner, H.-J.; Schütz, M. An efficient local coupled cluster method for accurate thermochemistry of large systems. J. Chem. Phys. 2011, 135, 144116. (33) Schwilk, M.; Usvyat, D.; Werner, H.-J. Communication: Improved pair approximations in local coupled-cluster methods. J. Chem. Phys. 2015, 142, 121102.

58 ACS Paragon Plus Environment

Page 58 of 74

Page 59 of 74

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Theory and Computation

(34) Köppl, C.; Werner, H.-J. On the use of Abelian point group symmetry in density-fitted local MP2 using various types of virtual orbitals. J. Chem. Phys. 2015, 142, 164108. (35) Kats, D.; Manby, F. R. Sparse tensor framework for implementation of general local correlation methods. J. Chem. Phys. 2013, 138, 144101. (36) Kats, D. Speeding up local correlation methods. J. Chem. Phys. 2014, 141, 244101. (37) Russ, N. J.; Crawford, T. D. Local correlation in coupled cluster calculations of molecular response properties. Chem. Phys. Lett. 2004, 400, 104–111. (38) Maslen, P. E.; Dutoi, A. D.; Lee, M. S.; Shao, Y.; Head-Gordon, M. Accurate local approximations to the triples correlation energy: formulation, implementation and tests of 5th-order scaling models. Mol. Phys. 2005, 103, 425–437. (39) DiStasio, R. A.; Jung, Y. S.; Head-Gordon, M. A resolution-of-the-identity implementation of the local triatomics-in-molecules model for second-order Møller-Plesset perturbation theory with application to alanine tetrapeptide conformational energies. J. Chem. Theory Comput. 2005, 1, 862–876. (40) Subotnik, J. E.; Head-Gordon, M. A local correlation model that yields intrinsically smooth potentialenergy surfaces. J. Chem. Phys. 2005, 123, 064108. (41) Lawler, K. V.; Parkhill, J. A.; Head-Gordon, M. Penalty functions for combining coupled-cluster and perturbation amplitudes in local correlation methods with optimized orbitals. Mol. Phys. 2008, 106, 2309–2324. (42) Schweizer, S.; Kussmann, J.; Doser, B.; Ochsenfeld, C. Linear-scaling Cholesky decomposition. J. Comp. Chem. 2008, 29, 1004–1010. (43) Schweizer, S.; Doser, B.; Ochsenfeld, C. An atomic orbital-based reformulation of energy gradients in second-order Møller-Plesset perturbation theory. J. Chem. Phys. 2008, 128, 154101. (44) Doser, B.; Lambrecht, D. S.; Kussmann, J.; Ochsenfeld, C. Linear-scaling atomic orbital-based second-order Møller–Plesset perturbation theory by rigorous integral screening criteria. J. Chem. Phys. 2009, 130, 064107.

59 ACS Paragon Plus Environment

Journal of Chemical Theory and Computation

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

(45) Doser, B.; Zienau, J.; Clin, L.; Lambrecht, D. S.; Ochsenfeld, C. A linear-scaling MP2 method for large molecules by rigorous integral-screening criteria. Z. Phys. Chem. 2010, 224, 397–412. (46) Kats, D.; Schütz, M. A multistate local coupled cluster CC2 response method based on the Laplace transform. J. Chem. Phys. 2009, 131, 124117. (47) Freundorfer, K.; Kats, D.; Korona, T.; Schütz, M. Local CC2 response method for triplet states based on Laplace transform: Excitation energies and first-order properties. J. Chem. Phys. 2010, 133, 244110. (48) Maschio, L. Local MP2 with density fitting for periodic systems: A parallel implementation. J. Chem. Theory Comput. 2011, 7, 2818–2830. (49) Neese, F.; Wennmohs, F.; Hansen, A. Efficient and accurate local approximations to coupled-electron pair approaches: An attempt to revive the pair natural orbital method. J. Chem. Phys. 2009, 130, 114108. (50) Neese, F.; Hansen, A.; Liakos, D. G. Efficient and accurate approximations to the local coupled cluster singles doubles method using a truncated pair natural orbital basis. J. Chem. Phys. 2009, 131, 064103. (51) Hansen, A.; Liakos, D. G.; Neese, F. Efficient and accurate local single reference correlation methods for high-spin open-shell molecules using pair natural orbitals. J. Chem. Phys. 2011, 135, 214102. (52) Liakos, D. G.; Hansen, A.; Neese, F. Weak molecular interactions studied with parallel implementations of the local pair natural orbital coupled pair and coupled cluster methods. J. Chem. Theory Comput. 2011, 7, 76–87. (53) Izsak, R.; Hansen, A.; Neese, F. The resolution of identity and chain of spheres approximations for the LPNO-CCSD singles Fock term. Mol. Phys. 2012, 110, 2413–2417. (54) Schmitz, G.; Helmich, B.; Hättig, C. A O(N 3 ) scaling PNO-MP2 method using a hybrid OSV-PNO approach with an iterative direct generation of OSVs. Mol. Phys. 2013, 111, 2463–2476.

60 ACS Paragon Plus Environment

Page 60 of 74

Page 61 of 74

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Theory and Computation

(55) Riplinger, C.; Neese, F. An efficient and near linear scaling pair natural orbital based local coupled cluster method. J. Chem. Phys. 2013, 138, 034106. (56) Riplinger, C.; Sandhoefer, B.; Hansen, A.; Neese, F. Natural triple excitations in local coupled cluster calculations with pair natural orbitals. J. Chem. Phys. 2013, 139, 134101. (57) Pinski, P.; Riplinger, C.; Valeev, E. F.; Neese, F. Sparse maps—A systematic infrastructure for reduced-scaling electronic structure methods. I. An efficient and simple linear scaling local MP2 method that uses an intermediate basis of pair natural orbitals. J. Chem. Phys. 2015, 143, 034108. (58) Riplinger, C.; Pinski, P.; Becker, U.; Valeev, E. F.; Neese, F. Sparse maps-A systematic infrastructure for reduced-scaling electronic structure methods. II. Linear scaling domain based pair natural orbital coupled cluster theory. J. Chem. Phys. 2016, 144, 024109. (59) Maurer, S. A.; Clin, L.; Ochsenfeld, C. Cholesky-decomposed density MP2 with density fitting: accurate MP2 and double-hybrid DFT energies for large systems. J. Chem. Phys. 2014, 140, 224112. (60) Schurkus, H. F.; Ochsenfeld, C. Communication: An effective linear-scaling atomic-orbital reformulation of the random-phase approximation using a contracted double-Laplace transformation. J. Chem. Phys. 2016, 144, 031101. (61) Werner, H.-J.; Knizia, G.; Krause, C.; Schwilk, M.; Dornbach, M. Scalable electron correlation methods I.: PNO-LMP2 with linear scaling in the molecular size and near-inverse-linear scaling in the number of processors. J. Chem. Theory Comput. 2015, 11, 484–507. (62) Schwilk, M.; Ma, Q.; Köppl, C.; Werner, H.-J. Scalable electron correlation methods. 3. Efficient and accurate parallel local coupled cluster with pair natural orbitals (PNO-LCCSD). J. Chem. Theory Comput. 2017, 13, 3650–3675. (63) Saebø, S.; Tong, W.; Pulay, P. Efficient elimination of basis set superposition errors by the local correlation method: Accurate ab initio studies of the water dimer. J. Chem. Phys. 1993, 98, 2170– 2175. (64) Kato, T. On the eigenfunctions of many-particle systems in quantum mechanics. Comm. Pure Appl. Math. 1957, 10, 151–177.

61 ACS Paragon Plus Environment

Journal of Chemical Theory and Computation

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

(65) Pack, R. T.; Byers Brown, W. Cusp Conditions for Molecular Wavefunctions. J. Chem. Phys. 1966, 45, 556–559. (66) Hylleraas, E. A. Neue Berechnung der Energie des Heliums im Grundzustande, sowie des tiefsten Terms von Ortho-Helium. Z. Phys. 1929, 54, 347–366. (67) Kutzelnigg, W. r12 -Dependent terms in the wave function as closed sums of partial wave amplitudes for large l. Theor. Chim. Acta 1985, 68, 445–469. (68) Klopper, W.; Kutzelnigg, W. Møller–Plesset calculations taking care of the correlation cusp. Chem. Phys. Lett. 1987, 134, 17–22. (69) Klopper, W.; Kutzelnigg, W. MP2-R12 calculations on the relative stability of carbocations. J. Phys. Chem. 1990, 94, 5625–5630. (70) Klopper, W.; Röhse, R.; Kutzelnigg, W. CID and CEPA calculations with linear r12 terms. Chem. Phys. Lett. 1991, 178, 455–461. (71) Kutzelnigg, W.; Klopper, W. Wave functions with terms linear in the interelectronic coordinates to take care of the correlation cusp. I. General theory. J. Chem. Phys. 1991, 94, 1985–2001. (72) Termath, V.; Klopper, W.; Kutzelnigg, W. Wave functions with terms linear in the interelectronic coordinates to take care of the correlation cusp. II. Second-order Møller-Plesset (MP2-R12) calculations on closed-shell atoms. J. Chem. Phys. 1991, 94, 2002–2019. (73) Klopper, W.; Kutzelnigg, W. Wave functions with terms linear in the interelectronic coordinates to take care of the correlation cusp. III. Second-order Møller-Plesset (MP2-R12) calculations on molecules of first row atoms. J. Chem. Phys. 1991, 94, 2020–2030. (74) Klopper, W. Orbital-invariant formulation of the MP2-R12 method. Chem. Phys. Lett. 1991, 186, 583–585. (75) Noga, J.; Kutzelnigg, W.; Klopper, W. CC-R12, a correlation cusp corrected coupled-cluster method with a pilot application to the Be2 potential curve. Chem. Phys. Lett. 1992, 199, 497–504.

62 ACS Paragon Plus Environment

Page 62 of 74

Page 63 of 74

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Theory and Computation

(76) Noga, J.; Kutzelnigg, W. Coupled cluster theory that takes care of the correlation cusp by inclusion of linear terms in the interelectronic coordinates. J. Chem. Phys. 1994, 101, 7738–7762. (77) Klopper, W.; Samson, C. C. M. Explicitly correlated second-order Møller-Plesset methods with auxiliary basis sets. J. Chem. Phys. 2002, 116, 6397–6410. (78) Manby, F. R. Density fitting in second-order linear-r12 Møller-Plesset perturbation theory. J. Chem. Phys. 2003, 119, 4607–4613. (79) Ten-no, S. Initiation of explicitly correlated Slater-type geminal theory. Chem. Phys. Lett. 2004, 398, 56–61. (80) Ten-no, S. Explicitly correlated second order perturbation theory: Introduction of a rational generator and numerical quadratures. J. Chem. Phys. 2004, 121, 117–129. (81) May, A. J.; Manby, F. R. An explicitly correlated second order Møller-Plesset theory using a frozen Gaussian geminal. J. Chem. Phys. 2004, 121, 4479–4485. (82) Valeev, E. F. Improving on the resolution of the identity in linear R12 ab initio theories. Chem. Phys. Lett. 2004, 395, 190–195. (83) Tew, D. P.; Klopper, W. New correlation factors for explicitly correlated electronic wave functions. J. Chem. Phys. 2005, 123, 074101. (84) Kedžuch, S.; Milko, M.; Noga, J. Alternative formulation of the matrix elements in MP2-R12 theory. Int. J. Quantum Chem. 2005, 105, 929–936. (85) Fliegl, H.; Klopper, W.; Hättig, C. Coupled-cluster theory with simplified linear-r12 corrections: The CCSD(R12) model. J. Chem. Phys. 2005, 122, 084107. (86) Fliegl, H.; Hättig, C.; Klopper, W. Inclusion of the (T) triples Correction into the linear-r12 corrected coupled-cluster model CCSD(R12). Int. J. Quantum Chem. 2006, 106, 2306–2317. (87) Manby, F. R.; Werner, H.-J.; Adler, T. B.; May, A. J. Explicitly correlated local second-order perturbation theory with a frozen geminal correlation factor. J. Chem. Phys. 2006, 124, 094103.

63 ACS Paragon Plus Environment

Journal of Chemical Theory and Computation

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

(88) Werner, H.-J.; Adler, T. B.; Manby, F. R. General orbital invariant MP2-F12 theory. J. Chem. Phys. 2007, 126, 164102. (89) Noga, J.; Kedžuch, S.; Šimunek, J. Second order explicitly correlated R12 theory revisited: A second quantization framework for treatment of the operators’ partitionings. J. Chem. Phys. 2007, 127, 034106. (90) Adler, T. B.; Knizia, G.; Werner, H.-J. A simple and efficient CCSD(T)-F12 approximation. J. Chem. Phys. 2007, 127, 221106. (91) Ten-no, S. A simple F12 geminal correction in multi-reference perturbation theory. Chem. Phys. Lett. 2007, 447, 175–179. (92) Tew, D. P.; Klopper, W.; Neiss, C.; Hättig, C. Quintuple-zeta quality coupled-cluster correlation energies with triple-zeta basis sets. Phys. Chem. Chem. Phys. 2007, 9, 1921–1930. (93) Knizia, G.; Werner, H.-J. Explicitly correlated RMP2 for high-spin open-shell reference states. J. Chem. Phys. 2008, 128, 154103. (94) Shiozaki, T.; Kamiya, M.; Hirata, S.; Valeev, E. F. Explicitly correlated coupled-cluster singles and doubles method based on complete diagrammatic equations. J. Chem. Phys. 2008, 129, 071101. (95) Shiozaki, T.; Kamiya, M.; Hirata, S.; Valeev, E. F. Equations of explicitly-correlated coupled-cluster methods. Phys. Chem. Chem. Phys. 2008, 10, 3358–3370. (96) Noga, J.; Kedžuch, S.; Šimunek, J.; Ten-no, S. Explicitly correlated coupled cluster F12 theory with single and double excitations. J. Chem. Phys. 2008, 128, 174103. (97) Tew, D. P.; Klopper, W.; Hättig, C. A diagonal orbital-invariant explicitly correlated coupled-cluster method. Chem. Phys. Lett. 2008, 452, 326–332. (98) Valeev, E. F. Coupled-cluster methods with perturbative inclusion of explicitly correlated terms: A preliminary investigation. Phys. Chem. Chem. Phys. 2008, 10, 106–113.

64 ACS Paragon Plus Environment

Page 64 of 74

Page 65 of 74

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Theory and Computation

(99) Valeev, E. F.; Crawford, T. D. Simple coupled-cluster singles and doubles method with perturbative inclusion of triples and explicitly correlated geminals: The CCSD(T)R12 model. J. Chem. Phys. 2008, 128, 244113. (100) Torheyden, M.; Valeev, E. F. Variational formulation of perturbative explicitly-correlated coupledcluster methods. Phys. Chem. Chem. Phys. 2008, 10, 3410–3420. (101) Bokhan, D.; Ten-no, S.; Noga, J. Implementation of the CCSD(T)-F12 method using cusp conditions. Phys. Chem. Chem. Phys. 2008, 10, 3320–3326. (102) Köhn, A.; Richings, G. W.; Tew, D. P. Implementation of the full explicitly correlated coupled-cluster singles and doubles model CCSD-F12 with optimally reduced auxiliary basis dependence. J. Chem. Phys. 2008, 129, 201103. (103) Knizia, G.; Adler, T. B.; Werner, H.-J. Simplified CCSD(T)-F12 methods: Theory and benchmarks. J. Chem. Phys. 2009, 130, 054104. (104) Shiozaki, T.; Kamiya, M.; Hirata, S.; Valeev, E. F. Higher-order explicitly correlated coupled-cluster methods. J. Chem. Phys. 2009, 130, 054101. (105) Bokhan, D.; Bernadotte, S.; Ten-no, S. Implementation of the CCSD(T)(F12) method using numerical quadratures. Chem. Phys. Lett. 2009, 469, 214–218. (106) Bischoff, F. A.; Wolfsegger, S.; Tew, D. P.; Klopper, W. Assessment of basis sets for F12 explicitlycorrelated molecular electronic-structure methods. Mol. Phys. 2009, 107, 963–975. (107) Hättig, C.; Tew, D. P.; Köhn, A. Accurate and efficient approximations to explicitly correlated coupled-cluster singles and doubles, CCSD-F12. J. Chem. Phys. 2010, 132, 231102. (108) Werner, H.-J.; Knizia, G.; Manby, F. R. Explicitly correlated coupled cluster methods with pairspecific geminals. Mol. Phys 2011, 109, 407–417. (109) Shiozaki, T.; Knizia, G.; Werner, H.-J. Explicitly correlated multireference configuration interaction: MRCI-F12. J. Chem. Phys. 2011, 134, 034113.

65 ACS Paragon Plus Environment

Journal of Chemical Theory and Computation

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

(110) Werner, H.-J.; Manby, F. R. Explicitly correlated second-order perturbation theory using density fitting and local approximations. J. Chem. Phys. 2006, 124, 054114. (111) Manby, F. R.; Werner, H.-J.; Adler, T. B.; May, A. J. Explicitly correlated local second-order perturbation theory with a frozen geminal correlation factor. J. Chem. Phys. 2006, 124, 094103. (112) Werner, H.-J. Eliminating the domain error in local explicitly correlated second-order Møller-Plesset perturbation theory. J. Chem. Phys. 2008, 129, 101103. (113) Adler, T. B.; Werner, H.-J. Local explicitly correlated coupled-cluster methods: Efficient removal of the basis set incompleteness and domain errors. J. Chem. Phys. 2009, 130, 241101. (114) Adler, T. B.; Werner, H.-J.; Manby, F. R. Local explicitly correlated second-order perturbation theory for the accurate treatment of large molecules. J. Chem. Phys. 2009, 130, 054106. (115) Adler, T. B.; Werner, H.-J. An explicitly correlated local coupled cluster method for calculations of large molecules close to the basis set limit. J. Chem. Phys. 2011, 135, 144117. (116) Krause, C.; Werner, H.-J. Comparison of explicitly correlated local coupled-cluster methods with various choices of virtual orbitals. Phys. Chem. Chem. Phys. 2012, 14, 7591–7604. (117) Ma, Q.; Werner, H.-J. Scalable electron correlation methods. 2. Parallel PNO-LMP2-F12 with near linear scaling in the molecular size. J. Chem. Theory Comput. 2015, 11, 5291–5304. (118) Tew, D. P.; Helmich, B.; Hättig, C. Local explicitly correlated second-order Moller-Plesset perturbation theory with pair natural orbitals. J. Chem. Phys. 2011, 135, 074107. (119) Hättig, C.; Tew, D. P.; Helmich, B. Local explicitly correlated second- and third-order Møller-Plesset perturbation theory with pair natural orbitals. J. Chem. Phys. 2012, 136, 204105. (120) Li, W. Linear scaling explicitly correlated MP2-F12 and ONIOM methods for the long-range interactions of the nanoscale clusters in methanol aqueous solutions. J. Chem. Phys. 2013, 138, 014106. (121) Wang, Y. M.; Hättig, C.; Reine, S.; Valeev, E.; Kjaergaard, T.; Kristensen, K. Explicitly correlated second-order Møller-Plesset perturbation theory in a Divide-Expand-Consolidate (DEC) context. J. Chem. Phys. 2016, 144, 204112.

66 ACS Paragon Plus Environment

Page 66 of 74

Page 67 of 74

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Theory and Computation

(122) Pavoševi´c, F.; Pinski, P.; Riplinger, C.; Neese, F.; Valeev, E. F. Sparse maps-A systematic infrastructure for reduced-scaling electronic structure methods. IV. Linear-scaling second-order explicitly correlated energy with pair natural orbitals. J. Chem. Phys. 2016, 144, 144109. (123) Marchetti, O.; Werner, H.-J. Accurate calculations of intermolecular interaction energies using explicitly correlated wave functions. Phys. Chem. Chem. Phys. 2008, 10, 3400. (124) Marchetti, O.; Werner, H.-J. Accurate calculations of intermolecular interaction energies using explicitly correlated coupled cluster wave functions and a dispersion weighted MP2 method. J. Phys. Chem. A 2009, 113, 11580–11585. (125) Sirianni, D. A.; Burns, L. A.; Sherrill, C. D. Comparison of explicitly correlated methods for computing high-accuracy benchmark energies for noncovalent interactions. J. Chem. Theory Comput. 2017, 13, 86–99. (126) Adler, T. B.; Werner, H.-J. An explicitly correlated local coupled cluster method for calculations of large molecules close to the basis set limit. J. Chem. Phys. 2011, 135, 144117. (127) Adler, T. B.; Werner, H.-J. Local explicitly correlated coupled-cluster methods: Efficient removal of the basis set incompleteness and domain errors. J. Chem. Phys. 2009, 130, 241101. (128) Schmitz, G.; Hättig, C. Perturbative triples correction for local pair natural orbital based explicitly correlated CCSD(F12*) using Laplace transformation techniques. J. Chem. Phys. 2016, 145, 234107. (129) Schmitz, G.; Hättig, C.; Tew, D. P. Explicitly correlated PNO-MP2 and PNO-CCSD and their application to the S66 set and large molecular systems. Phys. Chem. Chem. Phys. 2014, 16, 22167–22178. (130) Schmitz, G.; Hättig, C. Accuracy of explicitly correlated local PNO-CCSD(T). J. Chem. Theory Comput. 2017, 13, 2623–2633. (131) Pavoševi´c, F.; Peng, C.; Pinski, P.; Riplinger, C.; Neese, F.; Valeev, E. F. SparseMaps — A systematic infrastructure for reduced scaling electronic structure methods. V. Linear scaling explicitly correlated coupled-cluster method with pair natural orbitals. J. Chem. Phys. 2017, 146, 174108.

67 ACS Paragon Plus Environment

Journal of Chemical Theory and Computation

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

(132) Rendell, A. P.; Lee, T. J.; Lindh, R. Quantum chemistry on parallel computer architectures: Coupledcluster theory applied to the bending potential of fulminic acid. Chem. Phys. Lett. 1992, 194. (133) Kobayashi, R.; Rendell, A. P. A direct coupled cluster algorithm for massively parallel computers. Chem. Phys. Lett. 1997, 265. (134) Olson, R. M.; Bentz, J. L.; Kendall, R. A.; Schmidt, M. W.; Gordon, M. S. A novel approach to parallel coupled cluster calculations: Combining distributed and shared memory techniques for modern cluster based systems. J. Chem. Theory Comput. 2007, 3, 1312–1328. (135) Janowski, T.; Ford, A. R.; ; Pulay, P. Parallel calculation of coupled cluster singles and doubles wave functions using array files. J. Chem. Theory Comput. 2007, 3, 1368–1377. (136) Lotrich, V.; Flocke, N.; Ponton, M.; Yau, A. D.; Perera, A.; Deumens, E.; Bartlett, R. J. Parallel implementation of electronic structure energy, gradient, and Hessian calculations. J. Chem. Phys. 2008, 128, 194104. (137) Janowski, T.; Pulay, P. Efficient parallel implementation of the CCSD external exchange operator and the perturbative triples (T) energy calculation. J. Chem. Theory Comput. 2008, 4, 1585–1592. (138) Ku´s, T.; Lotrich, V. F.; Bartlett, R. J. Parallel implementation of the equation-of-motion coupledcluster singles and doubles method and application for radical adducts of cytosine. J. Chem. Phys. 2009, 130, 124122. (139) Asadchev, A.; Gordon, M. S. Fast and flexible coupled cluster implementation. J. Chem. Theory Comput. 2013, 9, 3385–3392. (140) Anisimov, V. M.; Bauer, G. H.; Chadalavada, K.; Olson, R. M.; Glenski, J. W.; Kramer, W. T. C.; Aprà, E.; Kowalski, K. Optimization of the coupled cluster implementation in NWChem on petascale parallel architectures. J. Chem. Theory Comput. 2014, 10, 4307–4316. (141) Solomonik, E.; Matthews, D.; Hammond, J. R.; Stanton, J. F.; Demmel, J. A massively parallel tensor contraction framework for coupled-cluster computations. J. Parallel Distrib. Comput. 2014, 74, 3176–3190.

68 ACS Paragon Plus Environment

Page 68 of 74

Page 69 of 74

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Theory and Computation

(142) Peng, C.; Calvin, J. A.; Pavoševi´c, F.; Zhang, J.; Valeev, E. F. Massively parallel implementation of explicitly correlated coupled-cluster singles and doubles using TiledArray framework. J. Phys. Chem. A 2016, 120, 10231–10244. (143) Knizia, G. Intrinsic atomic orbitals: An unbiased bridge between quantum theory and chemical concepts. J. Chem. Theory Comput. 2013, 9, 4834–4843. (144) Yousaf, K. E.; Peterson, K. A. Optimized auxiliary basis sets for explicitly correlated methods. J. Chem. Phys. 2008, 129, 184108. (145) Yousaf, K. E.; Peterson, K. A. Optimized complementarty auxiliary basis sets for explicitly correlated methods: aug-cc-pVnZ sets. Chem. Phys. Lett. 2009, 476, 303. (146) Weigend, F. A fully direct RI-HF algorithm: Implementation, optimized auxiliary basis sets, demonstration of accuracy and efficiency. Phys. Chem. Chem. Phys. 2002, 4, 4285–4291. (147) Werner, H.-J.; Adler, T. B.; Knizia, G.; Manby, F. R. In Recent Progress in Coupled Cluster Methods; ˇ Cársky, P., Paldus, J., Pittner, J., Eds.; Springer: Dordrecht, Heidelberg, London, New York, 2010; pp 573–619. (148) Werner, H.-J.; Knizia, G.; Adler, T. B.; Marchetti, O. Benchmark studies for explicitly correlated perturbation- and coupled cluster theories. Z. Phys. Chem. 2010, 224, 493–511. (149) Fedorov, A.; Batiste, L.; Couzijn, E. P. A.; Chen, P. Experimental and theoretical study of a gold(I) aminonitrene complex in the gas phase. ChemPhysChem 2010, 11, 1002–1005. (150) Huenerbein, R.; Schirmer, B.; Moellmann, J.; Grimme, S. Effects of London dispersion on the isomerization reactions of large organic molecules: a density functional benchmark study. Phys. Chem. Chem. Phys. 2010, 12, 6940–6948. (151) Peterson, K. A.; Adler, T. B.; Werner, H.-J. Systematically convergent basis sets for explicitly correlated wavefunctions: The atoms H, He, B-Ne, and Al-Ar. J. Chem. Phys. 2008, 128, 084102. (152) Friedrich, J.; Hänchen, J. Incremental CCSD(T)(F12*)|MP2: A black box method to obtain highly accurate reaction energies. J. Chem. Theory Comput. 2013, 9, 5381–5394.

69 ACS Paragon Plus Environment

Journal of Chemical Theory and Computation

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

(153) Pavoševi´c, F.; Neese, F.; Valeev, E. F. Geminal-spanning orbitals make explicitly correlated reducedscaling coupled-cluster methods robust, yet simple. J. Chem. Phys. 2014, 141, 054106. (154) Nieplocha, J.; Palmer, B.; Tipparaju, V.; Krishnan, M.; Trease, H.; Aprà, E. Advances, applications and performance of the global arrays shared memory programming toolkit. Int. J. High Perform C. 2006, 20, 203–231. (155) Karypis, G.; Kumar, V. A fast and high quality multilevel scheme for partitioning irregular graphs. SIAM J. Sci. Comput. 1998, 20, 359–392. (156) Patkowski, K. On the accuracy of explicitly correlated coupled-cluster interaction energies — have orbital results been beaten yet? J. Chem. Phys. 2012, 137, 034103. (157) Dunning, Jr., T. H. Gaussian basis sets for use in correlated molecular calculations. I. The atoms boron through neon and hydrogen. J. Chem. Phys. 1989, 90, 1007–1023. (158) Kendall, R. A.; Dunning, T. H.; Harrison, R. J. Electron affinities of the first-row atoms revisited. Systematic basis sets and wave functions. J. Chem. Phys. 1992, 96, 6796–6806. (159) Figgen, D.; Rauhut, G.; Dolg, M.; Stoll, H. Energy-consistent pseudopotentials for group 11 and 12 atoms: adjustment to multi-configuration Dirac–Hartree–Fock data. Chem. Phys. 2005, 311, 227– 244. (160) Peterson, K.; Puzzarini, C. Systematically convergent basis sets for transition metals. II. Pseudopotential-based correlation consistent basis sets for the group 11 (Cu, Ag, Au) and 12 (Zn, Cd, Hg) elements. Theor. Chem. Acc. 2005, 114, 283–296. (161) Peterson, K. A.; Figgen, D.; Dolg, M.; Stoll, H. Energy-consistent relativistic pseudopotentials and correlation consistent basis sets for the 4d elements Y–Pd. J. Chem. Phys. 2007, 126, 124101. (162) Figgen, D.; Peterson, K. A.; Dolg, M.; Stoll, H. Energy-consistent pseudopotentials and correlation consistent basis sets for the 5d elements Hf–Pt. J. Chem. Phys. 2009, 130, 164108. (163) Dunning, Jr., T.; Peterson, K. A.; Wilson, A. Gaussian basis sets for use in correlated molecular calculations. X. The atoms aluminum through argon revisited. J. Chem. Phys. 2001, 114, 9244–9253.

70 ACS Paragon Plus Environment

Page 70 of 74

Page 71 of 74

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Theory and Computation

(164) Weigend, F.; Köhn, A.; Hättig, C. Efficient use of the correlation consistent basis sets in resolution of the identity MP2 calculations. J. Chem. Phys. 2002, 116, 3175–3183. (165) Knizia, G.; Werner, H.-J. Explicitly correlated RMP2 for high-spin open-shell reference states. J. Chem. Phys. 2008, 128, 154103. (166) Werner, H.-J.; Knowles, P. J.; Knizia, G.; Manby, F. R.; Schütz, M.; and others, MOLPRO, development version 2017.1, a package of ab initio programs. see http://www.molpro.net. (167) Werner, H.-J.; Knowles, P. J.; Knizia, G.; Manby, F. R.; Schütz, M. Molpro: A general-purpose quantum chemistry program package. WIREs Comput. Mol. Sci. 2012, 2, 242–253. (168) Narancic, S.; Bach, A.; Chen, P. Simple fitting of energy-resolved reactive cross sections in threshold collision-induced dissociation (T-CID) experiments. J. Phys. Chem. A 2007, 111, 7006–7013. (169) Luo, S.; Zhao, Y.; Truhlar, D. G. Validation of electronic structure methods for isomerization reactions of large organic molecules. Phys. Chem. Chem. Phys. 2011, 13, 13683–13689. (170) Rayne, S.; Forest, K. A comparative examination of density functional performance against the ISOL24/11 isomerization energy benchmark. Comput. Theor. Chem. 2016, 1090, 147–152. (171) Zhao, Y.; Lynch, B. J.; Truhlar, D. G. Multi-coefficient extrapolated density functional theory for thermochemistry and thermochemical kinetics. Phys. Chem. Chem. Phys. 2005, 7, 43–52. (172) Liakos, D. G.; Neese, F. Is it possible to obtain coupled cluster quality energies at near density functional theory cost? Domain-based local pair natural orbital coupled cluster vs modern density functional theory. J. Chem. Theory Comput. 2015, 11, 4054–4063. ˇ (173) Jureˇcka, P.; Šponer, J.; Cerný, J.; Hobza, P. Benchmark database of accurate (MP2 and CCSD(T) complete basis set limit) interaction energies of small model complexes, DNA base pairs, and amino acid pairs. Phys. Chem. Chem. Phys. 2006, 8, 1985–1993. (174) Sedlak, R.; Janowski, T.; Pitonak, M.; Rezac, J.; Pulay, P.; Hobza, P. Accuracy of quantum chemical methods for large noncovalent complexes. J. Chem. Theory Comput. 2013, 9, 3364–3374.

71 ACS Paragon Plus Environment

Journal of Chemical Theory and Computation

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

(175) Brauer, B.; Kesharwani, M. K.; Martin, J. M. L. Some Observations on Counterpoise Corrections for Explicitly Correlated Calculations on Noncovalent Interactions. Journal Of Chemical Theory And Computation 2014, 10, 3791–3799. (176) Weymuth, T.; Couzijn, E. P. A.; Chen, P.; Reiher, M. New Benchmark Set of Transition-Metal Coordination Reactions for the Assessment of Density Functionals. J. Chem. Theory Comput. 2014, 10, 3092–3103. (177) Grimme, S. Semiempirical GGA-type density functional constructed with a long-range dispersion correction. J. Comput. Chem. 2006, 27, 1787–1799. (178) Grimme, S.; Antony, J.; Ehrlich, S.; Krieg, H. A consistent and accurate ab initio parametrization of density functional dispersion correction (DFT-D) for the 94 elements H-Pu. J. Chem. Phys. 2010, 132, 154104. (179) Becke, A. D. Density-functional exchange-energy approximation with correct asymptotic behavior. Phys. Rev. A 1988, 38, 3098–3100. (180) Perdew, J. P. Density-functional approximation for the correlation energy of the inhomogeneous electron gas. Phys. Rev. B 1986, 33, 8822–8824. (181) Kats, D. Communication: The distinguishable cluster approximation. II. The role of orbital relaxation. J. Chem. Phys. 2014, 141, 061101. (182) Kats, D.; Kreplin, D.; Werner, H.-J.; Manby, F. R. Accurate thermochemistry from explicitly correlated distinguishable cluster approximation. J. Chem. Phys. 2015, 142, 064111. (183) Janssen, C. L.; Nielsen, I. M. New diagnostics for coupled-cluster and Møller–Plesset perturbation theory. Chem. Phys. Lett. 1998, 290, 423–430. (184) Lee, T. J.; Taylor, P. R. A diagnostic for determining the quality of single-reference electron correlation methods. Int. J. Quantum Chem. 1989, 36, 199–207.

72 ACS Paragon Plus Environment

Page 72 of 74

Page 73 of 74

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Theory and Computation

(185) Datta, D.; Kossmann, S.; Neese, F. Analytic energy derivatives for the calculation of the first-order molecular properties using the domain-based local pair-natural orbital coupled-cluster theory. J. Chem. Phys. 2016, 145, 114101. (186) Frank, M. S.; Schmitz, G.; Hättig, C. The PNO–MP2 gradient and its application to molecular geometry optimisations. Mol. Phys. 2017, 115, 343–356. (187) Dornbach, M.; Werner, H.-J. To be published.

73 ACS Paragon Plus Environment

PNO-LCCSD-F12 / cc-pVTZ-F12

Journal of Chemical Theory and Page Computation 74 of 74 Computational time: < 3 hours 1 2 3 4

• 92 atoms ACS Paragon Plus Environment • 3345 basis functions • 80 cores on 4 nodes