Scalable Electron Correlation Methods. 3. Efficient and Accurate

Subscriber access provided by NEW YORK UNIV

Article

Scalable Electron Correlation Methods. 3. Efficient and accurate parallel local coupled cluster with pair natural orbitals (PNO-LCCSD) Max Schwilk, Qianli Ma, Christoph Köppl, and Hans-Joachim Werner J. Chem. Theory Comput., Just Accepted Manuscript • DOI: 10.1021/acs.jctc.7b00554 • Publication Date (Web): 29 Jun 2017 Downloaded from http://pubs.acs.org on June 30, 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 84

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. 3. Efficient and Accurate Parallel Local Coupled Cluster with Pair Natural Orbitals (PNO-LCCSD) Max Schwilk, Qianli Ma, Christoph Köppl, and Hans-Joachim Werner∗ Institut für Theoretische Chemie, Universität Stuttgart, Pfaffenwaldring 55, D-70569 Stuttgart, Germany E-mail: [email protected]

Abstract A well parallelized local singles and doubles coupled-cluster (LCCSD) method using pair natural virtual orbitals (PNOs) is presented. The PNOs are constructed using large domains of projected atomic orbitals (PAOs) and orbital specific virtual orbitals (OSVs). We introduce a hierarchy of close, weak, and distant pairs, based on pair energies evaluated by local Møller-Plesset perturbation theory (LMP2). In contrast to most previous implementations of LCCSD methods, the close and weak pairs are not approximated by LMP2, but treated by higher-order methods. This leads to greatly improved accuracy for large systems, in particular when long-range dispersion interactions are important. Close pair amplitudes are optimized using approximate LCCSD equations, in which slowly decaying terms that mutually cancel at long range are neglected. For weak pairs the same approximations are used, but in addition the non-linear terms are neglected (coupled electron pair approximation). Distant pairs are treated by spin-component scaled (SCS)-LMP2 using multipole approximations. For efficiency reasons, various projection approximations are also introduced. The impact of these

1 ACS Paragon Plus Environment


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

approximations on absolute and relative energies depends on the PNO domain sizes. The errors are found to be negligible, provided that sufficiently large PNO domains are used for close and weak pairs. For the selection of these domains the usual natural orbital occupation number criterion is found to be insufficient, and an additional energy criterion is used. For extended one-dimensional systems the computational effort of the method scales nearly linearly with the number of correlated electrons, but the linear scaling regime is usually not reached in real-life applications for three-dimensional systems. Nevertheless, due to the parallelization that is efficient up to about 100-200 CPU cores (dependent on the molecular size), accurate calculations for three-dimensional molecules with about 100 atoms and augmented triple-ζ basis sets (e.g. cc-pVTZ-F12) can be carried out in 1-3 hours of elapsed time (depending on the molecular structure and the number of CPU cores, and excluding the time for Hartree-Fock). The convergence of the results with respect to the thresholds and options that control the domain, pair, and projection approximations is extensively tested. Benchmark examples for several difficult and large cases are presented, which demonstrate that the errors of relative energies (e.g. reaction energies, barrier heights) caused by the pair and projection approximations can be reduced to below 1 kJ mol−1 . The remaining errors are mainly caused by the finite PAO domains. The larger these are made, the more intramolecular or intermolecular basis set superposition errors are introduced, and the canonical results are approached only very slowly. This problem is eliminated in the explicitly correlated variant of our method (PNO-LCCSD-F12), which will be described in a separate paper, along with a larger set of benchmark calculations.

1 Introduction The coupled-cluster method with singles and doubles amplitudes and a perturbative treatment of triples amplitudes [CCSD(T)] is considered the gold standard of quantum chemistry. Provided that the single-determinant Hartree-Fock wave function is a good zeroth-order approximation and large basis sets are used, CCSD(T) yields relative energies with chemical accuracy (< 1 kcal mol−1 ), and often the errors are even in the range of 1 kJ mol−1 . However, using conventional methods, such accurate calculations can only be carried out for very small molecules, since the computational 2 ACS Paragon Plus Environment

Page 2 of 84

Page 3 of 84

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 effort scales as O(Nel 7 ) and O(NAO ), where Nel is the number of correlated electrons, and NAO is the

number of basis functions (for fixed Nel ). This problem is particularly severe since the correlation energy converges very slowly with basis set size. The latter problem can be alleviated by explicit correlation methods, 1–43 which include terms in the wave function that depend explicitly on the inter-electronic distances ri j and thus allow to describe the wave function cusp for ri j → 0 accurately. This leads to a drastic reduction of the basis set incompleteness errors, and with modern F12 methods 13–43 near complete basis set (CBS) quality can be achieved with triple-ζ basis sets. There are several ways to reduce the scaling of the computational effort with molecular size, which exploit the short-range character of dynamical electron correlation in insulators using localized orbitals. In traditional local correlation methods 44–77,77–107 usually two kinds of approximations are introduced: The treatment of distant electron pairs is simplified or neglected, and the excitation space for each localized electron pair is restricted to a subspace (domain) of local virtual orbitals that are spatially close to the corresponding occupied orbitals. Asymptotically, this (formally) leads to linear scaling of all computational resources with the number of correlated electrons. A disadvantage of these methods is that they are complicated and difficult to implement, and the treatable molecular size is limited by the available computational resources (mainly memory and disk space). Nevertheless, the correlation energy of rather large molecules (100-200 atoms) can be very accurately computed using local coupled cluster methods, and the errors can be well controlled. The method described in the current paper belongs to this class. Local approximations have also been developed for multi-reference wave functions. 108–116 Another approach to reduce the scaling is to split the molecule or the occupied orbital space into smaller pieces, which are treated independently, mostly using conventional methods (although the use of local correlation methods is also possible). Advantages of such fragmentation methods 117–144 is the ease of parallel implementation and strict linear scaling with the molecular size in the asymptotic limit. However, since the fragments must strongly overlap, there is a lot of redundancy and the prefactor is high. Furthermore, in order to obtain accurate results, the fragments



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

must in some cases be very large, making canonical calculations with triple-ζ (or better) basis sets for the largest fragments very expensive or even impossible. A special way of assembling the correlation energy using a many-body expansion is used in the so-called incremental methods, 145–153 which also belong to the group of fragmentation methods. In the following, we will only consider traditional local correlation methods. For a long time mostly projected atomic orbitals (PAOs) have been used to span the domains, as originally proposed by Pulay. 44 Already 30 years ago Saebø and Pulay have shown that with small domains (typically all PAOs at two atoms for a bicentric bond and at one atom for lone pairs) 98-99% of the canonical correlation energy (for a given basis set) can be recovered. 45–49 The errors caused by the domain approximation (in the following “domain errors”) decrease with increasing basis set size. Distant pairs have been treated by local MP2 (LMP2) and multipole approximations. 55,56,83 Based on these approximations, integral-direct linear-scaling PAO-LMP2 and PAO-LCCSD(T) methods have been developed in our group. 57–63 Somewhat later, a significant improvement of the efficiency was achieved by the introduction of local density fitting approximations. 65,66,74 Furthermore, the combination of local and explicit correlation (F12) treatments lead to a drastic improvement of the accuracy. 78–81,154–166 This is not only caused by the reduction of basis set incompleteness errors, but also due to a strong reduction of the domain errors. The latter improvement is caused by implicit excitations to the space outside the domains in the explicitly correlated treatment. Since the F12 terms are fast decaying, near linear scaling can also be realized for explicitly correlated LMP2-F12 and LCCSD(T)-F12. 80,156,157 A disadvantage of the PAO based local methods is that the errors are not always easy to control. In most cases, the domains are chosen based on a selection of atoms for each pair of localized molecular orbital (LMO), and all PAOs arising from AOs at these centers are included. If the electronic structure changes strongly along a reaction path, it may be difficult to achieve a balanced treatment. Furthermore, steps on the potential energy surface (PES) can occur, 167 unless the domains are frozen, 70 or smoothening techniques 168 are used. Energy based domain selection schemes have also been proposed 69,77 but not much used in practice. Unfortunately, the domain


Page 4 of 84

Page 5 of 84

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


error decreases rather slowly with increasing domain size. It has been demonstrated that the small “primary domains” must be extended by including the PAOs at 2–3 shells of neighboring atoms in order to recover 99.9% of the canonical correlation energy. 69,79 For triple-ζ basis sets, the average domain sizes may then become larger than 600 per LMO pair. This leads to excessive CPU and disk requirements in calculations for large molecules. It should be noted, however, that this problem is very much alleviated in the local F12 methods. 78–81,154–157 The explicit correlation treatment typically reduces the domain errors by one order of magnitude, and often 99.9% of the canonical correlation energy can already be recovered with the small primary domains. It is well-known since the early days of computational chemistry that pair natural orbitals (PNOs) lead to the fastest convergence of the correlation energy with domain sizes. Meyer was the first to use PNO domains in CISD and CEPA methods. 169,170 He obtained about 90% of the correlation energies of H2 O and CH4 in the early seventies - a breakthrough at that time. Later, similar PNO methods were also implemented by Ahlrichs et al., 171 Staemmler, 172 and Taylor. 173,174 PNOs are pair-specific virtual orbitals that are obtained by diagonalizing external pair density matrices, which are computed from approximate (e.g. MP2) doubles amplitudes. The eigenvalues of the pair density matrices are the natural pair occupation numbers, and domains can be selected by omitting PNOs that belong to occupation numbers which are smaller than some threshold T PNO . This allows for a fine-grained error control. The PNO domain sizes are typically one order of magnitude smaller than the PAO domains needed to achieve similar accuracy. This leads to 2-3 orders of magnitude speedup in the iterative solution of the LCCSD equations. A further advantage of PNOs is that the domains decrease with the distance of the correlated electrons, facilitating the accurate treatment of the very numerous weak and distant pairs. 83 However, two major problems are that (i) the total number of PNOs increases with the number of pairs and can become huge in large molecules, and (ii) PNOs for different pairs are non-orthogonal. For a long time it was therefore believed that PNOs are unsuitable for large molecules, due to the very difficult and expensive integral transformations. In impressive and pioneering work, Neese and co-workers 96–104,107 have shown that these problems can be overcome by modern local density fitting techniques. Since



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

2009, they developed efficient PNO-LCCSD(T) methods for closed- and open-shell molecules, and demonstrated that these can be used for real chemistry applications for systems with > 100 atoms, often achieving chemical accuracy (1 kcal mol−1 ) for relative energies (as compared to canonical calculations with the same basis). 175 However, more recently the errors were found to be much larger when dispersion interactions are important, and then very tight thresholds were needed to obtain converged results. 107,176 This is due to the overestimation of dispersion energies by LMP2, which was used to approximate the weak and distant pairs. Originally, Neese et al. created the PNOs from semi-canonical MP2 amplitudes in the full virtual orbital basis, which scales as O(Nel5 ) [or O(Nel4 ) if distant pairs are neglected]. Later on, in the so-called domain-based local PNO (DLPNO) methods, 101–104 the initial semi-canonical amplitudes were created in large domains of PAOs, which leads asymptotically to linear scaling. More recently, PNO-based local correlation methods have also been developed in our group 78–83 and by Tew and Hättig et al. 158–163 Some of these methods used domains of orbital specific virtuals 177–180 (OSVs) to obtain the initial amplitude approximation, which makes it possible to reduce the scaling for the PNO generation to cubic. 161 In our methods, a sequence of transformations PAO → OSV → PNO is implemented, which leads to linear scaling and minimizes the CPU and memory demands for generating PNOs, without loosing significant accuracy. 78–82 Our current efforts are concentrated on developing new PNO-based low-order scaling correlation methods which are well parallelized. Ideally, the target are “scalable” implementations, i.e., linear scaling of the computational effort with molecular size and inverse linear scaling of the compute time with the number of CPU cores. Near linear scaling with molecular size has been achieved previously for PNO-LMP2 79,82 and PNO-LMP2-F12, 80 and these methods are also well parallelized. In the current work this is extended to PNO-LCCSD. With “well parallelized” we mean an implementation which runs well on multiple computer nodes with communication over a fast network (e.g. Infiniband) and typically a total of 100-200 CPU cores. This is to be distinguished from “moderately” parallel codes that work only on a single node with shared memory and a shared file system, or “massively” parallel programs that are designed for many thousands of


Page 6 of 84

Page 7 of 84

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


CPU cores. Due to the large communication demands in CCSD methods, massive parallelization is very difficult to achieve without introducing large redundancy in the computations. Even for our well parallelized PNO programs true inverse linear scaling with the number of cores or nodes could not be achieved, due to communication overheads and imperfect load balancing. True linear scaling of the computational effort with system size is also difficult to achieve in PNO-LCCSD. While all dominant parts of our current PNO-LCCSD method scale (in the asymptotic limit) nearly linearly with the molecular size, we decided for the sake of accuracy and robustness to keep a few terms which scale quadratically or cubically, though with a very small pre-factor. Efficient PNO-LCCSD methods require a number of approximations beyond the use of PNO domains. First of all, approximations for close, weak, and distant pairs are necessary in order to make calculations for large molecules feasible. Second, many summations in the residual equations have to be truncated in order to reduce the computational effort and to achieve low-order scaling. Third, various projection approximations as first introduced by Neese et al. 96,97,101 can be used to reduce the number of transformed integrals that need to be computed and stored. And finally, density fitting and screening techniques are used for efficient computation of the transformed integrals in the LMO/PNO basis. The correlated electron pairs can be classified according to their distance or their correlation energy contributions. In the current work we distinguish strong, close, weak, distant, and very distant pairs. Energy thresholds are generally preferred, and the different pair classes can then be defined based on LMP2 pair energies and energy thresholds. Distant and very distant pairs are treated by LMP2 and multipole approximations. 55,56,83 The impact of this approximation for PNO-LMP2 correlation energies has recently been investigated in detail by one of us. 83 It was found to be very small provided that sufficiently high multipole orders (up to dipole-octupole and quadrupole-quadrupole interactions) are included, and couplings between the distant pairs are not neglected. In the current work these couplings are included for distant pairs, but not for very distant pairs. In our previous PAO-LCCSD(T) as well as in the DLPNO-LCCSD(T) methods of Riplinger and Neese 101,104 weak pairs were also treated by LMP2 (or even semi-canonical LMP2). As al-



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

ready pointed out in earlier studies 42,181–183 this is not a good approximation for intermolecular interaction energies, since LMP2 significantly overestimates long-range dispersion contributions. Schütz et al. 181–183 therefore proposed more accurate approximate coupled-cluster treatments of intermolecular pairs. In a previous communication, 81 we demonstrated that the overestimation of intramolecular dispersion energies by LMP2 can also have a significant effect on reaction energies of large systems. We therefore proposed an approximate PNO-LCCD treatment of close and weak pairs, which is closely related to the approximations used by Schütz et al. 182 In fact, pair approximations based the same principles have already been proposed and tested 30 years ago by Saebø and Pulay in their seminal work on local MP4(SDQ) (4th-order Møller-Plesset perturbation theory with single, double, and quadruple excitation operators). 47 Most importantly, they showed that the expensive contributions of 4-external integrals cancel with 0-external and 2-external integral contributions at long range, and that the neglect of the cancelling terms leads to large savings without much loss of accuracy. This has been confirmed in our previous work 81 for very much larger systems than could be treated by Saebø and Pulay. In the current paper, we generalize this accurate weak and close pair treatment for PNO-LCCSD. We will also discuss and test the projection approximations. In their general form they have been introduced by Neese et al., 96,97,101 but to the best of our knowledge some of them have not been tested in detail so far. These approximations project integrals from one domain to another, thereby strongly reducing the number of non-redundant integrals that must be computed and stored. They become exact for complete PNO domains. Testing the local approximations is possible in various ways. For relatively small molecules, local and canonical calculations can be directly compared. With tighter and tighter thresholds, the local method should then converge exactly to the canonical result. Pair approximations can be tested by comparing local and canonical calculations for molecular dimers as a function of the intermolecular distance. If the intermolecular pairs are treated by pair approximations, the errors caused by these approximations can be seen directly. Unfortunately, such tests are not possible for large molecules, for which canonical calculations are not feasible any more. In this case one


Page 8 of 84

Page 9 of 84

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


can only study the dependence of absolute or relative energies on the truncation thresholds. In the current work we will use all these possibilities to establish the accuracy of our method. As already pointed out, explicitly correlated terms strongly reduce the basis set incompleteness and domain errors, and in our opinion they are indispensable for accurate calculations on large systems. However, the current work focusses only on the PNO-LCCSD method itself and the related domain, pair, and projection approximations. The F12 treatment is rather independent of these and will be described in a following paper. 184 Extensive benchmarks for the overall accuracy of our PNO-LCCSD-F12 method will be given therein as well. Clearly, connected triple excitations are also indispensable for accurate calculations. Linear scaling PAO, 59,60,62,85 OSV, 180 and PNO 102 implementations of the perturbative (T) treatment have been described in the past. Alternatively, Laplace transform methods 163,185 can be used. A new parallel implementation of the PNO-LCCSD(T)-F12 method is currently developed in our group and will be described in a future publication. The organization of this paper is as follows. First we will describe our benchmark systems and the relevant computational details. Then, we will introduce our notation and discuss some general aspects of local coupled-cluster theory. Subsequently, we will discuss the domain, pair, and projection approximations that are used in our method. For each approximation, we will present representative benchmark results. Finally, we will describe some technical details about our implementation and demonstrate the scaling with the molecular size and the number of processing cores for some model systems. A summary concludes the paper.

2 Benchmark systems and computational details In order to test our PNO-LCCSD and PNO-LCCSD-F12 methods, we have computed reaction energies, isomerization energies, barrier heights, and intermolecular interaction energies. On the one hand, the benchmark systems should be representative for real-life chemistry and large enough to test all kinds of local approximations. In particular, strong dispersion interactions, which can



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

