Sparse maps—A systematic infrastructure for reduced-scaling electronic structure methods. II. Linear scaling domain based pair natural orbital coupled cluster theory

The exact solution of the time-independent many-particle Schrödinger equation has been the “holy grail” of quantum chemical method development for many decades. It is well-known that wave function based ab initio methods are a systematic way to approach this goal.¹ For the significant subset of systems for which the mean field Hartree-Fock method provides a qualitatively valid starting point, there appears to be consensus that coupled cluster (CC) theory is the best way to tackle this problem. Coupled cluster theory, while not possessing the desirable property to be variational, is size consistent. The expansion of the cluster operator in terms of single-, double-, triple, and higher replacements leads to a series of coupled cluster approximations that rapidly converges towards the exact solution, the full configuration interaction method.² Furthermore, the theoretical apparatus to treat molecular properties and excited states is well-defined and well-understood.^3–5 A host of validation studies have shown that the best compromise between computational effort and achievable accuracy is offered by CC theory with single-, double-, and perturbative triple (CCSD(T)) excitations. This method, CCSD(T), has repeatedly been termed the gold standard of quantum chemistry.⁶ It is widely anticipated that this method is, in principle, able to deliver results of chemical accuracy (defined as within 1 kcal/mol of the exact solution) within its domain of applicability.^7,8 The main obstacle that has plagued the widespread application of CCSD(T) in computational quantum chemistry has been the steeply increasing computational effort with system size (O(N⁶) for CCSD and O(N⁷) for CCSD(T)). With such steep scaling, solving the problem with ever increasingly powerful hardware can only be successful to a very limited extent. Hence, approximations that exploit the locality of the dynamic electron correlation are a highly attractive approach to bring CCSD(T) in large scale use for every day computational chemistry applications. The fact that locality exists in the dynamic correlation has been realized early on. In fact, the dispersion interaction with a quickly decaying R⁻⁶ distance dependence is already the longest range surviving dynamic correlation contribution. However, systematically exploiting the locality of the electron correlation has been found to be surprisingly challenging. The origins of the problems in developing low-order scaling approximations to CCSD(T) are manifold. First of all, the complexity of the canonical CCSD(T) method is substantial and the technical complexity of local approximations raises it much higher, thus rendering method development a complex and time-consuming task. Second, and perhaps more important, the accuracy that must be achieved in order to avoid compromising the inherent accuracy of CCSD(T) is very high. Larger molecules have correlation energies of several tens of hartrees, i.e., several thousand kcal/mol. Hence, to achieve an absolute accuracy of about 1 kcal/mol, one needs to solve the CCSD(T) equations to an accuracy of, preferably, 99.99%. If the errors in the local approximations are systematic and smooth, one will profit from error cancellation in calculating chemically relevant energy differences. However, our experience indicates that error cancellation might at best improve the results by about one to two orders of magnitude. Therefore, around 99.9% of the total correlation energy is a practical goal for method development in local correlation theories. One should also mention that for properly constructed local correlation methods, the error is a size extensive quantity and will hence increase with system size. In many chemical processes, however, the changes in a given structure are essentially local in which case error cancellation will keep the error in relative energies reasonably constant with system size.

While it is clear from early results that 98% can be readily achieved,^9–12 an accuracy of 99.9% is much more difficult to reach and involves an exact (non-truncated) treatment of correlation effects that extend over a distance spanning four or more chemical bonds. Hence, it is very difficult to meet the third requirement—an early crossover of the computational effort between the canonical and local CCSD(T) implementations. Fourth, in order to be suitable for applications in computational chemistry in the same routine fashion that density functional theory or canonical CC theory currently are, the local approximations must be robust, reliable, and involve a minimum of heuristic elements.^13,14

Early local correlation theories were developed by Pulay and Saebø by assigning to each localized, occupied Hartree-Fock orbital a set (domain) of unoccupied orbitals composed of “nearby” projected atomic orbitals (PAOs); for each electron pair/triple/etc., the unoccupied space used to generate excited determinants is composed of unions of PAO domains for each constituent orbital.^15–18 Using this prescription, Werner, Schütz, and co-workers have subsequently demonstrated that linear scaling CCSD(T) (LCCSD(T)) on the basis of PAO domains can be realized.^9–12,19 However, the 99.9% accuracy goal is difficult to reach with these methods and requires PAO domains of some 20-40 atoms.²⁰ With at least triple-zeta basis sets, such domains contain more than 500 PAOs each which renders the PAO based LCCSD(T) method very expensive. Alternative local correlation schemes with demonstrated linear or low order scaling include the cluster in molecules approach by Li et al.,^21–23 its extensions by Kállay and co-workers,^24,25 and the incremental method^26–30 and the divide-expand-consolidate method by Jörgensen and co-workers.^31,32

