Scalable Electron Correlation Methods. 5. Parallel Perturbative Triples

Subscriber access provided by READING UNIV

Article

Scalable Electron Correlation Methods. 5. Parallel Perturbative Triples Correction for Explicitly Correlated Local Coupled Cluster with Pair Natural Orbitals Qianli Ma, and Hans-Joachim Werner J. Chem. Theory Comput., Just Accepted Manuscript • DOI: 10.1021/acs.jctc.7b01141 • Publication Date (Web): 06 Dec 2017 Downloaded from http://pubs.acs.org on December 12, 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 51 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. 5. Parallel Perturbative Triples Correction for Explicitly Correlated Local Coupled Cluster with Pair Natural Orbitals Qianli Ma and Hans-Joachim Werner∗ Institut für Theoretische Chemie, Universität Stuttgart, Pfaffenwaldring 55, D-70569 Stuttgart, Germany E-mail: [email protected]

Abstract A well-parallelized perturbative triples correction implementation for the pair natural orbital based coupled cluster method PNO-LCCSD(T)-F12 is presented. A composite approach is adopted in addressing the coupling due to off-diagonal Fock matrix elements, in which the local triples amplitudes are iteratively solved using small domains of triples natural orbitals, and a semicanonical (T0) domain correction with larger domains is applied to the reduce the domain errors. This treatment adds only about 20% to the computational cost of (T0) calculations with large domains, and the errors due to the use of small domains in the iterations are very small. In addition, a two-step triple list selection method is applied: First, an initial triple list is generated using LCCSD-F12 pair energy criteria and the (T0) triples energies are computed using small domains. Second, this list is reduced by neglecting triples with small energy contributions, and the final calculation with large domains is only carried out for the reduced list. The cost of the (T) calculation scales asymptotically linear with the molecular size and

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

shows excellent parallelization efficiency up to hundreds of CPU cores. The convergence of the (T) contribution to the relative energies of large molecular systems is carefully tested, and for most of the cases the results obtained with our default thresholds agree within ∼ 0.5 kcal mol−1 with those computed with very tight thresholds. For all tested molecular systems where canonical calculations are still feasible, the PNO-LCCSD(T)-F12 relative energies also agree within 0.5 kcal mol−1 with the canonical CCSD(T)-F12 results using the same F12 treatment. The (T) calculation generally takes 3–5 times the cost of the preceding PNO-LCCSD-F12 calculation, primarily due to the large number of triples required in obtaining accurate relative energies. We find that for large molecular systems the criteria used in previous local triples methods are insufficient, and a much larger number of triples are required than it was assumed so far. Still, using a commodity computer cluster, the PNO-LCCSD(T)-F12 calculation of molecules with ∼ 100 atoms using augmented triple-ζ basis sets can be carried out in a few hours of elapsed time.

1

Introduction

The coupled-cluster method with single and double excitations and a perturbative treatment of the triple excitations, CCSD(T), is widely recognized as the “gold standard” of quantum chemistry. For single reference cases, chemical accuracy can be achieved with the method, provided that basis-set errors are minimized using basis set extrapolation with large basis sets, or by the inclusion of explicitly correlated (F12) terms. 1–31 However, due to the steep scaling of the computational resources with the molecular size, applications of the CCSD(T) method with good basis sets are limited to molecules of 20–30 atoms. The steep scaling of conventional electron correlation methods can be alleviated by using local correlation methods, 32–98 which exploit the short-range character of dynamical electron correlation by using local orbitals and introducing local approximations. The pair natural orbital (PNO) based methods, which were known since the 1960s 99–106 and revived in 2009 by Neese and coworkers, 87,88 have been of interest to several groups including us. 67–70,72–74,87–95,107–110 These methods provide very fast convergence of the correlation energy

2


Page 2 of 51

Page 3 of 51 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


with the number of virtual orbitals per electron pair and enable calculations for large molecules with high accuracy and affordable cost. The combination of explicit correlation (F12) methods with local approximations has been particularly successful. 67–69,74,98,110–114 The F12 terms reduce not only the basis-set incompleteness errors, but also the errors caused by local domain approximations. In addition, the basis-set superposition errors (BSSE), which play an important role in the calculation of large molecules, are largely eliminated by the F12 terms. Recently, our group reported a well-parallelized PNO-LCCSD-F12 implementation. 72–74 The method employs fine-tuned pair and domain approximations and makes it possible to carry out accurate LCCSD-F12 calculations of three-dimensional molecules with more than 100 atoms in a few hours of elapsed time, using several nodes on a small computer cluster. Excellent agreement with the canonical results for small molecules has been demonstrated, and all local approximations have also been tested thoroughly for large molecular systems. However, the missing perturbative triples (T) correction prevented the method from being useful for real-life applications. There are several difficulties in implementing an efficient local (T) method. First, a proper choice of virtual orbitals is not obvious. In the past, projected atomic orbitals 47,48,50 (PAOs), orbital specific virtuals 115 (OSVs), and triples natural orbitals 93,110 (TNOs) have been employed. Among these choices, the TNOs provide the most compact virtual space for the triples amplitudes and are therefore employed in the present work. The TNOs are created in analogy to the PNOs. However, since there is no cheap way of estimating the “triple density matrices”, the pair density matrices for the three pairs associated to a given triple are averaged and then diagonalized to obtain the TNOs. Second, the number of triples 116 in large molecules can be huge, and to our best knowledge an efficient screening method does not yet exist. As a rough treatment, Schütz and Werner 47 used a pair-class criterion, which selects a subset of triples based on the number of “strong” pairs in the three pairs (i j, ik, and jk) associated to a triple i jk. Schmitz and Hättig 110 tested various ways of estimating the triples correction based on CCSD pair energies and proposed a purely empirical formula.

3



Third, the LMOs are not eigenvectors of the Fock matrix and therefore the non-iterative canonical (T) theory cannot be applied. This issue can be overcome by either solving the local (T) amplitude equations iteratively, 48 or using the Laplace transformation technique 110,117–120 to avoid the energy denominators. Ignoring the couplings by off-diagonal Fock matrix elements fi j , which leads to the (T0) or “semicanonical” approximation, normally causes only small errors in relative energies, but in certain cases the errors can reach 0.5 kcal mol−1 or even more. 120,121 Both approaches beyond (T0) require several times more computational time compared to a (T0) calculation. Due to the high cost for computing the (T) energy, efficient parallelization of the method over multiple computer nodes is highly desirable. In the present work we present an accurate (T) implementation on top of our PNO-LCCSDF12 method. We use the TNOs proposed by Riplinger et al. 93 and employ a composite approach in treating the off-diagonal Fock matrix coupling, in which small domains of TNOs are used for the iterative solution of the triples amplitudes and a (T0) domain correction obtained with larger domains is applied. This approach does not bare the very high cost of Laplace-transformed or iterative (T) calculations using large domains, yet the domain errors of the iterative correction are very small. We also implemented a two-step triple list screening process, in which an initial triple list is generated based on LCCSD-F12 pair classes (which in turn depend on the pair energies), and then reduced using the triple energies obtained from a small domain (T0) calculation. This screening provides a compact triple list and a systematic convergence behavior, which is very helpful for accurate calculations on large molecules. Similar to previous work in the series, our aim is to develop a method that runs efficiently on several computer nodes connected with a fast network. This is to be distinguished to moderately parallelized code designed for a single node or for large-scale supercomputers. The (T) program shows nearly ideal speedup with up to more than 300 CPU cores. This drastically reduces the time-to-solution and makes it possible to carry out accurate LCCSD(T)-F12 energy calculations of real-life molecules with ∼ 100 atoms using augmented triple-ζ basis sets within a few hours of elapsed time. It should be mentioned that our method requires rather large amounts of distributed memory, and several computer nodes may be

4


Page 4 of 51

Page 5 of 51 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


required to carry out accurate calculations for large molecules. Alternatively, shared disk files can be used in place of the distributed memory in the single-node case 73 to reduce the memory usage to 4–6 GB per CPU core for most calculations presented in this study. The (T) treatment presented in this paper can also be used for LCCSD(T) calculations without the F12 terms. However, it was found previously 74 that the F12 terms significantly reduce the basis set and domain errors, while adding only 25–40% computational time to a PNO-LCCSD calculation. Relative to a PNOLCCSD(T) calculation this fraction is further reduced to roughly 10%. We therefore only consider the F12 variant in the present work, and recommend using it in all future applications. Testing the accuracy of local (T) calculations is rather difficult. While for small molecules the results can be compared with canonical ones, for larger molecules the accuracy can only be tested by using tighter and tighter thresholds until convergence is reached. Previous works on local (T) methods 47,48,50,93,110,115 primarily used the first method of testing, but this can lead to error estimates that are too optimistic, since many errors caused by local approximations do not show up in small molecules. Furthermore, the total error of the correlation energy is an extensive property, and, as in all quantum chemical calculations, one relies on the fact that most of these errors cancel out if energy differences are computed. In this work we therefore focus on convergence tests of the local approximations on the relative energies for several large molecular systems that were found to be difficult cases in our previous studies. The paper is organized as follows: First, we describe the benchmark systems and computational details in Section 2, since these systems are used to test the approximations described in the following sections. Then, we describe briefly the local (T) theory and our working equations in Section 3. The two major local approximations, namely the truncations of virtual orbitals and of the triple list, are discussed in Sections 4 and 5, respectively. Technical details on the implementation are provided in Section 6, which is followed by some additional benchmark calculations in Section 7. A summary in Section 8 concludes the paper.

5



2 Benchmark Systems and Computational Details In order to investigate how the local approximations affect the perturbative triples correction in large molecules, we primarily used the same three benchmark systems that were employed previously in testing our PNO-LCCSD-F12 method. 73,74 All the three systems are difficult for local correlation methods, due to strong dispersion contributions to the reaction energies. These include the reaction energies of the “Auamin” and “Isomer4” reactions (see Figure SI of the supporting information), and the binding energy of the coronene dimer. The Auamin reaction is the dissociation of a gold(I)-aminonitrene complex (AuC41 H45 N4 P, for simplicity denoted “Auamin”) taken from ref 122. It represents a transition metal compound in which local approximations may behave differently from main group element compounds. The Isomer4 reaction is the fourth reaction from the ISOL24 test set, 123 and is the largest and most difficult case of the whole set, due to the drastic change of molecular structure in the reaction. The coronene dimer (in a parallel displaced stacked configuration, C2C2PD) is taken from the L7 test set. 124 It represents a molecular system with delocalized electronic structure and strong dispersion interactions, and has also been used as a benchmark system for the DLPNO-CCSD(T)-F12 method. 98 Furthermore, we present in Section 7 results for the “Androstendione reaction” and the “Testosterone reaction” used as benchmark reaction in the PAO-LCCSD(T)-F12 method developed earlier in our group. 113 In addition, we show the computed intermolecular interaction energies of the guanine–cytosine (G–C) dimers from the JSCH-2005 benchmark set, 125 which have been used in testing the OSV-L(T) method of Schütz et al. 115 Finally, for a comparison with the canonical CCSD(T)-F12 method, we used the FH benchmark set of Friedrich and Hänchen 126 and several dimers from the S22 set of intermolecular interactions. 125 The results for intermolecular interactions are not counterpoise (CP) corrected unless explicitly marked. The occupied molecular orbitals are localized using the intrinsic bond orbital (IBO) method. 127 Unless otherwise noted, the cc-pVTZ-F12 basis set (denoted VTZ-F12 in this work) of Peterson et al. 128 was used. An exception is the calculation on the coronene dimer, in which we used the aVTZ basis set, which consists of the cc-pVTZ basis sets 129 for H atoms, and the aug-cc-pVTZ 6