for example arise in large molecules between aromatic rings, may have a large effect and their importance increases with system size. On the other hand, the benchmark systems should not be too large, so that the many necessary calculations can be run in a reasonable time. The systems we have chosen comprise 60-174 atoms, and include the most difficult test cases that we have encountered so far. The following benchmark sets have been used [cf. Figure 1; 3-dimensional representations can be found in the supporting information (SI)]: 1. The FH benchmark set of Friedrich and Hähnchen, 151 which contains 104 small to medium size molecules and 55 reactions. This includes the “difficult cases” which have been omitted in some previous benchmarks. 2. The ISOL24 benchmark of Hünerbein et al. 186 This contains 24 real-life isomerization reactions with up to 81 atoms, which have been extensively used in the past to test the performance of DFT methods with numerous different functionals. 186–188 In many of these reactions the electronic structure changes significantly. In the current paper, we only use isomerization reaction 4, which was found to be the most difficult one, due to its strong dispersion contribution to the reaction energy. In the original paper of Hünerbein et al. 186 it was also found to be the most difficult one for density functional theory (DFT), with a dispersion correction of ≈40 kcal mol−1 . Results for the full ISOL24 set will be presented in our accompanying paper on PNO-LCCSD-F12. 184 3. The WCCR10 benchmark set of Weymuth, Couzijn, Chen, and Reiher, 189 which contains 10 ligand dissociation energies of large transition metal complexes for which gas-phase experimental values are available. In the current paper we will only present results for reaction 4 of this set, denoted as “WCCR10-4” in the following, which is by far the largest system of the whole benchmark set. Again, results for all 10 reactions will be presented in the PNO-LCCSD-F12 paper. 184 4. The dissociation of a gold(I)-aminonitrene complex (AuC41 H45 N4 P, for simplicity denoted Auamin, see Figure 1). This “Auamin” reaction is taken from ref 190. It plays an important 10 ACS Paragon Plus Environment

Page 10 of 84

Page 11 of 84

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


role in catalytic aziridination and insertion reactions. The Auamin molecule contains 92 atoms and has three phenyl and three mesityl groups. There are strong dispersion interactions between these groups and with the gold atom. The computed dissociation energy can be directly compared to an experimental gas-phase value 190 of 47.0 ± 2.7 kcal mol−1 , which has been obtained by subtracting the PW91/cc-pVTZ-PP zero-point correction 190 of −2.0 kcal/mol from the measured value. The Hartree-Fock value is computed to be 22.0 kcal mol−1 . Thus, the correlation contribution to the dissociation energy is large, and it is very sensitive to local approximations. The reaction therefore provides an excellent “difficult” benchmark system. 5. The “Androstendione reaction”, which is the last step in the synthesis of androstendione. This has also been used as a benchmark reaction in our previous work on local correlation methods. 79–81,157 6. An organocatalysis reaction, 191 which plays a role in a highly enantioselective organocatalysis. 7. The coronene dimer in a parallel displaced stacked configuration (C2C2PD) from the L7 test set of ref 192. This has also been used as benchmark system in the recent paper of Pavošević et al. on a DLPNO-CCSD-F12 method. 107 For consistency with the other papers in this series,, 79,80,184 all calculations have been carried out using the cc-pVTZ-F12 (VTZ-F12) basis sets 193 unless otherwise noted. These were specially optimized for explicitly correlated (F12) calculations. They contain diffuse functions, which are particularly important for F12 calculations (the s and p sets correspond to the aug-cc-pVQZ basis sets 194 ). Furthermore, these basis sets yield considerably better Hartree-Fock energies than the standard cc-pVTZ 195 or aug-cc-pVTZ 194 basis sets. For the gold and ruthenium atoms, the ECP60MDF and ECP28MDF effective core potentials, 196 respectively, with the associated augcc-pVTZ-PP basis sets 197 were employed. In some calculations of intermolecular interactions we used the standard aug-cc-pVTZ basis set. In these cases the non-augmented cc-pVTZ basis set was 11 ACS Paragon Plus Environment


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

Page 12 of 84

The Auamin reaction Mes

Mes N

N

N

N

Au N

Au N PPh3

N

+

PPh3

N

Mes Auamin

Mes

Mes=Mesityl

The Androstendione reaction O

O O

O

OH H

H H

H

H

HO

+ H

O

O

The Isomer4 reaction

O OH

The Organocatalysis reaction

S H t-Bu Me N Ph Ph

H

Ph

N

Ph

S

N

H t-Bu Me N Ph

H H N

O

H Ph Ph

t-Bu N

Ph

t-Bu Me N Ph

H

N

H N O

H

H Ph Ph

H H

Ph

S

N

H

N

t-Bu

Ph

N

O H Ph Ph

H H

t-Bu N H

H

The WCCR10-4 reaction PCy3 Cl Ru Cy3 P

PCy3 Cl Ru -PCy3

Cl PCy3

Cl

Cy3 P

Cy=Cyclohexane

Figure 1: The benchmark reactions


N

H

H N

H

Ph

S H

N

t-Bu Me N Ph Ph

H

H

N

H N O H Ph Ph

t-Bu N H H

Page 13 of 84

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


used for the hydrogen atoms, which reduces linear dependencies (in particular for the coronene dimer). In the following, this mixed basis is denoted aVTZ. Local density fitting (DF) approximations 65 were employed throughout this work for the evaluation of two-electron integrals, and the aug-cc-pVTZ/MP2FIT basis sets 198 were used as the DF auxiliary basis in all calculations. For computing the Fock matrix, we used the aug-cc-pVTZ/JKFIT basis sets, which were derived from the cc-pVTZ/JKFIT sets 199 by adding for each angular momentum another shell of diffuse functions. We note that the additional diffuse functions in the fitting basis sets have a relatively small effect in standard LCCSD calculations. However, they sometimes have a non-negligible effect in F12 calculations, and were used here for consistency with the F12 calculations in ref 184. Intrinsic bond orbitals 200 were used as LMOs. Inner-shell core orbitals were not correlated. The local fitting domains correspond to the IEXT=2 orbital domains (cf. section 4.1 and ref79), i.e. IDFDOM=2, RDFDOM=5 a0 . For the calculations with IEXT=3 larger DF domains (IDFDOM=3, RDFDOM=7 a0 ) were employed. Unless otherwise stated, all other parameters for local approximations are our program defaults values, as described in the following sections. All benchmark calculations were performed using the development version 201 of the Molpro quantum chemistry package. 202 The PNO-LCCSD method will be made available for public use in the near future.

3 Theory and implementation 3.1 Definitions and Notation The local closed-shell coupled-cluster singles and doubles (LCCSD) wave function is defined as ˆ

Ψ = eT |0i


(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 84

where |0i is the Hartree-Fock reference function and Tˆ the cluster operator Tˆ =

XX i

tai Eˆ ia +

a∈[i]

1 X X ij â ˆb T E E 2 i j a,b∈[i j] ab i j

(2)

† ij tai and T ab = T baji are the singles and doubles amplitudes, respectively, and Eˆ ia = ηαa † ηαi + ηβa ηβi

are the usual spin-summed one-electron excitation operators. [i] and [i j] denote orbital and pair domains of virtual orbitals, respectively. We will assume [i] = [ii], even though this restriction could also be lifted. If necessary, the type of virtual orbitals is given by a subscript: for example [i]PAO denotes an orbital-specific domain of PAOs, and [i j]PNO a pair-specific domain of PNOs. Here and throughout this paper we use the following index notation for the various orbital spaces: i, j, k, l denote occupied localized molecular orbitals (LMOs), and a, b, c, d denote virtual orbitals, belonging to the domains specified in the summations. In addition, the indices r, s refer to non-orthogonal projected atomic orbitals (PAOs), and µ, ν to AOs. Indices ai j , bi j , ci j , di j will denote virtual orbitals that are specific for a pair i j. In order to avoid excessive and redundant indexing, we will mostly omit the superscripts of the virtual orbital indices and use the convention that they are implied by the domain specifications in the summations. The orbitals in each domain are assumed to be orthonormal, but those belonging to different pairs are mutually non-orthogonal. The overlap matrices are defined as

S[i j,kl]

ab

= hai j |bkl i

(3)

We will in most cases express the PNO-LCCSD equations in matrix form. Matrices are written in bold face, and unless otherwise noted their elements refer to PNO domains given in square brackets. The definition of all matrices and vectors is given in the Appendix. Here we mention only a few, which are used in the following text. The singles and doubles amplitudes are considered


Page 15 of 84

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 vectors and matrices, respectively, i t[ii] a = tai ,

Ti[ij j,i j]

ab

a ∈ [i]

ij = T ab ,

(4)

a, b ∈ [i j]

(5)

The domain specification in square brackets is redundant here and could be omitted, but it is sometimes kept for clarity. Furthermore, we define integral matrices of 2-external integrals

j Ki[kl,mn]

= (akl i|bmn j)

(6)

ij J[kl,mn] ab = (akl bmn |i j)

(7)

ab

as well as contractions of 3-external and 4-external integrals with amplitudes h

K(T )[i j,i j] ij

i

ab

=

X

ij (ai j ci j |di j bi j )T cd

(8)

c,d∈[i j]

h i K(Ti j )(k) [mn]

a

=

X

ij (amn ci j |di j k)T cd

(9)

c,d∈[i j]

J(Ek j )[ik,i j]

ab

=

X

(aik bi j |kc j j )tcj

(10)

c∈[ j j]

K(Ek j )[ik,i j] ab =

X

(11)

(aik k|bi j c j j )tcj

c∈[ j j]

Some typical terms in the LCCSD residual equations are ij + Riabj = Kab

−

X kl

X

X

k

c,d∈[ik]PNO

αi j,kl

X

kj ik ki S ac (2T cd − T cd )Kdb kl S ac T cd S db + . . . ,

a, b ∈ [i j]

(12)

c,d∈[kl]

Here the indices a, b refer to PNOs for pair i j, while c, d refer to PNOs for pair ik or kl. In matrix



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

Page 16 of 84

form, this can be written as j Ri[ij j,i j] = Ki[ij j,i j] + S[i j,ik] 2Tik[ik,ik] − Tki[ik,ik] Kk[ik,i j] αi j,kl S[i j,kl] Tkl[kl,kl] S[kl,i j] + . . .

(13)

and then the ranges of the summations in the matrix multiplications are defined by the domains given in square brackets. Summations over repeated occupied indices (in this case k, l) are also implied. For brevity in later expressions, we define “projected” quantities, which are indicated by an overbar

¯ti[i j] = S[i j,ii] ti[ii]

(14)

¯ ij ¯i ¯ j† D [kl,mn] = t[kl] t[mn]

(15)

ij ¯ ij T [kl,mn] = S[kl,i j] T[i j,i j] S[i j,mn]

(16)

¯ ji ¯ ij ˜ ij T [kl,mn] = 2T[kl,mn] − T[kl,mn]

(17)

¯ ik = S[i j,ik] Tik , T¯ i j = Ti j , and ¯ti = ti . Note that S[kl,kl] is the unit matrix, and therefore T [ii] [ii] [i j,ik] [ik,ik] [i j,i j] [i j,i j] Thus, there is some redundancy in this notation, but we keep the quantities without overbars for clarity. Eq 13 can now be written most compactly as ˜ ik[i j,ik] Kk j − αi j,kl T¯ kl[i j,i j] + . . . Ri[ij j,i j] = Ki[ij j,i j] + T [k j,i j]

(18)

The equations are then formally analogous to the canonical case. The domains can be spanned by different equivalent orbital representations, and the above definitions are independent of the choice, provided that the orbitals within each domain are orthonormal. A PAO domain [i j]PAO comprises a subspace of non-orthogonal PAOs r, s. The orbitals in this domain can be orthogonalized and transformed to a pseudo-canonical (PC) representation


Page 17 of 84

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


[i j]PC-PAO by a transformation Wi j X

|ai =

a ∈ [i j]PC-PAO

ij |riWra ,

(19)

r∈[i j]PAO

so that S[i j,i j]PC-PAO = 1 and ij Fab =

X

ij ij Wra Frs W sb = δab ǫai j ,

a, b ∈ [i j]PC-PAO

(20)

r,s∈[i j]PAO

Similarly, OSV or PNO pair domains can be transformed to pseudo-canonical representations. The PNO-LMP2 and PNO-LCCSD calculations are carried out in the PC-PNO basis, and for brevity we will use in the following [i j] ≡ [i j]PC-PNO .

3.2 Dressed and undressed formulations of the PNO-LCCSD equations As shown in the Appendix, the PNO-LCCSD equations can be compactly described in terms of dressed orbitals |˜ii = |ii +

X

(21)

|ci tci

c∈[ii]

|¯aii i = |aii i −

X

h i k |ki ¯t[ii]

a

(22)

k

Using “dressed” integrals such as ˜ i j ]ab = (¯ai j˜i|b¯ i j ˜j) [K [i j,i j] h ij i kl ¯ mn ˜ ¯ j) K [kl,mn] ab = (a i|b h ij i = (akl b¯ mn |i ˜j) J¯ [kl,mn] ab

(23) (24) (25)

the doubles amplitude equations then take exactly the same form as for PNO-LCCD. An advantage of this formulation is the simple and compact form of the equations to be implemented. Furthermore, many summations over occupied indices are absorbed in the dressed orbitals, and therefore



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 84

truncations of these summations in large molecules can be avoided. However, a serious disadvantage is that the dressed integrals depend on the singles amplitudes and must therefore be recomputed in each LCCSD iteration. This is very expensive, particularly for the dressed 3-external and 4-external integrals, as well as for the dressed 2-external Coulomb integrals in eq 25. Another disadvantage of the dressed equations is that the dressed orbitals are less sparse than the original orbitals, and therefore linear scaling is more difficult to achieve. In the first stage of our work, we have implemented a partially dressed form of the PNOLCCSD equations, in which only the dressed 0–2-external integrals need to be recomputed in each iteration. This still leads to great simplifications, since for example eq 23 expands to 16 terms in the undressed form. Nevertheless, this approach is still rather expensive and not linearly scaling. In our final implementation all terms are fully expanded, and all required integrals are precomputed. The only exception is the dressed Fock matrix F˜ rs = hrs +

Xh

= Frs +

i ˜ − (rl|ls) ˜ 2(rs|ll)

l X X l

[2(rs|cl) − (rl|cs)] tcl

(26)

c∈[ll]

which is recomputed in each iteration. The terms involving this matrix describe orbital relaxation contributions and are slowly decaying. This will be further discussed in section 5. In order to achieve linear scaling (except for the dressed Fock terms), many summations have to be truncated, and domain as well as projection approximations have to be introduced. The principles of these approximations will be discussed in the following sections; full details are given in the Appendix, so that the method is exactly reproducible.


Page 19 of 84

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 Domain approximations 4.1 Definition of domains The construction of PAO, OSV, and PNO domains has been described in detail in ref 79, and therefore we will only briefly outline the procedure in this section. The thresholds and their default values are summarized in Table 1. For LMP2 the same thresholds as have been recommended in our previous papers 79,80 are used. The subsequent PNO-LCCSD calculation can be carried out with smaller domains, as will be explained in more detail in section 4.2. Table 1: Default values for the thresholds for domains and pair approximations used in PNO-LMP2 and PNO-LCCSD calculations Threshold T LMO IEXT REXT T OSV T PNO T en

PNO-LMP2 0.2 2 5 10−9 10−8 0.997

PNO-LCCSD 0.2 2 5 10−9 10−7 0.990

description primary PAO domains PAO domain extension (bonds) PAO domain extension (distance, in a0 ) OSV domains PNO occ. number threshold PNO energy threshold

a) In Eh ; these thresholds are applied using the converged LMP2 pair energies (large domains, without F12 correction) b) Using SCS-LMP2 distant pair energies computed with the iterative multipole approximation, p = 3, see ref 79 The LMOs and PAOs are represented in the AO basis by matrices Lµi and Pµr , respectively, with

P = 1 − LL† SAO

(27)

where [SAO ]µν = hµ|νi is the AO overlap matrix. The PAO orbital domains [i]PAO are determined by center lists, and all PAOs arising from AOs at these centers are included. In the first step, primary center lists are selected on the basis of IBO partial charges 200 (threshold T LMO ). These lists are then extended by adding all atoms which are connected by ≤IEXT bonds, or are within a distance REXT (in a0 ) to any atom in the primary domain. The extended PAO domains are orthogonalized 19 ACS Paragon Plus Environment


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

Page 20 of 84

and made pseudo-canonical. The eigenvalues are denoted ǫai . Subsequently, the integrals (ri|si), r, s ∈ [i]PC-PAO are computed and used to determine semi-canonical (SC)-PAO amplitudes ii T ab =

−(ai|bi) , a, b ∈ [i]PC-PAO ǫai + ǫbi − 2 fii

(28)

Next, OSV orbital domains [i]OSV are created by diagonalizing pair densities Dii = Tii Tii (or Tii directly; in this case the eigenvalues have to be squared to get the natural occupation number nia ). All OSVs with occupation numbers nia ≥ T OSV are kept and made pseudo-canonical. The transformations from the non-orthogonal PAO domains [i]PAO to the PC-OSV domains [i]PC-OSV are described by matrices Qi . The OSVs are used to compute semi-canonical distant pair energies EiOSV using multipole exj pansions for the integrals. The approximate pair energies are used to select distant pairs, using an energy threshold T dist (for details see section 5.1 and ref 83). For the remaining non-distant pairs OSV pair domains [i j]OSV = [i]PC-OSV ∪ [ j]PC-OSV are generated, orthogonalized, and made pseudocanonical. The transformation from the non-orthogonal domains [i j]OSV to the pseudo-canonical domains [i j]PC-OSV is denoted Wi j . Then the integrals (ai|b j), a, b ∈ [i j]PC-OSV are computed using local density fitting techniques as described in refs79,80. They are first evaluated in the PAO basis, and then transformed to the PC-OSV domains using the transformation matrices Qi and Wi j . The OSV integrals are used to compute semi-canonical OSV amplitudes Ti j , and from these OSV pair densities