In recent years, we have shown that pair natural orbitals (PNOs^33,34), which were extensively used in the early to mid-1970s,^35–47 form an ideal basis for local correlation treatments.^20,48–53 The tremendous compaction of the virtual space through the PNO expansion allows for the implementation of very accurate, robust, and highly efficient local correlation methods.^20,48–53 If the PNOs themselves are expanded in terms of PAOs, the domain based local pair natural orbital CCSD(T) (DLPNO-CCSD(T)) method emerges.^20,51 We have previously developed a near linear scaling version of this method and demonstrated its application to systems with more than 600 atoms and nearly 10 000 basis functions.^20,51 The method has been carefully benchmarked and found to be accurate and robust.^54–61 Hence, it is an extremely powerful tool for computational chemistry. In order to overcome the more stringent basis set requirements of correlated ab initio theories, several authors, including some of us, have demonstrated the possibility of explicitly correlated PNO calculations.^62–64 An alternative scheme in which orbital specific virtual orbitals (OSVs) are exploited instead of PNOs or PAOs has been developed and an efficient OSV-LCCSD(T) method has been achieved.^65,66

In the original DLPNO-CCSD(T) implementation, a few non-linear scaling steps remained that led to large disk space requirements.²⁰ In particular, the initial integral transformation from atomic orbitals to PNOs via the “density-fitting” (DF) or “resolution of the identity” (RIJ) approximation required up to hundreds of gigabytes of disk storage for the largest calculations. Hence, a fully linear scaling version of this method has remained an important method development target. In order to deal with the complexity of such a development, we have recently revisited the problem of algorithm design for local correlation implementations. In the first paper of the present series, we have reported a sparse map infrastructure that allows one to efficiently and straightforwardly develop algorithms that operate on a multitude of related incomplete index sets, arising from lists of AOs, PAOs, PNOs, auxiliary basis functions or atoms.⁵² Hence, using the sparse map infrastructure, linear scaling algorithms of considerable complexity can be developed efficiently. We have demonstrated this in detail in Part I of this series by developing the linear scaling DLPNO-MP2 algorithm.⁵² The algorithm has a very early crossover with canonical RI-MP2, yielding 99.9% of the correlation energy in a robust and systematic fashion. In the same paper, we have introduced the differential overlap integral (DOI) as a non-heuristic, yet computationally efficient, way to encode ‘multiplicative’ sparsity.⁵² Finally, we have improved the original pair prescreening algorithm that allows for the elimination of electron pairs with negligible contributions to the correlation energy before even the initial guess local MP2 calculation is entered.

In the present work, we have extended our use of the sparse map infrastructure to the much more elaborate case of the DLPNO-CCSD(T) method and have hence arrived at a DLPNO-CCSD(T) implementation that scales linearly in all essential computational steps, the integral transformations, the initial guess amplitude construction, the CCSD amplitude iterations, and the calculation of the triples contribution. The initial orbital localization as well as the preceding self-consistent field (SCF) procedure has not yet been linearized. Furthermore, the pair prescreening is formally quadratically scaling, albeit with a prefactor that is so small that it is irrelevant even for systems with up to a few thousand atoms. Linear scaling techniques for localization and Fock matrix construction are well known in the literature (localization: Refs. 67–69, Coulomb and exchange matrices: Refs. 70–74, overall SCF procedures: Refs. 75–81). The combination of these techniques with the linear scaling DLPNO-CCSD(T) method described in this paper is certainly worthwhile but outside the scope of this work.

In developing the linear scaling DLPNO-CCSD(T) method, utmost care was taken to not compromise on the inherent accuracy of the method. Since the linear scaling implementation will replace the original implementation with negligible differences in the calculated correlation energies, we have refrained from creating another acronym and in the future the term DLPNO-CCSD(T) will refer to its linear scaling implementation.