Page 6 of 51

Page 7 of 51 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 130 for other atoms. The diffuse polarization functions are important for an accurate discription of the dispersion interaction, 131 which makes the aug-cc-pVTZ basis set more appropriate than the VTZ-F12 basis set for the study of intermolecular interactions. We were not able to perform Hartree–Fock calculation of the coronene dimer with the full aug-cc-pVTZ or VTZ-F12 basis sets due to linear dependence problems. For the gold atom, the Stuttgart effective core potential ECP60MDF along with the aug-cc-pVTZ-PP basis set 132,133 was used. Local density fitting approximations were employed for the calculation of all two-electron integrals in the local correlation calculations. 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 134 by adding for each angular momentum another shell of diffuse functions in an even tempered manner. For all other integral types the aug-cc-pVTZ/MP2FIT basis set 135 was used. The aug-cc-pVTZ/JKFIT basis set was employed as RI basis for the F12 correction. The triple-ζ auxiliary basis sets were also used for double-ζ calculations. The F12b approximation and Ansatz 3*A with fixed amplitudes 1,2 and a geminal exponent 12,25 of 1 a−1 of the Hartree-Fock 0 were used in all F12 calculations. The CABS singles correction

energies was included in all calculations. Currently, this is computed in a separate program without local approximations, and the time for computing this correction is not included in the times presented in this paper. All benchmark calculations were performed using the development version 136 of Molpro. 137 Each node of our computer cluster has two 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 using all 20 cores per node, and one MPI process per core was created.

3 Theory In the present work we consider the closed-shell case and adapt the orbital-invariant local (T) theory developed by Schütz and Werner 47,48 to use TNOs. Here we only briefly outline the working

7



Page 8 of 51

equations, assuming that the TNOs are pseudocanonical, i.e., the external–external block of the Fock matrix for each triple is diagonal. In the following, indices i, j, k, l, ... denote occupied LMOs, ai j , bi j , ... denote PNOs, and ai jk , bi jk , ... denote TNOs. Unless otherwise indicated in the equation, we omit the domains of the TNOs and assume that

a, b, c, d ∈ [i jk]TNO

(1)

where [i jk]TNO represents a domain of pseudocanonical TNOs for the triple i jk. Using this notation, the perturbative triples (T) correction to LCCSD-F12 is given by X

X

i jk i jk ˜ i jk Wabc + Vabc T abc

(2)

  X  X  jk  lk    = Pˆ iabc (ia|bd)T¯ dcjk − (ia| jl)T¯ bc 

(3)

E (T ) =

(2 − δi j − δ jk )

i≥ j≥k

abc

with intermediates

i jk Wabc

d∈[i jk]TNO

l

and i jk Vabc =

1 ˆ i jk Pabc (ia| jb)t¯ck 2

(4)

where T¯ dcjk is a projected LCCSD-F12 doubles amplitude defined as T¯ dcjk =

X

[i jk, jk] jk [ jk,i jk] S da T ab S bc

(5)

a,b∈[ jk]PNO

[i jk, jk] with S da = hdi jk |a jk i being the overlap integral between a TNO and a PNO, and [ jk]PNO repre-

senting a domain of PNOs for the pair jk. t¯ck is a projected LCCSD-F12 singles amplitude defined similarly. Restricting the summation over d to d ∈ [i jk]TNO in eq 3 is a projection approximation.

8


Page 9 of 51 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


Since the LCCSD doubles amplitudes T dcjk are restricted to c, d ∈ [ jk]PNO , the exact form is X

i jk Wabc =

[i jk, jk] (ia|bd)T dejk S ec + ···

(6)

d,e∈[ jk]PNO

In the present method we assume that [i jk]TNO spans the three PNO domains [i j]PNO , [ik]PNO , and [ jk]PNO in order to avoid 3-external integrals of the type (iai jk |bi jk d jk ) and to exploit the permutation symmetry of the 3-external integrals (iai jk |bi jk di jk ). This treatment becomes exact in the limit of complete domains, and the errors should therefore decrease with increasing domain size. The contravariant amplitudes are defined as i jk i jk i jk i jk i jk i jk i jk T˜ abc = 4T abc − 2T acb − 2T bac + T bca + T cab − 2T cba

(7)

and the permutation operator is jk jk jik jki ji Pˆ iabc (·) = ·iabc + ·ikacbj + ·bac + ·bca + ·kicabj + ·kcba

(8)

i jk ik j i jk The triples amplitudes have the usual six-fold symmetry, i.e., T abc = T acb , etc. T abc must satisfy

the following system of linear amplitude equations: jk jk i jk i jk i jk Riabc = Diabc T abc − Zabc + Wabc =0

(9)

jk where Diabc is the energy denominator

jk Diabc = faa + fbb + fcc − fii − f j j − fkk

(10)

and i jk Zabc =

X l,k

i jl fkl T¯ abc +

X l, j

9

ilk f jl T¯ abc +

X l,i


l jk fil T¯ abc

(11)


i jl T¯ abc represents a projected triples amplitude defined similarly as in eq 5. i jk In the canonical case the Fock matrix is diagonal, which implies Zabc = 0, and the triples

amplitudes can be readily obtained non-iteratively from eqs 3 and 9. When local orbitals are used, the triples amplitudes need to be solved iteratively. The commonly used (T0) (or semicanonical) i jk ≈ 0. The (T0) approximation normally recovers about 97% of the approximation assumes Zabc

triples correction and the errors are largely canceled out in many cases when relative energies are considered. 47,48 Still, in some difficult cases this approximation can lead to errors of ∼ 0.5 kcal mol−1 or more. 120,121 The major disadvantage of the iterative solution of eq 9 is that the triples amplitudes have to be stored and updated in each iteration. Using PAOs or OSVs, this easily creates a storage bottleneck and leads to very high computational cost. As shown in Section 4, this bottleneck is largely removed when the TNOs are used as virtual orbitals, in particular since relatively small domains of TNOs are sufficient in the iterative solution of eq 9. Using the Jacobi method, 3–4 iterations are typically needed to converge the (T) correction to within 5×10−5 Eh for the molecules considered in this study. A blocked Gauss–Seidel method can significantly improve the speed of convergence, 48 but efficient parallelization of the method is extremely difficult. Schütz 48 also proposed a perturbative treatment of the internal–internal coupling called (T1). This is a good approximation that corrects for the majority of the error caused by the (T0) approximation. The (T1) working equations are identical to performing only one iteration in the solution of eq 9 using the Jacobi method and the (T0) amplitudes as the starting guess. However, in our present implementation, this does not lead to significant savings of CPU time, memory, or disk space compared to the full iterative method. An alternative approach of addressing the internalinternal coupling is the Laplace transformation technique, 117,118 and this has been implemented in some local or fragmentation-based (T) methods. 110,120 We will compare the Laplace-transformed (T) and iterative (T) from a cost perspective in Section 7. The (T) treatment does not include F12 terms. To improve the basis-set convergence of the triples energy, Knizia et al. 25 proposed a scaled triples treatment, in which the (T) correction is

10


Page 10 of 51

Page 11 of 51 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


scaled by the ratio of the LMP2-F12 and LMP2 correlation energies. Given that the F12 terms also reduce the domain errors in local methods, it would be interesting checking whether the scaled triples approach improves the convergence of the triples energy with respect to the domain sizes. By default our PNO-LCCSD-F12 program uses rather small PNO domains for the coupled-cluster calculation with an LMP2-F12 domain correction. 73,74 In this case, the scaled triples correction is computed by multiplying the total (T) correction of a molecule by the ratio of the LMP2-F12 and LMP2 correlation energies from the small domain calculations. Scaled triples are denoted (T*) in this work.

4 Selection of Virtual Orbitals In order to achieve linear scaling with the molecular size, two kinds of local approximations are applied: triples i jk with insignificant contribution to the (T) correction are neglected or treated with a less accurate approach, and the excitation space for each triple is restricted to a domain of local virtual orbitals. In this section we consider the latter approximation. The selection of the triples will be discussed in Section 5. In the present work we adopt the TNOs proposed by Riplinger et al. 93 Following the PNOLCCSD-F12 method, we use a transformation sequence PAO→OSV→TNO to reduce the domain sizes, and the PAO, OSV, and PNO domains are defined as described in detail in our previous work. 73,74 A PAO domain [i]PAO is determined by an atom list obtained by extending the primary domains (which include all atoms that have IBO partial charges ≥ 0.2) by including all atoms that are connected by ≤ IEXT bonds or within a distance REXT to any atoms in the primary domain. By default, we use IEXT = 2 and REXT = (2 × IEXT + 1) a0 . OSVs are then made for each LMO and truncated by an occupation number threshold T OSV (defaults to 10−9 ). With the default options, the OSV domains in the (T) calculation are the same as in the LCCSD-F12 calculation. For each triple i jk, we make a triple OSV domain [i jk]OSV , which is constructed from the union of the OSV domains for the three LMOs. The three pair density matrices Di j , Dik , and D jk in the

11



Page 12 of 51

pseudocanonical OSV triple domains are computed by transforming the PNO-LCCSD amplitudes to the OSV triple domains. These pair density matrices are averaged to obtain the “triple” density matrix 1 Di jk = (Di j + Dik + D jk ) 3

(12)

Diagonalizing Di jk yields the OSV→TNO transformation matrices with the eigenvalues being the TNO occupation numbers. The TNOs are truncated with an occupation number threshold T TNO and made pseudocanonical. Considering that the (T0) approximation gives ∼ 97% of the correlation and the resource usage of the iterative solution of the triples amplitudes scales cubically with the TNO domain sizes, we use a composite approach of running the iterative (T) calculations with smaller domains and applying a (T0) domain correction obtained with large domains. The final (T) contribution is computed as (T) (T0) (T0) E (T) = Esmall − Esmall + Elarge

(13)

(T) (T0) where two different TNO occupation number thresholds, T TNO and T TNO , are used for making the

small and large TNO domains, respectively. This approach is similar to the one we used for the PNO-LCCSD method, in which a small domain LCCSD calculation is performed and an LMP2 domain correction is applied. In practice, the (T0) calculations of the large and small domains can be done in the same loop over the triples, in which the two-electron integrals are first transformed from the PAO basis to the large-domain TNO basis. After the (T0) contribution in the large domain is computed, the integrals are transformed to the small-domain TNO basis for further calculations. The latter integral transformation is cheap and the extra cost of the small domain (T0) calculation is relatively small. In our PNO-LCCSD-F12 method, the PNOs are truncated using two criteria. In addition to CC the standard occupation number criterion with a threshold T PNO (defaults to 10−7 ), we include CC additional PNOs until the semicanonical PNO-LMP2 pair energy reaches a certain fraction (T en ,

defaults to 0.99) of the semicanonical OSV-LMP2 pair energy. The latter criterion ensures that suf-

12


Page 13 of 51