Di j =

1 ˜ i j† i j ˜ i j i j† T T +T T , 1 + δi j

ij ij ij T˜ ab = 2T ab − T ba

(29)

which are diagonalized to yield the OSV → PNO transformations Vi j . The eigenvalues of Di j are the pair natural occupation numbers niaj . PNO domains for non-distant pairs are selected using two thresholds T PNO and T en . The PNOs for a pair i j are always processed in the order of decreasing occupation numbers niaj . A PNO ai j is included if either the pair occupation number niaj is greater or equal to the threshold T PNO ; or the 20 ACS Paragon Plus Environment

Page 21 of 84

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


semi-canonical PNO pair energy EiPNO is still smaller than T en × EiSC-OSV . As will be demonstrated j j further below, the latter criterion is important to ensure that the PNO domains of weak pairs are sufficiently large. Finally, the PNO domains are made pseudo-canonical, using a transformation matrix Ui j . The whole process can be summarized as follows: Qi

P

AO −→ [i]PAO −→ [i]PC-OSV

(30)

i j OSV = [i]PC-OSV ∪ [ j]PC-OSV

(31)

Ui j Vi j Wi j i j OSV −→ [i j]PC-OSV −→ [i j]PNO −→ [i j]PC-PNO

(32)

In each step of the transformation sequence PAO → OSV → PNO the domains are reduced. For the VTZ-F12 basis, 193 typical average pair domain sizes of large molecules with default settings are 600, 300, 30. The intermediate OSV step can optionally be skipped (as in the DLPNO methods of Riplinger and Neese 101,104 ), but the advantage of it is that the PNO generation, which scales cubically with the domain size, is significantly faster. Also the transformation matrices Vi j are smaller, which is important to minimize the local memory demands of each processor in parallel calculations. Timings for some typical cases can be found in ref 82. In the following, the dependence of the domain sizes on the PNO thresholds will be demonstrated for the so-called Auamin molecule, which is shown in Figure 1. Figure 2 shows the LMP2 pair energies as a function of the distances Ri j = |rr − r j |, where ri and r j are the charge centers of the LMOs i and j, respectively. Each of the of 7503 points represents one pair. As expected, at long range the pair energies decay approximately as R−6 i j . Also shown are the numbers of strong, close, weak, and distant pairs, using default thresholds. One can see that there are 4295 weak and distant pairs. The individual correlation energies of these pairs are small, but add up to a significant contribution to the total correlation energy. As will be shown later, these pairs also have a significant effect on the reaction energy. Figure 3 show the PNO domain sizes as a function of Ri j . Again, each point represents one pair. The first row shows the domain sizes for T PNO = 10−7 (left) and T PNO = 10−8 (right) without the



10−1

Pair energy −Eij / Eh

10−2 10−3

1602 strong pairs (LCCSD-F12)

10−4

1606 close pairs (approx. LCCSD)

10−5

2127 weak pairs (CEPA) 10−6 2168 distant pairs (SCS-LMP2) multipole approximation

10−7 10−8

1

2

4 8 Orbital distance Rij / a0

16

Figure 2: Pair energies of the Auamin molecule as a function of the orbital distance Ri j . The red line corresponds to R−6 decay.

200 150

TPNO = 10−7 Ten = 0

TPNO = 10−8 Ten = 0

TPNO = 10−7 Ten = 0.99

TPNO = 10−8 Ten = 0.99

TPNO = 10−7 Ten = 0.997

TPNO = 10−8 Ten = 0.997

100 50 0 PNO pair domain size

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

200 150 100 50 0 200 150 100 50 0

0

5

10 15 20 25 0 5 10 15 20 25 Orbital distance Rij

Figure 3: Auamin PNO domain sizes for various combinations of the thresholds T PNO and T en


Page 22 of 84

Page 23 of 84

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


energy criterion (T en = 0). In these calculations all pairs have been fully optimized by PNO-LMP2, without multipole approximations. The domains obtained with T PNO = 10−8 are nearly twice as large as the ones with T PNO = 10−7 and the maximum domain size increases from about 120 to 220. The domain sizes quickly decrease with increasing distance Ri j (decreasing pair energies), and many distant pairs vanish (zero domain size). Also the domains of weak pairs are rather small. This leads to an underestimation of the weak and distant pair energies. In the following rows of the figure the energy completeness criterion is switched on (see also Figure S1 in the SI for additional values of T en ). The distant pair domains now remain finite, though very small. As has been shown in ref 79, these very small domains are sufficient to recover a very large fraction of the long-range dispersion energies. Up to a threshold of about T PNO = 0.95 not much happens for strong and close pairs, but the domains of weak pairs begin slowly to grow (see Figure S2 of SI). For a value of 0.99 (second row of Figure 3) many close and weak pair domains grow significantly, while the strong and distant pair domains are still hardly affected. And for T PNO = 0.997 the weak and close pair domains become even larger than most strong pair domains (which are still not much affected). Overall, the average domain size increases almost by a factor of 2 when increasing T en from 0.9 to 0.997. These findings show that even with an occupation number threshold of 10−8 the close and weak pair energies are significantly underestimated, unless the energy criterion is used. This affects the long-range interactions, and can lead to significant errors of the reaction energy. However, the convergence of the absolute and relative energies with respect to the PNO domain thresholds can be much improved by domain corrections, as shown in the next section.

4.2 Domain corrections The PNO domain sizes significantly affect the computational cost of PNO-LCCSD calculations, in particular the memory or disk space needed to store 4-external integrals (ab|cd), 3-external integrals (ab|ck), and the very many PNO-overlap matrices S[i j,kl] . On the other hand, as will be demonstrated below, very tight thresholds are needed to converge the relative energies to within 23 ACS Paragon Plus Environment


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 24 of 84

∼ 1 kJ mol−1 . To a large extent, this is probably due to intramolecular basis set superposition effects (BSSE) on the reaction energies, which can only be recovered by using large domains. However, accurate results can also be obtained by first performing a LMP2 calculation with tight PNO domain thresholds, and then the LCCSD calculation with smaller PNO domains. The difference of the LMP2 correlation energies obtained with large and small domains is then added to the LCCSD correlation energy. For this, we determine the large and small domains simultaneously using the semi-canonical OSV (or PAO) amplitudes. The initial small PNO domains are a subset of the larger ones. However, this is no longer true once the transformation to PC-PNOs is carried out. This must be done separately for the large and small domains, and the eigenvectors of the j j Fock matrices for pair i j are the columns of unitary matrices Uilarge and Uismall , respectively. The j† j , where transformation from the large to the small pseudo-canonical domains is then Uilarge Uismall j only the vectors of Uilarge that belong to the small domains are used. After carrying out the LMP2

calculation with the large domains, the integrals, amplitudes, as well as the OSV-PNO transformation matrices are transformed from the large to the small domains. The overlap matrices are recomputed, and the LMP2 is repeated with the smaller domains. The extra cost is quite small, since only very few iterations are needed. Then the following domain corrections are evaluated: sc sc sc ∆ELMP2(large) = EOSV-LMP2 − EPNO-LMP2(large)

(33)

∆ELMP2(small) = ELMP2(large) − ELMP2(small) sc +∆ELMP2(large)

(34)

sc ∆ELMP2(large) is a semi-canonical domain correction for the large domains, which amounts to the dif-

ference of the semi-canonical OSV-LMP2 and PNO-LMP2 energies. Similar domain corrections are used in the DLPNO methods of Riplinger and Neese. 101,104 In our program they are obtained for free, since the semi-canonical energies are computed anyway in the domain selection process (cf. section 4.1). If the correction ∆ELMP2(small) is added to ELMP2(small) , the domain-corrected energy


Page 25 of 84

of the large PNO-LMP2 is exactly recovered. The correction ∆ELMP2(small) is also added to the final LCCSD energy, which is computed with the small PNO domains. 66

Isomer4 reaction

Auamin reaction 44 Reaction energy / kcal mol−1

64 Reaction energy / kcal mol−1

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


62 60 58 56 54 −6 10

PNO-LCCSD (0.00) PNO-LCCSD (0.90) PNO-LCCSD (0.99) PNO-LCCSD (0.00) + ∆E(LMP2) PNO-LCCSD (0.90) + ∆E(LMP2) PNO-LCCSD (0.99) + ∆E(LMP2)

10−7

10−8

10−9

42

40

38

36 10−6

TPNO

PNO-LCCSD (0.00) PNO-LCCSD (0.90) PNO-LCCSD (0.99) PNO-LCCSD (0.00) + ∆E(LMP2) PNO-LCCSD (0.90) + ∆E(LMP2) PNO-LCCSD (0.99) + ∆E(LMP2)

10−7

10−8

10−9

TPNO

Figure 4: The PNO-LCCSD reaction energies of the Isomer4 (left) and Auamin (right) reactions as a function of the PNO threshold T PNO with and without LMP2 domain corrections and for various thresholds T en (given in parenthesis in the caption). Default PAO domains (IEXT=2) are used. The effect of the PNO domain thresholds and domain corrections on the PNO-LCCSD reaction energies of the Auamin and Isomer4 reactions are shown in Figure 4 as a function of T PNO , using various thresholds T en . It is found that without the energy criterion (T en = 0) the convergence is rather slow and similar for both molecules. For a threshold T PNO = 10−7 only about 99.5% of the estimated correlation energy for full domains is recovered (see Figure S2 of SI), and the domain error of the reaction energy of Auamin is still larger than 10 kJ mol−1 . With increasing energy threshold the convergence is strongly improved. With T en = 0.99 the domain corrected values are very close to the estimated limit, even for the largest threshold, T PNO = 10−6 . We have therefore chosen T PNO = 10−7 and T en = 0.99 as default values for our PNO-LCCSD program. These corrections are very much in the spirit of composite correlation approaches, as used for example in G3 203 or W3 204 theories. Here the corrections affect only domain errors, and no recomputation of any integrals is required. It is also possible to carry out the PNO-LCCSD calculation using smaller PAO domains and/or smaller basis sets. Then the correction ELCCSD(small) −ELMP2(small) is added to an accurate PNO-LMP2 energy. We found that this also works well, but then two inde25 ACS Paragon Plus Environment


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

pendent calculations have to be carried out. In Figure S2 of the SI we also compare the above approach with the simpler semi-canonical SC domain corrections ∆ELMP2(small) as used in the DLPNO-CCSD program of Riplinger and Neese.

In this case the correction is computed as in eq 33, but using the small PNO domains, which SC are used in the PNO-LCCSD calculation. The correction ∆ELMP2(small) is somewhat smaller than

∆ELMP2(small) and underestimates the domain error of PNO-LCCSD, while ∆ELMP2(small) overestimates it, in particular for small PNO domains. However, for tight thresholds (T PNO ≤ 3 × 10−8 ) the full LMP2 correction yields in both cases faster convergence to the limit. The overshooting of the full LMP2 correction for very small domains is particularly strong for Isomer4. This is probably due to the strong overestimation of long-range dispersion interactions by LMP2, which have a large effect on the isomerization energy. The semi-canonical correction then appears to be better, but this is probably due to error compensation. Obviously, the domain correction and overall accuracy also depends on the PAO domains. In fact, the finite PAO domain sizes are probably causing the largest errors in PNO-LMP2 and PNOLCCSD calculations. This is demonstrated for our two test reactions in Figure 5. The convergence of the relative energies with increasing PAO domain sizes is frustratingly slow. Most likely, this is due to intramolecular BSSE effects, which are partly removed in local correlation methods. It is well known that this improves the accuracy of intermolecular interaction energies 57,205 [without counterpoise (CP) correction], but unfortunately this is not always true for reaction energies. As demonstrated in the figure, this problem is almost entirely removed if F12 terms are included, which strongly reduce the basis set incompleteness errors and consequently the BSSE. Furthermore, they also strongly reduce the domain errors. The importance of BSSE effects will be discussed more extensively and for more examples in our forthcoming paper on PNO-LCCSD-F12. 184 Finally we note that the convergence of the results with the PNO domain thresholds is also much improved when F12 terms are included. In this case it is again possible to apply domain corrections to the LCCSD-F12 energy, using LMP2-F12 instead of LMP2. But these corrections are very much smaller than the ones shown in Figure 4. This will be demonstrated in our paper on


Page 26 of 84

Page 27 of 84

82

62

Isomer4 reaction Reaction energy / kcal mol−1

Reaction energy / kcal mol−1

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


PNO-LMP2-F12

80

78 PNO-LMP2 76

74

1

2 IEXT

3

Auamin reaction PNO-LMP2-F12

60

58 PNO-LMP2 56 1

2 IEXT

3

Figure 5: The PNO-LMP2 and PNO-LMP2-F12 reaction energies of the Isomer4 (left) and Auamin (right) reactions as a function of the PAO domain parameter IEXT, using REXT = (2 × IEXT + 1) a0 . Default PNO domains are used. PNO-LCCSD-F12. 184

5 Pair Approximations The computational cost of local coupled cluster calculations depends strongly on the number of pairs that are fully treated by LCCSD. Pulay 44 has therefore proposed to exploit the fast decay of the pair energies as function of the distances Ri j (cf. Figure 2), and to treat “weak” pairs with small energies (or large Ri j ) more approximately by LMP2. This has been adopted by our group in early PAO-LCCSD methods, 42,50,61 and also by Neese and co-workers in their DLPNO-CCSD methods. 101,104 However, in the course of the current work we found that this approximation can sometimes lead to large errors, since LMP2 strongly overestimates long-range dispersion interactions. To demonstrate this, we consider the reaction energies of the Auamin and Isomer4 reactions as a function of a threshold T cutpairs . Pairs with energies below this threshold are considered as weak pairs and are either neglected, or approximated by LMP2 or SCS-LMP2 (this includes distant pairs). The remaining “strong pairs” are fully treated by PNO-LCCSD. The results are shown in Figure 6. It is seen that the many weak pairs with individually small energy contributions have a large overall effect on the reaction energies. If they are entirely neglected, the reaction energies



90 80

Isomer4 reaction SCS-LMP2 weak pair correction

70 60 50

approx. LCCSD weak pairs neglected

40 30 20 −3 10

10−4 10−5 Tcutpairs / Eh

Auamin reaction

60

LMP2 weak pair correction



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

10−6

50

Page 28 of 84

LMP2 weak pair correction SCS-LMP2 weak pair correction

40

approx. LCCSD

30 20

weak pairs neglected

10 −3 10

10−4 10−5 Tcutpairs / Eh

10−6

Figure 6: The PNO-LCCSD reaction energies of the Isomer4 (left) and Auamin (right) reactions as a function of the threshold T cutpairs . Pairs with energies below this threshold are either neglected (black), or approximated by LMP2 (violet), SCS-LMP2 (blue), or by our new close and weak pair approximations (red). In the latter case T close = T cutpairs , T weak = T close /10. All other pairs are fully treated by LCCSD. of both systems are strongly underestimated, since the weak pair correlation contributions to the reactant molecules are larger than those of the products. On the other hand, if the weak pairs are treated by LMP2, their correlation contribution is overestimated, leading to too large reaction energies. Better convergence is obtained with spin-component-scaled (SCS)-LMP2. 206 But even this approximation converges rather slowly with increasing T cutpairs and is not sufficiently accurate for pairs with Ei j > 10−6 Eh . A very similar qualitative behavior has been found for several other reactions. We have therefore introduced a much more accurate hierarchy of pair approximations, which treat strong, close, weak, distant and very distant pairs at different levels. A summary of these classes and the corresponding default thresholds is shown in Table 2. The individual approximations are described in the following.

5.1 Distant pair approximation The integrals (ai|b j) for distant pairs are approximated by dipolar multipole approximations, 55,56,83 using pseudo-canonical OSV domains a ∈ [i]PC-OSV , b ∈ [ j]PC-OSV . Exchange and charge transfer contributions are neglected (NOEX approximation 83 ). In the current work, we include all contri28 ACS Paragon Plus Environment

Page 29 of 84

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


Table 2: Definition of pair classes, thresholds, and their default values Pair class Threshold Default valuea (Eh ) Treatment strong pairs Ei j > T close full LCCSD close pairs T close ≥ Ei j > T weak T close = 10−4 approximate LCCSD weak pairs T weak ≥ Ei j > T dist T weak = 10−5 linearized LCCSD distant pairs T dist ≥ Ei j > T vdist T dist = 10−6 iterative LMP2 multipole approximationb −7 Very distant pairs T vdist ≥ Ei j T vdist = 10 non-iterative SC-LMP2 multipole approximationb a) The thresholds T close and T weak are applied using the converged LMP2 pair energies, while T dist and T vdist are applied using semi-canonical LMP2 pair energies obtained with the non-iterative multipole approximation, p = 3. b) Using SCS-LMP2 distant pair energies computed with the iterative multipole approximation, p = 3 (see ref 79).

butions up order p = 3, which means all contributions to the integrals that decay with R−(p+2) ij or slower. These are the dipole–dipole (p = 1), dipole–quadrupole (p = 2), and quadrupolequadrupole as well as dipole–octupole interactions (p = 3). For details and explicit expressions we refer to ref 83. The approximated integrals are used to calculate semi-canonical distant pair endist ergies Eidist j (non-iterative multipole approximation). Very distant pairs with energies E i j ≤ T vdist