Coupled cluster theory and its formulation in terms of pair natural orbitals is well documented in the literature.^{20,48,49,51,52} Hence, in the present work, we discuss the relevant algorithmic aspects that lead to a linear scaling implementation of the DLPNO-CCSD(T) method.

The algorithm for linear scaling DLPNO-CCSD(T) is outlined in Scheme 1. Most of the changes to the original DLPNO-CCSD(T) implementation concern the early stages of the calculation in which the lists of significant electron pairs, the correlation domains, and PNOs are constructed and all necessary integral transformations take place. Therefore, this part of the algorithm will be described in some detail, while the remaining steps, most importantly the solution of the coupled cluster amplitude equations and the calculation of the triple excitation correction, are identical to the previously described algorithms.^20,51 As before, the accuracy of the entire DLPNO-CCSD(T) calculation is controlled by three thresholds: T_CutPNO, T_CutPairs, and a third threshold T_CutDO. The thresholds T_CutPNO and T_CutPairs control the length of the PNO expansion and the number of electron pairs that are included in the CCSD iterations, respectively. They are identical to the previous implementation and have been carefully investigated.^52,54 The third threshold, T_CutDO, controls the domain sizes. It has been investigated in detail in Ref. 52 and will be further described below (for additional calibration data see also Figures S1–S3 in the supplementary material⁸²). In order to reach full linear scaling, a number of additional truncations are unavoidable. However, we have followed the strategy to bind the values of these thresholds to the values of the three main thresholds in order to keep the procedure as transparent as possible to the user. In general, sensible defaults are chosen, such that the DLPNO-CCSD(T) method can be used in a black-box fashion.⁵⁴ 

For clarity, Table I summarizes labels used in this manuscript. In this section, we will make use of the results of Part I of this series.⁵² There the formal concept of a sparse map was introduced. In a nutshell, a sparse map is a multimap, a data structure that maps each index in the domain (input) index set to a set of indices in the image (output) index set. For example, a sparse map occupied orbitals to PAOs is denoted as $L ( i → μ ̃ A )$⁠. The subscript “A” indicates that all functions attached to a given atom are included. A sparse map can be implemented in a variety of ways, the simplest of which is to represent it as a list of lists; composition of operations on sparse maps is straightforward. Although sparse map is a concept that arises naturally in describing sparse matrix algebra, sparse maps allow for an easy functional-style composition of algorithms on sparse tensor data as well as avoiding the exponential explosion of metadata in describing the sparse structure of tensors. For more details and operations that can be performed on sparse maps, the reader is referred to Ref. 52.

The first three steps of the algorithm (localization of the occupied orbitals in the reference determinant, construction of the Fock matrix, and determination of PAOs) are identical to the previous implementation. The first new step is the calculation of the DOI integrals. In Part I of this series,⁵² we have shown that the integrals

are excellent substitutes for the Schwartz prescreening integrals $SPI i k ≡ f i g k | f i g k$ (here, f_i, g_k are square integrable one-electron functions, i.e., orbitals or basis functions). The DOI integrals can be calculated efficiently in a linear scaling fashion using numerical integration techniques.⁵² They are used extensively in order to exploit the multiplicative sparsity defined in Part I of this series.⁵² In particular, the DOI integrals between occupied, localized MOs i and PAOs $μ ̃$ are used to construct local orbital domains in step 5 of the algorithm. This replaces the previous heuristic domain definition that was only based on the spatial extent of the occupied orbital i. Unlike previous authors, we normalize the PAOs to unity (see discussion in Ref. 52). The domain definition using the DOI has the key advantage that both the occupied and the unoccupied space are taken into account in the definition of the domains. In this way, it is ensured that all PAOs that have a significant differential overlap with the occupied orbital i are part of the correlation domain for this orbital. Obviously, a significant differential overlap between a given MO and a given PAO is a pre-requisite for significantly large exchange integrals $i μ ̃ | j ν ̃$ and related integrals that dominate the calculation of the correlation energy. The construction of the domains is controlled by the single threshold T_CutDO (note that the much smaller fitting domains are still based on the threshold T_CutMKN as discussed in Ref. 52). The domains have been carefully investigated in Part I of this series and have been found to be of comparable size and accuracy to the heuristically constructed ones.⁵² However, we do believe that the present definition is to be preferred over the heuristic construction since it is more rigorous and does not contain ad hoc elements.