ficient numbers of PNOs are included for weaker pairs. This significantly improves the accuracy of the PNO-LCCSD-F12 calculations. Unfortunately, it is not straightforward to use such a completeness criterion for the TNO domains, and for weak triples the number of TNOs selected with a moderate T TNO can drop to zero in large molecules. To ensure that such triples are still included in our calculations, we restrict the minimal number of TNOs per triple to be nine, based on our previous observations that three PNOs per LMO are normally sufficient to describe the correlation energy of a distant LMO pair. 71 5 (T) correction / 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


4 3 2

(T) (T), TTNO varying

3

(T0) TTNO

1

1 0 10−6

Isomer4 reaction 10−7

−4

(T) varying

10−8 TTNO

10−9

10−10

0 10−6

(T0) TTNO varying

10−8 TTNO

(T0) TTNO varying

(T) (T0)

−5

Auamin reaction 10−7

coronene dimer (CP)

−3

(T0)

2

(T0)

(T)

−2

(T) (T), TTNO varying

10−9

10−10

−6

10−6

(T)

(T), TTNO varying 10−7

10−8 TTNO

10−9

10−10

(T) Figure 1: (T) corrections to reaction energies as a function of the TNO domain thresholds T TNO (T0) (T0) (T) and T TNO . The black lines show the results obtained by varying T TNO with T TNO = 10−7 , and the (T) (T0) red lines show the results obtained by varying T TNO with T TNO = 10−9 .

(T) The computed reaction energies for the three test systems with various choices of T TNO and (T0) T TNO are plotted in Figure 1. The dashed black lines show the (T0) corrections as a function of (T0) T TNO . We see that the convergence is smooth, although a very tight threshold of 10−9 is necessary

to converge the reaction energies to ∼ 0.1 kcal mol−1 . The solid black lines show the iterative (T) (T) correction when T TNO is fixed at 10−7 . We see that (T0) underestimates the reaction energy of

the Isomer4 reaction by ∼ 0.5 kcal mol−1 , but the errors are smaller for the other two benchmark systems. This is in agreement with the findings of previous studies that (T0) normally behaves well for relative energies but is insufficient in certain cases. The red lines in Figure 1 show the (T) (T0) reaction energies obtained by varying T TNO while keeping T TNO = 10−9 fixed. We see that the (T) (T) reaction energies are insensitive to T TNO and the errors are very small with T TNO = 10−7 , even

for the Isomer4 reaction. The storage requirement and computational cost for the iterations scale

13



(T) cubically with the the TNO domain size and increases very fast with reduced T TNO . We were (T) therefore unable to perform the calculation with T TNO = 10−9 for the coronene dimer. However, (T) the binding energies should already be well converged with T TNO = 10−8 , since those obtained with (T) (T) T TNO = 10−7 and T TNO = 10−8 differ only by 0.03 kcal mol−1 .

(T) correction / 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 51

4.5 4.0

3.5 CC Ten

= 0.99

3.0

3.5

2.5

1.5 Isomer4 reaction

1.5 10−7

10−8 CC TPNO

10−9

1.0 −6 10

CC = 0 Ten

−5.0

2.0

CC = 0 Ten

2.0 10−6

−4.5

2.5

3.0

coronene dimer (CP)

−4.0

CC = 0.99 Ten

−5.5

CC = 0 Ten

Auamin reaction 10−7 CC TPNO

10−8

−6.0

10−9

10−6

CC = 0.99 Ten

10−7 CC TPNO

10−8

10−9

CC Figure 2: (T) corrections to reaction energies as a function of the PNO domain threshold T PNO . The solid and dashed lines show the (T) and (T*) corrections, respectively. Default parameters for all other local approximations are used.

CC The accuracy of the (T) correction also depends heavily on the PNO domain threshold T PNO

used in the LCCSD-F12 calculations, which affects the accuracy of the LCCSD singles and doubles amplitudes. The LCCSD amplitudes in turn determine the triples energy. This is illustrated in CC Figure 2. The solid black lines in the figure show the (T) correction with respect to T PNO without

the energy completeness criterion. We can see the convergence is rather slow, and the results CC CC obtained with T PNO = 10−7 and T PNO = 10−9 for the three reactions differ 0.6–0.9 kcal mol−1 .

These errors can be compared with the difference between the LCCSD reaction energies obtained with these thresholds, which range from 2.0–3.4 kcal mol−1 . While the errors in the LCCSD energies can be largely eliminated using either the F12 variant or a semicanonical LMP2 domain correction, these treatments do not improve the accuracy of the LCCSD amplitudes and of the (T) energies. The PNO domain errors are significantly reduced by applying the energy completeness criterion in the domain selection, as shown by the red lines in Figure 2. Still, our default PNO thresholds lead to errors of 0.2–0.3 kcal mol−1 in the (T) corrections for the three reactions relative CC to T PNO = 10−9 results. These errors are about twice as large as those in the LCCSD-F12 reaction

14


Page 15 of 51

energies and are a major part of the total errors in the (T) corrections. The dashed lines in Figure 2 show the F12 scaled triples correction, which unfortunately does not help in reducing the PNO domain errors. 4.5 (T) correction / 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


3.5 (T*) IEXT=IEXT T

4.0

−4.5 (T*) IEXT=IEXT T

3.0

3.5

2.5

(T) IEXT=2

(T) IEXT=IEXT T 3.0

2.0

2.5

1.5

2.0

Isomer4 reaction 0

1

2 IEXT T

3

1.0

(T) IEXT=IEXT T

−5.0

(T) IEXT=2

(T) IEXT=2 −5.5

(T) IEXT=IEXT T

−6.0 Auamin reaction 0

1

2 IEXT T

3

−6.5

(T*) IEXT=IEXT T coronene dimer (CP) 0

1

2

3

IEXT T

Figure 3: (T) corrections to reaction energies as a function of the PAO domain parameter IEXT_T. The black line shows the results with IEXT = 2 (default in LCCSD-F12). The red and blue lines show the (T) and (T*) corrections, respectively, obtained by varying IEXT and IEXT_T simultaneously. The dashed red line in the leftmost panel shows the (T) corrections (IEXT = IEXT_T) obtained with the smaller VDZ-F12 basis set.

Next we consider the effect of the PAO domain sizes on the (T) corrections. The results are shown in Figure 3. Here we introduce new parameters IEXT_T and REXT_T [again using by default REXT_T = (2 × IEXT_T + 1) a0 ] to control the PAO domain sizes, which have the same meaning as IEXT and REXT, except that they only affect the triples calculations. The PAO domains in LCCSD-F12 calculations are still controlled by IEXT, independent of IEXT_T. When the two parameters are different, the TNOs are generated directly from domains of PAOs, without using OSVs as intermediates. We observe from Figure 3 that using larger PAO domains for (T) than in LCCSD-F12 is unnecessary. When using IEXT = 2 for the LCCSD calculation, IEXT_T = 1 in the (T) calculation already leads to rather small errors as compared with calculations with larger PAO domains in (T), e.g. 0.2 kcal mol−1 for the Isomer4 reaction. On the other hand, the (T) correction is sensitive to IEXT, which controls the PAO domain size in the LCCSD calculation, and the convergence is slow. As pointed out in our previous work, 73,74 this slow convergence is at least partly caused by basis set superposition errors, which are largely removed in local calculations with small PAO domains, but reintroduced with increasing the PAO domains. While the F12 terms 15



largely eliminate this effect in LCCSD calculations, the quality of the amplitudes is not improved and the slow convergence with the PAO domain sizes still appears in the (T) correction. The PAO domain errors and the effects of BSSE are larger when a smaller basis set is used. In that case, the dependence of the triples energy with respect to the PAO domain size is stronger, as shown by the dashed red line in Figure 3. Since the (T) correction is relatively small compared to the LCCSD correlation energy, the overall errors from using domains of PAOs should still be a small fraction of 1 kcal mol−1 in most cases. The scaled triples treatment slightly improves the convergence with IEXT for the Isomer4 reaction but barely helps in the other two systems. It should be noted that the truncations of PAOs, PNOs, and TNOs domains all lead to underestimated (T) contributions. Triple screening, as presented in the following section, will also lead to errors in the same direction. Therefore, in order to achieve chemical accuracy (∼ 1 kcal mol−1 ), all thresholds need to be chosen tight so that the errors from each individual approximation are a small fraction of 1 kcal mol−1 . From combined cost and accuracy considerations, we choose IEXT (T0) (T) = IEXT_T = 2, T TNO = 10−9 , and T TNO = 10−7 as the defaults, and these thresholds are used in the

following tests of the triple list truncation.

5

Selection of the Triples

Selecting a proper list of triples to be included in the (T) calculation is critical to the efficiency and the scalability of a local (T) method. To our best knowledge, there is no cheap and reliable method for estimating the triples energy corrections. In the PAO-LCCSD(T) method of Schütz and Werner, 47,48,62 an option “triptyp” is used to control the list of triples. For a triple i jk to be included in the list, the three related pairs, i.e., i j, ik, and jk, shall have at most “triptyp” close pairs and the rest (if any) must be strong pairs. Various tests have shown that triptyp = 2 leads to very small errors. Riplinger et al. 93 further recommend including only triples where at least two of the three triples are strong, which corresponds to triptyp=1. We will show below that this choice can lead to large errors unless very tight pair selection criteria are used.

16


Page 16 of 51

Page 17 of 51 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


Figure 4: Correspondence of (T0) triple energies from small (IEXT_T = 0, T TNO = 10−7 Eh ) and large domain (IEXT_T = 2, T TNO = 10−9 Eh ) calculations. The colors of the points indicate the number of strong pairs (LCCSD-F12 pair correlation energy greater than 10−4 Eh ) among i j, ik, and jk for a given triple i jk. Triples with (T0) contributions less than 10−10 Eh from the small domain calculation are excluded. We show in Figure 4 in different colors the energies of triples classified by the strong pair (T) criterion. Here we use an LCCSD-F12 pair energy threshold T close to determine strong pairs.

Similar figures for the Isomer4 reactant and product are shown in Figure S2 of the Supporting (T) Information. By default we use T close = 10−4 Eh , which is the same as the close pair threshold in the

previous work 73,74 and the DLPNO-CCSD(T) method with the “NormalPNO” preset of options. 138 But note that here we use LCCSD-F12 pair energies, not PNO-LMP2 energies as employed to determine the pair classes in the PNO-LCCSD-F12 method. We see that there is some correlation between the triples energies and the number of strong pairs associated with each triple. However, this does not provide an ideal triple screening method. With triptyp = 1, many triples with energy contributions in the 10−5 Eh range are screened out, which may leads to rather large errors (see below). On the other hand, with triptyp = 2, a huge number of triples with contributions ≤ 10−8 Eh are included and the calculation can be become rather expensive. In order to generate a compact yet accurate triple list, we use a two-step triple selection method. First, an initial list is generated using the pair class criterion with triptyp = 2. The rest of the triples

17



are considered “distant”, and their contributions to the correlation energy are neglected. A (T0) calculation using TNOs generated from the small “primary” PAO domains and truncated using a loose T TNO threshold of 10−7 is then performed. Using the resulting triple energies, the nondistant triples are categorized as “strong” or “weak” with an energy threshold T trip . The large domain (T0) calculation and iterations are then performed only for the strong triples, while the energy contributions of the weak triples are simply taken from the small domain calculation. Figure 4 also shows the correlation of the (T0) contributions from the small domain and large domain calculations. We see that the small domain calculation typically underestimates the (T0) correction, but still provides very good energy estimates for the triple screening.


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 51