are dropped, and their semi-canonical pair energies are added later to the final LMP2 and LCCSD correlation energies. For pairs with T dist ≥ Eidist j > T vdist two different approximations are possible. Either they are treated non-iteratively by the semi-canonical approximation as the very distant pairs (in this case there is no difference between distant and very distant pairs). Or their amplitudes and pair energies are optimized iteratively together with all other pairs, using the approximate integrals. ij For this purpose, semi-canonical amplitudes T ab , a ∈ [i]PC-OSV , b ∈ [ j]PC-OSV are used to generate

distant pair PNOs (dPNOs, cf. ref 83). Due to the latter restriction of the indices a and b, the distant pair density matrices Di j are block diagonal, and separate PNOs are obtained for i and j. A j j dist completeness criterion T PNO = 0.997 is used to determine the dPNO domains [i]idPNO and [ j]idPNO . j j dist All PNOs a ∈ [i]idPNO and b ∈ [ j]idPNO are kept which are needed to satisfy EiPNO ≥ T PNO × Eidist j j ,

where Eidist j are the semi-canonical distant pair energies in the PC-OSV domains. The dPNOs a and j b are added stepwise in the order of decreasing energy contributions. Note that the domains [i]idPNO j and [ j]idPNO are pair specific. Typically, each of them contains only 3–4 dPNOs. The integrals are

transformed to the dPNO basis and stored for later use in the iterative LMP2. Due to the very small domains, the storage requirements for these integrals are small, despite the large number of distant 29 ACS Paragon Plus Environment


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

pairs. It has been demonstrated in ref 83 that the iterative multipole approximation is much more accurate than the non-iterative one, and it is therefore used throughout this paper. The dominant part of the remaining error stems from the neglect of charge transfer and exchange terms (NOEX approximation), and not from the multipole approximation. More significant errors can arise if only the dipole-dipole approximation (p = 1) is used as in other works. 101,103 As shown in Figure 6, for distant pairs the SCS-LMP2 approximation is more accurate than LMP2. Therefore, SCS-LMP2 distant pair energies are used throughout this work (but note that the distant pair selection uses LMP2 distant pair energies). Due to the NOEX approximation, the SCSLMP2 distant pair energies are simply obtained by multiplying the LMP2 ones by s = ( fk + f⊥ )/2, where fk and f⊥ are Grimme’s scaling factors 206 for parallel and antiparallel spin contributions, respectively. Using the recommended values of fk = 6/5, f⊥ = 1/3 a scaling factor of s = 23/30 = 0.7666 is obtained. It should be noted that this is closely related to the so-called scaled opposite spin (SOS)-MP2 approximation, 207 in which case the scaling factor would be s = cSOS /2 = 0.65.

5.2 Weak, close, and strong pair approximations One can classify the contribution to the LCCSD correlation energy of each individual term in the residuals Ri j and ri according to its decay with respect to the distances Rik , R jk , Ril , R jl , Ri j , and Rkl where k and/or l denote summation indices in the residual equations (cf. Appendix). The decay of a particular term with respect to these distances can be different, and the distance which causes the fastest decay is used to classify the term and will in the following discussion be simply denoted RP . Apart from the decay properties of a term, the expected magnitude of its contribution is an important consideration. A demonstration how to elaborate the contributions of the individual terms of the residuals Ri j and ri to the correlation energy is given in the SI. The long-range decay of a term with respect to R can be either exponential (if dependent on the overlap of two LMOs), or with R−n (1 ≤ n ≤ 6) (if dependent on multipole interactions). It turns out that all terms in the residuals that decay with R−1 , R−2 , R−4 , and R−5 in the summation 30 ACS Paragon Plus Environment

Page 30 of 84

Page 31 of 84

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


indices cancel at long range and only three types of decay occur: fast exponential decay, medium fast R−6 decay, and slow R−3 decay. The cancellations are related to particle-hole symmetries: in a diagrammatic representation, the terms that cancel differ only in the direction of the arrows in the related diagrams. This has already been worked out for many LCCSD terms by Schütz et al., 181,182 and the resulting approximations for PNO-LCCD have been described in detail in ref81. Here they are extended and refined for the LCCSD case. With the exception of one high order term in the singles residual (the second term in eq 65, entering only in MP6), the slowly decaying R−3 terms all involve single excitations and can be fully included through the singles Fock-like contribution F˜ (cf. eq 26). 61,74,100,182 This matrix is computed in each iteration in the AO basis and then transformed to the LMO/PNO basis. This proceeds in close analogy to the computation of the two-electron contributions in a Fock matrix, and local density fitting approximations can be used to achieve linear scaling for the exchange contribution. 208 The Coulomb contribution scales quadratically, and the transformation of F˜ from the AO into the PAO basis even scales cubically. Earlier attempts to approximate these contributions in a linear scaling fashion showed that this is a very sensitive approximation: in some case the errors were small and it worked well, but in others the LCCSD iterations did not even converge. 61,74 We therefore decided to sacrifice linear scaling for accuracy and robustness. The additional effort in each LCCSD iteration is less than for a Hartree-Fock iteration. The general strategy for the pair approximations used in our method is summarized in Table 3. In the first column of this table the pair class of the residuals Ri j is given, while the remaining Table 3: Truncation of summations in the PNO-LCCSD equations (see text) pair classes of P included in the summations linear terms non-linear terms ij −RP −6 −RP R e decay RP decay e decay R−6 P decay strong strong,close strong,close,weak strong,close strong,close,weak close strong,close strong,close,weak strong strong,close a) weak strong,close strong,close,weak — —a) a) not included for weak pair residuals columns refer to the classes of index pairs P in the summations of these residuals. As outlined 31 ACS Paragon Plus Environment


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

above, the decay of each individual term depends on a distance RP , where P is a pair index ik, il, jk, jl, or kl, and k, l are the summation indices. The summations in the residuals are restricted by the condition that the index P belongs to the listed pair classes. For the selection of the pair classes we use the converged PNO-LMP2 pair energies, see Table 2. If the LCCSD is carried out with smaller PNO domains, then the pair energies of the PNO-LMP2 calculation with large domains (cf. eq 34) are used to determine the pair classes. The strong pair residuals are based on the full LCCSD equations, but the summations in the terms are restricted to certain pair classes as indicated in Table 3. The singles residuals ri are treated similarly to the strong pair residuals Rii , which is justified by considerations given in the SI. All details of the summation restrictions in the strong pair and singles residuals are given in the Appendix, those of the close and weak pair residuals are given in the SI. The general goal of the close and weak pair approximations is to simplify the full LCCSD equations without sacrificing much accuracy. The close pair residuals include all terms that formally contain contributions to the MP4 correlation energy with exception of the terms whose summed contributions decay as R−9 i j , as described in detail in ref81 and references therein. Two higher order groups of terms are also included into the close pair residuals, even though they only contribute in the order of MP5 or higher to the correlation energy. One of the groups contain dressed Fock matrices. Their magnitudes can be significant while the computational overhead is minimal. Also, we discovered that the contributions containing the 3-external integral contractions K(Ek j ) (cf. eq 11) and their related 1-external counterparts as defined in eq 77 give significant contributions to the close pair residuals and are therefore included. This can be rationalized by the fact that if the orbital k is located spatially in between the orbitals i and j, it acts as a “bridging” orbital that slows down the decay of the overall contribution with the distance Ri j . In the weak pair residuals we use the same approximations as for close pairs, and in addition all non-linear terms in the LCCSD equations are neglected, leading to the coupled electron pair approximation (CEPA). 169,170 The detailed summation restrictions of the close and weak pair residual can be found in the SI. If T weak ≤ T dist the weak pair approximation is not used.


Page 32 of 84

Page 33 of 84

Figure 7 demonstrates the accuracy of these approximations for the Auamin and Isomer4 reactions on a finer scale than in Figure 6. The figure shows the reaction energies as a function of the threshold T close . In order to minimize the number of independent parameters, we have chosen T weak = T close /10. In addition, the results without weak pair approximation (T weak = 0) are also shown. For the largest threshold (T close = 10−3 Eh ), the errors of the Isomer4 and Auamin reaction 43.0

Isomer4 reaction

Auamin reaction


63.6 Reaction energy / kcal mol−1

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


63.4

63.2

Tweak = Tclose /10 Tweak = 0

63.0

42.8

42.6 Tweak = Tclose /10 Tweak = 0 42.4

42.2

62.8 10−3

10−4 10−5 Tclose / Eh

10−6

42.0 10−3

10−4 10−5 Tclose / Eh

10−6

Figure 7: The reaction energies of the Isomer4 (left) and Auamin (right) reactions as a function of the pair selection thresholds T close and T weak with default domains and default projection approximations. The red curves correspond to those in Figure 6. energies caused by the close pair approximation amount to 0.8 and 0.6 kcal mol−1 , and the additional weak pair errors are 0.03 and 0.2 kcal mol−1 , respectively. The total errors quickly decrease to negligible values for smaller thresholds T close . For our default threshold T close = 10−4 Eh , they are reduced to 0.03 and 0.1 kcal mol−1 , respectively. This should be compared with the much larger errors that arise if the close and weak pairs are treated by LMP2 or SCS-LMP2 as shown in Figure 6. The improved treatment of close and weak pairs is particularly important if strong dispersion interactions are present. An extreme example is the interaction energy of the coronene dimer, which will be discussed in Section 9.



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

5.3 A case study: the benzene dimer As a further test and demonstration of the improved pair approximations, we apply in this section our method to the stacked and T-shaped forms of the benzene dimer, with geometries taken from the S22 benchmark set. 209 For these calculations we have chosen the aVTZ orbital basis set. Note that these calculations do not include F12 terms and triples corrections. The total interaction energies may therefore have large errors, but here we are only interested in the additional errors caused by the pair approximations. Figure 8 shows the computed PNO-LCCSD interaction energies of the stacked and T-shaped structures as a function of the pair selection threshold T close , using different approximations for the pairs with energies below this threshold. The reference is a calculation without pair approximations, which corresponds to a threshold T close = 1 µEh . If only the close-pair CCSD approximation is applied (red curves), the errors for the largest threshold (1 mEh ) amount only to 0.07 kcal mol−1 for both structures. If only the weak pair CEPA approximation is used (violet curves, T weak = T close ), the errors increase to about 0.29 kcal mol−1 , still being quite similar for both structures. However, if the pairs with energies below the threshold of 1 mEh are treated by LMP2 (black curves, no multipole approximations), the errors increase drastically to 3.35 and 1.48 kcal mol−1 for the stacked and T-shaped structures, respectively. This is due to the well known overshooting of the MP2 approximation for dispersion interactions. These errors are reduced to 1.63/0.62 kcal mol−1 by the SCS-MP2 variant (green curves). The blue curves show the results for the case that both the close and weak pair approximations are used, with T close = T weak /10. In this case the errors become negligibly small for T close ≤ 0.3 mEh . The open symbols (placed arbitrarily at the threshold of 1 mEh ) correspond to results obtained with all intramolecular pairs fully treated by PNO-LCCSD, while all intermolecular pairs are treated by the respective approximation. A general problem associated with the assignment of pairs to different classes is that this assignment depends on the geometry, which can lead to small steps on potential energy surfaces. This problem can be avoided by selecting the pair classes at one distance and then keep them fixed for all other distances using a restart option. To test and demonstrate this effect, we have com34 ACS Paragon Plus Environment

Page 34 of 84

0 Interaction energy / kcal mol−1

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


Benzene dimer (stacked)

−1

−2

−3 close pair CCSD approx. CEPA (Tweak = Tclose ) SCS-LMP2 (Tdist = Tclose ) LMP2 (Tdist = Tclose ) Tweak = Tclose /10

−4

−5 10−3

10−4 10−5 Pair threshold Tclose / Eh

−1.5 Interaction energy / kcal mol−1

Page 35 of 84

10−6

Benzene dimer (T-shaped)

−2.0

−2.5

−3.0 close pair CCSD approx. CEPA (Tweak = Tclose ) SCS-LMP2 (Tdist = Tclose ) LMP2 (Tdist = Tclose ) Tweak = Tclose /10

−3.5

−4.0 10−3

10−4 10−5 Pair threshold Tclose / Eh

10−6

Figure 8: Dependence of the interaction energy of the benzene dimer on the close and weak pair thresholds using different pair approximations (see text). Upper panel: stacked structure; lower panel: T-shaped structure. The open symbols (arbitrarily placed at the threshold 10−3 Eh ) show the interaction energies if only the intermolecular pairs are treated by the respective approximations.



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

puted the interaction energy of the benzene dimer for the sandwich structure (without horizontal displacement of the two molecules) as a function of the distance between the two molecular planes. The upper panel of Figure 9 shows the computed PNO-LCCSD potential energy curves with (blue) and without (red) close pair approximations. In the latter case, the close pair list has been selected at the longest distance (8 a0 ) using T close = 3 × 10−4 Eh and then kept fixed. To keep things simple, the weak and distant pair approximations have been disabled (T weak = 0). This means that all intermolecular pairs (and some others) are treated by the close pair approximation. The error caused by this approximation is shown in the lower panel of Figure 9 (blue curve). As expected, it increases with decreasing R, but amounts to only 0.11 kcal mol−1 at R = 3.5 a0 . An alternative approach is to select the close pairs at a short distance, e.g. R = 3.5 a0 . In this case some intermolecular pairs are strong pairs, and consequently the error is much smaller (red curve). In fact, this significant improvement is caused by only 9 additional pairs; in the latter calculation there are 207 strong pairs and in the former case 198. In both cases the potentials are smooth. However, if the close pairs are selected at each individual distance R, small steps occur on the potential where pairs change their class. Using the chosen threshold T close = 3 × 10−4 Eh this happens near the minimum of the potential. The steps are small (of the order 0.02 kcal mol−1 ), but could cause trouble when fitting the potential. It should be noted that with our default threshold T close = 1 × 10−4 Eh the steps are an order of magnitude smaller and hardly visible on the scale of the figure (see Figure S4 in the SI). Small discontinuities can also be caused by variations of the PAO or PNO domain sizes as a function of geometry, but these are not visible on the scale of the figure. Microscopically smooth potentials could be obtained by freezing the domains, as described earlier. 70


Page 36 of 84

1.0 Interaction energy / kcal mol−1

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


Benzene dimer (stacked)

0.5 0.0

−0.5

−1.0 all pairs strong intermolecular pairs close

−1.5 3

4

5

6

7

8

R / a0 0.05 Error of close-pair approx. / kcal mol−1

Page 37 of 84

Benzene dimer (T-shaped)

0.00

−0.05

pair selection at R = 3.5 a0 pair selection at R = 8.0 a0 pair selection at each R

−0.10

3

4

5

6

7

8

R / a0

Figure 9: Intermolecular PNO-LCCSD interaction energy of the benzene dimer sandwich structure as a function of the distance R between the benzene planes. The two molecules are parallel to each other and not horizontally displaced. Upper panel: potential without pair approximations (red) and with close pair approximation for intermolecular pairs (blue); lower panel: errors relative to the PNO-LCCSD (strong) potential caused by the close pair approximation for T close = 3 × 10−4 , T weak = 0. Red curve: close pair list selected at R = 3.5 a0 and then frozen; blue curve: close pair list determined at R = 8.0 a0 and then frozen; black curve: close pairs selected at each R.



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 38 of 84

6 Projection Approximations A major problem in PNO-LCCSD is the large number of integral blocks that have to be computed and stored. A typical term in the PNO-LCCSD residuals is (summations over k, l are implied) ¯ ik[i j,ik] Kkl[ik, jl] T ¯ lj Ri[ij j,i j] = . . . T [ jl,i j]

(35)

Since the PNOs are different for each pair, one has to compute and store a huge number of integral blocks Kkl[ik, jl] . In order to reduce the number of these blocks, Neese et al. 96,97,101 proposed projection approximations, in which an orbital in a particular domain, e.g. |aik i, is projected onto another domain, e.g. [kl]. This leads to the expansion

|aik i ≈

X

|ckl ihckl |aik i =

X c∈[kl]

c∈[kl]

|ckl i S[kl,ik] ca

(36)

It is easily shown that this also minimizes the mean square deviation haik − a˜ ik |aik − a˜ ik i. This leads to the approximation j Ki[kl,mn] ≈ S[kl,i j] Ki[ij j,i j] S[i j,mn]

(37)

The term in eq (35) can then be written as ¯ ik[i j,kl] Kkl[kl,kl] T ¯ lj Ri[ij j,i j] ≈ . . . T [kl,i j]

(38)

The projection of the amplitudes (or integrals, whatever is more convenient) is done on the fly. The big advantage is that only the Kkl[kl,kl] blocks have to be precomputed and stored. The number of these integrals is drastically smaller than the number of all integrals Kkl[ik,l j] , but more matrix multiplications and additional S-matrix blocks are needed. Overall, the reduction of data is very significant, and the CPU time is comparable (due to savings in the integral transformation). The projection approximation might look rather crude since the domain [kl] only partly spans the do38 ACS Paragon Plus Environment

Page 39 of 84

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


mains [ik] and [ jl]. However, we have implemented this as an optional approximation (denoted project_k) for all contributions of the integrals Kkl in non-linear terms of the LCCSD equations and found that the errors are surprisingly small (see below). We also tested a slightly different projection approximation for the contributions of the Coulomb integrals Jk j in the intermediates Yi jk and Zi jk : ¯ ki[i j,ki] Jk j Ri[ij j,i j] = . . . T [ki,i j] ¯ ki[i j,i j] Jk j ≈ ...T [i j,i j]

(39)