The linear scaling DLPNO-CCSD(T) method uses the same electron pair prescreening that was reported in Part I of this series⁵² and differs from the original procedure in Ref. 20. The refined prescreening procedure is controlled by two thresholds T_CutPre and T_CutDOij that are both bound to the value of the pair cut-off T_CutPairs. It is more conservative and further reduces the already small errors introduced by the original procedure. An alternative procedure that improves upon the accuracy of the LMP2 estimate itself has recently been discussed by Usyvat, Schütz, and co-workers.^83–85 Such techniques are certainly helpful for keeping the number of pairs that have to be treated by CCSD to a minimum. They hence serve a complimentary purpose to the pair prescreening, which aims at identifying the negligible pairs efficiently rather than trying to estimate any pair energy to high accuracy.

Since the refined pair prescreening eliminates fewer electron pairs, and for the sake of robustness, initially will retain electron pairs that will finally be found to be negligible, we also decided to redesign the initial guess semicanonical local MP2 (SC-LMP2^36,37) treatment. Without this refinement, the initial SC-LMP2 calculation could dominate the overall computational effort, which is clearly undesirable. The refined procedure guarantees that the SC-LMP2 correction for the PNO errors does not create errors that are larger than those arising from the coupled cluster part of the calculation, while at the same time, the computational effort of the SC-LMP2 step and the PNO construction remains small (detailed numbers are reported below). The refined procedure does not change the scaling of the overall procedure.

In linear scaling DLPNO-CCSD(T), we initially perform a crude SC-LMP2 calculation in small orbital and fitting domains. The small orbital domains are constructed by considering the sparse maps of the DOIs, scaling the differential overlap cutoff T_CutDO by a factor of 2 (T_{CutDO_Crude}). Similarly, the cutoff T_CutMKN for the small fitting domains is scaled by a factor of 100 (T_{CutMKN_Crude}). These default values ensure that the SC-LMP2 energy recovers more than 99.9% of the value without scaling.

Based on the initial small domains, an initial linear scaling integral transformation⁵² from AO basis three-index integrals $μ ν | K$ to $i μ ̃ | K$ is performed, where μ, ν, κ, τ label atomic orbitals and K, L refer to auxiliary basis functions, i, j, k, l denote localized occupied orbitals, and $μ ̃ , ν ̃$ refers to PAOs.⁵² In order to arrive at a linear scaling integral transformation, it is inevitable to exploit the additive sparsity that is present in the molecular orbitals (⁠$ψ i r = ∑ μ c μ i μ r$⁠) and the PAOs (⁠$μ ̃ r = ∑ μ P ̃ μ μ ̃ μ r$⁠) by screening the expansion coefficients c and $P ̃$⁠. While screening of the MO coefficients has been found to be unproblematic, screening of $P ̃$ must be done carefully in order to avoid numerical instabilities. We refer to Part I of this series⁵² as well as the discussion of Schwilk et al.,⁸³ for further information. Throughout this work, our default consists of cutting coefficients that are smaller than 10⁻⁴ in absolute value (cutoff T_CutC). This, in our experience, leads to numerically very stable results that have minor impact (typically below 1 μE_h) on the computed energies. However, one should realize that with such conservative thresholds, significant savings due to this approximation only occur for large or extended one-dimensional systems. In fact, in their PNO-LMP2 method, Schwilk et al. have chosen to not cut the PAO vectors by default.⁸³ 

Following the local integral transformation, the local exchange operators $K ij μ ̃ , ν ̃$ are constructed and used to calculate SC-LMP2 amplitudes (defined as the uniterated amplitudes $T ij μ ̃ , ν ̃ = K ij μ ̃ , ν ̃ / ε μ ̃ + ε ν ̃ − f ii − f jj$⁠), where $ε μ ̃$ and $ε ν ̃$ are energies of quasi-canonical, non-redundant PAOs as described in Refs. 20, 48, 49, and 52 and f_ii and f_jj are diagonal elements of the Fock matrix in terms of localized internal orbitals. Using T^ij and K^ij, one computes an estimate of the pair correlation energy ε_ij. These initial pair energies are used to divide the linear scaling number of electron pairs that survived the prescreening into three classes: (1) electron pairs with a pair correlation energy >T_CutPairs that will enter the coupled cluster iterations, (2) electron pairs with a pair correlation energy <T_CutPairs and >T_{CutPairs_MP2}. These electron pairs will be treated together with the pairs of type (1) in a second, more accurate SC-LMP2 step to be described below. (3) Pairs with a pair correlation energy <T_{CutPairs_MP2} will not be considered further. However, since the crude SC-LMP2 energies are still more accurate than the pair energy estimates provided by the prescreening procedure, the energies of these pairs will be summed up and together with the estimated sum of the pairs rejected by the prescreening procedure will provide an estimate for the error introduced by the local approximations. It should be pointed out that with the chosen default thresholds, these corrections are very small. For example, in our largest example, C₃₀₀H₆₀₂, where 97.3% of all pairs are skipped during prescreening or crude MP2, these corrections add up to only 0.05% of the correlation energy.