(T) triptyp=2,3

4

(T0) triptyp=2,3 (T) triptyp=1

3

2

1 10−5

10−8

(T0) triptyp=2,3

10−9

(T0) triptyp=1

−1 10−5

10−7 Ttrip / Eh

10−8

(T0) triptyp=1 (T0) triptyp=2,3 (T) triptyp=2,3

−6

Auamin reaction 10−6

(T) triptyp=1

−2 −4

(T) triptyp=1

0

Isomer4 reaction 10−7 Ttrip / Eh

2 1

(T0) triptyp=1

10−6

(T) triptyp=2,3

3

10−9

coronene dimer (CP) 10−5

10−6

10−7 Ttrip / Eh

10−8

10−9

Figure 5: (T) contributions to the reaction energies as a function of T trip for different choices of triptyp. The lines for triptyp = 2 (solid lines) and triptyp = 3 (dashed lines) almost fully overlap and are therefore plotted in the same color. We show in Figure 5 the (T0) and (T) energy corrections to the reaction energies as a function for T trip . It can be clearly seen that triptyp = 2 is the proper choice for making the initial triple list. At the limit of small T trip , using triptyp = 1 leads to very large errors, ranging from ∼ 1 kcal mol−1 for the Isomer4 reaction to more than 3 kcal mol−1 for the coronene dimer. On the other hand, including triples where none of the three associated pairs is strong (triptyp = 3) has a negligible effect. Depending on the reaction, the convergence of the (T0) correction with respect to T trip can be rather slow, and a tight threshold of 10−8 Eh may be needed to converge the relative energy to ∼ 0.1 kcal mol−1 . We found this rather disappointing, considering the fact that the results include the (T0)

18


Page 19 of 51

correction computed with small domains for weak triples, and the variation of the (T0) energy as a function of T trip purely reflects the domain errors of these weak triples. On the other hand, the difference between (T) and (T0) does not change much for thresholds below 10−7 Eh . In the current work we have decided to use by default T trip = 10−8 Eh , to keep the accuracy as high as possible at still affordable cost. Rather significant savings are possible with a threshold of T trip = 10−7 Eh , if the somewhat larger errors can be tolerated. It should be noted, however, that due to the slow convergence of the (T0) energy with the TNO domain sizes the prescreening with small domains replacements does PSfrag not save much computational time. At slightly larger expense (typically 20%), one could do

the large domain (T0) calculation for all triptyp = 2 triples, and select only a smaller number of triples and smaller domains for the iterations. This will be explored in future work. percent of (T) correction of the best calc.

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


100 (T) triptyp=2

99 98

(T) triptyp=1

97

(T0) triptyp=2

96

(T0) triptyp=1

95 94 10−5

Auamin molecule 10−6

10−7

10−8

10−9

Ttrip / Eh

Figure 6: Percentage of the (T0) and (T) correlation energy contributions relative to the results of (T) calculations with triptyp = 3, T trip = 10−9 for the Auamin molecule. The VTZ-F12 basis set and default settings for all domain approximations are used. In several previous studies on local (T) corrections, the percentage of the canonical (T) energy correction recovered was used as the criterion in tuning the local approximations. 47,48,93 This can be quite misleading for calculations of large molecules. As an example, we show in Figure 6 the percentage of the (T) energy contribution for the Auamin molecule relative to the calculation with a very large triple list (triptyp = 3, T trip = 10−9 ). This can be compared with the middle panel of Figure 5. We see that with T trip = 10−9 Eh , the triptyp = 1 calculation yields 97.8% of the (T) correction to the correlation energy compared with the best calculation, but only 10% of the 19



Page 20 of 51

Sfrag replacements

triples contribution to the reaction energy! In our previous work 73 it has been shown that longrange interactions strongly affect the dissociation energy of Auamin. As already pointed out by Schütz, 48 the inclusion of triples where one or two pairs are weak are important to describe longrange Axilrod–Teller interactions. With triptyp=1, these triples are missing, and this is probably the main reason for the error. Since for the Auamin molecule the (T) contribution to the correlation energy is ∼ 0.45 Eh , one percent of the (T) correlation amounts to ∼ 2.8 kcal mol−1 , and thus the additional 2% of the triples energy recovered with triptyp=2 can indeed have a large effect on the dissociation energy. In fact, the same triple list truncation thresholds lead to less errors for the two dissociated fragments of Auamin, which amounts to only 1.5% and 1.0% of the correlation energy, respectively. 4 (T) correction / 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

3

(T) (T0)

3

coronene dimer (CP)

−3

(T0)

2

−4 1

2 Isomer4 reaction 10−2

−2

(T)

10−3

(T) Tclose

10−4

10−5

0 10−2

/ Eh

−5 Auamin reaction 10−3

(T) Tclose

10−4

10−5

/ Eh

−6

10−2

(T0) (T) 10−3

(T) Tclose

10−4

10−5

/ Eh

(T) Figure 7: (T) contributions to the reaction energies as a function of T close using triptyp = 2. Default parameters are used for all other local approximations.

Clearly, whether the initial triple list is good enough also depends on the definition of strong (T) pairs. We show in Figure 7 the (T) correction as a function of T close . It can be seen that a threshold

of 10−4 Eh is necessary to converge the (T) correction. This is a rather tight threshold, since it will include most pairs i j with orbital distances between the charge centers of the LMOs i and j up to 7 a0 . Triples where at least one of the three associated pairs is distant in the initial triple list are neglected. To test this approximation, we have performed calculations in which the distant pair threshold T dist in LCCSD-F12 was decreased from 10−6 Eh to 10−7 Eh , which leads to larger initial triple lists. It was found that this change affects the (T) contribution to the reaction energies of all three systems by less than 0.03 kcal mol−1 . 20


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

PSfrag replacements


0.50 (T) correction / kcal mol−1

Page 21 of 51

0.25 (T) w/o weak triple amplitudes

0.00 −0.25 −0.50 −0.75

(T0)

−1.00

(T) with weak triple amplitudes

−1.25 10−5

10−6

10−7 Ttrip / Eh

10−8

10−9

Figure 8: (T) contributions to the binding energy of the stacked benzene dimer as a function of T trip . The red and blue lines show the results obtained with and without the inclusion of weak triple amplitudes, respectively, in the strong triple residuals. The black line shows the results from (T0) calculations. The aug-cc-pVTZ basis set and default settings for all domain approximations are used. We note that when using an energy threshold, it is critical to save the semicanonical triples amplitudes of the weak triples in the prescreen process and include them in the computation of the (T) residuals of strong triples (eqs 9 and 11). This treatment is an analogy of the “keepcls = 1” option used in the PAO-LCCSD method, 49,113 in which LMP2 amplitudes of weak pairs are included in the strong pair LCCSD residuals to reduce the error caused by the pair approximations. The effect of the weak triple amplitudes is illustrated in Figure 8 for the CP corrected binding energy of the benzene dimer. In this example, the (T0) approximation leads to very small errors and can be used as a reference of the triple truncation errors. The convergence with respect to T trip is not smooth and errors are huge if one ignores the contributions of weak triple amplitudes in the strong triple residuals. This issue is resolved by including the fixed semicanonical weak triples amplitudes. Then the convergence with respect to T trip is smooth and is very similar to that of (T0) calculations. Given that the TNO domains are small for weak triples, this treatment does not significantly increase the storage requirements of the triples amplitudes and is applied in all calculations presented in this work. On the other hand, we find that considering the off-diagonal coupling of strong pairs only, i.e., assuming fi j , 0 in eq 11 only if i j is a strong pair, leads to less than 0.01 kcal mol−1 errors for all three test systems, and this is exploited by default to reduce the

21



communication cost in the iterations. Isomer4 reaction

4.10


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 22 of 51

Auamin reaction

2.85

−5.50 −5.55

4.05

2.80

4.00

2.75

−5.65

2.70

−5.70

2.65

−5.75

3.95 3.90 10−3

10−4

10−5

Tclose / Eh

10−6

10−3

−5.60

10−4

10−5

Tclose / Eh

10−6

10−3

Coronene dimer (CP) 10−4

10−5 Tclose / Eh

10−6

Figure 9: (T) contributions to the reaction energies as a function of T close . Default parameters are used for all other local approximations.

We now consider the effect of the pair approximation in the LCCSD-F12 calculation on the (T) correction, and the results are shown in Figure 9. In the LCCSD-F12 method, we catagorize the nondistant pairs to strong, close, and weak pairs and use different treatments from full to linearized LCCSD. By default, the weak pair threshold is T weak = T close /10. From Figure 9 one can see that the (T) correction does not change significantly with T close , indicating that the approximate LCCSD treatments for close and weak pairs provide sufficiently accurate amplitudes for the (T) calculations. With the default T close = 10−4 Eh , the errors on the (T) correction from close and weak pair approximations are less than 0.05 kcal mol−1 for all three reactions. We note that the pair approximations in LCCSD-F12 calculations lead to an overestimation of the correlation energy, both in LCCSD and indirectly through the amplitudes in the (T) correction. This provides an error cancellation of the domain errors. The overestimated pair energies also lead to a slightly larger triple list, and therefore reduce the triple list truncation errors. However, the effect of these error cancellations is rather minor and the primary sources of error with our defaults are the domain approximations.

22


Page 23 of 51 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 Parallelization and Integral Transformation Similar to our previous PNO-based methods, 68,69,73,74 the (T) program is parallelized using MPI and the GlobalArrays (GA) toolkit. 139 The latter provides efficient remote memory access interfaces and can simplify the programming of many procedures. In our current implementation a (T) calculation consists of three steps: the triples screening using minimal TNO domains, the (T0) calculation using large domains, and the iterative solution of the triples amplitudes. We have adopted different parallelization schemes for the (T0) calculation and the iterations. The triples screening shares the same program with the (T0) calculation. In a typical (T0) calculation, the major part of the computational time is spent in the calculation of the 2-electron integrals, in particular the 3-external integrals in eq 3. We used the local density fitting technique as described in ref 74, using the union of the fitting domains of the three LMOs for each triple as the fitting domain. At the beginning of the (T0) calculation all required 3-index PAO integrals are computed and stored in GAs. Even though most of these integrals are already computed in the preceding LCCSD-F12 calculation, they are recomputed in order to minimize the GA usage during the LCCSD-F12 iterations, in which the GA memory is needed to accommodate the large amounts of PNO overlap matrices. For each triple, we transform the 3-index PAO integrals to the TNO basis before assembling the integrals and computing the (T0) correction. Considering the large number of triples involved, we use a dynamic parallelization scheme in the (T0) calculations. The triples are distributed in batches on a first-come first-serve basis, and each processor takes as many triples as scratch memory allows to perform the PAO→PNO integral transformation, while buffering the PAO integrals so that a block of integrals is never retrieved from GA twice within the batch. The 3-index PNO integrals are then computed for the batch of triples and dumped to local disks. This is followed by a loop over the triples in the batch to compute the (T0) correction, loading the PNO integrals from disk. The maximum number of triples in a batch is also limited depending on the number of triples remaining and the number of CPU cores to reduce the waiting time. Also a fallback algorithm is implemented for triples involving large PAO or TNO domains, for which buffering the PAO integrals in local memory 23