The advantage of this approximation is that in the density fitting less integrals of the (A|rs) type are required. This greatly simplifies and speeds up the evaluation of the integrals Jk j , and the memory required to hold the 3-index integrals is significantly reduced. However, since these terms are linear in the doubles amplitudes, the error is larger than that caused by applying the projection approximation only to the non-linear terms. We implemented an option project_j to control this projection, which can take the values off (no projection), weak (only terms in weak pair residuals projected), close (terms in weak and close pair residuals projected), or all (terms in all residuals projected). The quality of this approximation depends on the PNO domain sizes. The errors are particularly sensitive to the energy criterion T en , which mainly affects the close and weak pair domains. This can be rationalized by the fact that without the energy criterion the close and weak pair domains [i j] decrease with Ri j (see Figure 3), and then do not overlap sufficiently with the pair domains [ik]. However, with larger domains (e.g. T en ≥ 0.99), the errors become very small even if the approximation is used for all pairs. This effect is demonstrated in Figure 10 for the Auamin reaction. The figure also shows that the errors are much smaller than the domain correction described in Section 4.2. In order to reduce the number of 3- and 4-external integrals we introduce further projection approximations. First of all, we approximate products of singles amplitudes by eq 15 in all terms where they occur (see Appendix). We tested a few difficult molecules in the FH test set with the



43

Auamin reaction with LMP2 domain correction


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 40 of 84

42 project j=off project j=weak project j=all 41

without domain correction 40

0.90

0.92

0.94

0.96

0.98

1.00

Ten

Figure 10: Dependence of the Auamin reaction energy on the energy threshold T en and project_j. IEXT=2, T PNO = 10−7 (see text). The dashed line corresponds to the result with large domains, T PNO = 10−8 , T en = 0.997. Results are shown with and without the domain correction described in Section 4.2, eq 34. dressed form of PNO-LCCSD equations and found the errors caused by this projection less than our default energy convergence threshold of 1 µEh . The linear contribution of the singles amplitudes to the doubles residual can use the approximation

K(Ei j )[i j,i j] ≈

¯ i j )[i j,i j] K(E

ab

=

X

(ai j i|bi j ci j ) ¯t[ij j] c

(40)

c∈[i j]

where the contractions of 3-external integrals with singles amplitudes in eq 11 are approximated using the definition of the projected singles in eq 14. This should be a rather good approximation, since in the PAO basis the singles domain [ j]PAO is contained in the pair domain [i j]PAO (which is the union of [i]PAO and [ j]PAO ). It is not guaranteed that [ j]PNO is contained in [i j]PNO , but the error will become smaller with decreasing PNO threshold and zero if T PNO = 0. This approximation is controlled by the option project_s, whose possible values are the same as project_j. When this projection is enabled for a pair ik, we also approximate a term in the singles residuals ri[ii] with ˜ ik[i j,i j] )k[ii] ≈ S[ii,ik] K(T ˜ ik[ik,ik] )k[ik] ri[ii] + = K(T 40 ACS Paragon Plus Environment

(41)

Page 41 of 84

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 otherwise the projection would not reduce the number of required 3-external integrals. For the terms involving the J(Ek j ) contractions we use ¯ ik[i j,ik] J(Ek j )[ik,i j] ≈ T ¯ ik[i j,i j] J(E ¯ k j )[i j,i j] Gi[ij j,i j] + = T

(42)

where in analogy to eq 40 X ¯ kj J(E )[ik,i j] ab = (aik bi j |kci j ) ¯t[ij j] c

(43)

c∈[i j]

This has the advantage that only integrals (ai j bi j |ci j k) involving a single PNO domain are needed. The index k is restricted by the condition that it needs to be a strong pair with either i or j (see Appendix). The effect of these approximations is small. They can be enabled using option project_je=on. For consistency, when project_je=on, we also project the 1-external contribution in eq 75 by lk j ¯l† klk[ik]j ¯tl† [i j] ≈ S[ik,i j] k[i j] t[i j]

(44)

Similar approximations can be applied for the corresponding exchange contributions Gi[ij j,i j] + = T¯ ik[i j,ik] K(Ek j )[ik,i j]

(45)

In this case three different projection approximations have been implemented and tested: (i) No projection, i.e. use eq 45 (project_ke=0). The number of the required 3-external integrals (aik k|b j j ci j ) is large, especially since they are also included into the close pair residuals. Their storage can be avoided by recomputing the contractions K(Ek j )[ik,i j] in each iteration j by a partial integral transformation, similar to the integrals Kk[ik,i j] . However, this is rather

expensive in terms of CPU time.



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 42 of 84

(ii) Only project [ j j] to [i j] as defined in eq 40 (project_ke=1): ¯ ik[i j,ik] K(E¯ k j )[ik,i j] Gi[ij j,i j] + ≈ T

(46)

The effect of this approximation is very small, but the number of 3-external integrals (aik k|bi j ci j ) in these terms still quite large. (iii) In addition, project both [ j] and [ik] to [i j] (project_ke=2): ¯ k j )[i j,i j] Gi[ij j,i j] + ≈ T¯ ik[i j,i j] K(E

(47)

This leads to the smallest number of integrals and has also been used by Riplinger and Neese. 104 The integral set is the same as needed for the contributions of the projected J(Ek j ) integrals, cf. eq 42. Table 4: Statistical analysis (in kcal mol−1 ) of the errors caused by projection approximations on the reaction energies of the FH set without pair approximationsa

j off off all all all a

b c d

k on on on on on

project_ je ke on 1 on 2 on 1 on 2 on 2

s off off off off all

T en = 0.99 (default) MAXb RMSDc MADd 0.01 0.00 0.00 0.09 0.03 0.02 0.07 0.02 0.02 0.16 0.03 0.02 0.15 0.03 0.02

b

MAX 0.01 0.10 0.16 0.14 0.26

T en = 0.90 RMSDc 0.00 0.03 0.05 0.04 0.06

MADd 0.00 0.02 0.04 0.02 0.04

Default domains, except for T en , which is shown in the table, have been used. The reference values have been obtained with all projections turned off, except for the one shown in eq 76. Maximum deviation from the reference values. Root-mean-square deviation. Mean absolute deviation. In Table 4 we show the error statistics of various projection approximations of the FH test set

with the VTZ-F12 basis set and two values of T en . In these tests pair approximations were turned off and the reference values were computed with the dressed PNO-LCCSD equations, enabling only the single unavoidable projection approximation shown in eq 76. This should have a neg42 ACS Paragon Plus Environment

Page 43 of 84

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


ligible effect. We can see from the first row of Table 4 that the combined errors of project_k, project_je, and project_ke=1 are negligible for both domain sizes. The errors of other projections are somewhat larger and overall increase when smaller PNO domains are used. Still, the RMSD projection errors are well within 0.1 kcal mol−1 in all cases. Table 5: Effect of projection approximations on reaction energies of some large systems (in kcal mol−1 ) using default domains Isomer4a

Aureacta

s off weak

62.84 62.84

42.34 42.33

Coronene dimerb 14.09 14.06

1 1 1

weak weak weak

62.84 62.82 62.82

42.32 42.29 42.25

14.04 14.01 13.95

on on on

0 1 2

weak weak weak

62.81 62.82 62.91

42.24 42.25 42.27

14.00 13.95 13.94

on on

2 2

close all

62.89 62.84

42.25 42.20

13.92 13.93

project_ k je ke on off 0 on on 1

weak close all

on on on

on on on

all all all

on on on

all all

on on

j off off

a b

Basis VTZ-F12 Basis aVTZ, without counterpoise correction

In Table 5 the computed reaction energies for some large molecular systems are shown, using our default domain sizes and various combinations of projections. We found that the approximations are nearly additive, and therefore only a selection of the many possible combinations is presented, in which the various projections are successively turned on. Furthermore, they are rather independent of the pair approximations, and we therefore only show the results for our default values T close = 10−4 Eh , T weak = 10−5 Eh . The effect of all projection approximations on the relative energies is found to be very small, even using the most aggressive projections project_j=all and project_ke=2. The latter approximation changes the decay behavior of the integrals, since |cik i contains contributions close to k, but this may not be the case for |ci j i. We had therefore anticipated larger errors for this approximation. In fact, the errors of the absolute correlation energies are of the order of 2 mEh , which is 43 ACS Paragon Plus Environment


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 44 of 84

somewhat larger than for other projection approximations. Note that in some cases there are some small error compensations between project_s and project_ke. Table 6: Default values for the various options affecting integral projections Option default value project_j all project_k on project_je on project_ke 2 project_s weak

Description j kj Jk[ik,i j] ≈ J[i j,i j] in all pair residuals i j j Ki[kl,mn] ≈ Ki[ij j,i j] in all nonlinear terms ¯ k j )[i j,i j] J(Ek j )[ik,i j] ≈ J(E ¯ k j )[i j,i j] K(Ek j )[ik,i j] ≈ K(E ¯ i j )[i j,i j] and K(Ei j )[i j,i j] ≈ K(E ˜ ik )k ˜ ik )k ≈ S[ii,ik] K(T K(T [ik,ik] [ii] [ik,ik] [ik]

Our default choice of projections are summarized in Table 6. Given that the projection errors are usually small with our choices of T PNO and T en , we decided to use by default the rather aggressive projections project_j=all and project_ke=2 to reduce the computational cost. The performance impact of project_s is rather small and we use weak as the default. We do note, however, the error of project_j increases rather noticeably when small PNO domains or doubleζ basis sets are used. This is illustrated in Figure S3 in the SI. However, the projection errors in these cases are still likely to be less significant than the domain and basis set errors. Table 7: Comparison of reaction energies (in kcal mol−1 ) with default parameters to those obtained with very tight thresholds and less projection approximations Large domains, reduced Number of projectionsb , T close / Eh a a Reaction atoms basis functions 10−4 10−5 Isomer4 reaction 81 2543 62.66 62.59 Coronene dimerd 72 2544 14.03 14.05 Auamin reaction 92 3345 42.28 42.13 Androstendione reaction 61 2113 4.58 4.52 Organocatalysis ∆Er 106 3772 −20.66 −20.62 Organocatalysis ∆ET S 1 106 3772 13.84 13.87 Organocatalysis ∆ET S 2 106 3772 10.73 10.77 WCCR10-4 reactione 174 5168 46.42 46.19 a) For the largest molecule of the reaction. b) IEXT=2, T PNO = 10−8 ,T en = 0.997, project_ke=1, project_j=close. c) IEXT=2, T PNO = 10−7 ,T en = 0.990, default projection approximations see Table 6. d) Basis set: aVTZ, without CP correction. e) Basis set: VTZ-F12, cc-pVTZ for hydrogen atoms.

Default domains and projectionsc , T close / Eh 10−4 10−5 62.91 62.87 13.94 14.00 42.26 42.18 4.56 4.51 −20.70 −20.67 13.88 13.90 10.81 10.83 46.26 46.13

In order to demonstrate the overall effect of the domain, pair, and projection approximations, 44 ACS Paragon Plus Environment

Page 45 of 84

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


we show in Table 7 a comparison of the reaction energies of some large systems obtained with the default values (cf. Tables 1, 2, 6) and calculations with large PNO domains and without the most critical projection approximation project_ke=2. Furthermore, results with T close = 10−4 Eh (default) and T close = 10−5 Eh are compared; the latter should be converged with respect to the pair approximations. A similar table showing the errors relative to these reference values can be found in the supporting information We believe that for the calculations with large domains and T close = 10−5 Eh the remaining errors caused by the PNO, pair, and projection approximations are negligible (< 0.1 kcal mol−1 ). The deviations of the results obtained with default options from these reference values are found to be always smaller than 0.2 kcal mol−1 , except for the isomer4 reaction, where it amounts to 0.31 kcal mol−1 . This demonstrates the outstanding accuracy of our method. The overall accuracy of the results (relative to canonical CCSD for a given basis set) is probably most affected by the finite PAO domain sizes. This problem is largely removed in the explicitly correlated PNO-LCCSD-F12 method, which will be presented in a forthcoming publication. 184

7 Parallelization Our PNO-LCCSD method is parallelized with the help of MPI and the Global Arrays (GA) toolkit. 210 The majority of the time of a PNO-LCCSD calculation is spent in two steps, the evaluation and the transformations of the 2-electron integrals, and the iterative solution of the LCCSD equations. We used different parallelization strategies for these steps. The 2-electron integrals in the PNO-LCCSD method are computed using an efficient local density fitting method similar to that used in our PNO-LMP2 method. 79 In this step, we use primarily a dynamic parallelization method, where the tasks are allocated to the processors on a first come first serve basis. We first evaluate all required 3-index PAO/LMO integral sets and store them in GAs. This step is distributed over blocks of fitting functions. We then parallelize over pairs of LMOs and assemble and transform the integrals to PNOs. Our integral transformation routines are designed



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 explicit correlated terms in mind, and we will describe the integral transformations in greater detail in a forthcoming paper describing our PNO-LCCSD-F12 method. 184 In the PNO-LCCSD iterations we use a static parallelization model, in which pairs of LMOs are distributed to processors before the iterations. This treatment allows us to store the amplitudes and some integrals locally in memory or on disk, reducing the communication demands. However, it is highly non-trivial to obtain a good load balancing with this model, in particular since the sharing of data among the pairs distributed to a processor should be maximized as well. Presently we perform the pair distribution using the metis graph partitioning program. 211 We first construct a graph with vertices representing LMO pairs. The vertices are given weights equal to the estimated cost of computing the doubles residuals Ri j of the pair, and are connected with edges weighted by the number of atomic centers that are common to the primary PAO domains of the pairs (i.e. vertices). The graph is then partitioned into parts with nearly equal vertex weights while minimizing the weights of cut edges. The PNO-LCCSD method requires the temporary storage and communication of large amounts of data. While most types of quantities are stored in distributed memory by default as in our PNO-LMP2 and PNO-LMP2-F12 methods, we use more flexible algorithms to handle difficult cases. For example, we store the 3- and 4-external integrals in processor-specific scratch files, and communicate and sort these integrals once before the iterations. This treatment avoids excessive memory demands for keeping these integrals in distributed memory. As another example, the PNO overlap matrices are accessed very frequently in the iterations, and we choose to store them in node-specific shared memory, which provides a good compromise between data redundancy and access latency. Such flexible options of handling large amounts of data are critical to the performance and portability of our program.


Page 46 of 84

Page 47 of 84

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


Table 8: Elapsed time (in second) of various steps and resource usages of PNO-LCCSD calculations with the VTZ-F12 basis set and default options molecule Isomer4 reactant Isomer4 product Coronene dimerb Coronene dimerb AuAmin AuAmin AuAmin (gly)32 b WCCR10-4 reactante a b c d e f

number of AO functions nodesa 2543 2 2543 2 2544 2 2544 4 3345 3 3345 4 3345 6 7306 2 5168 10f

LMP2 total 160 72 342 202 289 238 192 520 451

integrals 4683 1110 11548 6667 8136 6179 4379 3160 8746

LCCSD iterations 682 293 2184 1260 1589 1269 1056 2676 2607

total 5365 1403 13732 7927 9725 7448 5435 5836 11353

resource usage GAc diskd 17.9 30.3 8.6 17.2 25.0 133.8 29.0 73.7 36.1 44.2 36.1 36.1 36.1 24.6 17.7 47.6 93.6 16.2

20 CPU cores per node are used unless otherwise stated. aVTZ basis set. global, in gigawords (1 GW = 8 GB). maximum disk usage on each node, in GB. Basis set: VTZ-F12, cc-pVTZ for hydrogen atoms. 18 CPU cores per node.

8 General Performance and Scalability 8.1 Computational Cost In Table 8 we present timing examples for some molecules using different numbers of nodes. Each of our computer nodes has 20 CPU cores (Xeon E5-2680 v2 @ 2.80 GHz), 256 GB of memory, and 2 TB of solid state disk space. The nodes are connected by a QDR-Infiniband network. We show in the table the elapsed times of both the integral transformations and the PNO-LCCSD iterations. The cost of preceding PNO-LMP2 calculations, which has been demonstrated previously, 79,82 is also shown separately and not included in the PNO-LCCSD timings. For all the calculations shown, the majority of the time is spent in the evaluation and transformation of the integrals. This is a natural consequence of the large total number of PNOs as compared to the number of AOs, and is in line with the findings of Riplinger et al., 104 as well as with our previous works. 79,80 We will provide more detailed timings for the integral evaluation and transformation in a subsequent paper. 184 Within the iterations, a substantial amount of CPU time is spent on the evaluation of the quadratic doubles amplitudes contractions. This is on the one hand due to additional contractions with overlap matrices in the case of the option project_k 47 ACS Paragon Plus Environment


(switched on by default, cf. section 6) and on the other hand due to the fact that many of these quadratic LCCD terms exhibit a slow r−6 decay in at least one of the summation index pairs ik or jl, and that therefore the summations cannot be restricted aggressively without loosing significant accuracy. The computation times for both the integrals as well as for the iterations depend strongly on the molecular structure. This is very clearly seen for the isomer4 reaction, where the reactant is very compact and the product very extended. Thus, the calculation of the product is almost 4 times faster. This is partly due to the much larger number of distant and very distant pairs in the product (2197 vs. 1186 in the reactant), and partly to smaller domains. Even more expensive was the calculation for the coronene dimer. This is not only due to the compact structure, but also to the relatively large domains of the conjugated π-orbitals.

8.2 Scaling with the molecular size 4000

Elapsed time / s

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

Page 48 of 84

1500

1500 other integral steps

3000 PNO-LCCSD (total)

1000

1000 dressed Fock matrices

2000

integrals 500

500

1000 iterations 0

4

8

12

16

20

24

28

32

0

other iteration steps

3-index PAO integrals 4

8

12

16

20

24

28

32

0

4

8

12

16

20

24

28

Number of glycine units n