After the crude SC-LMP2 step, it is determined which electron pairs are treated by CCSD (strong pairs) and which will remain at the SC-LMP2 level (weak pairs). Note that the computation of the triples correction will also make use of these pairs, in the same way as described in Ref. 51. For the accurate SC-LMP2 calculation, the final, large orbital, and fitting domains are constructed using the default values T_CutDO = 10⁻² and T_CutMKN = 10⁻³. After this step is completed, all sparse maps can be determined and the second, more accurate linear scaling integral transformation is performed using the MO and PAO truncations as outlined before.

It should be noted that our program has the option to provide fully iterated LMP2 amplitudes (including the coupling terms from the off-diagonal Fock matrix elements) in place of the SC-LMP2 amplitudes described above. The question whether the more elaborate fully iterated LMP2 offers advantages over the initially chosen SC-LMP2 approach was studied in some detail in Ref. 54. The benchmark data assembled there indicate that the accuracy gain from fully iterated amplitudes is very modest such that the increased computational effort cannot be justified for default truncation thresholds. Only for the most accurate PNO calculations (referred to as “TightPNO” in Ref. 54), iterated LMP2 amplitudes are used as a default in our program.

From the LMP2 or SC-LMP2 guess amplitudes, approximate pair densities $D ij = T ij T ̃ i j + + T ij+ T ̃ ij$ are constructed as in the original DLPNO-CCSD(T) method,^20,51 followed by the PNO construction and the estimation of the PNO errors in the same way as described in Refs. 20 and 51.

The remainder of the calculation proceeds as in the original DLPNO-CCSD(T) method. A technical complication concerns the fact that previously complete integral lists $i μ ̃ | K ̄$ were available, while the linear scaling integral transformation only provides a linear scaling list of integrals that must be further processed by the integral transformations to create integrals over PNOs. These transformations hence require the matching of several sparse index sets. They are greatly facilitated by the sparse map infrastructure described in Part I of this series.⁵² 

The only other significant change to the previous DLPNO-CCSD(T) algorithm concerns the calculation of the singles-singles interaction contribution to the CCSD residual. It is well known that these contributions can be expressed in terms of a “singles Fock matrix” $G t 1$ that can be computed in AO basis using a variety of approximations.⁵³ If done this way, this step is of the same level of complexity as the calculation of a Fock matrix during the SCF procedure. It could most certainly be made linear scaling by employing the techniques for linear scaling SCF procedures cited in the introduction. However, the prefactor for this construction is high. In our experience, this step can become computationally expensive for calculations on large molecules with large and accurate basis sets. In the new implementation, the necessary integrals with one- and three-external labels, respectively, are calculated in the PNO basis using the same MO/PAO integrals that enter the remaining calculation. These integrals are generated in a linear scaling fashion. This leads to large speedups and linear scaling for the entire sigma-vector construction. For example, for C₂₀₀H₄₀₂, the largest alkane chain that we compute with the nonlinear code, a speedup of 34 is observed with the linear version (620 min vs. 18 min) for the terms involving the singles Fock matrix.

A number of calculations have been performed to evaluate the effect of the new domain construction and to define sensible default values for the newly introduced thresholds. The results are summarized in the supplementary material (Tables S2 and S3).⁸² All truncation parameters had their default values defined in Ref. 54 as well as above (for an overview see Table S1⁸²). The only change concerns the default value for cutting off the triples natural orbitals.⁵¹ Given that the computation of the triples correction is computationally considerably less expensive than the PNO integral transformation and the solution of the CCSD equations, a tight threshold T_CutTNO = 10⁻⁹ is affordable and further improves the fraction of recovered correlation energy to typically >99.9%.