is not possible. This parallelization strategy significantly reduces the communication cost of the integral transformations and leads to excellent parallelization efficiency in the (T0) calculation (see Section 7). Surprisingly, we found that using the small primary fitting domains in assembling the integrals required by the (T) calculations leads only to negligible errors (≤ 0.05% of the triples correction to the correlation energy or ≤ 0.01 kcal mol−1 for the relative energies we tested). With the default parameters, this choice implies that the fitting domains contain two layers of atoms less than the PAO domains. We also observed similar small effects of the fitting domains in our PNO-LMP2F12 method (see Figure 10 of ref 69). The effect of the fitting domains on the PNO-LCCSD-F12 method is more involved, given that various projection approximations are made. In the present work, we use the extended fitting domains for the PNO-LCCSD-F12 calculation and then the primary fitting domains in the (T) calculation. This can reduce the computational time of the (T0) calculation by a factor of three as compared to using the extended fitting domains. In the iterations we use a static parallelization scheme, in which each processor is responsible for computing the residuals of triples for which it has computed the (T0) correction. After the (T0) calculation of each triple, the PAO→TNO transformation matrices, (T0) amplitudes, as well as i jk i jk the intermediates Wabc and Vabc are dumped to local disk. At the beginning of the iterations, GAs

are created to hold the transformation matrices and the amplitudes, while the other intermediates are kept on local disk. The static parallelization avoids the communication of these intermediates entirely. To avoid having multiple sets of triples amplitudes in GAs, within an iteration, each processor computes the updated triples amplitudes and dumps them to local disk, and these amplitudes are uploaded to the GAs at the end of the iteration. Even though the cost of the (T0) calculations is not directly related to the cost of the iterations, the static parallelization seems to provide a reasonable efficiency. We note that the (T) amplitudes are independent of the CCSD singles amplitudes i jk i jk intermediates. Vabc can be computed after the iterations for the (T) energy to avoid or the Vabc

their storage entirely. Alternatively, if the (T0) energy is also desired, 2-external TNO integrals i jk can be stored in the (T0) calculation and Vabc be recomputed after the iterations. These are not yet

24


Page 24 of 51

Page 25 of 51 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


implemented in our current program. As opposed to the PNO-LCCSD program, where we store the PNO overlap integrals in nodespecific shared memory, all the overlap integrals that occur in the (T) calculations are computed on the fly. To reduce the cost, in the beginning of the (T0) calculation for a triple i jk we create a united PAO domain that spans the PAO domains of all pairs that occur in eqs 3 and 4. The overlap integrals hr|ai jk i are then computed for all the PAOs r within the united domain and kept in memory. When a block of overlap integrals hbmn |ai jk i is needed, we only have to perform the second-half PAO→PNO transformation. A similar technique is used in the iterations, except that the required overlap integrals are between two different TNO domains. The computational cost of the (T0) correction can potentially be reduced by precomputing the projected doubles amplitudes used in eq 3 at the cost of extra memory requirements, given that a same block of projected of amplitudes may be used in different terms for a triple i jk. However, the savings were found to be quite small in the PAO case 48 and are expected to be even smaller when PNOs and TNOs are used.

7

Illustrative Calculations

7.1 Scalability As an illustration of the scalability with the molecular size we show in Figure 10 the elapsed time for the (T) calculations as a function of the molecular size for linear glycine polypeptides (gly)n . The near linear scaling of the preceding PNO-LCCSD-F12 calculations for these systems was already demonstrated in our previous work 74 and is not repeated here. We can observe that all major steps of the (T) calculations show near linear scaling with the molecular size, although a slight rise of the computational time for (T0) is observed for the two largest molecules. This is due to memory limits of the computer. In the current implementation we have the full PAO overlap and Fock matrices, as well as lists of integral blocks stored in local memory for each processor. Thus, for a same number of nodes, the amount of available memory for integral transformation becomes smaller with increasing molecular size. 25


PSfrag replacements


(gly)n / aVTZ

PSfrag replacements elapsed time / s

6000

4000

(T0) calculations

(T) correction total

2000

iterations triple screening

0 4

8

12 24 16 20 number of glycine units n

28

32

Figure 10: Elapsed times of various steps in the (T) calculations for glycine polypeptides (gly)n using 80 CPU cores on four nodes. The aVTZ basis set and default parameters as described in the text were employed. Three iterations for solving the (T) amplitudes were performed for each calculation. speedup relative to 60 processors

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 26 of 51

LCCSD(T)-F12 speedup

(T) correction speedup 5

5 (T0) calculations 4

4

(T) correction total

3

ideal scaling

ideal scaling

3 LCCSD(T)-F12 total

2

2 iterations 1 60

120 180 240 number of processors

1 60

300

120 180 240 number of processors

300

Figure 11: Speedup of the (T) calculations for the Auamin molecule relative to 60 cores (three nodes) using default options and the VTZ-F12 basis set. Figure 11 shows the speedup with respect to 60 cores (3 nodes) for the (T) calculations of the Auamin molecule, using our default options and the VTZ-F12 basis set. Three nodes are needed to have enough memory for computing the LCCSD-F12 energy (unless the slower disk option is used). We see that the most expensive step, namely the (T0) calculation, shows almost perfect scaling between 120 and 320 CPU cores. The distributed memory required for the calculation is roughly independent of the number of cores. Therefore, when more nodes are used, more scratch memory is available for each processor. This leads to a more efficient batching in the integral transformation and compensates the increased communication cost. When the number of nodes is 26


Page 27 of 51 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


small, the amount of scratch memory available can be too small to carry out the normal integral transformation for some triples, and then the use of our fallback algorithms leads to a drop of efficiency. As the number of nodes increases, the distributed memory used by each node will eventually become negligible and then the parallelization efficiency will drop as the communication cost increases. The parallelization efficiency of the iterations is less ideal. This is expected, given that most i jk work of the iterations is to compute the Zabc quantities (eq 11), which requires frequent fetching

of triples amplitudes from distributed memory. Also we used a static parallelization algorithm in this step to reduce the distributed memory usage, and the relative waiting time in this step also increases, from 10% for 60 cores to 30% for 320 cores. It is possible that the parallelization efficiency of the iterations could be further improved by buffering the triples amplitudes. However, this is a relative fast step with our default settings, and the overall scaling of the (T) calculation is quite satisfactory and noticeably better than of the preceding PNO-LCCSD-F12 calculations. The overall scaling with the number of processors of the PNO-LCCSD(T)-F12 method is plotted in the right panel of Figure 11. The parallel efficiency of the PNO-LCCSD-F12 method can potentially be further improved using the integral transformation strategies described in Section 6, provided that projection approximations are applied to all 2-external Coulomb integrals and 3-external integrals. Although not specifically optimized, the PNO-CCSD(T)-F12 program can run on a single computer node using MPI files instead of GAs. Parallelization efficiency with this option is demonstrated and compared to that obtained using GAs in Figure S3 in the Supporting Information.

7.2 FH Test Set For a comparison with canonical results we used the FH test set. 126 We show in Table 1 the error statistics of the PNO-CCSD(T)-F12 results against two sets of reference values. The first set is from our own canonical CCSD(T)-F12b calculations using Ansatz 3*A and fixed amplitudes, as in our local methods. A comparison to this set only shows the errors from the local approximations. With the default settings, the PNO-LCCSD(T)-F12 method gives an RMS error of 0.18 kcal mol−1 . 27



Table 1: Error Statistics (kcal mol−1 ) for the Reaction Energies of the FH Test Set Using the VTZF12 Basis Set vs CCSD(T)-F12/3*A vs Friedrich reference method MAXa RMSb MADc MAXa RMSb MADc canonical calculations CCSD(T)-F12/3*A 0.00 0.00 0.00 0.44 0.19 0.15 CCSD(T*)-F12/3*A 0.30 0.11 0.08 0.62 0.22 0.17 local calculations with default parameters LCCSD(T0)-F12 0.85 0.28 0.21 1.28 0.39 0.29 LCCSD(T1)-F12 0.49 0.19 0.15 0.93 0.30 0.22 LCCSD(T)-F12 0.42 0.18 0.13 0.85 0.28 0.21 LCCSD(T*)-F12 0.45d 0.18d 0.14d 0.68 0.26 0.20 local calculations with large domainse LCCSD(T0)-F12 0.70 0.21 0.16 1.13 0.34 0.26 LCCSD(T1)-F12 0.30 0.11 0.08 0.74 0.24 0.18 LCCSD(T)-F12 0.22 0.10 0.07 0.65 0.23 0.17 d d d LCCSD(T*)-F12 0.25 0.10 0.07 0.58 0.23 0.18 a b c d e

Maximum deviation. Root-mean-square deviation. Mean absolute deviation. Using CCSD(T*)-F12/3*A as the reference. CC = 10−8 , T CC = 0.997, T (T0) = 10−10 , Computed with T PNO en TNO (T) T TNO = 10−8 , and T trip = 0.

28


Page 28 of 51

Page 29 of 51 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


As expected, this is larger than the local RMS error excluding the (T) correction, which amounts to 0.08 kcal mol−1 . 74 With larger PNO and TNO domains, the RMS error is reduced to 0.10 kcal mol−1 . The (T0) approximation introduces an additional error of ∼ 0.1 kcal mol−1 for both domain sizes, which is already largely eliminated using the (T1) approximation. For most reactions in this set, the errors of (T0) relative to (T) become dominant only when large domains are used. A comparison with the canonical (T*) results shows that using scaled triples correction does not affect the local errors noticeably. The other set of reference values (denoted “Friedrich reference”) is taken from ref 140. The author used a composite scheme to determine the reference values, in which the CCSD(F12*)/ccpVTZ-F12 results were corrected by the HF and MP2-F12 results from cc-pVQZ-F12 calculations, and extrapolated (T) corrections from cc-pVDZ-F12 and cc-pVTZ-F12 calculations were applied. These calculations also used a different RI basis set and a different F12 approximations than our calculations. While it is difficult to estimate the errors of these values relative to the real CCSD(T)/CBS values, the difference between this set and our canonical results should provide a reasonable estimation of the intrinsic errors of the CCSD(T)-F12/cc-pVTZ-F12 calculations. From Table 1 we see that the two reference sets have an RMS deviation of 0.19 kcal mol−1 , which is very comparable to the local errors of our PNO-LCCSD(T)-F12 method with the default settings. Schmitz and Hättig 110 have reported errors of their PNO-CCSD(F12*)(T) results against the same reference and obtained an RMS error of 1.08 kJ mol−1 (0.26 kcal mol−1 ), using PNO and TNO thresholds of 10−8 . This is very similar to the results we obtained. Their calculations used PNOs and TNOs made from canonical MOs without pair or triple list approximations, and reflect purely errors from the PNO/TNO truncations. Also a different F12 treatment, the Laplace transformed (T), and different auxiliary basis sets were used.

7.3 Intermolecular Interactions In this section we consider intermolecular interaction energies. We applied the PNO-LCCSD(T)F12 method to selected dimers from the S22 benchmark set, 125 for which canonical CCSD(T)-F12 29