Figure 11: Elapsed time of various steps in PNO-LCCSD calculations (excluding the preceeding PNO-LMP2 calculations) of glycine polypeptides (gly)n . The calculations were computed with 80 cores on 4 nodes using the aVTZ basis set and default parameters for local approximations as described in the text. The dashed line in the middle panel shows the elapsed times in evaluating the 3-index PAO integrals using 60 cores on 3 nodes. In this section we demonstrate the scaling of the PNO-LCCSD method with respect to the number of correlated electrons for linear glycine polypeptides (gly)n , using our default options for the local approximations. These molecules are favorable test cases for local correlation methods and can reveal the asymptotic scaling behavior of the method. Similar scaling behavior should in 48 ACS Paragon Plus Environment

32

Page 49 of 84

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


principle hold for three-dimensional systems, but these usually have significantly larger numbers of electron pairs in each class, and the asymptotic scaling behavior shows up very late, in particular when good (at least augmented triple-ζ) basis sets and reasonable thresholds are used. Investigating the scaling of such systems with the PNO-LCCSD method is beyond the capability of our computer hardware. It is also quite irrelevant, since for most real-life applications it is more important that the method is fast and accurate for systems with up to 150–200 atoms. In Figure 11 the elapsed time of various steps of PNO-LCCSD calculations with respect to the number of glycine units are shown. These calculations were performed using our default options, which are found to give converged results even for difficult cases (cf. Table 7). We observe from the left panel of Figure 11 that the overall scaling of the PNO-LCCSD method is near-linear, with slightly nonlinear scaling behaviors both in the integral transformations and the iterations. We were able to identify the major contributors to the nonlinear scaling and they are illustrated in the middle and the right panels of Figure 11. For the integral transformations, the evaluation and transformation of 3-index PAO integrals scales slightly worse than the other steps. This is because the average number of AOs contributing to each LMO or PAO grows with the molecule size, even for this 1-dimensional system. 79 This growing number of AOs also leads to an increased memory usage, forcing us to use a slower fallback algorithm that uses less memory for some of the integral blocks. We believe that this algorithm switching is the primary reason for the significantly increased time of computing 3-index integrals for (gly)32 , given that such a jump shows up earlier when less scratch memory is available (dashed line in the middle panel of Figure 11). The jump disappears when smaller PAO domains (IEXT=1) are used. In each iteration, the dressed Fock matrix needs to be computed and transformed. The calculation of the exchange part of the Fock matrix scales linearly and employs local density fitting techniques similar to those described in ref 208. The Coulomb part is computed using standard density fitting and scales quadratically. Consequently, the elapsed time for the dressed Fock matrix scales between linearly and quadratically. This is illustrated in the right panel of Figure 11. We note, however, this problem is of less importance for real-life molecules with comparable number



of basis functions, where increased numbers of LMO pairs in each class make other steps more expensive. Typically, the dressed Fock matrices then take only ∼ 1/3 of the time spent in the iterations, and the latter is then ∼ 1/5 of the total cost of a PNO-LCCSD calculation.

8.3 Demonstration of Scaling with Number of Processors coronene dimer / aVTZ speedup 4 Speedup relative to 2 nodes

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

integrals PNO-LCCSD total 3 iterations 2

1

with disk files (on a faster computer) with disk files

0

2

4

6 8 Number of nodes

10

12

Figure 12: Speedups with the number of nodes (relative to 2 nodes) for PNO-LCCSD calculations of the coronene dimer using the default options and the aVTZ basis set. Each node uses 20 cores and 256 GB of memory. Two additional points are shown for single node calculations using disk files instead of distributed memory. We show in Figure 12 the speedup with respect to 2 nodes for PNO-LCCSD calculations of the coronene dimer using our default options and the aVTZ basis set. Overall, our aim of developing a “well-parallelized” program is achieved. Relative to the timings obtained with 40 cores, we observed a speedup of 2.5 when using 120 cores (83% efficiency) and 4.0 when using 240 cores (67% efficiency). The parallel efficiency is found to be comparable to those in our previous tests with the PNO-LMP2 program for an alanine helix (ala)32 . 79 An example for the speedups obtained on a single computer node is given in Figure S5 of the SI, where we see a speedup of 8.9 using 20 cores (relative to using 2 cores, 89% efficiency). The integral transformations require a large amount of communication, which contributes to its non-ideal parallel efficiency. The communication bottleneck in our present program is that 50 ACS Paragon Plus Environment

Page 50 of 84

Page 51 of 84

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


of the 3-index integrals with 2 external (PAO) labels. These integrals are expensive to compute and redundant computation is undesired for our applications. Also buffering of such integrals is hardly possible on our machines due to memory limits. Presently, we designed our program to load these integrals only once per pair from GAs and assemble Jk j integral blocks as well as the 3and 4-external integral blocks using them. The most time-consuming steps of the integral transformations are parallelized over pairs of LMOs. In these calculations we have 1965 strong pairs and a dynamic distribution of these pairs over a few hundred processors is expected to yield good efficiency. However, the computation of 3-index integrals is parallelized over blocks of fitting functions, and with a block size of 32 we have only 192 blocks for this molecule. The lack of tasks to be distributed contributes to the drop of parallel efficiency when a large number of nodes are used. The uneven load balancing due to static parallelization over pairs in the iterations is another major contribution to the reduced parallel efficiency when a large number of nodes are used. However, the static pair distribution significantly reduces the amount of communication, in particularly for loading the PNO overlap matrices, which are stored redundantly in a node-specific way. In some cases we observed a multifold increase of the time spent in the iterations when these matrices were non-redundantly stored in a global array and loaded on the fly. We therefore believe that the static parallelization is a good compromise in this step. The static parallelization also allows us to store many integrals on local disks, which reduces the memory requirements significantly. Due to the imperfect parallel efficiency, it is most cost effective to use a small number of nodes. The minimum number of nodes on which a calculation can run is determined by the available memory, given that the disk usage is usually moderate. In most of our calculations 4 GB per core were reserved as local scratch memory, and the remaining memory is used for GAs. This means that for example at least 3 nodes on our cluster (20 cores and 256 GB of memory on each) were needed for an Auamin calculation, which overall required 289 GB of GA space (cf. Table 8). In our current implementation the memory cannot be dynamically shared between GAs and local scratch memory, and memory pools for both need to be allocated prior to a calculation.



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

Alternatively, our program can also run on a single node using exclusively shared disk files to handle the large amounts of data. This is a common strategy used by “moderately” parallelized quantum chemistry programs. With this option, most calculations with default thresholds shown in this work can be performed on one computer node using 4–6 GB of memory per CPU core. It should be noted, however, that our program is not yet optimized for this environment and that the actual performance may vary significantly depending on the available memory and disk configurations. This is demonstrated in Figure 12, where a calculation on one node of our cluster is much slower than on a newer machine, which has Xeon E5-2660 v3 @ 2.60 GHz processors, 512 GB of memory, and faster disks. This processor is 20–30% faster than those in our cluster, but the main speedup is probably due to the larger memory, which allows more efficient I/O buffering.

9 Summary We have presented an efficient parallel PNO-LCCSD method, which asymptotically scales nearly linear with the number of correlated electrons. All local approximations, namely the domain, pair, and projection approximations have been described in detail and extensively tested. The general ansatz of our method is similar to the one used in the DLPNO-CCSD method of Riplinger and Neese. 101,102,104,107 However, there are major improvements of the accuracy and efficiency. The main differences are: (i) A much more accurate approximate coupled-cluster treatment of close and weak pairs has been introduced and tested in detail. This eliminates the problems caused by the LMP2 weak pair approximations, which are used in the DLPNO-CCSD(T) method as well as in our previous PAO-LCCSD(T) and PAO-LCCSD(T)-F12 methods. 74,157 Usually, LMP2 strongly overestimates the weak pair energy contributions, which can lead to errors of several kcal mol−1 in reaction energies or intermolecular interaction energies of large systems. With our new method, these errors are reduced to a few tenths of a kcal mol−1 or less. (ii) Larger PNO domains are used for close, weak, and distant pairs. This mainly results from an 52 ACS Paragon Plus Environment

Page 52 of 84

Page 53 of 84

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


additional energy-based completeness criterion for the domain selection. The larger domains increase the computational cost, but further improve the accuracy of long-range interactions. We have demonstrated that they also reduce the impact of projection approximations. (iii) A more accurate distant pair LMP2 approximation is employed, which includes higher orders in the multipole approximation of the integrals, 83 and fully optimizes the distant pair amplitudes (rather than using a semi-canonical approximation). This roughly halves the errors caused by the distant pair approximation (relative to a full LMP2 treatment of the distant pairs). (iv) Our program is designed to work in parallel on multiple computer nodes. If multiple nodes are used, as many data as possible is kept in distributed high-speed memory. The use of a shared file system by different nodes is avoided. In order to reduce the memory requirements, the 3-external and 4-external integrals which are local to each node are (optionally) written to local disks. This means that the I/O performance scales with the number of nodes. Optionally, also all data can be written to disk, and this makes it possible to carry out very large calculations even on a single computer node. (v) Some projection approximations, which are necessary to reduce the number of PNO integral blocks, are different than in the DLPNO-CCSD program. In the current work, all these approximations were for the first time tested in detail. It was found that most of them cause surprisingly small errors, in particular for relative energies. The errors depend on the PNO domain sizes, and can become significant if too small domains are used. The significance of all these improvements has been demonstrated for several large systems. Most impressively, they can perhaps be seen in calculations for the interaction energy of the coronene dimer, which is determined by strong dispersion interactions. Pavošević et al. 107 very recently published DLPNO-CCSD(T)-F12 calculations for this system, and reported a computation time for DLPNO-CCSD of 130768 sec (36 h) on 16 cores Intel Xeon [email protected] GHz, using their “tight” settings and the VDZ-F12 basis set. They also showed that due to the LMP2 weak pair 53 ACS Paragon Plus Environment


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

approximation this calculation still overestimates the binding energy by 5.3 kcal mol−1 , relative to an even much more expensive calculation with “verytight” options, in which all non-distant pairs were fully treated by LCCSD, without LMP2 weak pair approximations. We have demonstrated that the errors caused by the new close and weak pair approximations in our method only amount to less than 0.1 kcal mol−1 (using default options). Nevertheless, with the same basis sets, our PNO-LCCSD calculation for the coronene dimer took only 7672 sec (2.1 h) elapsed time on 16 cores Intel Xeon E5-2660 v3 @ 2.60 GHz (using shared files to store the integrals, without the preceding DF-HF calculation). Except for the slightly different clock speeds, the two mentioned processors are similar, but of course the machines may differ in other aspects (memory, disk), and therefore such comparisons on different machines can give only a rough estimate of the relative performance. According to the results presented in this paper, the main limitation of the accuracy of our PNOLCCSD method relative to canonical CCSD with the same basis is due to the finite PAO domains. Probably related to intra- or intermolecular BSSE effects, the convergence of the correlation energy as well as of relative energies with respect to the PAO domain size is sometimes frustratingly slow. This problem is almost entirely eliminated in our corresponding PNO-LCCSD-F12 method, which will be described in a following paper. 184 In that work we will present further extensive benchmark calculations for the overall accuracy of the PNO-LCCSD-F12 method, including direct comparisons with canonical CCSD-F12 calculations for medium size systems. Despite drastic improvements of the accuracy due to the reduction of basis set superposition and domain errors, the computation time for PNO-LCCSD-F12 is (with our program) typically only 20–30% longer than for a corresponding PNO-LCCSD calculation. The exclusive use of the explicitly correlated variant is therefore strongly recommended. Our next important next goal is the development of an efficient and accurate perturbative triples energy correction and the extension to open shell cases. Work in these directions is in progress.


Page 54 of 84

Page 55 of 84

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


Acknowledgement This work has been funded by the ERC Advanced Grant 320723 (ASES).

Associated Content 3D structures of the investigated molecular systems, additional results, and further elaborations on our implemented coupled cluster equations and theory can be found in the Supporting Information. The Supporting Information is available free of charge via the Internet at http://pubs.acs.org

References (1) Kutzelnigg, W. r12 -Dependent terms in the wave function as closed sums of partial wave amplitudes for large l. Theor. Chim. Acta 1985, 68, 445–469. (2) Klopper, W.; Kutzelnigg, W. Møller-Plesset calculations taking care of the correlation CUSP. Chem. Phys. Lett. 1987, 134, 17–22. (3) Klopper, W.; Kutzelnigg, W. MP2-R12 calculations on the relative stability of carbocations. J. Phys. Chem. 1990, 94, 5625–5630. (4) Klopper, W.; Röhse, R.; Kutzelnigg, W. CID and CEPA calculations with linear r12 terms. Chem. Phys. Lett. 1991, 178, 455–461. (5) Kutzelnigg, W.; Klopper, W. Wave functions with terms linear in the interelectronic coordinates to take care of the correlation cusp. I. General Theory. J. Chem. Phys. 1991, 94, 1985–2001. (6) Termath, V.; Klopper, W.; Kutzelnigg, W. Wave functions with terms linear in the interelectronic coordinates to take care of the correlation cusp. II. Second-order Møller-Plesset (MP2-R12) calculations on closed-shell atoms. J. Chem. Phys. 1991, 94, 2002–2019. 55 ACS Paragon Plus Environment


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

(7) Klopper, W.; Kutzelnigg, W. Wave functions with terms linear in the interelectronic coordinates to take care of the correlation cusp. III. Second-order Møller-Plesset (MP2-R12) calculations on molecules of first row atoms. J. Chem. Phys. 1991, 94, 2020–2030. (8) Klopper, W. Orbital-invariant formulation of the MP2-R12 method. Chem. Phys. Lett. 1991, 186, 583–585. (9) Noga, J.; Kutzelnigg, W.; Klopper, W. CC-R12, a correlation cusp corrected coupled-cluster method with a pilot application to the Be2 potential curve. Chem. Phys. Lett. 1992, 199, 497–504. (10) Noga, J.; Kutzelnigg, W. Coupled cluster theory that takes care of the correlation cusp by inclusion of linear terms in the interelectronic coordinates. J. Chem. Phys. 1994, 101, 7738– 7762. (11) Klopper, W.; Samson, C. C. M. Explicitly correlated second-order Møller-Plesset methods with auxiliary basis sets. J. Chem. Phys. 2002, 116, 6397–6410. (12) Manby, F. R. Density fitting in second-order linear-r12 Møller-Plesset perturbation theory. J. Chem. Phys. 2003, 119, 4607–4613. (13) Ten-no, S. Initiation of explicitly correlated Slater-type geminal theory. Chem. Phys. Lett. 2004, 398, 56–61. (14) Ten-no, S. Explicitly correlated second order perturbation theory: Introduction of a rational generator and numerical quadratures. J. Chem. Phys. 2004, 121, 117–129. (15) 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. (16) Valeev, E. F. Improving on the resolution of the identity in linear R12 ab initio theories. Chem. Phys. Lett. 2004, 395, 190–195.


Page 56 of 84

Page 57 of 84

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


(17) Tew, D. P.; Klopper, W. New correlation factors for explicitly correlated electronic wave functions. J. Chem. Phys. 2005, 123, 074101. (18) Kedžuch, S.; Milko, M.; Noga, J. Alternative formulation of the matrix elements in MP2R12 theory. Int. J. Quantum Chem. 2005, 105, 929–936. (19) 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. (20) Fliegl, H.; Hättig, C.; Klopper, W. Inclusion of the (T) triples Correction into the linearr12 corrected coupled-cluster model CCSD(R12). Int. J. Quantum Chem. 2006, 106, 2306– 2317. (21) Manby, F. R.; Werner, H.-J.; Adler, T. B.; May, A. J. Explicitly correlated local secondorder perturbation theory with a frozen geminal correlation factor. J. Chem. Phys. 2006, 124, 094103. (22) Werner, H.-J.; Adler, T. B.; Manby, F. R. General orbital invariant MP2-F12 theory. J. Chem. Phys. 2007, 126, 164102. (23) 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. (24) Adler, T. B.; Knizia, G.; Werner, H.-J. A simple and efficient CCSD(T)-F12 approximation. J. Chem. Phys. 2007, 127, 221106. (25) Ten-no, S. A simple F12 geminal correction in multi-reference perturbation theory. Chem. Phys. Lett. 2007, 447, 175–179. (26) 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.



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

(27) Knizia, G.; Werner, H.-J. Explicitly correlated RMP2 for high-spin open-shell reference states. J. Chem. Phys. 2008, 128, 154103. (28) 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. (29) Shiozaki, T.; Kamiya, M.; Hirata, S.; Valeev, E. F. Equations of explicitly-correlated coupled-cluster methods. Phys. Chem. Chem. Phys. 2008, 10, 3358–3370. (30) 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. (31) Tew, D. P.; Klopper, W.; Hättig, C. A diagonal orbital-invariant explicitly correlated coupled-cluster method. Chem. Phys. Lett. 2008, 452, 326–332. (32) Valeev, E. F. Coupled-cluster methods with perturbative inclusion of explicitly correlated terms: a preliminary investigation. Phys. Chem. Chem. Phys. 2008, 10, 106–113. (33) 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. (34) Torheyden, M.; Valeev, E. F. Variational formulation of perturbative explicitly-correlated coupled-cluster methods. Phys. Chem. Chem. Phys. 2008, 10, 3410–3420. (35) 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. (36) 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.


Page 58 of 84

Page 59 of 84

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