With these default values, the DLPNO-CCSD(T) correlation energies of a number of medium-sized molecules and associated reactions were computed and compared to the semicanonical results. Table S2⁸² shows that the average recovered correlation energy is 99.91%, compared to 99.82% using the nonlinear methodology with the old domain construction, thus decreasing the error in correlation energy by a factor of two. Due to the increase in accuracy also the prediction of relative energies has improved. The mean absolute deviation for reaction energies from the canonical reference energies (see Table S3⁸²) is only 0.36 kcal/mol compared to 0.53 kcal/mol with the old methodology, while the maximum deviation has improved from −2.15 kcal/mol to 1.81 kcal/mol.

The scaling behavior was investigated on a set of linear alkane chains. Figure 1 shows a comparison of the timings for the old nonlinear DLPNO-CCSD(T) with the new linear method. In order to put the results into perspective, the time needed to compute the reference determinant using the efficient RIJCOSX (chain-of-spheres) approximation^86,87 is also shown. It is evident from the results that the computational cost of the new implementation scales linearly with respect to system size and improves greatly over the previous DLPNO-CCSD(T) implementation that already allowed for calculations of unprecedented size. The new linear methodology is almost a factor of 9 faster than the nonlinear methodology for C₂₀₀H₄₀₂, the largest alkane chain that we computed with the nonlinear code. The effective exponent obtained by finite difference in the range from 100 to 900 atoms is 1.20. Furthermore, the time taken by the coupled cluster step is never much larger than the time taken for the SCF calculation. In fact, owing to the linear scaling of the coupled cluster calculation versus quadratic to cubic scaling for the SCF procedure, the SCF calculation will eventually become rate limiting. At this point, linear scaling SCF techniques become effective in order to avoid bottlenecks arising from the SCF procedure itself. However, for the more commonly met three-dimensional molecules, the point where linear scaling SCF procedures become a necessity occurs much later. In fact, for most of today’s applications with systems of up to 200 atoms it is, in our experience, more important to improve the prefactor of the SCF procedure rather than the scaling. Nevertheless, as mentioned in the introduction, the combination of linear scaling DLPNO-CCSD(T) with linear scaling SCF procedures is certainly a worthwhile goal for future developments.

For linear alkane chains, the crossover between linear scaling DLPNO-CCSD(T) and quadratic to cubic scaling SCF appears to occur around 350 atoms. Based on the wall-clock times reported in Figure 1, it is clear that coupled cluster calculations with more than 10 000 basis functions (triple-zeta basis set) and 1000 atoms can be performed in a routine fashion with the new DLPNO-CCSD(T) code.

The contributions for the different parts of the algorithm were investigated in more detail. Figure 2 shows detailed timings for all major contributions (that sum up to >99.8% of the entire DLPNO computing time) on a per atom basis. Linear scaling is achieved when the computing time per atom approaches an asymptotic value. The most expensive part of the DLPNO-CCSD(T) calculation for large systems, the 3-index integral transformation, reaches an asymptotic value of about 40 s per atom at around 500 atoms. For all other linear scaling parts of the algorithm, the asymptotic value is reached much earlier at about 100 atoms. This late onset of linear scaling for the integral transformation is mainly caused by the long delocalization tails of the PAOs. This subject clearly warrants further investigation.

Only the construction of a few matrices in the global PAO space (e.g., overlap and Fock matrix), the localization of occupied MOs, and the calculation of the differential overlap integrals show nonlinear scaling. For linear carbon chains with more than 1000 atoms, these will become significant and will deteriorate the scaling behavior. However, for more realistic 3D systems like crambin (644 atoms), these nonlinear scaling routines take less than 4% of the entire computation time (vide infra) and can thus be regarded as negligible for up to 1000 atoms.

As discussed above, the 3-index integral transformation exploits the additive sparsity that is present in the molecular orbitals and in the PAOs in order to achieve linear scaling.⁵² This means that not all integrals are calculated, i.e., some are screened. Figure 3 shows how this screening efficiency changes with the molecular size and the number of computed shell sets per atom depending on the molecular size.