Table 2: Comparison of Canonical CCSD(T)-F12 and Local PNO-LCCSD(T)-F12 Results for Intermolecular Interaction Energies (kcal mol−1 )a moleculeb benzene dimer (stacked) benzene dimer (stacked) CP benzene dimer (T-shaped) benzene dimer (T-shaped) CP uracil dimer (H-bonded) uracil dimer (H-bonded) CP uracil dimer (stacked) uracil dimer (stacked) CP pyrazine dimer pyrazine dimer CP a b c d

CCSD-F12 can. loc.c loc.d −1.21 −1.11 −1.15 −1.12 −1.06 −1.09 −2.04 −1.99 −2.01 −1.96 −1.93 −1.95 −19.73 −19.63 −19.74 −19.65 −19.60 −19.65 −7.91 −7.79 −7.85 −7.81 −7.71 −7.76 −2.67 −2.57 −2.62 −2.62 −2.53 −2.57

CCSD(T)-F12 can. (T) can. (T*) loc. (T0)c loc. (T)c loc. (T*)c loc. (T0)d loc. (T)d loc. (T*)d −2.71 −2.73 −2.30 −2.33 −2.41 −2.43 −2.46 −2.54 −2.50 −2.59 −2.23 −2.26 −2.36 −2.33 −2.36 −2.46 −2.83 −2.81 −2.58 −2.60 −2.63 −2.66 −2.68 −2.71 −2.65 −2.70 −2.50 −2.52 −2.57 −2.56 −2.58 −2.63 −20.80 −20.85 −20.34 −20.42 −20.50 −20.53 −20.62 −20.70 −20.60 −20.72 −20.30 −20.35 −20.47 −20.41 −20.47 −20.60 −9.90 −9.97 −9.32 −9.36 −9.50 −9.48 −9.53 −9.67 −9.62 −9.78 −9.20 −9.24 −9.40 −9.36 −9.41 −9.57 −4.25 −4.29 −3.83 −3.86 −3.95 −3.95 −4.00 −4.09 −4.08 −4.19 −3.76 −3.79 −3.90 −3.88 −3.93 −4.04

All calculations used the F12b approximation with Ansatz 3*A and fixed amplitudes. Basis set: aug-cc-pVTZ. Rows with “CP” show the counterpoise corrected results. With default settings. CC = 10−8 , T CC = 0.997, T (T0) = 10−10 , T (T) = 10−8 , and T With tight PNO and TNO thresholds and larger triple lists: T PNO trip = 0. en TNO TNO

calculations are still viable with the aug-cc-pVTZ basis. The results are shown in Table 2. We first consider the CP corrected results. In general we see that the (T) corrections contribute significantly to the intermolecular interaction energies, and good agreement between the local and canonical results is achieved. With the default thresholds, the errors are within 0.4 kcal mol−1 for all five systems. These errors are reduced by nearly a half using larger PNO and TNO domains. Overall, the errors of the (T) correction are several times larger than the errors from the LCCSDF12 calculations. Since these are relatively small molecules and the triple list truncation error is expected to be negligible (see Figure 8), the primary source of errors must be the domain approximations. Again, we see that the scaled triples approach has hardly any effect on these errors. For the interaction energies the (T0) approximation works very well and leads to no more than 0.05 kcal mol−1 deviations from the iteratively solved (T) corrections. This is largely due to the error cancellation between the (T) contributions in the dimer and the monomers, as it can be expected that the occupied LMOs of the dimer are localized on individual monomers and are very similar to the LMOs obtained in the monomer calculations. In contrast, for the Isomer4 reaction, where the structures of the reactant and the product are very different, the (T0) approximation leads to significantly larger errors. Without the CP correction, the agreement between local and canonical results is noticeably worse. As discussed extensively in the previous papers in this series, 73,74 this is largely contributed 30


Page 30 of 51

Page 31 of 51 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


to BSSE in the canonical calculations. It can be seen in Table 2 that the CP correction is less for local methods, and there is some error cancellation between BSSE and domain approximations. In general, when computing intermolecular interaction energies using local methods, more accurate results are likely obtained without CP correction. Only when very large domains are used, the CP correction can bring local results closer to canonical ones. This observation is consistent with the findings of other studies of noncovalent interactions, 141,142 where the authors recommend to use CP corrections only for large basis sets.

7.4 Performance for Large Molecules Table 3: Reaction Energies (kcal mol−1 ) Obtained from PNO-LCCSD(T)-F12 Calculations F12b, default settings F12b, with large domainsa basis HF + CABS LCCSD LCCSD(T0) LCCSD(T) LCCSD(T*) LCCSD LCCSD(T0) LCCSD(T) LCCSD(T*) VDZ-F12 19.07 63.90 66.56 67.03 67.46 63.92 66.97 67.57 68.05 VTZ-F12 18.69 63.59 66.99 67.53 67.73 63.43 67.13 67.80 67.99 coronene dimer VDZ-F12 15.68 −12.31 −17.23 −17.35 −17.87 −12.55 −18.00 −18.17 −18.76 aVTZ 15.75 −12.76 −18.26 −18.39 −18.69 −12.87 −18.91 −19.09 −19.41 coronene dimer (CP) VDZ-F12 15.98 −11.95 −16.78 −16.90 −17.50 −12.16 −17.51 −17.67 −18.35 aVTZ 15.97 −12.31 −17.72 −17.84 −18.22 −12.67 −18.62 −18.81 −19.21 Auamin VDZ-F12 22.01 44.38 46.40 46.60 47.14 44.58 46.95 47.21 47.81 VTZ-F12 22.03 44.57 47.05 47.26 47.55 44.67 47.43 47.70 47.99 androstendione VDZ-F12 −1.30 4.51 4.61 4.83 4.92 4.60 4.81 5.07 5.17 VTZ-F12 −1.38 4.45 4.66 4.90 4.94 4.50 4.79 5.07 5.11 testosterone VDZ-F12 −6.28 −5.72 −5.46 −5.48 −5.16 −5.74 −5.53 −5.56 −5.25 VTZ-F12 −6.27 −6.17 −5.87 −5.89 −5.81 −6.22 −5.95 −5.99 −5.91 G–C/Watson–Crick (CP) aVTZ −24.58 −29.87 −30.81 −30.87 −31.03 −29.98 −31.00 −31.07 −31.24 G–C/stacked (CP) aVTZ −5.47 −15.99 −17.92 −17.98 −18.19 −16.07 −18.15 −18.23 −18.45 reaction Isomer4

a

CC = 10−8 , T CC = 0.997, T (T0) = 10−10 , T (T) = 10−8 , and T Computed with T PNO trip = 0. en TNO TNO

In this section we present PNO-LCCSD(T)-F12 results for some large molecular systems, which are described in Section 2 and earlier work. The results of calculations with default settings and very tight PNO/TNO thresholds are shown in Table 3. In triple-ζ calculations, the results of the default and large domain calculations agree to within 0.5 kcal mol−1 , with the exception of the coronene dimer, where the difference amounts to 0.7 kcal mol−1 (1.0 kcal mol−1 if CP corrected). This is largely due to the delocalized electronic structure and the very large (T) contribution to the binding energy. In agreement with our findings for other benchmark systems, the (T) correction is the major source of error. However, as shown in the convergence tests provided in Section 4, a significant portion of these errors are inherited from the LCCSD-F12 calculations, where the errors 31



in the energies are corrected by the F12 terms, but the errors in the LCCSD amplitudes remain. We note that the primary source of error of the large-domain results as compared to the canonical CCSD(T)-F12 ones (using the same F12 treatment) is the PAO domain error. This should be less than 1 kcal mol−1 and is related to the BSSE that affects the canonical calculations as well (see Section 4). As expected, the (T) results are more sensitive to the basis set compared with the LCCSD-F12 results due to the lack of the F12 treatment. In most cases, the basis set errors of (T) and LCCSDF12 add up. However, this is not always the case, since the 3*A approximation that we used for the F12 treatment is known to slightly overshoot the correlation energy. The domain errors in double-ζ calculations are in the same range of, but somewhat larger than those in triple-ζ calculations. We do not see any systematic error cancellation between the basis-set errors and the local errors. The scaled triples do improve the basis-set convergence in many cases. The LCCSD(T) results for the Auamin reaction, 47.26 and 47.70 kcal mol−1 with the default and large domains, lies well within the error bounds of the zero-point-energy corrected experimental value of 47.0 ± 2.7 kcal mol−1 . Recently, Nagy and Kállay 120 reported a reaction energy of ∼ 198 kJ mol−1 (∼ 47.3 kcal mol−1 ) from fragmentation based LNO-CCSD(T) calculations using the cc-pVTZ basis set, which is also in good agreement with our results and the experimental values. It should be noted, however, that the F12 treatment contributes ∼ 2.3 kcal mol−1 if the VTZ-F12 basis set is used and probably even more with the cc-pVTZ basis set. Thus, the good agreement is partly due to significant error compensations. From the results presented in Section 7.3, we expect that the binding energy of the coronene dimer computed without CP correction is closer to the CP-corrected canonical result. Using IEXT = 2, our large-domain aVTZ calculation gives a binding energy of −19.1 kcal mol−1 . With IEXT=3 the binding energy would be ∼ 0.3 kcal mol−1 larger, as shown in Figure 3 and Table 8 of ref 74. Adding this correction to the results of the large-domain calculation yields a binding energy is −19.4 kcal mol−1 , which lies in between the best DLPNO-CCSD(T)-F12/VDZ-F12 result of Pavošević et al. 98 (−19.14 kcal mol−1 ) and a hybrid QCISD(T) result of Sedlak et al. 124 (−24.36

32


Page 32 of 51

Page 33 of 51 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


kcal mol−1 ). Despite the good agreement of the final results, the (T) correction of Pavošević et al., −5.25 kcal mol−1 , is close to our VDZ-F12 value with large domains (−5.41 kcal mol−1 ) but not to our aVTZ result (−6.22 kcal mol−1 with IEXT = 2). Hence there is still ∼ 1 kcal mol−1 difference between the two local CCSD-F12 results. Our calculations used larger PNO domains which is important to the accuracy. Other sources for the disagreement might be different PAO domains, and the different F12 approximations with different RI basis sets. It is known that our iterative F12 treatment slightly reduces the triples energy. 25 We believe that our result is the best one presently available for this system, given the larger basis set and the larger PNO domains used, and the carefully studied convergence behavior with respect to all local approximations. The results for the androstendione reaction are presented for a comparison with the previous generation of PAO-LCCSD(T)-F12 methods 143 in Molpro. Both programs use similar theoretical frameworks and should yield the same correlation energy when local approximations are not applied. The best result reported in ref 113 is 21.6 kJ mol−1 (5.07 kcal mol−1 ), which happens to agree perfectly with our best current result. However, a closer inspection shows that the former calculation yielded a (T0) correction of 0.10 kcal mol−1 , which is only a small fraction of our current value of 0.29 kcal mol−1 . In addition, our iterative (T) correction is significantly larger (0.57 kcal mol−1 ). The PAO-LCCSD(T)-F12 program used a more aggressive triple list truncation, and even the best result reported in ref 113 is based on Rclose = 5 a0 (orbital distances between strong (T) pairs are at most 5 a0 ), which roughly corresponds to T close = 10−3 Eh in our current implementa-