(37) Knizia, G.; Adler, T. B.; Werner, H.-J. Simplified CCSD(T)-F12 methods: Theory and benchmarks. J. Chem. Phys. 2009, 130, 054104. (38) Shiozaki, T.; Kamiya, M.; Hirata, S.; Valeev, E. F. Higher-order explicitly correlated coupled-cluster methods. J. Chem. Phys. 2009, 130, 054101. (39) Bokhan, D.; Bernadotte, S.; Ten-no, S. Implementation of the CCSD(T)(F12) method using numerical quadratures. Chem. Phys. Lett. 2009, 469, 214–218. (40) Bischoff, F. A.; Wolfsegger, S.; Tew, D. P.; Klopper, W. Assessment of basis sets for F12 explicitly-correlated molecular electronic-structure methods. Mol. Phys. 2009, 107, 963– 975. (41) 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. (42) Werner, H.-J.; Knizia, G.; Manby, F. R. Explicitly correlated coupled cluster methods with pair-specific geminals. Mol. Phys 2011, 109, 407–417. (43) Shiozaki, T.; Knizia, G.; Werner, H.-J. Explicitly correlated multireference configuration interaction: MRCI-F12. J. Chem. Phys. 2011, 134, 034113. (44) Pulay, P. Localizability of dynamic electron correlation. Chem. Phys. Lett. 1983, 100, 151– 154. (45) Saebø, S.; Pulay, P. Local Configuration-Interaction: An efficient approach for larger molecules. Chem. Phys. Lett. 1985, 113, 13–18. (46) Pulay, P.; Saebø, S. Orbital-invariant formulation and second-order gradient evaluation in Møller-Plesset perturbation theory. Theor. Chim. Acta 1986, 69, 357–368. (47) 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. 59 ACS Paragon Plus Environment


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

(48) Saebø, S.; Pulay, P. The local correlation treatment. II. Implementation and tests. J. Chem. Phys. 1988, 88, 1884–1890. (49) Saebø, S.; Pulay, P. Local treatment of electron correlation. Annu. Rev. Phys. Chem. 1993, 44, 213–236. (50) Hampel, C.; Werner, H.-J. Local treatment of electron correlation in coupled cluster theory. J. Chem. Phys. 1996, 104, 6286–6297. (51) Maslen, P. E.; Head-Gordon, M. Non-iterative local second order Moller-Plesset theory. Chem. Phys. Lett. 1998, 283, 102–180. (52) 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. (53) 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. (54) 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. (55) Hetzer, G.; Pulay, P.; Werner, H.-J. Multipole approximation of distant pair energies in local MP2 calculations. Chem. Phys. Lett. 1998, 290, 143–149. (56) 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. (57) 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. (58) 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. 60 ACS Paragon Plus Environment

Page 60 of 84

Page 61 of 84

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


(59) Schütz, M.; Werner, H.-J. Local perturbative triples correction (T) with linear cost scaling. Chem. Phys. Lett. 2000, 318, 370–378. (60) 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. (61) 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. (62) 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. (63) Schütz, M. A new, fast, semi-direct implementation of Linear Scaling Local Coupled Cluster Theory. Phys. Chem. Chem. Phys. 2002, 4, 3941–3947. (64) 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. (65) 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. (66) 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. (67) Auer, A.; Nooijen, M. Dynamically screened local correlation method using enveloping localized orbitals. J. Chem. Phys. 2006, 125, 024104. (68) 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. (69) 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. 61 ACS Paragon Plus Environment


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

(70) Mata, R.; Werner, H.-J. Calculation of smooth potential energy surfaces using local electron correlation methods. J. Chem. Phys. 2006, 125, 184110. (71) Mata, R.; Werner, H.-J. Local correlation methods with a natural localized molecular orbital basis. Mol. Phys. 2007, 105, 2753–2761. (72) 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. (73) Mata, R.; Werner, H.-J.; Schütz, M. Correlation regions within a localized molecular orbital approach. J. Chem. Phys. 2008, 128, 144106. (74) Werner, H.-J.; Schütz, M. An efficient local coupled cluster method for accurate thermochemistry of large systems. J. Chem. Phys. 2011, 135, 144116. (75) Kats, D.; Manby, F. R. Sparse tensor framework for implementation of general local correlation methods. J. Chem. Phys. 2013, 138, 144101. (76) Kats, D. Speeding up local correlation methods. J. Chem. Phys. 2014, 141, 244101. (77) Kats, D. Speeding up local correlation methods: System-inherent domains. J. Chem. Phys. 2016, 145, 014103. (78) 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. (79) 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. (80) Ma, Q.; Werner, H.-J. Scalable electron correlation methods. II. Parallel PNO-LMP2-F12 with near linear scaling in the molecular size. J. Chem. Theory Comput. 2015, 11, 5291– 5304. 62 ACS Paragon Plus Environment

Page 62 of 84

Page 63 of 84

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


(81) Schwilk, M.; Usvyat, D.; Werner, H.-J. Communication: Improved pair approximations in local coupled-cluster methods. J. Chem. Phys. 2015, 142, 121102. (82) 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. (83) Werner, H.-J. Communication: Multipole approximations of distant pair energies in local correlation methods with pair natural orbitals. J. Chem. Phys. 2016, 145, 201101. (84) Russ, N. J.; Crawford, T. D. Local correlation in coupled cluster calculations of molecular response properties. Chem. Phys. Lett. 2004, 400, 104–111. (85) 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. (86) 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. (87) Subotnik, J. E.; Head-Gordon, M. A local correlation model that yields intrinsically smooth potential-energy surfaces. J. Chem. Phys. 2005, 123, 064108. (88) Lawler, K. V.; Parkhill, J. A.; Head-Gordon, M. Penalty functions for combining coupledcluster and perturbation amplitudes in local correlation methods with optimized orbitals. Mol. Phys. 2008, 106, 2309–2324. (89) Schweizer, S.; Kussmann, J.; Doser, B.; Ochsenfeld, C. Linear-scaling Cholesky decomposition. J. Comp. Chem. 2008, 29, 1004–1010. (90) Schweizer, S.; Doser, B.; Ochsenfeld, C. An atomic orbital-based reformulation of energy



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

gradients in second-order Møller-Plesset perturbation theory. J. Chem. Phys. 2008, 128, 154101. (91) Doser, B.; Lambrecht, D. S.; Kussmann, J.; Ochsenfeld, C. Linear-scaling atomic orbitalbased second-order Møller–Plesset perturbation theory by rigorous integral screening criteria. J. Chem. Phys. 2009, 130, 064107. (92) 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. (93) Kats, D.; Schütz, M. A multistate local coupled cluster CC2 response method based on the Laplace transform. J. Chem. Phys. 2009, 131, 124117. (94) 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. (95) Maschio, L. Local MP2 with density fitting for periodic systems: A parallel implementation. J. Chem. Theory Comput. 2011, 7, 2818–2830. (96) 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. (97) 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. (98) 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. 64 ACS Paragon Plus Environment

Page 64 of 84

Page 65 of 84

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


(99) 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. (100) 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. (101) Riplinger, C.; Neese, F. An efficient and near linear scaling pair natural orbital based local coupled cluster method. J. Chem. Phys. 2013, 138, 034106. (102) 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. (103) 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. (104) 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. (105) 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. (106) 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. (107) 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 65 ACS Paragon Plus Environment


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

explicitly correlated coupled-cluster method with pair natural orbitals. J. Chem. Phys. 2017, 146, 174108. (108) Walter, D.; Venkatnathan, A.; Carter, E. Local correlation in the virtual space in multireference singles and doubles configuration interaction. J. Chem. Phys. 2003, 118, 8127–8139. (109) Chwee, T. S.; Szilva, A. B.; Lindh, R.; Carter, E. A. Linear scaling multireference singles and doubles configuration interaction. J. Chem. Phys. 2008, 128, 224106. (110) Chwee, T. S.; Carter, E. A. Density fitting of two-electron integrals in local multireference single and double excitation configuration interaction calculations. Mol. Phys. 2010, 108, 2519–2526. (111) Chwee, T. S.; Carter, E. A. Cholesky decomposition within local multireference singles and doubles configuration interaction. J. Chem. Phys. 2010, 132, 074104. (112) Krisiloff, D. B.; Carter, E. A. Approximately size extensive local multireference singles and doubles configuration interaction. Phys. Chem. Chem. Phys. 2012, 14, 7710–7717. (113) Krisiloff, D. B.; Krauter, C. M.; Ricci, F. J.; Carter, E. A. Density fitting and cholesky decomposition of the two-electron integrals in local multireference configuration interaction theory. J. Chem. Theory Comput. 2015, 11, 5242–5251. (114) Demel, O.; Pittner, J.; Neese, F. A local pair natural orbital-based multireference Mukherjee’s coupled cluster method. J. Chem. Theory Comput. 2015, 11, 3104–3114. (115) Guo, Y.; Sivalingam, K.; Valeev, E. F.; Neese, F. Sparse maps-A systematic infrastructure for reduced-scaling electronic structure methods. III. Linear-scaling multireference domainbased pair natural orbital N-electron valence perturbation theory. J. Chem. Phys. 2016, 144, 094111. (116) Menezes, F.; Kats, D.; Werner, H.-J. Local complete active space second-order perturbation theory using pair natural orbitals (PNO-CASPT2). J. Chem. Phys. 2016, 145, 124115. 66 ACS Paragon Plus Environment

Page 66 of 84

Page 67 of 84

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


(117) Fedorov, D. G.; Kitaura, K. Coupled-cluster theory based upon the fragment molecularorbital method. J. Chem. Phys. 2005, 123, 134103. (118) Fedorov, A.; Couzijn, E. P. A.; Nagornova, N. S.; Boyarkin, O. V.; Rizzo, T. R.; Chen, P. Structure and bonding of isoleptic coinage metal (Cu, Ag, Au) dimethylaminonitrenes in the gas phase. J. Am. Chem. Soc. 2010, 132, 13789–13798. (119) Hughes, T. F.; Flocke, N.; Bartlett, R. J. Natural linear-scaled coupled-cluster theory with local transferable triple excitations: Applications to peptides. J. Phys. Chem. A 2008, 112, 5994–6003. (120) Gordon, M. S.; Mullin, J. M.; Pruitt, S. R.; Roskop, L. B.; Slipchenko, L. V.; Boatz, J. A. Accurate methods for large molecular systems. J. Phys. Chem. B 2009, 113, 9646–9663. (121) Gordon, M. S.; Fedorov, D. G.; Pruitt, S. R.; Slipchenko, L. V. Fragmentation methods: A route to accurate calculations on large systems. Chem. Rev. 2012, 112, 632–672. (122) Findlater, A. D.; Zahariev, F.; Gordon, M. S. Combined fragment molecular orbital cluster in molecule approach to massively parallel electron correlation calculations for large systems. J. Phys. Chem. A 2015, 119, 3587–3593. (123) Ziółkowski, M.; Jansìk, B.; Kjærgaard, T.; Jørgensen, P. Linear scaling coupled cluster method with correlation energy based error control. J. Chem. Phys. 2010, 133, 014107. (124) Kristensen, K.; Ziółkowski, M.; Jansìk, B.; Kjærgaard, T.; Jørgensen, P. A locality analysis of the divide-expand-consolidate coupled cluster amplitude equations. J. Chem. Theory Comput. 2011, 7, 1677–1694. (125) Høyvik, I.-M.; Kristensen, K.; Jansìk, B.; Jørgensen, P. The divide-expand-consolidate family of coupled cluster methods: Numerical illustrations using second order Møller-Plesset perturbation theory. J. Chem. Phys. 2012, 136, 014105.



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) Kristensen, K.; Jørgensen, P.; Jansìk, B.; Kjærgaard, T.; Reine, S. Molecular gradient for second-order Møller-Plesset perturbation theory using the divide-expand-consolidate (DEC) scheme. J. Chem. Phys. 2012, 137, 114102. (127) Kristensen, K.; Kjærgaard, T.; Høyvik, I.-M.; Ettenhuber, P.; Jørgensen, P.; Jansik, B.; Reine, S.; Jakowski, J. The divide-expand-consolidate MP2 scheme goes massively parallel. Mol. Phys. 2013, 111, 1196–1210. (128) Eriksen, J. J.; Baudin, P.; Ettenhuber, P.; Kristensen, K.; Kjærgaard, T.; Jørgensen, P. Linearscaling coupled cluster with perturbative triple excitations: The divide-expand-consolidate CCSD(T) model. J. Chem. Theory Comput. 2015, 11, 2984–2993. (129) Baudin, P.; Ettenhuber, P.; Reine, S.; Kristensen, K.; Kjaergaard, T. Efficient linear-scaling second-order Møller-Plesset perturbation theory: The divide-expand-consolidate RI-MP2 model. J. Chem. Phys. 2016, 144, 054102. (130) Bykov, D.; Kristensen, K.; Kjaergaard, T. The molecular gradient using the divide-expandconsolidate resolution of the identity second-order Møller-Plesset perturbation theory: The DEC-RI-MP2 gradient. J. Chem. Phys. 2016, 145, 024106. (131) Ettenhuber, P.; Baudin, P.; Kjaergaard, T.; Jorgensen, P.; Kristensen, K. Orbital spaces in the divide-expand-consolidate coupled cluster method. J. Chem. Phys. 2016, 144, 164116. (132) Li, W.; Piecuch, P.; Gour, J. R.; Li, S. Local correlation calculations using standard and renormalized coupled-cluster approaches. J. Chem. Phys. 2009, 131, 114109. (133) Li, W.; Piecuch, P. Improved design of orbital domains within the cluster-in-molecule local correlation framework: Single-environment cluster-in-molecule ansatz and its application to local coupled-cluster approach with singles and doubles. J. Phys. Chem. A 2010, 114, 8644–8657.


Page 68 of 84

Page 69 of 84

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


(134) Li, W.; Piecuch, P. Multilevel extension of the cluster-in-molecule local correlation methodology: Merging coupled-cluster and Møller-Plesset perturbation theories. J. Phys. Chem. A 2010, 114, 6721–6727. (135) Li, W.; Guo, Y.; Li, S. A refined cluster-in-molecule local correlation approach for predicting the relative energies of large systems. Phys. Chem. Chem. Phys. 2012, 14, 7854–7862. (136) Guo, Y.; Li, W.; Li, S. Improved cluster-in-molecule local correlation approach for electron correlation calculation of large systems. J. Phys. Chem. A 2014, 118, 8996–9004. (137) Li, W.; Chen, C.; Zhao, D.; Li, S. LSQC: Low scaling quantum chemistry program. Int. J. Quantum Chem. 2015, 115, 641–646. (138) Li, W.; Ni, Z.; Li, S. Cluster-in-molecule local correlation method for post-Hartree-Fock calculations of large systems. Mol. Phys. 2016, 114, 1447–1460. (139) Rolik, Z.; Kallay, M. A general-order local coupled-cluster method based on the cluster-inmolecule approach. J. Chem. Phys. 2011, 135, 104111. (140) Rolik, Z.; Szegedy, L.; Ladjánszki, I.; Ladóczki, B.; Kallay, M. An efficient linear-scaling CCSD(T) method based on local natural orbitals. J. Chem. Phys. 2013, 139, 094105. (141) Kallay, M. Linear-scaling implementation of the direct random-phase approximation. J. Chem. Phys. 2015, 142, 204105. (142) Nagy, P. R.; Samu, G.; Kallay, M. An integral-direct linear-scaling second-order MøllerPlesset approach. J. Chem. Theory Comput. 2016, 12, 4897–4914. (143) Saha, A.; Raghavachari, K. Dimers of Dimers (DOD): A new fragment-based method applied to large water clusters. J. Chem. Theory Comput. 2014, 10, 58–67. (144) Yoshikawa, T.; Nakai, H. Linear-scaling self-consistent field calculations based on divideand-conquer method using resolution-of-identity approximation on graphical processing units. J. Comput. Chem. 2015, 36, 164–170. 69 ACS Paragon Plus Environment


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

(145) Stoll, H. On the correlation energy of graphite. J. Chem. Phys. 1992, 97, 8449–8454. (146) Friedrich, J.; Hanrath, M.; Dolg, M. Fully automated implementation of the incremental scheme: Application to CCSD energies for hydrocarbons and transition metal compounds. J. Chem. Phys. 2007, 126, 154110. (147) Friedrich, J.; Hanrath, M.; Dolg, M. Evaluation of incremental correlation energies for openshell systems: Application to the intermediates of the 4-exo cyclization, Arduengo carbenes and an anionic water cluster. J. Phys. Chem. A 2008, 112, 8762–8766. (148) Friedrich, J.; Dolg, M. Implementation and performance of a domain-specific basis set incremental approach for correlation energies: Applications to hydrocarbons and a glycine oligomer. J. Chem. Phys. 2008, 129, 244105. (149) Friedrich, J.; Coriani, S.; Helgaker, T.; Dolg, M. Implementation of the incremental scheme for one-electron first-order properties in coupled-cluster theory. J. Chem. Phys. 2009, 131, 154102. (150) Friedrich, J.; Dolg, M. Fully automated incremental evaluation of MP2 and CCSD(T) energies: Application to water clusters. J. Chem. Theory Comput. 2009, 5, 287–294. (151) 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. (152) Meitei, O. R.; Hesselmann, A. Molecular energies from an incremental fragmentation method. J. Chem. Phys. 2016, 144, 084109. (153) Hoyau, S.; Maynau, D.; Malrieu, J.-P. A regionally contracted multireference configuration interaction method: General theory and results of an incremental version. J. Chem. Phys. 2011, 134, 054125. (154) Werner, H.-J. Eliminating the domain error in local explicitly correlated second-order Møller-Plesset perturbation theory. J. Chem. Phys. 2008, 129, 101103. 70 ACS Paragon Plus Environment

Page 70 of 84

Page 71 of 84

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


(155) 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. (156) Adler, T. B.; Werner, H.-J.; Manby, F. R. Local explicitly correlated second-order perturbation theory for the accurate treatment of large molecules. J. Chem. Phys. 2009, 130, 054106. (157) 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. (158) 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. (159) Helmich, B.; Hättig, C. Local pair natural orbitals for excited states. J. Chem. Phys. 2011, 135, 214106. (160) 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. (161) 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. (162) Schmitz, G.; Hättig, C.; Tew, D. P. Explicitly correlated PNO-MP2 and PNO-CCSD and their application to the S66 set and large molecular systems. Phys. Chem. Chem. Phys. 2014, 16, 22167–22178. (163) 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.



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