In the nonlinear DLPNO code, the integral transformation producing the 2-external 3-index integral integrals was already linearized in terms of integral storage, i.e., the disk space needed to store the 3-index integrals increased linearly with system size. However, by increasing the molecular sizes beyond a few hundred atoms, we quickly hit the storage wall for the 1-external integrals. With the present implementation, all 3-index integrals have linear disk space requirements. Figure 4 shows the disk space needed to store the 1-external and 2-internal 3-index integrals. For both types of integrals, a linear increase with system size is observed. On the other hand, the nonlinear code needs almost 1 TB for the 1-external 3-index integrals with 600 atoms, and almost 3 TB for 900 atoms.

The previously reported calculation on crambin together with the def2-SV(P) basis set represented the first CCSD(T) level calculation on an entire protein. While this calculation took several weeks on one processor core and consumed 1.3 TB of disk space, the present code finishes this calculation in less than 4 days elapsed time and consumes only 0.3 TB of disk space. It is even efficient enough to complete an improved calculation with the def2-TZVP basis set (12 075 basis functions), which can be performed much more efficiently (two weeks on 4 processor cores and using less than 700 GB of disk space) with the present code, as is summarized in Table II. We even performed a DLPNO-CCSD(T) calculation with more than 20 000 basis functions and more than 1000 atoms (C₃₅₀H₉₀₂ def2-TZVPP). The contributions from the different parts of the algorithm are similar to those of the smaller carbon chain calculations, indicating that no further major bottlenecks arise for this huge calculation.

In the present paper, a linear scaling DLPNO-CCSD(T) method was developed. It is based on a substantial redesign of the previous DLPNO-CCSD(T) implementation of Refs. 20 and 51 with the help of the sparse map infrastructure of Ref. 52. The latter programming infrastructure allows for an efficient implementation of complicated operations on sparse index sets. Furthermore, several algorithmic refinements have been incorporated in the linear scaling DLPNO-CCSD(T) procedure: (a) the use of the differential overlap integrals for non-heuristic construction of domains (exploitation of multiplicative sparsity in general terms), (b) the refined prescreening procedure of Ref. 52, (c) a two-level local (and optionally semicanonical) MP2 procedure that ensures an efficient and accurate initial guess and PNO construction, (d) linear scaling integral transformation described in Ref. 52 has been used throughout (this requires additional thresholds for the truncation of MO and PAO vectors that have been carefully calibrated), and (e) the direct construction of the singles-singles coupling terms from transformed integrals rather than from an AO basis Fock-like matrix.

The algorithmic improvements lead to a method that is linear scaling in all major computational steps. Only steps that individually account for less than 2% of the computational time even for calculations of very large three-dimensional systems remain non-linear scaling (mostly the orbital localization, the grid setup for the calculation of the differential overlap integrals and formation of the initial Fock matrix). Every linear scaling implementation has such steps and hence we feel that it is justified to refer to our implementation as a linear scaling DLPNO-CCSD(T) method. Should the need arise to linearize the mentioned steps, this could be readily accomplished.

The numerical tests show that for large three-dimensional systems, the new method is up to 7 times faster and consumes up to 4 times less disk space than the previous implementation. The linear scaling DLPNO-CCSD(T) method does not sacrifice any accuracy relative to the original implementation. If anything, the algorithmic improvements described herein improve the already high accuracy of the method.

All available benchmark and chemical application studies show that DLPNO-CCSD(T) is a reliable, robust, and accurate black box method for computational chemistry.^54–61 Thus, the present work represents an important step forward towards the availability of large-scale coupled cluster calculations for routine applications in computational chemistry. In fact, the computational requirements for DLPNO-MP2 and DLPNO-CCSD(T) are not grossly different. Hence, for any system for which a DLPNO-MP2 calculation can be performed, DLPNO-CCSD(T) should also be possible. Clearly, given the superior accuracy of coupled cluster theory, the latter is the preferred choice. The implementation will be made publically available in the next release of the ORCA⁵⁸ quantum chemistry package.

F.N., P.P., and C.R. gratefully acknowledge financial support of this work by the Max Planck Society, the DFG SPP 1601, and the cluster of excellence (RESOLV, University of Bochum). P.P. gratefully acknowledges funding from Fonds der Chemischen Industrie. E.F.V. acknowledges support by the U.S. National Science Foundation (Grant Nos. CHE-1362655 and ACI-1047696) and by Camille and Henry Dreyfus Foundation.