tion and may not lead to converged relative energies (see Figure 7). Also the PAO domains used in the previous calculations are very small, due to the limitation of computational power at the time. These lead to too small a (T) correction, which is partly canceled by the LMP2 treatment of the close and weak pairs. The latter approximation leads to an overestimation of the LCCSD correlation energy and of the (T) correction. The OSV based (T) implementation of Schütz et al. 115 likely suffers from similar problems. The authors reported (T) corrections for the CP-corrected binding energies of Watson–Crick and stacked G–C dimers to be −0.48 and −1.61 kcal mol−1 , respectively. These are significantly less than our current best results of −1.09 and −2.16 kcal mol−1 using the

33



same basis sets. Table 4: Elapsed Time and GA Usage of PNO-LCCSD(T)-F12 Calculations Using the VTZ-F12 Basis Set and Default Options number of elapsed time / min molecule atoms AO functions nodes LCCSD-F12 triple screen (T0) iterations Auamin 92 3345 16 65.1 16.4 71.6 17.1 Auamin 92 3345 4 152.0 54.3 313.3 43.9 Auaminc 92 3345 1 319.1 146.7 766.9 65.7 (gly)32 227 7306 4 89.0 15.9 77.1 20.4 Isomer4 reactant 81 2543 2 122.0 27.9 255.7 17.4 Isomer4 product 81 2543 2 36.4 8.4 57.0 7.5 coronene dimerd 72 2544 4 169.0 60.2 485.2 94.2 a b c d

GA usage / gigaworda totalb LCCSD-F12 (T0) iterations 170.1 37.5 35.7 12.1 563.5 37.5 35.7 12.1 1298.4 not applicable 202.5 18.8 11.6 9.3 423.0 18.8 14.5 5.2 109.3 9.0 6.5 2.8 808.7 29.2 21.9 17.8

1 gigaword = 8 × 109 bytes. The maximum GA usage of the three steps determines the GA requirement of the whole calculation. Total elapsed time of the PNO-LCCSD(T)-F12 calculations, excluding that of the preceding HF + CABS calculations. Using a single faster computer node with two Xeon E5-2660 v3 @ 2.6 GHz processors (20 physical cores in total), 512 GB memory, and faster solid state hard drives. MPI files are used instead of GAs in this calculation. Using the aVTZ basis set.

In Table 4 we show the cost of the LCCSD(T)-F12 calculations for some large molecules with VTZ-F12 basis sets and default options. Typically, the (T) calculation is about three times as expensive as the preceding LCCSD-F12 calculation. This ratio decreases if the molecule has a less dense electronic structure, since then more triples can be screened to achieve the same level of accuracy. The ratio also decreases when more computer nodes are used since the (T0) calculation in our present implementation is very well parallelized. For the coronene dimer, the (T) calculation is more than four times as expensive as the LCCSD-F12 calculation due to the large number of triples (> 80000) required to achieve chemical accuracy. Pavošević et al. 114 also reported that computing the (T) energy takes most of the computational time in a DLPNO-CCSD(T)-F12 calculation of the coronene dimer using the “TightPNO” options. In all calculations presented in Table 4, most of the computational time is spent in the (T0) calculation, and the (T0) calculation is dominated by the PAO→PNO integral transformation, which typically takes 80% of the computational time of (T0). We also show the GA usage of the (T) calculations in Table 4. It can be seen that with our defaults, the (T) calculation does not require more distributed memory than in the LCCSD-F12 calculation. We also noticed that the local scratch memory requirements in (T) calculations does not exceed those of the preceding LCCSD-F12 calculations, although the performance of the (T) calculations benefits more from an increased amount of scratch memory. Therefore the minimal number of nodes required for a certain PNO-LCCSD(T)-F12 calculation is the same as in the 34


Page 34 of 51

Page 35 of 51 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-F12 calculation, as discussed in ref 74. With the default threshold (T) T TNO = 10−7 the amount of storage required for the iterations, primarily for the triples amplitudes,

is only a fraction of the requirements of other quantities used in the PNO-LCCSD-F12 calculations, such as the 3-index integrals and the PNO overlap matrices. The Laplace transformation technique is an alternative method for treating the off-diagonal couplings in the local (T) calculations. This technique avoids the storage of the triples amplitudes that is troublesome in iterative local (T) calculations when large domains of virtual orbitals are used. However, it comes with a big performance penalty as separate (T0)-like calculations need to be performed for 3–4 quadrature points for obtaining converged energies. Probably more important is that in the calculations for each quadrature point the permutation symmetry of the three-external integrals is lifted due to the introduction of Laplace modified orbitals. In our local methods this symmetry cannot be easily restored and the GA requirements would be nearly doubled compared to those in a (T0) calculation with our local DF integral transformation strategy. In this work, we (T ) have shown that using T TNO = 10−7 in the iterative (T) approach leads to very small errors when

a (T0) domain correction is applied. In this case, the cost of the iterative (T) method should be much less than that of the Laplace transformed (T) method. The Laplace transformation approach might be advantageous for very accurate calculations, when very large TNO domains are needed to describe the off-diagonal couplings in (T). However, such calculation would only be meaningful when based on a very accurate LCCSD-F12 calculation, and the latter would require huge amounts of disk and GA space on its own.

8

Summary

In the present work the PNO-LCCSD-F12 method previously developed in our group is extended to include a perturbative triples correction. The (T) correction is computed with TNOs as proposed by Riplinger et al. 93 We used a hybrid treatment with an iterative solution of the local (T) amplitude equations using small TNO domains, and correct the domain errors using large domains and the

35



semicanonical (T0) approximation. This treatment significantly reduces the computational cost over large domain iterative (T) calculations, yet the errors introduced are found to be very small, even for systems where the (T0) approximation does not work well. Local approximations used in the present method are carefully tested for a few large molecular systems, which we found to be difficult cases for local correlation methods in our previous studies. Overall, the (T) correction introduces somewhat larger errors than expected for LCCSD-F12. This is inevitable given that all major sources of errors seem to underestimate the triples correction. Furthermore, the F12 treatment, which is critical for reducing the domain and basis set errors for LCCSD, does not apply to the (T) calculation. As opposed to the conclusions from previous studies, we believe that very tight thresholds are needed to achieve chemical (∼ 1 kcal mol−1 ) accuracy. With these tight thresholds, the cost of the (T) calculation is 3–5 times the cost of the preceding LCCSD-F12 calculation. However, for all the systems that we tested, the (T) calculations did not require more disk space or memory. The cost can be significantly reduced by using smaller PAO domains and a slightly more aggressive triple list truncation threshold. This will cause an additional error of ∼ 0.5 kcal mol−1 , which may still be accurate enough for many applications. The current work is closely related to the local perturbative triples correction methods developed by Riplinger et al., 93 and by Schmitz and Hättig. 110 All programs use similar methods for generating the TNOs. The major differences between our work and the others are: First, in our method the triples amplitudes are iteratively solved with small TNO domains, which provides improved accuracy over the (T0) approximation and does not require as high computational resources as Laplace-transformed (T) treatments. Second, we introduced a two-step triple list screening process and have carefully studied the effect of triple list truncations on relative energies on large molecular systems. This showed that a much larger number of triples is needed to obtain converged results than assumed in previous works. Finally, our implementation is very well parallelized over multiple computer nodes, which provides a short time-to-solution even for large-domain calculations of large molecules. It has been shown that the major source of the remaining errors of the (T) energy are the PNO

36


Page 36 of 51

Page 37 of 51 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


and TNO domain approximations. Clearly, the TNOs are less suitable for a compact wave function representation than the PNOs in doubles treatments, and improved virtual orbitals for describing the triples would be highly desirable. Whether these can be found is still unclear. Presently we are exploring various possibilities of improving the accuracy and efficiency of the current (T) method. We are also trying to further improve the parallel efficiency of our PNOLCCSD-F12 program, based on the algorithms developed in the present work. Further high quality benchmark calculations will be presented once these improvements are completed. It is also very desirable to extend the current method for the open-shell cases, and work in this direction is in progress. Despite possible further improvements, the present PNO-LCCSD(T)-F12 method is a highly optimized parallel implementation of the local coupled-cluster method and opens up the possibility of studying large molecular systems with so far unprecedented accuracy and efficiency.

Acknowledgement This research has been supported by European Research Council (ERC) Advanced Grant 320723 (ASES).

Supporting Information Available

Benchmark reactions, triple energy plots of the Isomer4

reactant and product, and a plot demonstrating the speedup of the PNO-LCCSD(T)-F12 program on a single computer node. (pno-triples-si.pdf)

References (1) Ten-no, S. Initiation of explicitly correlated Slater-type geminal theory. Chem. Phys. Lett. 2004, 398, 56–61. (2) Ten-no, S. Explicitly correlated second order perturbation theory: Introduction of a rational generator and numerical quadratures. J. Chem. Phys. 2004, 121, 117–129.

37



(3) 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. (4) Valeev, E. F. Improving on the resolution of the identity in linear R12 ab initio theories. Chem. Phys. Lett. 2004, 395, 190–195. (5) Tew, D. P.; Klopper, W. New correlation factors for explicitly correlated electronic wave functions. J. Chem. Phys. 2005, 123, 074101. (6) 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. (7) 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. (8) 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. (9) 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. (10) Werner, H.-J.; Adler, T. B.; Manby, F. R. General orbital invariant MP2-F12 theory. J. Chem. Phys. 2007, 126, 164102. (11) 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. (12) Adler, T. B.; Knizia, G.; Werner, H.-J. A simple and efficient CCSD(T)-F12 approximation. J. Chem. Phys. 2007, 127, 221106. (13) Ten-no, S. A simple F12 geminal correction in multi-reference perturbation theory. Chem. Phys. Lett. 2007, 447, 175–179. (14) 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.

38


Page 38 of 51

Page 39 of 51 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


(15) Knizia, G.; Werner, H.-J. Explicitly correlated RMP2 for high-spin open-shell reference states. J. Chem. Phys. 2008, 128, 154103. (16) 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. (17) Shiozaki, T.; Kamiya, M.; Hirata, S.; Valeev, E. F. Equations of explicitly-correlated coupled-cluster methods. Phys. Chem. Chem. Phys. 2008, 10, 3358–3370. (18) 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. (19) Tew, D. P.; Klopper, W.; Hättig, C. A diagonal orbital-invariant explicitly correlated coupled-cluster method. Chem. Phys. Lett. 2008, 452, 326–332. (20) Valeev, E. F. Coupled-cluster methods with perturbative inclusion of explicitly correlated terms: A preliminary investigation. Phys. Chem. Chem. Phys. 2008, 10, 106–113. (21) 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. (22) Torheyden, M.; Valeev, E. F. Variational formulation of perturbative explicitly-correlated coupledcluster methods. Phys. Chem. Chem. Phys. 2008, 10, 3410–3420. (23) 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. (24) 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. (25) Knizia, G.; Adler, T. B.; Werner, H.-J. Simplified CCSD(T)-F12 methods: Theory and benchmarks. J. Chem. Phys. 2009, 130, 054104.

39