(164) Li, W. Linear scaling explicitly correlated MP2-F12 and ONIOM methods for the longrange interactions of the nanoscale clusters in methanol aqueous solutions. J. Chem. Phys. 2013, 138, 014106. (165) Wang, Y. M.; Hättig, C.; Reine, S.; Valeev, E.; Kjaergaard, T.; Kristensen, K. Explicitly correlated second-order Møller-Plesset perturbation theory in a Divide-Expand-Consolidate (DEC) context. J. Chem. Phys. 2016, 144, 204112. (166) 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 secondorder explicitly correlated energy with pair natural orbitals. J. Chem. Phys. 2016, 144, 144109. (167) Russ, N. J.; Crawford, T. D. Potential energy surface discontinuities in local correlation methods. J. Chem. Phys. 2004, 121, 691–696. (168) Subotnik, J. E.; Sodth, A.; Head-Gordon, M. A near linear-scaling smooth local coupled cluster algorithm for electronic structure. J. Chem. Phys. 2006, 125, 074116. (169) Meyer, W. Ionization energies of water from PNO-CI calculations. Int. J. Quantum Chem. 1971, 5, 341–348. (170) 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. (171) 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.


Page 72 of 84

Page 73 of 84

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


(172) Staemmler, V.; Jaquet, R. CEPA calculations on open-shell molecules. I. Outline of the method. Theor. Chim. Acta. 1981, 59, 487–500. (173) 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. (174) Taylor, P. R. A rapidly convergent CI expansion based on several reference configurations, using optimized correlating orbitals. J. Chem. Phys. 1981, 74, 1256–1270. (175) Liakos, D. G.; Neese, F. Is it possible to obtain coupled cluster quality energies at near density functional theory cost? Domain-based local pair natural orbital coupled cluster vs modern density functional theory. J. Chem. Theory Comput. 2015, 11, 4054–4063. (176) 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. (177) Yang, J.; Kurashige, Y.; Manby, F. R. Tensor factorizations of local second-order Møller– Plesset theory. J. Chem. Phys. 2011, 134, 044123. (178) Kurashige, Y.; Yang, J.; Chan, G. K. L.; Manby, F. R. Optimization of orbital-specific virtuals in local Møller-Plesset perturbation theory. J. Chem. Phys. 2012, 136, 124106. (179) Yang, J.; Chan, G. K.-L.; Manby, F. R.; Schütz, M.; Werner, H.-J. The orbital-specific-virtual local coupled cluster singles and doubles method. J. Chem. Phys. 2012, 136, 144105. (180) 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. (181) Masur, O.; Usvyat, D.; Schütz, M. Efficient and accurate treatment of weak pairs in local CCSD(T) calculations. J. Chem. Phys. 2013, 139, 164116.



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

(182) Schütz, M.; Masur, O.; Usvyat, D. Efficient and accurate treatment of weak pairs in local CCSD(T) calculations. II. Beyond the ring approximation. J. Chem. Phys. 2014, 140, 244107. (183) Schütz, M.; Maschio, L.; Karttunen, A. J.; Usvyat, D. Exfoliation Energy of Black Phosphorus Revisited: A Coupled Cluster Benchmark. J. Phys. Chem. Lett. 2017, 8, 1290–1294. (184) Ma, Q.; Schwilk, M.; Köppl, C.; Werner, H.-J. To be published. (185) Janowski, T.; Pulay, P. Efficient parallel implementation of the CCSD external exchange operator and the perturbative triples (T) energy calculation. J. Chem. Theory Comput. 2008, 4, 1585–1592. (186) 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. (187) Luo, S.; Zhao, Y.; Truhlar, D. G. Validation of electronic structure methods for isomerization reactions of large organic molecules. Phys. Chem. Chem. Phys. 2011, 13, 13683–13689. (188) Rayne, S.; Forest, K. A comparative examination of density functional performance against the ISOL24/11 isomerization energy benchmark. Comput. Theor. Chem. 2016, 1090, 147– 152. (189) Weymuth, T.; Couzijn, E. P. A.; Chen, P.; Reiher, M. New Benchmark Set of TransitionMetal Coordination Reactions for the Assessment of Density Functionals. J. Chem. Theory Comput. 2014, 10, 3092–3103. (190) 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. (191) Zuend, S. J.; Jacobsen, E. N. Mechanism of amido-thiourea catalyzed enantioselective imine


Page 74 of 84

Page 75 of 84

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


hydrocyanation: Transition state stabilization via multiple non-covalent interactions. J. Am. Chem. Soc. 2009, 131, 15358. (192) 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. (193) 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. (194) 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. (195) 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. (196) 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. (197) 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. (198) 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. (199) Weigend, F. A fully direct RI-HF algorithm: Implementation, optimised auxiliary basis sets, demonstration of accuracy and efficiency. Phys. Chem. Chem. Phys. 2002, 4, 4285–4291. (200) Knizia, G. Intrinsic atomic orbitals: an unbiased bridge between quantum theory and chemical concepts. J. Chem. Theory Comput. 2013, 9, 4834–4843. 75 ACS Paragon Plus Environment


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

(201) 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. (202) Werner, H.-J.; Knowles, P. J.; Knizia, G.; Manby, F. R.; Schütz, M. Molpro: A generalpurpose quantum chemistry program package. WIREs Comput. Mol. Sci. 2011, 2, 242–253. (203) Curtiss, L. A.; Raghavachari, K.; Redfern, P. C.; Rassolov, V.; Pople, J. A. Gaussian-3 (G3) theory for molecules containing first and second-row atoms. J. Chem. Phys. 1998, 109, 7764–7776. (204) Boese, A. D.; Oren, M.; Atasoylu, O.; Martin, J. M. L.; Kállay, M.; Gauss, J. W3 theory: Robust computational thermochemistry in the kJ/mol accuracy range. J. Chem. Phys. 2004, 120, 4129–4141. (205) Saebø, S.; Tong, W.; Pulay, P. Efficient elimination of basis set superposition errors by the local correlation method: Accurate ab initio studies of the water dimer. J. Chem. Phys. 1993, 98, 2170–2175. (206) Grimme, S. Improved second-order Møller-Plesset perturbation theory by separate scaling of parallel- and antiparallel-spin pair correlation energies. J. Chem. Phys. 2003, 118, 9095– 9102. (207) Jung, Y.; Lochan, R. C.; Dutoi, A. D.; Head-Gordon, M. Scaled opposite-spin second order Møller–Plesset correlation energy: An economical electronic structure method. J. Chem. Phys. 2004, 121, 9793–9802. (208) Köppl, C.; Werner, H.-J. Parallel and low-order scaling implementation of Hartree-Fock exchange using local density fitting. J. Chem. Theory Comput. 2016, 12, 3122–3134. ˇ (209) 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. 76 ACS Paragon Plus Environment

Page 76 of 84

Page 77 of 84

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


(210) 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. (211) Karypis, G.; Kumar, V. A fast and high quality multilevel scheme for partitioning irregular graphs. SIAM J. Sci. Comput. 1998, 20, 359–392.

A

The linked PNO-LCCSD equations with pair approximations

The PNO-LCCSD equations can be compactly described in terms of dressed orbitals |˜ii = |ii +

X

|ci tci ii

(48)

c

|¯aii i = |aii i −

X

|ki t¯akii

(49)

k

In practice, the dressed orbitals are generated in every iteration by subsequent backtransformation i of the singles amplitudes t[ii] from PNO domain space aii to PAO domain space ri :

h i ¯ ii ti tri i = Qi W i ii aii ra

(50)

¯ i j = Wi j Vi j Ui j (cf. eqs 30 – 32). Dressed PAOs P¯ µr and LMOs L˜ µi are then computed in where W the AO basis: L˜ µi = Lµi + Pµri tri i

(51)

P¯ µr = Pµr − Lµk S rsk tksk

(52)

We introduce the following dressed integrals and integral contractions ˜ i j ]ab = (¯ai j˜i|b¯ i j ˜j) [K [i j,i j] 77 ACS Paragon Plus Environment

(53)


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

h

¯ ij K [kl,mn]

h

j J¯ i[kl,mn]

i

i

= (akl i|b¯ mn ˜j)

(54)

= (akl b¯ mn |i ˜j)

(55)

ab

ab

i h ¯ i j )[i j,i j] = (¯ai j ci j |di j b¯ i j )T i ijj i j K(T c d

(56)

ab

h

¯ i j )(k) K(T [mn]

i

F˜ rs = hrs +

Page 78 of 84

= (¯amn ci j |di j k)T ci ijj di j

(57)

Xh i ˜ − (rl|ls) ˜ 2(rs|ll)

(58)

a

l

h

F¯ [i j,i j]

i

ab

= F˜ ai j b¯ i j

(59)

h i i = F˜ ˜i¯akl f¯[kl]

(60)

h i i f˜[kl] = F˜ akl i

(61)

a

a

h

˜

jk ki[mn]

i

a

˜ = (amn i| jk)

(62)

f˜i j = F˜ ˜i j

(63)

˜˜ Kkli j = (˜ik| ˜jl)

(64)

In the following equations, we will use underbraces to indicate the condition with which the term is included. We will use {i} s to represent a domain of LMOs k where ik is strong, and similarly {i}c for a domain of k’s where ik is strong or close. Given that full PNO-LCCSD calculations are not performed for distant pairs, all summations are always implicitly restricted so that no pairs of LMOs that appear in the equations as amplitudes or integral indices are distant. In terms of the quantities defined above, the LCCSD doubles residual takes the same form as the LCCD equations: i k ¯ T˜ ik[ik,ik] )(k) +K( ri = f¯[ii] +T˜ ik[ii,ik] f˜[ik] | {z } | {z [ii] } k∈{i}w

k∈{i}w

˜

˜ kl[ii,kl] klk[kl]i −T | {z }

k∈{i}c ∩{l}w ,l∈{i}w ∩{k}w


(65)

Page 79 of 84

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


˜ i j + K(T ¯ i j )[i j,i j] + Gi j + G ji† Ri j =K [i j,i j] [i j,i j] [i j,i j] ¯ kl[i j,i j] +αi j,kl T | {z }

(66)

(k∈{i} s ∧l∈{ j} s )∨(k∈{ j} s ∧l∈{i} s ) ∧k∈{l} s ∧l∈{k} s

with 1 ¯ ki jk i jk Gi[ij j,i j] =Ti j Xi[ij j,i j] + T˜ ik[i j,ik] Yi[ik,i j] − T[i j,ik] Z[ik,i j] 2 jk † ¯ ik − (T¯ ki[i j,ik] Zi[ik,i j] ) − βi jk T[i j,i j]

(67)

The intermediates are defined as ˜˜

βi jk =

αi j,kl = Kkli j + tr(Ti j Klk[i j,i j] )

(68)

f˜jk |{z}

(69)

k∈{ j}w ∩{i}w

Xi[ij j,i j] = F¯ [i j,i j]

+ tr(Lkl[ jl, jl] Tl j ) | {z }

k∈{ j}c ∩{i}w ∩{l}w ,l∈{ j}w ∩{i}w ∩{k}w

¯ lk −Lkl[i j,kl] T | {z [kl,i}j]

(70)

k∈{ j}c ∩{i}w ∩{l}w ,l∈{ j}w ∩{i}w ∩{k}w

1 kl ˜ l j 1 ¯k j jk ¯ kj Yi[ik,i j] = K[ik,i j] − J[ik,i j] + L[ik, jl] T[ jl,i j] 2 | 4 {z }

(71)

k∈{i}w ∩{l}w ,l∈{ j}w ∩{k}w

1 lk ˜ jl jk ¯k j Zi[ik,i j] = J[ik,i j] − K[ik, jl] T[ jl,i j] | 2 {z }

(72)

k∈{i}c ∪{ j}c ,l∈{ j}c ∪{i}c

In practice the density fitting evaluation of some dressed integrals requires 3-index integrals with two external labels to be computed in each iteration, which is prohibitively expensive. We 79 ACS Paragon Plus Environment


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 80 of 84

therefore expand these dressed integral (contractions) by ¯ T˜ ik )(k) = K(T˜ ik )(k) − tr(Llk[ik,ik] Tki[ik,ik] )¯tl[ii] K( [ii] [ii] | {z }

(73)

l∈{i} s ∩{k}w ,k∈{i}w ∩{l}w

¯ i j )[i j,i j] =K(Ti j )[i j,i j] K(T

¯ kl[i j,i j] + tr(Klk[i j,i j] Ti j )D | {z }

same restriction as for αi j,kl in eq 66

h i k† † i j (k) ¯k† ¯ − K(T ji )(k) t [i j] [i j] − K(T )[i j] t[i j] | {z }

(74)

k∈({i} s ∪{ j} s )∩({i}c ∩{ j}c )

lk j ¯l† j kj lk kj ¯ jl J¯ k[ik,i j] = J[ik,i j] +J(E )[ik,i j] − k[ik] t[i j] − K[ik, j j] D[ j j,i j] {z } |{z} | k∈{ j}c ∩{i}w

(75)

k,l∈({i} s ∪{ j} s )∩({i}c ∩{ j}c )∧k∈{l}w ∧l∈{k}w

where we introduced additional pair approximations on individual terms depending on their orders and asymptotic behaviors. Our most accurate implementation of the residuals uses this formulation. We refer to this implementation as the “dressed” version of the residuals. We notice that the last term in eq 75 contributes only to MP7 in its lowest order but requires a significant number of 2-external integrals. We therefore introduce a projection approximation ¯ jl ≈ Klk[ik, jl] D ¯ jl Klk[ik, j j] D [ j j,i j] [ jl,i j]

(76)

and the integrals required for this term are identical to those required in the last terms of eqs 71 and 72. With this form of the residuals, it is necessary to recompute the dressed Fock matrix as well as ¯ k j , kkl˜i , and K ˜i ˜j in each iteration. They can be computed similarly to ˜ i j, K the dressed integrals K kl the standard integrals by a partial integral transformation. Since the dressed orbitals are less sparse than the standard LMOs or PAOs, low-order scaling is difficult to achieve. To avoid this problem as well as reducing the computational cost, the default in our program is to expand the dressed 80 ACS Paragon Plus Environment

Page 81 of 84

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


integrals entirely: ¯ k j = Kk j K [ik,i j] [ik,i j] |{z}

k∈{i}w ∩{ j}w

+K(Ek j )[ik,i j] −

X

kkl[ik]j ¯tl† [i j] −

X

(77)

l

l

|

¯ jl Kkl[ik, jl] D [ jl,i j]

{z

}

k∈{i}c ∩{ j}c ,l∈{ j} s ∩{i}c ∩{k}c

˜

klk[kl]i = klki [kl]

(78)

+Klk[kl,ii] ti[ii] | {z }

k∈{i}c ∩{l}w ,l∈{i}w ∩{k}w

˜˜

j† lki kl j lk ¯ ij Kkli j = Kkli j + ti† [ii] k[ii] + t[ j j] k[ j j] − tr(D[i j,i j] K[i j,i j] )

(79)

˜ i j =Ki j + K(D ¯ i j ) + G′i j + G′ ji† K [i j,i j] [i j,i j] [i j,i j] [i j,i j] [i j,i j] X ˜˜ ¯ kl[i j,i j] + Kkli j D

(80)

kl

|

{z

}

same restriction as for αi j,kl in eq 66

where we used an additional intermediate G′i[i j,ij j] =K(Ei j )[i j,i j]

−

X

kik[i j]j ¯tk† [i j]

−

¯ jk Kik[i j, jk] D [ jk,i j]

k

k

|

X

{z

}

|

{z

}

k∈({i}c ∪{ j}c )∩({i}w ∩{ j}w ) k∈({i} s ∪{ j} s )∩({i}c ∩{ j}c )

−

X k

|

¯ ik[ik,i j] − J[ijkj,ik] D

X

¯ i j )(k) ¯tk† K(D [i j,i j] [i j] [i j]

{z k

k∈({i} s ∪{ j} s )∩({i}c ∩{ j}c )

(81)

}

We introduced additional pair and projection approximations in these integral expansions. The additional projection approximations are all in terms at least quadratic in the singles amplitudes and therefore are of MP6 or higher order. The errors introduced by these approximations are



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

discussed in Section 6. We refer to the residuals with these approximated dressed integrals as the “undressed” version. The only mandatory dressed integrals used are the Fock matrices. The “undressed” equations require many additional 3-external integrals as well as slightly more blocks of 2-external integrals, which leads to higher disk requirements. This problem can be alleviated by introducing additional projection approximations, which have been discussed in Section 6. We note that, with full domains and all pairs included, both the “dressed” and the “undressed” versions are equivalent to the canonical ones, and it has been tested that in each case exactly the same results are obtained. The summation restrictions given above are the ones used in the strong pair residual. The detailed summation restrictions for the close and weak pair residuals are given in the Supporting Information.


Page 82 of 84

Page 83 of 84

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


PNO-LCCSD / cc-pVTZ-F12

. 92 atoms . 3345 basis functions . 120 cores on 6 nodes Compute time 3.2 minutes PNO-LMP2 1.5 hours PNO-LCCSD

Figure 13: For Table of Contents Only


PNO-LCCSD / cc-pVTZ-F12

Journal of Chemical Theory and Page Computation 84 of 84

1 2 3 4

. 92 atoms . 3345 basis functions . 120 cores on 6 nodes

ACS Paragon PlusCompute Environment time 3.2 minutes PNO-LMP2 1.5 hours PNO-LCCSD

Scalable Electron Correlation Methods. 3. Efficient and Accurate

Recommend Documents