(26) Shiozaki, T.; Kamiya, M.; Hirata, S.; Valeev, E. F. Higher-order explicitly correlated coupled-cluster methods. J. Chem. Phys. 2009, 130, 054101. (27) Bokhan, D.; Bernadotte, S.; Ten-no, S. Implementation of the CCSD(T)(F12) method using numerical quadratures. Chem. Phys. Lett. 2009, 469, 214–218. (28) 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. (29) 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. (30) Werner, H.-J.; Knizia, G.; Manby, F. R. Explicitly correlated coupled cluster methods with pairspecific geminals. Mol. Phys 2011, 109, 407–417. (31) Shiozaki, T.; Knizia, G.; Werner, H.-J. Explicitly correlated multireference configuration interaction: MRCI-F12. J. Chem. Phys. 2011, 134, 034113. (32) Pulay, P. Localizability of dynamic electron correlation. Chem. Phys. Lett. 1983, 100, 151–154. (33) Saebø, S.; Pulay, P. Local Configuration-Interaction: An efficient approach for larger molecules. Chem. Phys. Lett. 1985, 113, 13–18. (34) Pulay, P.; Saebø, S. Orbital-invariant formulation and second-order gradient evaluation in MøllerPlesset perturbation theory. Theor. Chim. Acta 1986, 69, 357–368. (35) 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. (36) Saebø, S.; Pulay, P. The local correlation treatment. II. Implementation and tests. J. Chem. Phys. 1988, 88, 1884–1890. (37) Saebø, S.; Pulay, P. Local treatment of electron correlation. Annu. Rev. Phys. Chem. 1993, 44, 213– 236.

40


Page 40 of 51

Page 41 of 51 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


(38) Hampel, C.; Werner, H.-J. Local treatment of electron correlation in coupled cluster theory. J. Chem. Phys. 1996, 104, 6286–6297. (39) Maslen, P. E.; Head-Gordon, M. Non-iterative local second order Moller-Plesset theory. Chem. Phys. Lett. 1998, 283, 102–180. (40) 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. (41) 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. (42) 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. (43) Hetzer, G.; Pulay, P.; Werner, H.-J. Multipole approximation of distant pair energies in local MP2 calculations. Chem. Phys. Lett. 1998, 290, 143–149. (44) 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. (45) 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. (46) 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. (47) Schütz, M.; Werner, H.-J. Local perturbative triples correction (T) with linear cost scaling. Chem. Phys. Lett. 2000, 318, 370–378. (48) 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. (49) 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.

41



(50) 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. (51) Schütz, M. A new, fast, semi-direct implementation of linear scaling local coupled cluster theory. Phys. Chem. Chem. Phys. 2002, 4, 3941–3947. (52) 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. (53) 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. (54) 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. (55) Auer, A.; Nooijen, M. Dynamically screened local correlation method using enveloping localized orbitals. J. Chem. Phys. 2006, 125, 024104. (56) 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. (57) 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. (58) Mata, R.; Werner, H.-J. Calculation of smooth potential energy surfaces using local electron correlation methods. J. Chem. Phys. 2006, 125, 184110. (59) Mata, R.; Werner, H.-J. Local correlation methods with a natural localized molecular orbital basis. Mol. Phys. 2007, 105, 2753–2761. (60) 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. (61) Mata, R.; Werner, H.-J.; Schütz, M. Correlation regions within a localized molecular orbital approach. J. Chem. Phys. 2008, 128, 144106.

42


Page 42 of 51

Page 43 of 51 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


(62) Werner, H.-J.; Schütz, M. An efficient local coupled cluster method for accurate thermochemistry of large systems. J. Chem. Phys. 2011, 135, 144116. (63) Kats, D.; Usvyat, D.; Schütz, M. On the use of the laplace transform in local correlation methods. Phys. Chem. Chem. Phys. 2008, 10, 3430–3439. (64) Kats, D.; Manby, F. R. Sparse tensor framework for implementation of general local correlation methods. J. Chem. Phys. 2013, 138, 144101. (65) Kats, D. Speeding up local correlation methods. J. Chem. Phys. 2014, 141, 244101. (66) Kats, D. Speeding up local correlation methods: System-inherent domains. J. Chem. Phys. 2016, 145, 014103. (67) 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. (68) 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. (69) 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. (70) 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. (71) Werner, H.-J. Communication: Multipole approximations of distant pair energies in local correlation methods with pair natural orbitals. J. Chem. Phys. 2016, 145, 201101. (72) Schwilk, M.; Usvyat, D.; Werner, H.-J. Communication: Improved pair approximations in local coupled-cluster methods. J. Chem. Phys. 2015, 142, 121102. (73) 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.

43



(74) Ma, Q.; Schwilk, M.; Köppl, C.; Werner, H.-J. Scalable electron correlation methods. 4. Parallel explicitly correlated local coupled cluster with pair natural orbitals (PNO-LCCSD-F12). J. Chem. Theory Comput. 2017, 13, 4871–4896. (75) Russ, N. J.; Crawford, T. D. Local correlation in coupled cluster calculations of molecular response properties. Chem. Phys. Lett. 2004, 400, 104–111. (76) 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. (77) 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. (78) Subotnik, J. E.; Head-Gordon, M. A local correlation model that yields intrinsically smooth potentialenergy surfaces. J. Chem. Phys. 2005, 123, 064108. (79) 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. (80) Schweizer, S.; Kussmann, J.; Doser, B.; Ochsenfeld, C. Linear-scaling Cholesky decomposition. J. Comp. Chem. 2008, 29, 1004–1010. (81) 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. (82) 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. (83) 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.

44


Page 44 of 51

Page 45 of 51 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


(84) Kats, D.; Schütz, M. A multistate local coupled cluster CC2 response method based on the Laplace transform. J. Chem. Phys. 2009, 131, 124117. (85) 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. (86) Maschio, L. Local MP2 with density fitting for periodic systems: A parallel implementation. J. Chem. Theory Comput. 2011, 7, 2818–2830. (87) 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. (88) 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. (89) 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. (90) 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. (91) 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. (92) Riplinger, C.; Neese, F. An efficient and near linear scaling pair natural orbital based local coupled cluster method. J. Chem. Phys. 2013, 138, 034106. (93) 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.

45



(94) 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. (95) 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. (96) 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. (97) 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. (98) Pavošević, 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. (99) Edmiston, C.; Krauss, M. Configuration-interaction calculation of H3 and H2 . J. Chem. Phys. 1965, 42, 1119–1120. (100) Meyer, W. Ionization energies of water from PNO-CI calculations. Int. J. Quantum Chem. 1971, 5, 341–348. (101) Meyer, W. PNO-CI Studies of electron correlation effects. I. Configuration expansion by means of nonorthogonal orbitals, and application to the ground state and ionized states of methane. J. Chem. Phys. 1973, 58, 1017–1035. (102) Ahlrichs, R.; Driessler, F.; Lischka, H.; Staemmler, V.; Kutzelnigg, W. PNO-CI (pair natural orbital configuration interaction) and CEPA-PNO (coupled electron pair approximation with pair natural orbitals) calculations of molecular systems. II. The molecules BeH2 , BH, BH3 , CH4 , CH−3 , NH3 (planar and pyramidal), H2 O, OH+3 , HF and the Ne atom. J. Chem. Phys. 1975, 62, 1235–1247.

46


Page 46 of 51

Page 47 of 51 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


(103) Taylor, P. R.; Bacskay, G.; Hush, N.; Hurley, A. The coupled-pair approximation in a basis of independent-pair natural orbitals. Chem. Phys. Lett. 1976, 41, 444–449. (104) Taylor, P. R. A rapidly convergent CI expansion based on several reference configurations, using optimized correlating orbitals. J. Chem. Phys. 1981, 74, 1256–1270. (105) Staemmler, V.; Jaquet, R. CEPA calculations on open-shell molecules. I. Outline of the method. Theor. Chim. Acta. 1981, 59, 487–500. (106) Fink, R.; Staemmler, V. A multi-configuration reference CEPA method based on pair natural orbitals. Theor. Chem. Acc. 1993, 87, 129–145. (107) 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. (108) 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. (109) 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. (110) 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. (111) Werner, H.-J. Eliminating the domain error in local explicitly correlated second-order Møller-Plesset perturbation theory. J. Chem. Phys. 2008, 129, 101103. (112) 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. (113) 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. (114) Pavošević, 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.

47



(115) Schütz, M.; Yang, J.; Chan, G. K.-L.; Manby, F. R.; Werner, H.-J. The orbital-specific virtual local triples correction: OSV-L(T). J. Chem. Phys. 2013, 138, 054109. (116) In the current work a “triple” always means a selection of three local molecular orbitals (LMOs), and a “triple list” represents a subset of all possible triples. A “triple energy” is the energy contribution of a single triple, and the “triples energy” denotes the total triples energy. (117) Häser, M.; Almlöf, J. Laplace transform techniques in Møller–Plesset perturbation theory. J. Chem. Phys. 1992, 96, 489–494. (118) Almlöf, J. Elimination of energy denominators in Møller–Plesset perturbation theory by a Laplace transform approach. Chem. Phys. Lett. 1991, 181, 8449–8454. (119) Constans, P.; Ayala, P. Y.; Scuseria, G. E. Scaling reduction of the perturbative triples correction (T) to coupled cluster theory via Laplace transform formalism. J. Chem. Phys. 2000, 113, 10451. (120) Nagy, P. R.; Kállay, M. Optimization of the linear-scaling local natural orbital CCSD(T) method: Redundancy-free triples correction using Laplace transform. J. Chem. Phys. 2017, 146, 214106. (121) Schmitz, G.; Hättig, C. Accuracy of explicitly correlated local PNO-CCSD(T). J. Chem. Theory Comput. 2017, 13, 2623–2633. (122) 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. (123) 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. (124) 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. ˇ (125) 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.

48


Page 48 of 51

Page 49 of 51 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


(126) 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. (127) Knizia, G. Intrinsic atomic orbitals: An unbiased bridge between quantum theory and chemical concepts. J. Chem. Theory Comput. 2013, 9, 4834–4843. (128) 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. (129) 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. (130) 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. (131) 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. (132) 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. (133) 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. (134) 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. (135) 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. (136) 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.

49



(137) 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. (138) Liakos, D. G.; Sparta, M.; Kesharwani, M. K.; Martin, J. M. L.; Neese, F. Exploring the accuracy limits of local pair natural orbital coupled-cluster theory. J. Chem. Theory Comput. 2015, 11, 1525– 1539. (139) 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. (140) Friedrich, J. Efficient calculation of accurate reaction energies—Assessment of different models in electronic structure theory. J. Chem. Theory Comput. 2015, 11, 3596–3609. (141) Burns, L. A.; Marshall, M. S.; Sherrill, C. D. Comparing counterpoise-corrected, uncorrected, and averaged binding energies for benchmarking noncovalent interactions. J. Chem. Theory Comput. 2014, 10, 49–57. (142) Brauer, B.; Kesharwani, M. K.; Martin, J. M. L. Some Observations on Counterpoise Corrections for Explicitly Correlated Calculations on Noncovalent Interactions. J. Chem. Theory Comput. 2014, 10, 3791–3799. (143) 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.

50


Page 50 of 51

PNO-LCCSD(T)-F12 / cc-pVTZ-F12 Page Journal 51 of of 51Chemical Theory and Computation • 92 atoms • 3345 basis functions 1 Computational 2 ACS Paragon Plus Environmenttime: 3 9.4 h 80 cores on 4 nodes 2.8 h 320 cores on 16 nodes 4

Scalable Electron Correlation Methods. 5. Parallel Perturbative Triples

Recommend Documents