Pursuing the basis set limit of CCSD(T) non-covalent interaction energies for medium-sized complexes: case study on the S66 compilation

Our recently developed conventional, explicitly correlated, and local natural orbital (LNO) based coupled cluster with single, double, and perturbative triple excitations [CCSD(T)] methods are utilised to closely approach the basis set limit of interaction energies between medium-sized molecular complexes. The study is partly motivated by our previous finding that state-of-the-art diffusion Monte Carlo and CCSD(T) interaction energies are not in complete agreement even for some of the dimers in the popular S66 compilation [Al-Hamdani and Nagy et al. Nat. Commun. (2021) 12:3927]. We improve upon the previous S66 CCSD(T) reference by using up to the quadruple-ζ (QZ) level CCSD(T) contributions, utilising our recent (T+) approach and suggesting a more robust explicitly correlated CCSD reference. Our efficient LNO-CCSD(T) method allowed us to also compute tightly converged, quintuple-ζ level interaction energies for the S66 set. Moreover, the new CCSD(T)/QZ references enabled the so far most rigorous benchmark of the LNO and other popular local correlation approaches on the S66 set. The accuracy of local approximations was found in the few 10 to few 100 cal/mol interaction energy error range, assuming sufficiently tight settings, such as the Tight and the Normal settings of LNO-CCSD(T) threshold combinations, respectively. GRAPHICAL ABSTRACT


Introduction
Non-covalent intermolecular interactions are ubiquitous and play a key role in a large number of (bio)chemical systems. The precise and efficient modelling of intermolecular interactions is still one of the great challenges CONTACT  of computational modelling since accurate description is required for a large number of individually small but cumulatively important interaction components. Since electron correlation effects have to be taken into account due to the required sub-kcal/mol accuracy, a number of wave function based methods have been made available for intermolecular interaction computations [1][2][3][4]. Among others, the application of the symmetryadapted perturbation theory (SAPT) [5][6][7], diffusion Monte Carlo (DMC) [8][9][10][11], random phase approximation (RPA) [12,13], and coupled-cluster (CC) methods [14], in particular the gold standard CC model with single and double excitations (CCSD), augmented with perturbative triples correction [CCSD(T)] [15] have been particularly successful in this context.
For the performance assessment and comparison of the various computational methods, a number of test compilations have been suggested including molecular complexes from a few atoms to extensive supramolecular or biomolecular examples representative of all major types of non-covalent interactions [16][17][18][19][20][21][22][23]. Among those the S66 list of Řezáč, Riley, and Hobza [17] has become one of the standard benchmark sets. The S66 compilation contains 66 dimers composed from 14 monomer molecules and represents a wide range of interaction patterns occurring in organic and biomolecules, such as hydrogen bonding or dispersion interactions between both aliphatic and aromatic compounds. The largest dimers of S66 reach up to 34 atoms and 116 electrons, which pose a considerable difficulty for computationally demanding methods, such as conventional CCSD(T). While the previous CCSD(T) benchmarks of S66 appear to be well-established [17,[24][25][26], we found recently that highly-converged CCSD(T) and DMC interaction energies (IEs) are not completely consistent for a few π-π stacking complexes among the most challenging S66 dimers [27]. The deviation of at least 0.1 kcal/mol between the CCSD(T) and DMC IEs was outside of the methods' uncertainty estimates. This finding provided an additional motivation to improve on the basis set convergence of the previous S66 IEs, to explore and eliminate basis set incompleteness as a potential source of deviation between DMC and CCSD(T) IEs.
Since the evaluation of conventional CCSD(T) energies scale with the seventh power of the system size, the largest molecules accessible with sufficiently large atomic orbital (AO) basis sets are also in the range of 20-40 atoms, even with recent massively parallel implementations [28][29][30][31][32][33][34][35][36][37][38]. Though the use of convergent basis set hierarchies [39], complete basis set (CBS) extrapolation techniques [40], and counterpoise (CP) corrections [41,42] is invaluable, basis sets including diffuse and high angular momentum orbitals are still required to sufficiently approach the CBS limit of the CC IEs. Therefore, the original S66 IEs reference [17] and its revisited variants [24][25][26] had to rely on composite or focal point type approaches in combination with second-order Møller-Plesset (MP2) based basis set corrections [43,44]. However, the well-known challenges of MP2, e.g. in accurately describing aromatic interactions provide a reason for minimising the role of MP2-based corrections when aiming at improved CCSD(T)/CBS estimates.
To that end, the basis set used for the CCSD(T) part, and thus the computational effort, has to be increased. To accelerate such extensive CCSD(T) computations, we combined the t1-transformed CCSD expressions and the density-fitting (DF) approximation [45] with our integral-direct and efficiently parallelised CCSD and (T) algorithms into a highly-optimised, parallel, and dataeconomic CCSD(T) code, scalable up to over 30 atoms representing the largest quadruple-ζ level CCSD(T) computation ever presented [37]. Combined with the rapid development in computational hardware, our new CCSD(T) code enabled us here to update the original [17] double-ζ and the latest [26] triple-ζ CCSD(T) components of the S66 IEs, presented 11 and 5 years ago, respectively, by quadruple-ζ level CCSD(T) contributions.
Since an affordable explicitly correlated treatment of the (T) corrections is not yet available, in practise MP2-based scaling approaches are employed to the (T) terms [52,62]. Most recently, we showed that the first, not size-consistent variant of Knizia et al. [52], that is, the (T * ) approach, can be effectively replaced by a similar, but size-consistent scheme termed (T+) [62]. While the (T * ) and (T+) results are often similar when the basis set convergence patterns match for the dimer and the monomers, (T * ) can exhibit considerable sizeconsistency errors when the monomers of a dimer are markedly different. We showed that this issue is effectively resolved by the (T+) alternative [62], and thus we will employ here the resulting CCSD(F12 * )(T+) method for the S66 dimers too.
To further increase the efficiency of CC methods, a number of successful local correlation approaches have been introduced exploiting the relatively quick distance decay or local nature of electron correlation . Currently, strategies that combine these local approximations with natural orbital [93][94][95][96] (NO)-based compressions stand out [87,89,91,97]. Here, we will utilise and benchmark some of these state-of-the-art local CCSD(T) approaches, including the pair NO (PNO) based method of Hättig, Tew, and Schmitz [70,97,98], the PNO-L method of Werner and Ma [89,99], the domainbased local PNO (DLPNO) method of Neese, Valeev, and co-workers [87,100,101], as well as the local natural orbital (LNO) approach developed by Nagy and Kállay [36,91].
The efficiency of these local CCSD(T) methods has already been utilised to investigate the non-covalent interactions in large molecular systems beyond the capabilities of conventional CCSD(T) [4,20,27,87,89,91,[102][103][104][105][106][107][108][109][110][111]. For instance, we reported some of the largest local CCSD(T) computations for the up to 101atom complexes of the L7 compilation [18] and a 132atom supramolecular fullerene dimer using very tightly converged LNO-CCSD(T) with augmented quintupleζ basis sets [27]. Such challenging complexes with extended π -π interactions limit the applicability range of all local CCSD(T) codes aiming at robust accuracy. Systems with more rapidly decaying long-range interactions and smaller volumetric density of electrons, such as biomolecules can be treated more effectively. We illustrated this recently by computing the quadruple-ζ level LNO-CCSD(T) interaction energy of a 1023-atom protein-ligand complex [91].
Clearly, these large-scale applications are far from the range of the at most 20-30-atom molecules and triple-ζ basis sets, where the accuracy of various local CCSD(T) approaches are usually documented. Therefore, it is important to show that local approximations have well-controlled and homogeneous accuracy for a wide range of molecules, across multiple large basis set choices, CBS extrapolation, and CP corrections, etc. Thus, we will also utilise the new, conventional CCSD(T) results to assess the LNO and other, PNO-based local correlation schemes in a complexity range that was not accessible before. Currently, more and more of such benchmarks would be needed considering the increasing popularity and wide applicability domain of local CCSD(T) approaches. Moreover, comparisons treating different local CCSD(T) methods on an equal footing [27,91,110,[112][113][114][115] remain scarce. Additionally, by utilising the unique capabilities of our LNO-CCSD(T) implementation [36,91], we approach the CBS limit even closer by reporting very well converged, quintuple-ζ level LNO-CCSD(T) IEs for the first time for the entire S66 compilation.

Methods: algorithmic considerations and development
A brief summary is given of our previous and present algorithmic developments and programming efforts regarding the relevant CCSD(T) methods allowing for a more detailed analysis of the current capabilities and limitations.

Conventional and explicitly correlated CCSD(T)
In conventional or explicitly correlated CCSD(T), both the computational time and the data storage bottlenecks are often relevant since these scale with the seventhand fourth power of the number of correlated occupied (n o ) and virtual (n v ) molecular orbitals (MOs), respectively. For a single CCSD iteration and the (T) parts, respectively, the rate determining steps require about 2(n 2 o n 4 v /4 + 4n 3 o n 3 v ) and 2n 3 o n 3 v (n o + n v ) floating point operations (FLOP). For instance, an extensive computation with n o = 30 and n v = 1000, the FLOP requirement of 12 CCSD iterations and the (T) step amounts to about 96 and 55.6 petaFLOP (10 15 FLOP), respectively. It has become possible to meet this high demand with the latest parallel CCSD(T) codes exploiting the appearance of many-core CPUs with teraFLOPS, i.e. 10 12 FLOP per second capabilities. However, conventional CCSD(T) applications reaching the 1000-orbital range are still scarce [28,[30][31][32][33][34][35][36][37].
Using our new integral-direct and highly optimised CCSD(T) code and 224 CPU cores for 3 days, we reported the largest CCSD(T) computation so far at the quadrupleζ level involving 31 atoms and 1569 orbitals [37]. Consequently, the outstanding peak-performance utilisation and parallel-efficiency of our CCSD(T) code made the next step accessible in the basis set hierarchy for the entire S66 compilation. Namely, using the highly-reliable haug-cc-pVQZ (heavy-augmented correlation consistent valence quadruple-ζ ) basis set, the largest S66 dimers containing 36 atoms and 1520 orbitals are now within our reach.
When using such a large number of orbitals, we also have to consider the data storage bottlenecks of CCSD(T). Many CCSD and (T) implementations require the storage of at least a part of the four-external and threeexternal four-centre electron repulsion integral (ERI) tensors, respectively. With the above dimensions of n o = 30 and n v = 1000, the corresponding arrays alone would require 1863 and 224 GB storage space, which is far from practical without an integral direct implementation. Consequently, various algorithms have been implemented to deal with these data size bottlenecks [28, 33-35, 37, 45, 116]. We introduced the first completely integral-direct approach eliminating the storage of all four-centre ERI tensors with at least two-external indices. Both of our CCSD and (T) algorithms affordably assemble these four-centre ERIs when needed from threecentre ERIs exploiting DF and t1-transformed expressions, making the memory and disk input/output (I/O) demand affordable even well above 1000 orbitals [37].
Considering the slow convergence of the CC correlation energies with the AO basis set, the use of F12based CC methods is highly advantageous, at least at the CCSD level. Multiple variants of explicitly correlated CCSD approaches [52,55,117,118] exhibit improved accuracy per cost performance compared to conventional CCSD [57]. Recently, we also combined the CCSD(F12 * ) formulation of Hättig and co-workers [55,117] with our efficient DF-CCSD(T) algorithm discussed above [62], although its data-storage and parallel scaling capabilities are not yet fully optimised. Nevertheless, dozens of CCSD(F12 * )(T+) computations presented here with the haug-cc-pVTZ-F12 basis [119] exceeded the 1000 AO mark, illustrating that the implementation is already highly capable. Since the latest DF-and parallelisationbased advancements are not yet available for F12-based CCSD codes, except for the CCSD(2) F12 implementation of Valeev et al. [33], our S66 computations with 1184 AOs are close to the capability limits of current F12-based CCSD codes without additional approximations.
Unfortunately, current F12-based CC methods, in practise, cannot directly accelerate the basis set convergence of the (T) contribution. A frequently employed heuristic approach called (T * ) scales the conventional, that is, F12 contribution-free (T) correction, (E (T) ) by the ratio of the MP2-F12 and MP2 correlation energies (E MP2−F12 and E MP2 , respectively) [52] as which is, however, not size-consistent. Inspired by the practicality of (T * ), we recently proposed a similar scaling approach denoted by (T+), which is, however, sizeconsistent [62]. In order to introduce the (T+) approach, let us rewrite the CCSD(T) correlation energy as where one of the occupied summation indices, i was separated to define orbital-specific correlation energy contributions.
is size-consistent and also improves upon the basis set convergence of the (T) energy contribution [62]. Sometimes a size consistent (sc) variants of (T * ), that is (T * sc ) and (Tc sc ), are also utilised, where the E MP2−F12 E MP2 or the analogous E CCSD(F12 * ) E CCSD scaling factors obtained for the dimer are also used for all monomers [25,26].

Local natural orbital CC methods
It is clear from the above analysis that it is very challenging to extend the limits of the seventh-power-scaling CCSD(T) with parallelisation and F12 methods alone. Therefore, the aim of the numerous approximate local correlation approaches is to accelerate CC methods while retaining their reliability. It is common to consider orbital or orbital pair specific correlation energy contributions, such as in Equation (2), in order to introduce orbital (pair) specific approximations [36,71,78,81,87,89,91,120,121].
Let us briefly summarise the main ideas of the LNO family of local correlation approaches [36,[90][91][92][122][123][124][125][126] in order to provide some background for the introduction of a methodological development proposed here. All of our local correlation methods, starting from the MP2 level through CCSD(T) and CCSDT(Q) up to general order CC, employ orbital specific approximations utilising the following energy expression: Here, I and J denote localised MOs (LMOs) and all LMO pairs are assigned to either the strongly correlated or the distant orbital pair group. The individually small distant pair contributions, δE MP2 IJ , are accounted for by an efficient, approximate MP2-level expression [124]. Then domain approximations are employed for the dominant, strong pair contributions. For each LMO I, an orbital domain, that is a set of atomic orbitals, strong pair LMOs, domain specific virtual MOs, etc., is defined to evaluate their correlation energy contribution at the MP2 level [36,91]. In a subsequent step, domain-or LMO-specific LNOs are constructed to compress the orbital spaces of the domain and to evaluate the δE LNO−M I correlation energy contribution of the domain at the M = CCSD(T), CCSDT(Q), etc. level of theory. Finally, the LNO compression approximation is corrected at the MP2 level in each domain via a E MP2 As discussed before in detail [36,91,124,125], this procedure involves a number of different approximations, for which we introduced a hierarchical set of composite thresholds: Loose, Normal, Tight, very Tight, etc. [91]. Note that the composite thresholds combine a set of settings which individually form a systematically convergent series for each approximation. Although the systematic convergence property is not necessarily inherited by the composite thresholds in all cases, in practise, they usually still provide a monotonically convergent energy series which is also suitable for extrapolation toward the local approximation free correlation energy [91]. To that end, we proposed to use a rather cautious extrapolation expression assuming only that monotonic convergence occurs in the threshold series [91]. Namely, an energy estimate is formed from the two tightest available results supposing that the subsequent step in the series will be smaller than the last. For instance, the extrapolation using Normal and Tight settings will give a Normal-Tight (N-T) result of where an LNO error estimate is also defined via the ±(E Tight − E Normal )/2 term [91]. We note that the expression in Equation (5) can be identically reformulated as E N−T = E Normal + 1.5(E Tight − E Normal ), which formula was recently adopted by Bistoni et al. [101] to extrapolate DLPNO-CCSD(T) results in a very similar manner.
We have successfully employed this strategy to obtain systematically improved LNO-CCSD(T) interaction energies for a number of studies. For instance, complexes of alkali and alkaline earth metal ions representative of their non-covalent interaction in protein environments were investigated up to the very Tight convergence level in LNO-CCSD(T) combined with augmented quintuple-ζ basis sets [107][108][109]. Reaching this level of accuracy with a few tenths of a kcal/mol uncertainties compared to CCSD(T)/CBS in both the local correlation and the basis set error estimates has become relatively affordable for systems including 50-100 atoms. Our recent algorithmic developments also made well converged local CCSD(T)/CBS estimates accessible up to more than 100 atoms even for challenging complexes with extended, delocalised π-electron systems [27,91]. Although with a considerable computational effort, even very Tight and aug-cc-pV5Z level LNO-CCSD(T) computations were feasible for molecular complexes of up to 130 atoms, including the L7 set of complexes and an extensive fullerene complex [18,27]. The challenging complexes of that study highlight the technical difficulties occurring with unusually extended, delocalised π -systems. Namely, their slowly decaying long-range interactions lead to a large number of important orbital pair interactions and thus to a considerable decrease in the efficiency of any local correlation approach [27,89,106]. Finally, in the cases free of large π -systems, such as protein-ligand complexes, our LNO-CCSD(T) implementation scales up to 1000 atoms [91].
This extensive exploration of intermolecular interactions revealed that the previous LNO algorithms did not provide completely balanced error measures for systems with and without ghost AOs used for the counterpoise (CP) correction of the basis set superposition error (BSSE). As expected, the BSSE is even more problematic for large molecular complexes with extensive diffuse AO basis sets than for the well-documented case of smallto medium sized complexes. The source of this unbalanced accuracy is in the construction of the LMO-specific orbital domains. Here, first those virtual MOs are gathered into the domains which are expected to be important for the correlation of the domain's LMOs. These virtual MOs reside on atoms whose AOs contribute considerably to the expansion of the domain's LMO. As pointed out in our previous work [36], a drawback of this initial virtual MO list surfaces if there are ghost atoms in the system holding AOs but not nuclei or electrons. Then, there could be some ghost atoms where the expansion coefficients of all LMOs are especially small and the virtual MOs residing on these ghost atoms are not added to the initial virtual MO list of the LMO domain.
To overcome this, we extended the initial virtual MO definition via a second virtual MO selection criterion [36]. In brief, we identify potentially important ghost atoms not only by the size of the LMO coefficients but also using a three-centre AO ERI estimate. First, a fixed threshold was determined for identifying the additional set of domain atoms and the corresponding orbitals, which was suitable for Normal threshold settings and up to medium-sized molecular complexes. In the present work, we update this selection by making it LNO threshold and target accuracy dependent as well as somewhat tighter than before. In this way, the accuracy is increased for large complexes with many ghost atoms and long-range interactions with a moderate additional computational effort. For the sake of reproducibility, the updated algorithm extending the one explained in Ref. [36] employs the (1-T EDo )10 4 E h threshold for the screening integrals. Thus, this threshold is now tied to the LMO completeness threshold T EDo . For instance, the new threshold amounts to 1, 0.5, and 0.1 E h with Normal, Tight, and very Tight settings, instead of the previous value of 2 E h used irrespective of the target accuracy. Let us note that this update is turned on by default in the latest, 2022 release of Mrcc for the default localcc = 2021 set of options [127,128]. Therefore, here we provide updated LNO-CCSD(T) benchmarks with these settings in Section 4.3.

Computational details
All new DF-based conventional [37], explicitly correlated [62], and LNO-based local [36,91] CCSD(T) computations have been performed using the 2022 release of the Mrcc quantum chemistry programme suite [127,128]. For the explicitly correlated CCSD(T) computations the CCSD(F12 * ) approach of Hättig and co-workers [55] is combined with our (T+) scaled-triples corrections [62]. The core electrons were frozen in all correlation calculations.
Correlation consistent basis sets augmented with diffuse functions only for the non-hydrogen, i.e. heavy atoms, that is, haug-cc-pVXZ [39,129] (X = D, T, Q, 5) were employed, which are abbreviated throughout the manuscript as haXZ. As the AO basis set of the CCSD(F12 * )(T+) computations, the correlation consistent X-tuple-ζ (haug-)cc-pVXZ-F12 (X = D, T, Q), or haXZ-F12 in brief, basis sets developed for explicitly correlated calculations by Peterson et al. [130] were employed together with the corresponding cc-pVXZ-F12-OPTRI complementary auxiliary basis (CABS) bases of Yousaf and Peterson [131,132]. The diffuse functions extending the original cc-pVXZ-F12 bases were adopted from the work of Martin and co-workers [119]. Following the recommendation of Shaw and Hill, their cc-pVDZ-F12-OPTRI+ basis set was employed instead of the original 'OPTRI' variant at the double-ζ level [133]. All Hartree-Fock (HF) IEs corresponding to the CCSD(F12 * )(T+) results contain the CABS correction (HF+CABS).
Results utilising (ha)XZ(-F12) and (ha)(X + 1)Z(-F12) basis sets as well as complete basis set (CBS) extrapolation will be denoted by (ha)(X,X + 1)Z(-F12). For the extrapolation of HF energies, the two-point formula suggested by Karton and Martin is used with the recommended parameters [137,138]. Conventional, LNO-, and DLPNO-CCSD(T) correlation energies were extrapolated employing the formula of Helgaker and co-workers with the exponent of 3 [40]. CCSD(F12 * ) correlation energies using the (ha)XZ-F12 and (ha)(X + 1)Z-F12 basis sets are extrapolated via the formula of Ref. [40] with the exponents of 3.1458 for X = D, as well as 4.3448 and 4.6324 for X = T without and with the diffuse AOs, respectively, following the recommendation of Martin and co-workers [26].
To characterise the performance of the methods, the mean signed error (MSE), the mean absolute error (MAE), the root mean square deviation (RMSD), the standard deviation (STD), and the maximum absolute (MAX) error of the computed results will be applied.
Regarding the computational cost, obtaining all DF-CCSD(T)/haXZ (X = T, Q) and the CCSD(F12 * )(T+)/ haXZ-F12 (X = D, T) references for the entire S66 set, including CP corrections, required the execution of about 30 exaFLOP that is, 3·10 19 operations. Using a current many-core CPU of about 1 teraFLOP performance, this amounts to about 350,000 core hours or about a node year effort, which is rather affordable compared to the highperformance computing grants and resources available today.

Conventional and F12-based CCSD(T): convergence to the CBS limit
Let us start with an analysis of one characteristic example, illustrating the convergence of the CCSD(T) correlation energies very close to the CBS limit, within cca. 1% or 0.01 kcal/mol remaining basis set uncertainty. Here, we select the dimer No. 20 of S66, the acetic acid dimer, which is ideal for this illustration for a number of reasons.
The dimer contains 16 atoms, which is around the size limit accessible with basis sets one cardinal number higher than our reference computations. That is, highly demanding but very well converged CCSD(T)/ha5Z and CCSD(F12 * )(T+)/haQZ-F12 computations are within the limit of feasibility (1456 and 1328 AOs, respectively). Additionally, the (AcOH) 2 IE exhibits one of the slowest convergence to the CBS limit of CCSD(T) and exhibits a double H-bond motif that may be challenging for certain DMC algorithms [11].
The most common approach to analyse the performance of CCSD(T)/XZ and CCSD(F12 * )/XZ-F12 is to compare their convergence with respect to the cardinal number X. Here, we also take into account the computational cost in the analysis by considering the about one quarter to one third higher demand of CCSD(F12 * ) with respect to CCSD with the same basis set, and more importantly, the marked size difference of the (ha)XZ and (ha)XZ-F12 basis sets even with the same X value. In other words, we are interested in the accuracy achieved compared to the invested number of operations or computer time taking into consideration that CCSD(F12 * )(T)/XZ-F12 computations are more demanding than CCSD(T)/XZ with the same X. To that end, Figures 1 and 2 plot the CCSD and (T) correlation energy contributions, respectively, to the IE of (AcOH) 2 as the function of the operation count requirement. For the sake of simplicity, we transform the operation counts to compatational time estimates assuming a 2 teraFLOPS, that is 10 12 FLOP per second, (cca. 40 core) CPU and 50% peak performance utilisation, which efficiency is representative of our CCSD(T) code in Mrcc [37]. The time estimates on the horizontal axes take into account the number of operations required for the rate determining fifth-and sixth-power-scaling terms for CCSD and the sixth-and seventh-power-scaling terms for (T). For instance, the large, ha5Z computations requiring about 10 5.2 and 10 5.1 teraFLOP for the CCSD and (T) parts, respectively, can be translated to about 1.9 and 1.4 days of wall time at a useful performance of 1 teraFLOPS, which is in accord with our measurements. Finally, in Figures 1 and 2, we also take into consideration the cost of CBS extrapolation, that is, the operation counts for the (X−1)Z and XZ computations, which are required to obtain the CBS(X−1,X) extrapolated values, are added together on the x axes. Consequently, the CBS(X−1,X) extrapolated results in the second panels of Figures 1 and 2 are slightly more demanding than the corresponding XZ results. The plotted data can be found in Tables S1  and S2   coloured vertical blocks collect results obtained with the same X cardinal number. Considering first the conventional CCSD results (red 'CCSD noCP' and orange 'CCSD CP' curves of Figure 1 with square symbols), we find a more rapid convergence for the raw IE contributions than with CP corrections. This pattern is not the most common case for IE computations as CP corrections usually improve the results, but this trend is characteristic of the first, hydrogen bonding dominated third of the S66 compilation (see Figure 3 below). While the CCSD/haDZ and CP corrected CCSD/haTZ results did not fit into Figure 1, the raw and the CP corrected CCSD results tend to about −2.872 and −2.875 kcal/mol, respectively, at the CCSD/ha(Q,5)Z level. Compared to that, the CCSD(F12 * ) results, both with XZ-F12 (triangles) and haXZ-F12 (circles), and both without (blue and purple) and with CP (dark and light green) are already fairly accurate at the (ha)DZ-F12 level. However, considering the convergence pattern of the F12-based IEs, the good (ha)DZ-F12 performance is partly the result of error compensation. Nevertheless, (CP corrected) CCSD(F12 * )/ha(D,T)Z-F12 or better results can be considered converged and accurate within 1%. The agreement of raw and CP corrected CCSD(F12 * )/(T,Q)Z-F12 within 0.025 kcal/mol and of raw and CP corrected CCSD(F12 * )/ha(T,Q)Z-F12 within 0.004 kcal/mol is again convincing. Additionally, the CCSD(F12 * )/ha(T,Q)Z-F12 results without and with CP of −2.866 and −2.862 kcal/mol agrees with the CCSD/ha(Q,5)Z results within 0.01 kcal/mol or 0.34% remaining uncertainty.
Considering the computational cost, it increases by about an order of magnitude with one step increase in the cardinal number. A much less emphasised aspect is that CCSD(F12 * )/XZ-F12 and CCSD(F12 * )/haXZ-F12 computations are about 3-5 and 10-12 times more demanding than CCSD/haXZ in this case. Since CCSD(F12 * )/(h)aXZ-F12 could be an order of magnitude more costly than CCSD/XZ, it is not fully informative when these kinds of cost analyses are simplified to only the X cardinal number value. Nevertheless, the additional costs of CCSD(F12 * )/XZ-F12 computations certainly return, whereas it may depend on the particular application and target accuracy if the investment of the twice more operations with the additional diffuse functions in haXZ-F12 is worthwhile.
Continuing with the (T) contributions in Figure 2, the CP corrected (orange) and uncorrected (red) (T)/XZ results (square point type) behave similarly to the CCSD ones: the raw curve is fairly flat while the CP corrected results converge slower but in a satisfactory way. The two sets of results both approach −0.895 kcal/mol within 0.001 kcal/mol at the (T)/ha(Q,5)Z level, but their agreement is already excellent at the ha(D,T)Z and ha(T,Q)Z levels. However, we note that the CP corrected (T)/haDZ result of −0.686 kcal/mol (not plotted), exhibiting 23% error, does not appear to be highly reliable for use in extrapolation. In accord with numerical experience in the literature, the triples corrections using the CCSD(F12 * ) amplitudes, even the scaled (T+) variant showed here, converge more slowly to the CBS limit. Since the basis set incompleteness error of the triples corrections is not decreased directly via the explicitly correlated treatment, the BSSE leads to overbinding, especially with the raw (T+) results. The CP correction uniformly improves the (T+) results, which are again within 0.01-0.02 kcal/mol of the (T)/ha(Q,5)Z values starting from the (T+)/(ha)TZ-F12 level. Compared to the case of CCSD, the cost increase when using XZ-F12 and haXZ-F12 bases instead of XZ is milder as only about 1.5 to 3 times more operations are required for the (T) part with the (ha)XZ-F12 bases. Therefore, (ha)TZ-F12 level (T+) results, especially with CP corrections, are competitive, but (T)/ha(T,Q)Z appears to be somewhat more robust in this respect.
The comparison of the CCSD and (T) trends in Figures 1 and 2 is also instructive. Since the (T) term does not benefit directly from the F12 treatment, at least for this example, it is similarly hard to converge the CCSD and (T) IE contributions to close to 0.01 kcal/mol accuracy. The comparison of the CCSD and (T) expenses also show a trend that could be unexpected considering only the higher-order scaling of (T) and the corresponding common assumption of its higher cost. However, here, we find the opposite cost relation, partly caused by the 12 CCSD iterations, only 24 correlated electrons, and the better parallel scaling of (T). Namely, using the same basis set, the (T) term requires more operations than CCSD only with haDZ and haTZ; with the other basis sets, the CCSD or CCSD(F12 * ) part is more costly than (T) by 10-220%. This adds an interesting perspective to the previous best reference of Ref. [26] labelled SIL-VER, where CCSD(F12 * )/aTZ-F12 was combined with (T)/ha(D,T)Z. Considering the required operations and not just the basis set cardinal numbers, it is interesting to see that obtaining (T)/haQZ for (AcOH) 2 is only about 40% more demanding than CCSD(F12 * )/aTZ. While the extra cost of (T) is higher with more correlated electrons, still, this analysis uncovers an opportunity to improve upon the previous best S66 reference with an affordable computational effort. Naturally, the (T)/haQZ over CCSD(F12 * )/aTZ cost ratio is higher for the dimers with more electrons but remains lower than 3.4 even for the largest uracil dimers.

Conventional and F12-based CCSD(T): new reference for the S66 set
Utilizing our recent, well-parallelised, and integral direct codes [37,62] we obtained DF-based CCSD(T) results up to ha(T,Q)Z and CCSD(F12 * )(T+) results up to ha(D,T)Z-F12 for the entire S66 list. Considering the numerical experience of related studies for smaller or even S66-sized dimers [17,25,26,[42][43][44][58][59][60][61], it is still not obvious which set of benchmarks is better from the following perspectives. For instance, the trends noted for (AcOH) 2 are not in accord with the common recommendations of using CP corrected or half-CP corrected results. The application of CBS extrapolation should also be investigated carefully, both for CCSD(F12 * ) and especially for (T+). It is not obvious which set of results are more reliable: the CCSD/ha(T,Q)Z or the CCSD(F12 * )/ha(D,T)Z-F12 results and the (T)/ha(T,Q)Z or the scaled (T+)/ha(D,T)Z-F12 results. Finally, it is unclear, whether the current standard of small/medium basis CCSD(T) and MP2-based CBS correction type composite approaches should be considered here, e.g. due to the well-documented overestimation of MP2-based π -π stacking IEs.
First, one may look at the difference of the CP corrected and uncorrected results as a measure of the basis set convergence. This is a helpful indicator, but cannot determine alone whether the raw, CP, or half-CP results are the most accurate. Therefore, we also analysed the basis set convergence patterns of these possibilities. However, we could not find a clear connection between the absolute accuracy and measures of the level of convergence, such as the size of the last step taken between, e.g. the haDZ-F12 and haTZ-F12 results. The reason behind this is illustrated by the smaller steps of CCSD(F12 * )/XZ-F12 compared to CCSD(F12 * )/haXZ-F12 in Figure 1. Therefore, a somewhat more complicated route is chosen, that is, we looked at the convergence trends for all S66 dimers individually in order to find the more reliable reference.
The correlation energy contributions of the CCSD and (T) results are collected in Figures 3 and 4, respectively. The dispersion dominated dimers are separated via gray shade in the middle from the hydrogen-bonded (left) and mixed (right) dimers. Conventional and mostly F12-based CCSD deviations from the CP corrected CCSD(F12 * )/ha(D,T)Z-F12 reference are gathered on the top and bottom panels of Figure 3, respectively. The first group of dimers with H-bonding character on the left, including (AcOH) 2 highlighted in Section 4.1, exhibits different trends compared to the dispersion dominated (middle) and mixed (right) dimers. Interestingly, the raw results (solid symbols) overbind only above the haDZ level, while the raw CCSD/haDZ IEs are underestimated, which do not appear to be dominated by BSSE. For the dispersion and mixed groups, CCSD/haXZ IEs with X = D, T, Q exhibit a satisfactory trend much more suitable for CBS extrapolation. In contrast to that, the CP results (hollow symbols) are overcorrected as usual but exhibit a more reliable convergence pattern for the entire S66 set and, with a few exceptions, smaller errors than their raw counterparts.
Consequently, without CP, CCSD/ha(D,T)Z results (light orange pentagons) are often even worse than or comparable to raw CCSD/haTZ. However, with CP, the CCSD/ha(D,T)Z values are about 2-3 times better than even the CP corrected CCSD/haQZ (dark green diamonds). The raw and CP corrected CCSD/ha(T,Q)Z errors (green circles on the bottom panel) are more comparable to the performance of CCSD(F12 * ). The ha(T,Q)Z extrapolation overshoots in both cases, but much less with CP than without it, thus a half-CP type averaging of raw and CP results is not helpful for CCSD/ha(T,Q)Z. The CCSD(F12 * )/haTZ-F12 calculations, both with and without CP (blue squares on the bottom panel), are more efficient and accurate than CCSD/ha(T,Q)Z, while CCSD(F12 * )/haDZ-F12 is comparable to CCSD/haQZ at a much smaller cost. Except for a few H-bonded dimers, the raw and CP results approach the reference consistently from below and above, respectively. The outliers again limit the accuracy of the raw and half-CP variants of CCSD(F12 * )/ha(D,T)Z-F12 due to their dependence on raw CCSD(F12 * )/haDZ-F12. In contrast, the raw and CP corrected CCSD(F12 * )/haTZ-F12 results are arranged fairly symmetrically on the two sides of our reference. Therefore, the half-CP CCSD(F12 * )/haTZ-F12 combination, which lacks only the diffuse AOs on H atoms with respect to the half-CP CCSD(F12 * )/aTZ-F12 reference choice of Ref. [26], is within 0.01 (0.02) kcal/mol RMSD (MAX error) of the CP CCSD(F12 * )/ha(D,T)Z-F12 reference selected in the present work. Our reference choice is supported by the much more homogeneous convergence trend of CCSD(F12 * )/haDZ-F12 and CCSD(F12 * )/haTZ-F12 with CP correction than without CP and thus ensures that the less well-behaved raw CCSD(F12 * )/haDZ-F12 results are not employed in the final reference.
The case of the (T) corrections, collected in Figure 4, is different since (i) the (T) contribution to the IE is 3-5 times smaller than that of CCSD, (ii) the (T) term does not directly benefit from the F12 treatment, but (iii) it is often assumed that (T) converges to its CBS limit faster than CCSD. Indeed, compared to the CP (T)/ha(T,Q)Z reference, the ratio of CCSD and (T) IE contributions and their deviations from the reference appear to be similar, leading to considerably smaller absolute errors in the (T) term. Similar to the case of CCSD, raw (T) errors with haDZ (red triangles) and haTZ (blue squares) have much less homogeneous behaviour than their CP corrected counterparts. However, for the (T) corrections, the dispersion dominated complexes are the outliers with larger negative errors at the haDZ and haTZ levels. Even though the CP corrected haXZ (X = D, T, Q) errors are, on the average, larger than the corresponding uncorrected ones, their homogeneous sign and convergence trend is much more suitable for extrapolation. Since mostly both the raw and the CP corrected results underestimate the reference (without CBS extrapolation), the use of half-CP is not advised. Considering the extrapolated values, the raw and CP corrected (T)/ha(D,T)Z results (light orange and magenta pentagons) agree relatively well, except for the most complicated π-π stacking dimers containing benzene, pyridine, and uracil monomers (no. [24][25][26][27][28][29]. Thus, especially for the dispersion dominated complexes in the gray shaded area, the use of the half-CP (T)/ha(D,T)Z reference of the previous best S66 benchmark [26] is not the best choice in the light of our new reference numbers.
While CP (T)/ha(D,T)Z is fairly convincing with only about twice as large errors as raw (T)/ha(T,Q)Z (purple circles of the bottom panel), the homogeneous convergence pattern of the CP corrected results suggests that CP (T)/ha(T,Q)Z is the best choice of these three options. Nevertheless, the average absolute deviation of raw (T)/ha(T,Q)Z from the CP (T)/ha(T,Q)Z reference is 0.003 kcal/mol, thus their average, that is, half-CP (T)/ha(T,Q)Z could be just as good for most purposes. Since the (T+)/haXZ-F12 results are found to be better than the unscaled (T)/haXZ-F12 ones [both obtained with CCSD(F12 * ) amplitudes], here, we focus on the performance of (T+) and report the corresponding (T)/haXZ-F12 results in Figure S1 of the SM. Interestingly, the (T+) results, especially with CP, are fairly accurate even with the smallest haDZ-F12 basis set (red triangles of the bottom panel). The raw (T+)/haDZ-F12 and (T+)/haTZ-F12 results are surprisingly close to each other, while CP (T+)/haTZ-F12 offers substantial improvement over CP (T+)/haDZ-F12. In fact, at least compared to the CP (T)/ha(T,Q)Z reference, CP (T+)/haTZ-F12 is the best choice if haQZ level results are not affordable. Let us also note the remarkable agreement of the raw and CP corrected (T+)/ha(D,T)Z-F12 results (light orange and green pentagons) with each other, but not as much with the CP (T)/ha(T,Q)Z reference. However, raw and CP (T)/ha(T,Q)Z agree with each other twice as well and appear to be more reliable. Moreover, we know much less about the behaviour of the very recent (T+) method [62]. For instance, it is unclear whether the (T+) correlation energies are subjectable to CBS extrapolation and if so, what would be the suitable exponent of X. This question will be a subject of a future study. Until then, we recommend to use CP (T+)/haTZ-   the basis set convergence. Indeed, in Table 1 the maximum 'raw-CP' values are considerably (2-6 times) higher than the MAX error of the best performing CCSD variant for a given basis set. This also suggests that CP CCSD(F12 * )/ha(D,T)Z-F12 is converged considerably better than its 0.061 kcal/mol maximum 'raw-CP' measure. The largest 'raw-CP' value corresponds to (AcOH) 2 (closely followed by the uracil dimer and the uracil-AcOH dimer). However, we know from the analysis of Section 4.1 that CP CCSD(F12 * )/ha(D,T)Z-F12 differs from the CBS limit by about 0.003 kcal/mol, and thus the major part of the 0.061 kcal/mol 'raw-CP' value comes from the error of the raw number. We also note that the CCSD(F12 * )/aTZ-F12 values used as reference in Ref. [26] (last row of Table 1) closely agree with our (DF-based) CCSD(F12 * )/haTZ-F12 results, indicating that the additional p-and d-type diffuse AOs on the hydrogen atoms of aTZ-F12 contribute negligibly to the IEs. Additionally, the 'raw-CP' value of CCSD(F12 * )/(h)aTZ-F12 is about twice as large as that of CCSD(F12 * )/ha(D,T)Z-F12, thus up to 0.02 kcal/mol improvement can be expected from updating the S66 reference to CP CCSD(F12 * )/ha(D,T)Z-F12.
Considering the corresponding statistical measures for the (T) contributions in Table 2, the hierarchy of the raw, CP, and half-CP results is much less clear cut compared to the case of CCSD. One commonality is that full CP appears to be best of all types of CBS extrapolated results, while raw is mostly but not always better than half-CP and CP without CBS extrapolation. The behaviour of the 'raw-CP' measure also differs for (T): here, it is often very close to or the same as the sum of the MAX raw and MAX CP errors, while in a few cases 'raw-CP' does not appear to be larger than the MAX errors [cf. the (T+)/ha(D,T)Z-F12 case]. Nevertheless, the at most 0.008 kcal/mol MAX deviation of the raw and CP corrected (T)/ha(T,Q)Z results is very convincing and suggests that our new (T) reference is accurate within cca. 0.01 kcal/mol. This estimation holds for the case of (AcOH) 2 , where the (T)/ha(T,Q)Z IE components are within 0.008 and 0.005 kcal/mol of the (T)/ha(Q,5)Z values with and without CP, respectively.
The best converged result sets, that is, the CP (T)/ha(T,Q)Z reference, CP (T)/ha(D,T)Z-F12, and CP (T+)/haTZ-F12 are only about twice as close to each other as the analogous differences for CCSD in Table 1, indicating that converging the (T) contribution is similarly important, already for medium-sized complexes. The half-CP (T)/ha(D,T)Z variant employed in Ref. [26], with 0.011 and 0.050 kcal/mol RMSD and MAX errors, respectively, seems to be less reliable than its CP corrected analogue.
It is also worth noting the performance of the sizeconsistent (T+) scaling approach compared to the previously employed (T * ) and (T * sc ) scaling alternatives. In agreement with the observations of Ref. [62], the difference of the (T+) and (T * ) results is the largest for dimers with markedly different monomers, e.g. for the complexes of benzene or pyridine formed with water, MeOH, or AcOH. This size-consistency related difference can amount up to 0.1 kcal/mol with haDZ-F12 and remains too large for our purposes (up to 0.06 and 0.04 kcal/mol) even at the haTZ-F12 and ha(D,T)Z-F12 levels. Even if, for the dimers with markedly different monomers, the size-consistent (T * sc ) version is considerably better than (T * ), when all complexes are considered, the maximum deviations of (T+) and (T * sc ) are even higher than those of (T+) and (T * ), namely 0.2, 0.05, and 0.05 kcal/mol at the haDZ-F12, haTZ-F12, and ha(D,T)Z-F12 levels, respectively. The worst performers for (T * sc ) are the π -π stacking dimers containing benzene, pyridine, and uracil monomers (no. [24][25][26][27][28][29], where the MP2-based scaling factors obtained for the dimers are not adequate for the monomers. The performance of (T * ) and (T * sc ) found here is in accord with the findings of Ref. [25] for the S66x8 set, where such high size-consistency errors were also characterised as unacceptable.
In summary, we suggest combining CP corrected CCSD(F12 * )/ha(D,T)Z-F12 and (T)/ha(T,Q)Z to update the CCSD(T) reference IEs for the S66 set. This corresponds to the level of theory for the (T) term labelled as GOLD in the previous best converged S66 benchmark of Ref. [26]. The CCSD IEs employed here should be more converged than the SILVER level of Ref. [26], but we cannot test if their accuracy reaches that of the GOLD level of CCSD(F12 * )/cc-pVQZ-F12 suggested in Ref. [26]. Let us note in passing, that the largest CCSD(F12 * )/cc-pVQZ-F12 computation for S66 would involve 1664 AOs, which should be feasible in the near future upon some memory optimisation due in our CCSD(F12 * ) code. Regarding the mean-field contribution, the half-CP HF+CABS/aQZ-F12 results of Ref. [26] are highly converged and suitable for combination with our new correlation energy contributions.
One remaining question is whether to add an MP2 based basis set correction ( MP2) to the CCSD results. Previous numerical experience shows that less converged CCSD results benefit from the MP2 correction. However, it is not possible to perform similar analysis comparing to better references at the size of the S66 set and at the convergence level of CCSD(F12 * )/ha(D,T)Z-F12. On the one hand, higher level of basis set convergence can be achieved for MP2 than for CCSD. For instance, the MP2-F12/a(T,Q)Z-F12 results of Ref. [26] exhibit for the 'raw-CP' convergence measure smaller RMSD and MAX values of 0.001 and 0.004 kcal/mol, respectively. This indicates higher level of basis set convergence for MP2-F12 at the a(T,Q)Z-F12 level than obtained for CCSD(F12 * ) using only the ha(D,T)Z-F12 basis set level. For comparison, CP MP2-F12/ha(D,T)Z-F12 deviates from MP2-F12/a(T,Q)Z-F12 by 0.013 and 0.040 kcal/mol RMSD and MAX error values. One counterargument to using MP2-based basis set corrections is that the MP2/CBS results are quite far from the CCSD/CBS ones. Namely, the relation of MP2-F12 and CCSD(F12 * ) IEs at the ha(D,T)Z-F12 level ranges from 10% underestimated [e.g. for (AcOH) 2 ] to 30% overestimated (e.g. for π-π stacking dimers), and it is little known about how these errors could influence the MP2-based CBS correction in combination with such well converged CCSD data.
For a more detailed understanding, the differences of CCSD and MP2 IEs contributions are collected in Table 3 compared to our best converged CP [CCSD(F12 * ) − MP2-F12]/ha(D,T)Z-F12 reference. In brief, the results without explicit correlation do not exhibit strictly monotonic convergence with increasing basis set size, while the F12-based results show convincing trends. However, compared to the basis set completeness level of the CCSD results in Table 1, the [CCSD(F12 * ) − MP2-F12]/ha(D,T)Z-F12 differences appear to be only about three times better converged. Therefore, the inaccuracies in MP2 can also affect the MP2-based CBS corrections, and we cannot expect that such MP2 corrected CCSD results match the convergence level of MP2-F12/a(T,Q)Z-F12. All in all, adding the MP2 correction to our final reference model seems to be beneficial, leading to: half-CP HF+CABS/aQZ-F12 + CP MP2/a(T,Q)Z-F12 + CP [CCSD(F12 * ) − MP2-F12]/ha(D,T)Z-F12 + CP (T)/ha(T,Q)Z. This composite scheme involves higher level CCSD and (T) components than the one termed '14-carat GOLD' in Ref. [26]. In order to simplify the discussion below we adopt the phrase of Ref. [26] and label our new reference as 14k-GOLD. The resulting 14k-GOLD reference IEs are given in Table 4.
Using the 14k-GOLD CCSD(T) reference, it is worthwhile analysing the performance of the lower-cost results collected in Table 5 both without and with MP2based CBS corrections (left and right halves of Table 5, respectively). The RMSD measures of Table 5 show that the use of the MP2 CBS corrections is beneficial with smaller basis sets up to haTZ, while the haQZ and ha(T,Q)Z results do not benefit from the correction. This trend appears to be reversed for the CCSD(F12 * )(T+) results, that is, MP2 corrected haTZ-F12 results are somewhat improved, while the uncorrected CCSD(F12 * )(T+)/haDZ-F12 IEs are preferred. The picture is much less clear when comparing the raw, CP, and half-CP as well as their MP2 (un)corrected versions with the same basis set. Depending on the AO basis set and the use of F12 approaches, basically any of the six CP and MP2 correction combinations can come up on top, and it is very hard to disentangle the compensation between the various sources of errors.
From the combinations of Table 5, CP CCSD(F12 * ) (T+)/haTZ-F12 + MP2 and CP CCSD(F12 * )(T+)/ haDZ-F12 stand out as particularly accurate lower-cost alternatives. These only differ in a few aspects from the SILVER and BRONZE level composite approaches of Ref. [26] (cf. caption of Table 6). For a more detailed comparison, statistical measures for previous CCSD(T) benchmarks, including the SILVER, STERLING, and BRONZE composite schemes of Ref. [26], as well as   the original (orig.) [17] and revised (rev.) [24] references of Řezáč and Hobza are collected in Table 6 with respect to the 14k-GOLD reference proposed here. Apparently, all but maybe the original set of IEs are found highly reliable, being within 0.075 kcal/mol MAX error and 0.033 kcal/mol RMSD of the new 14k-GOLD results. Considering the level of convergence for the individual terms in these composite schemes, such as the CCSD(T)/ha(D,T)Z part of "S66 revised" or the CCSD(F12 * )/(Tc sc )/DZ-F12 part of BRONZE, they benefit from some compensation of errors as also noted in Ref. [26]. These findings highlight the particularly successful design of previous composite schemes developed at the time when higher-level references were available only for the smallest few dimers in S66 [17,24]. To assess the progress made via 14k-GOLD, let us also look at the difference of the raw and CP corrected results corresponding to the energy components of the composite schemes as a measure of their basis set convergence level (see last row of Table 6). The respective values corresponding to 14k-GOLD are 0.008 and 0.024 kcal/mol RMSD and MAX difference, indicating that the 0.1 kcal/mol interaction energy gold-standard is indeed achieved. Compared to that, the same measures for SILVER are 0.023 and 0.056 kcal/mol RMSD and MAX difference, that is, they are about 2-3 times less converged according to this measure. The maximum of the |raw − CP| measure increases considerably up to 0.16, 0.47, 0.43, and 0.27 kcal/mol for the STERLING, BRONZE, S66 revisited, and S66 original benchmarks. Thus, the new 14k-GOLD reference exhibits a marked improvement regarding this basis set completeness measure, while it is proven to be very challenging to considerably improve upon the accuracy of previous composite Table 5. RMSD with raw, CP, and half-CP (CP/2) and without and with the MP2 basis set correction for CCSD(T) and CCSD(F12 * )(T+) correlation energy contributions to the IEs compared to the here updated 14k-GOLD S66 reference [in kcal/mol].   energy expression based IEs. As expected, the largest difference between SILVER and our new reference is found for the dimers with the slowest basis set convergence, namely π-π stacked complexes (e.g. uracil 2 and uracilpyridine) and H-bonded complexes [e.g. (AcOH) 2 ]. Having CCSD(F12 * )/haQZ-F12 and (T)/ha(Q,5)Z references at hand for (AcOH) 2 , we can be very confident that the new CCSD and (T) reference IE components are converged within 0.003 and 0.005 kcal/mol, respectively, while it would require an enormous effort to similarly assess the accuracy for the larger π-π stacked dimers. Besides the analysis of IE errors, it is also important to consider the relative errors as the S66 IEs spread across a wide range from −1.3 to −19.7 kcal/mol. In brief, we inspect the |raw−CP| RMSD and MAX error basis set convergence measures relative to the 14k-GOLD reference. These are found for the individual energy components to be as follows: 0.036% and 0.103% for MP2-F12/a(T,Q)Z-F12; 0.207% and 0.708% for [CCSD(F12 * )−MP2-F12]/a(D,T)Z-F12; 0.109% and 0.286% for (T)/ha(T,Q)Z. That is, the corresponding relative RMSD and MAX uncertainties are 0.266% and 0.784% for the total CCSD(T) IEs. Thus, this robust measure indicates that the IEs are converged with at most 0.8% relative accuracy, while the same measure is still accurate for the SILVER IEs being about twice as large (1.5%) as for 14k-GOLD.

Local natural orbital and other local CCSD(T) approaches
Local correlation approaches, such as the LNO-CCSD(T) as well as PNO-based schemes considered here, offer significant cost reduction for CCSD(T). For instance, they facilitate the increase of the basis set used for mediumsized molecules in the size range of S66 or allow for computations for much larger complexes than accessible via conventional CCSD(T). However, all of these local correlation methods introduce a set of approximations which have to be carefully assessed. Thus, first we exploit the possibility to test the accuracy of various local approximations against the new, large basis set references since so far local CCSD(T) benchmarks against exact CCSD(T) were not possible in the size range of S66 up to haQZ or ha(T,Q)Z [89,91,98,100]. Second, we utilise wellconverged LNO-CCSD(T) IEs up to ha(Q,5)Z to assess the basis set convergence level at the ha(T,Q)Z level used also for the (T) part of the new 14k-GOLD S66 reference IEs.
First, it is worthwhile utilising the CCSD(T)/ha5Z IEs of (AcOH) 2   Loose to very Tight (vTight) threshold settings as well as the corresponding Loose-Normal (L-N), etc. extrapolated values. The right panels show DLPNO-CCSD(T1) IEs with the NormalPNO and TightPNO settings as recommended by the developers [87,101]. Reference DF-CCSD(T) IEs, indicated via dashed horizontal lines in Figure 5, represent the local approximation free limit of both LNO-CCSD(T) and DLPNO-CCSD(T1) IEs due to the consistent use of the DF approximation in all three approaches.
Considering that (AcOH) 2 is one of the most complicated dimers of the S66 set in terms of the basis set requirement, the rate of convergence of the LNO-CCSD(T) results and their balanced accuracy across the basis sets is highly satisfactory. The Loose LNO-CCSD(T) results are only useful for quick but crude exploration, thus Normal, Tight, or vTight settings are required to reach, respectively, 0.3-0.4, 0.1-0.2, or even 0.02-0.05 kcal/mol deviations from the reference both with and without CP. In this case the LNO extrapolation consistently improves the IEs by about a whole step in the LNO setting hierarchy for all basis sets. For instance, the quality of the L-N and Normal-Tight (N-T) extrapolated LNO-CCSD(T) results match that obtained with the more demanding Tight and vTight settings, respectively. More importantly, the error bars assigned to the extrapolated IEs according to Equation (5) estimate the size of the next step in the hierarchy well, that is, Tight (vTight) results are indeed within the range defined by the error bars of L-N (N-T). As the L-N results are still about 0.2-0.3 kcal/mol away from the exact reference, their respective error bars slightly underestimate this distance toward the local approximation free limit. Therefore, for high accuracy production computations N-T or better settings are recommended, while Loose thresholds are useful to obtain an additional set of results cost-effectively to monitor the convergence pattern. It is interesting to note that the case of raw (AcOH) 2 IEs does not follow the trend observed in our previous LNO-CCSD(T)-based investigations of IEs [27,91,[107][108][109], that is, the basis set error is usually larger than the error of the local approximations. Instead, the CP corrected (AcOH) 2 IEs on the bottom of Figure 5 exhibit the more typical case of larger basis set incompleteness than LNO errors.
For the DLPNO-CCSD(T1) results on the right panels of Figure 5, the homogeneity of the IE errors with various basis sets is similarly satisfactory. At least for the case of (AcOH) 2 , the NormalPNO settings are a bit less accurate than Loose LNO-CCSD(T). In general, NormalPNO DLPNO-CCSD(T1) performs better than Loose LNO-CCSD(T) for the entire S66 set as shown below, e.g. in Figure 6. Moreover, in a large step of improvement, TightPNO DLPNO-CCSD(T1) IEs for (AcOH) 2 become comparable to Tight or L-N LNO-CCSD(T) for all basis sets evaluated.
Let us continue with a similar convergence analysis for the entire S66 set. Tables 7 and 8   correlation and interaction energies were found to be practically identical to those without CP corrections for both the LNO and the DLPNO approach.
Let us start with the MAE and MAX measures of the relative correlation energy errors collected in Table 7 with the haTZ, haQZ, and ha(T,Q)Z basis set levels. The top third of Table 7 shows the systematic convergence of the LNO-CCSD(T) correlation energies from Loose up to very very Tight (vvTight) settings. For instance, Loose LNO-CCSD(T)/haTZ correlation energies are already at the 0.045% relative error level, which decreases, in a highly regular manner, by about a factor of 2, in each step throughout the tightest settings. The MAX LNO-CCSD(T)/haTZ errors are similarly consistent being about 3 times higher than the corresponding MAE, that is, 99.9% accuracy is provided for all S66 dimers already with the Normal settings. Both error measures increase notably by about a factor of 2 when moving from haTZ to haQZ and then to ha(T,Q)Z, but the Tight results still maintain the at least 99.9% accuracy even at the ha(T,Q)Z level. Interestingly, LNO extrapolated results, collected in the middle third of Table 7 exhibit much more homogeneous relative correlation energy errors across the three basis set levels than those without extrapolation. While the LNO-CCSD(T) correlation energies show minor improvements upon LNO extrapolation with the haTZ basis, the LNO extrapolation drastically improves the haQZ and ha(T,Q)Z correlation energy errors in Table 7 to a level comparable to the accuracy of LNO-CCSD(T)/haTZ. Since this is the first possibility to perform such an analysis on the LNO extrapolation for sufficiently large molecules at the haQZ basis set level, much more data are required to check the generality of the above trends.
Finally, let us compare the CP corrected LNO-CCSD(T) results using the updated domain construction algorithm of Section 2.2. As noted, the CP uncorrected and corrected interaction energies are almost identical. Namely, Loose MAE values for raw and CP differ at most by only 0.004%, which gradually disappears as the LNO thresholds are tightened. Interestingly, the CP corrected monomer correlation energies are barely affected by the updated domain construction algorithm, changes in the above relative correlation energy error statistics are negligible. The largest improvement is observed for the water monomer of the water-pyridine dimer (no. 18), where the change is 0.02% with Loose settings. The average absolute change for the entire set with Loose settings is only 0.002% and decreases even further when the LNO thresholds are tightened. This trend can be explained considering the relatively small size and the relatively short, equilibrium distance of the monomers. While the results with the increased domain sizes are encouraging from the perspective of the previous algorithm [36,91], nonnegligible changes may occur for larger complexes (e.g. for the L7 set [18,27]) or for more distant monomers, such as those in the S66x8 compilation [17].
The analogous DLPNO-CCSD(T) relative correlation energy errors are also collected in the bottom third of Table 7 with both NormalPNO and TightPNO settings, as well as with both the DLPNO-(T0) and the DLPNO-(T1) treatment of the triples corrections [87,100]. It is important to highlight that the DLPNOapproximation free limit of the DLPNO-CCSD(T0) correlation energies would be a DF-CCSD(T0) set of reference values, which is not identical to the DF-CCSD(T) reference used here. Nevertheless, looking at the DLPNO-CCSD(T0) results too reveals that the differences of the (T0) and (T1) correlation energies are about 0.2% across all basis sets and DLPNO thresholds settings considered. Interestingly, the DLPNO-CCSD(T1) results improve with increasing basis set size. While TightPNO DLPNO-CCSD(T1)/haTZ is comparable to Loose LNO-CCSD(T)/haTZ, NormalPNO DLPNO-CCSD(T1)/haQZ is closer to Normal LNO-CCSD(T)/haQZ, which trend may be related to the compensation between basis set and DLPNO errors noted Table 7. MAE and MAX relative correlation energy errors (in %) with respect to the DF-CCSD(T) reference for various LNO-and DLPNO-CCSD(T) settings without CP corrections. Top and middle thirds: LNO-CCSD(T) using various composite threshold sets without and with LNO extrapolation, respectively. Bottom third: DLPNO-CCSD(T) using NormalPNO and TightPNO settings as well as the (T0) and (T1) triples correction algorithms.  previously [140]. Considering both the MAE and MAX relative correlation energy measures and all basis sets, TightPNO DLPNO-CCSD(T1) practically reaches 99.9% accuracy for all species in S66. This is, however, slightly worse than the convergence level of the more costefficient Loose-Normal LNO-CCSD(T) correlation energies, indicating that a smaller compensation of errors is sufficient in LNO correlation energies to provide comparable quality IEs. Finally, since not many benchmarks are available in the literature dedicated to assess the performance of the DLPNO method with CP corrections, it is valuable to note on the near identical performance of all DLPNO variants with and without CP corrections for all considered basis set choices.
To continue with the analogous analysis for IEs, the corresponding MAE and MAX IE error measures with respect to the DF-CCSD(T) reference are collected in Table 8 for the same LNO and DLPNO settings and basis sets. For the sake of brevity, we again show only the CP uncorrected statistics because the CP corrected ones are practically identical for both local correlation approaches and for all three levels of basis set.
Regarding the LNO-CCSD(T) IEs, we again find systematic convergence and a factor of 2, sometimes even a factor of 3, improvement in the MAE values for all three levels of basis set. The MAX errors, being consistently about 3-4 times larger than the MAE values, thus follow a similar trend. The Normal (Tight) LNO settings provide 0.1-0.2 (sub 0.1) kcal/mol IE accuracy on the average, but Tight or vTight settings are needed for similar accuracy regarding the MAX errors. Comparing the haTZ, haQZ, and the ha(T,Q)Z results, we find a slight decrease in the quality of the IEs with increasing the basis set size. However, this effect is much milder than what is observed for the basis set dependence of the LNO-CCSD(T) correlation energies above. In other words, the larger correlation energy errors compensate very well when forming the interaction energies, indicating that the IEs may converge somewhat faster with the basis set size than the correlation energies. Such an effect also contributes to the trend that, except for the L-N variant, the LNO extrapolation does not considerably improve the IE error statistics. An additional factor to consider is that so close to the converged IEs, the natural orbital approximation may not dominate the overall LNO error. Thus, other sources of errors with possibly different signs could become comparable, which, in turn, can affect the performance of the LNO extrapolation. Nevertheless, such cases can be identified by observing the convergence pattern of the IEs of individual dimers.
The analogous data for the DLPNO-CCSD(T) variants are collected in the bottom third of Table 8. The performance of all four DLPNO-CCSD(T) variants, just as for the correlation energies, is consistent for all three basis set levels. Based on the MAE and MAX IE error measures alone, both the (T0) and the (T1) variants with NormalPNO settings are about 15-45% more accurate than Loose LNO-CCSD(T), while the TightPNO variants are comparable to or similarly outperform Normal LNO-CCSD(T) with the haTZ and the larger basis sets, respectively.
Besides considering only the MAE and MAX error measures, it is useful to look at the complete error distributions too. Since the IE error trends are fairly consistent for all three basis sets, and with and without CP for both the LNO and DLPNO methods, the IE error distributions are only displayed for the case of haTZ in Figure  6. Violin curves of Figure 6 show the distribution of the signed IE errors, where the height of the curve indicates the frequency of the signed IE errors corresponding to an error value on the vertical axis. On the right side of the error distribution curves, the horizontal lines of the boxes indicate the lower, median, and upper quartiles, respectively, encompassing the middle 50% of the data points. Whiskers extend to the most distant data point whose IE error value lies within 1.5 times of the interquartile range (IQR), the latter of which is the difference of the lower and upper quartiles. Outliers beyond the whiskers, if any, are represented by dots.
One apparent feature of Figure 6 is the high symmetry of the LNO-CCSD(T) error curves around the 0 kcal/mol error line and the correspondingly small mean of the signed errors. Clearly, most of the Loose and about half of the Normal LNO-CCSD(T) IEs show errors above 0.1 kcal/mol, but the two middle quartiles of Normal LNO-CCSD(T) are almost completely within the gray lines indicating the ±0.1 kcal/mol error mark. Moreover, most of the Tight and all of the vTight results fall in the ±0.1 kcal/mol error range. It is interesting to have a closer look at the dots representing the outliers. For both Normal and Tight LNO-CCSD(T), the IE of the most challenging H-bonded (π -π stacking) dimers are consistently underestimated (overestimated) [e.g. H-bonded complexes of AcOH, AcNH2, and uracil, as well as the ππ stacking dimers of uracil, benzene, and pyridine]. The latter trend was also observed for more extended π -π dimers investigated in Ref. [27].
Considering the DLPNO-CCSD(T) results, the middle and right panels of Figure 6  Having the expected accuracy of the local approximations established at least up to the ha(T,Q)Z level of basis set, let us exploit their efficiency to approach the basis set limit of the IEs even closer, at the ha(Q,5)Z level. One complication is that the AO basis set for some of the dimers exhibit considerable linear dependency in the ha5Z basis. To our knowledge, only the LNO approach has features implemented to treat the quasi-redundancy of such large basis sets [91], thus the ha5Z results are only computed with the LNO-CCSD(T) method. For comparison, we also gather recently reported and well-converged PNO-based explicitly correlated CCSD(T) results with up to haQZ-F12 basis sets [98,99].
First, we consider the difference of the raw and CP corrected LNO-CCSD(T) results as a measure of the basis set convergence. While the average |raw−CP| value converges convincingly as 0.74, 0.26, and 0.09 kcal/mol with the haTZ, haQZ, and ha5Z basis sets, respectively, the ha(T,Q)Z and ha(Q,5)Z |raw−CP| results tend to similar values in the 0.04-0.05 kcal/mol range using even vvTight LNO settings. Compared to that, the average LNO-CCSD(T) error in the |raw−CP| measure is about 0.01 kcal/mol, at least at the ha(T,Q)Z level where the reference is known. This is in accord with the 0.017 kcal/mol MAE value in Table 8 obtained for the vvTight LNO-CCSD(T)/ha(T,Q)Z IEs. The corresponding MAX error being 0.07 kcal/mol indicates that the LNO error is probably somewhat larger than the basis set uncertainty for the CBS extrapolated results, but some compensation of LNO errors can be expected in the |raw−CP| measures. We can also compare the LNO-CCSD(T)/ha(Q,5)Z results to the new 14k-GOLD reference, keeping in mind that the latter has still up to a few tens of cal/mol uncertainty in its basis set convergence. The corresponding MAE and MAX error values are collected in Table 9 both without and with CP corrections. Again, the convergence toward the LNO approximation free limit is very convincing as the Tight MAE values already drop below 0.06 kcal/mol, while vTight and vvTight are very close to each other in the 0.02-0.04 kcal/mol MAE range. Interestingly, vTight LNO-CCSD(T)/ha(Q,5)Z results with CP are somewhat closer to the 14k-GOLD reference than with vvTight settings, indicating that CP corrected LNO-CCSD(T)/ha(Q,5)Z could tend to a slightly different CBS limit than the 14k-GOLD values obtained with smaller AO basis sets.
To inspect this in more detail, distributions of the signed deviation of LNO-CCSD(T)/ha(Q,5)Z with respect to the 14k-GOLD reference are plotted on the violin curves of Figure 7. The corresponding boxes and whiskers on the right side of the distributions again indicate the middle two quartiles and the 1.5 IQR. A notable difference can be observed if we compare the convergence pattern of the LNO-CCSD(T)/ha(Q,5)Z deviations with respect to the 14k-GOLD reference to the case of LNO-CCSD(T) errors against the exact DF-CCSD(T) reference in Figure 6. For instance, the widths of the distributions do not drop as significantly when stepping from Tight to vTight as against the exact reference, which indicates that the remaining basis set incompleteness error of 14k-GOLD may interfere with the systematically decreasing LNO errors when tending from Tight toward vTight and vvTight settings. This observation holds for both the raw and the CP corrected LNO-CCSD(T)/ha(Q,5)Z results on the two left panels  Figure 6 for details of the violin and the boxwhisker plot notation.
of Figure 7. However, raw LNO-CCSD(T)/ha(Q,5)Z IEs tend to an LNO approximation free limit that underestimates 14k-GOLD, while their CP corrected counterparts settle at the opposite side overestimating 14k-GOLD. Since 14k-GOLD is also constructed using CP corrected energy components, CP corrected vvTight LNO-CCSD(T)/ha(Q,5)Z may reveal here that 14k-GOLD is not attractive enough, at least on the scale of few tens of cal/mol. On the other hand, the maximum LNO uncertainty is probably somewhat larger than a few tens of cal/mol, especially for the π-π stacking dimers, which are the outliers among the CP corrected vvTight results. Therefore, it remains very challenging to separate the comparable sized local and basis set incompleteness errors in the difference of LNO-CCSD(T)/ha(Q,5)Z and 14k-GOLD.
Additional information may be obtained in this regard by inspecting alternative local correlation schemes, such as the PNO-LCCSD(T * )-F12b/haQZ-F12 results of Ma and Werner [99] as well as the PNO-CCSD(F12 * )(T)/aTZ IEs reported by Schmitz and Hättig [98]. Both sets of results utilise the explicitly correlated treatment at the CCSD level and fairly tight PNO truncation thresholds (T = 10 −8 ). The PNO-LCCSD(T * )-F12b results benefit from the large, haQZ-F12 basis set and CP correction, while the non CP corrected PNO-CCSD(F12 * )(T)/aTZ IEs take advantage of the excellent performance of CCSD(F12 * )(T) noted above in Sect. 4.2. Inspecting Table 9 and Figure 7, we find the PNO-LCCSD(T * )-F12b results with its default domain construction options about halfway between our CP corrected Normal and Tight LNO-CCSD(T)/ha(Q,5)Z results. The performance of PNO-CCSD(F12 * )(T)/aTZ is close to Tight LNO-CCSD(T)/ha(Q,5)Z with an even smaller number of outliers. The PNO-LCCSD(T * )-F12b results clearly benefit from the tight domain settings being even a bit closer to 14k-GOLD than the tighter LNO variants. Interestingly, the tight PNO-LCCSD(T * )-F12b/haQZ-F12 results, on the average, slightly underestimate the IEs of 14k-GOLD despite both having CP corrections.
All in all, the investigated local CCSD(T) approaches all have excellent qualities as well as some remaining drawbacks, which prevent us from drawing decisive conclusions. For instance, the LNO results are very tightly converged and reach the ha(Q,5)Z level, but do not benefit from the advantages of the F12 treatment. While the CCSD part of tight PNO-LCCSD(T * )-F12b/haQZ-F12 could be very well converged, the scaled (T * ) treatment still has uncertainties up to few tens of cal/mol, as pointed out here as well as in Ref. [99]. At this point, because of the agreement of several different sources, we can be confident that the uncertainty of our CCSD(T)/CBS estimates does not exceed 0.05 kcal/mol, but it will be quite challenging to shrink this uncertainty window further in a more decisive manner.
Finally, it is important to analyse the same data from the perspective of relative IE errors. This is a challenging property, especially for the dispersion dominated dimers. For these dimers, a large portion of the repulsive HF and attractive correlation energy contributions to the IEs partly cancel each other, and thus the total IE is often considerably smaller than the separate HF and CC components. This amplifies the errors caused in the correlation energy contributions, when they are compared to the smaller total IEs. Such relative IE errors, evaluated against the exact DF-CCSD(T) reference, are collected in Table 10.
Here, we find the MAE of Normal LNO-CCSD(T) at the 3-4% range, while Tight and vTight settings improve this metric to about 1% and 0.5%, respectively. Compared Table 10. MAE and MAX relative errors of total IEs (in %) with respect to the DF-CCSD(T) reference for various LNO-and DLPNO-CCSD(T) settings without CP corrections. to that, NormalPNO and TightPNO DLPNO-CCSD(T1) MAEs are at the 5-6% and 2% range, respectively. Thus, the comparison of LNO-CCSD(T) and DLPNO-CCSD(T) relative errors shows trends similar to the case of the absolute IE errors. However, the MAX errors being about 4-5 times larger than the MAE values reveal an important aspect. Namely, large relative errors can occur also for cases with fairly acceptable absolute errors if the total IE is also small. For instance, the largest relative error of 2.2% for vTight LNO-CCSD(T)/haTZ corresponds to the 0.04 kcal/mol IE error obtained for the neopentane dimer due to its about −1.94 kcal/mol total IE. Nevertheless, the MAX LNO-CCSD(T) errors reaching 2-3% and about 1.5% with vTight and vvTight settings, respectively, are fairly promising as it is feasible to obtain LNO-CCSD(T) results with such well-converged settings even above 100 atoms [27,91].

Conclusions and outlook
In this report, we utilised our state-of-the-art DF-based conventional [37], explicitly correlated [62], and local natural orbital (LNO) [36,91] CCSD(T) approaches to push forward the understanding of the basis set convergence of CCSD(T) interaction energies in medium-sized molecular complexes. These recent methodological and algorithmic developments allowed us for the first time to compute augmented quadruple-and quintuple-ζ level DF-CCSD(T) and LNO-CCSD(T) interaction energies (IEs), respectively, for the entire S66 molecular complex compilation. A comprehensive analysis of previous references and our new benchmarks suggest updating the S66 reference using the CCSD(F12 * )/ha(D,T)Z-F12based CCSD contributions combined with (T)/ha(T,Q)Z triples corrections. Following the terminology of Martin and co-workers [26], the new CCSD(T)/CBS estimate is labelled by (14-karat) 14k-GOLD. The new reference exhibits about 0.01 (0.05) kcal/mol average (maximum) basis set convergence uncertainty, while the relative accuracy of the correlation energies appear to fall below the maximum 1% mark for the first time. The CCSD(F12 * )(T+)/haTZ-F12 interaction energies, relying on our recent size-consistent (T+) method [62], offer a close second best option with only at most 0.016 kcal/mol deviation from 14k-GOLD without the need for quadruple-ζ level computations. While shrinking these uncertainties even further at the moment poses considerable challenges for conventional CCSD(T) codes, there is a reason for optimism. For instance, the largest CCSD(T) computation five years ago in the previous S66 benchmark work [26] contained 1200 atomic orbitals for the uracil dimer, which required 8 days runtime on 96 cores, 1.5 TB memory, and 18 TB SSD disk space. Compared to that, our largest CCSD(T) computation so far with 31% more orbitals required only 3 days on 224 cores and 0.33 TB memory using our recent implementation in the Mrcc package [37]. Employing the approximate but much more efficient LNO-CCSD(T) approach, we were able to estimate the S66 IEs at the LNO-CCSD(T)/ha(Q,5)Z level. Considering also the local approximation error in the latter, the uncertainty level of about 0.05 kcal/mol of the best CCSD(T) estimates has been corroborated using multiple independent approaches including previous benchmarks [26,99], 14k-GOLD, and tightly converged LNO-CCSD(T)/ha(Q,5)Z. Such highly converged results show that the deviation of DMC and CCSD(T) IEs at the 0.3±0.2 kcal/mol scale, e.g.for the benzene or uracil dimers [27], cannot solely be explained by the basis set incompleteness error of CCSD(T) IEs.
The new DF-CCSD(T)/haQZ references also allowed us to benchmark our LNO and the commonly employed DLPNO local correlation approaches against exact references for the first time with such a large basis set at the system size of the S66 dimers. Importantly, both the LNO and the DLPNO approaches were found to maintain their accuracy levels established up to the haTZ level also with the larger haQZ basis set, upon basis set extrapolation, and in combination with counterpoise corrections. For instance, LNO-CCSD(T) IEs with Tight settings exhibit about 0.05 (0.2) kcal/mol average (maximum) deviations with respect to the exact reference, while the TightPNO DLPNO-CCSD(T1) IE errors are about twice as high. If needed, more tightly converged LNO-CCSD(T) can reliably approach the local approximation free CCSD(T) results matching the about 0.01 (0.05) kcal/mol average (maximum) uncertainty of the new basis set incompleteness estimate. Due to its efficient implementation with small memory and disk I/O requirements, LNO-CCSD(T) is thus an excellent choice to generate large data-sets with highly-accurate references used routinely for the parametrization or assessment of DFT methods [141,142] and for the training step in machine learning approaches [21][22][23][143][144][145][146][147][148][149][150][151].
It is important to note that there is a large gap between the size and complexity of the S66 dimers and those of the systems that are in the applicability range of local CCSD(T) approaches as, for example, our LNO-CCSD(T) code was already shown to scale up to a 1000atom protein-ligand complex with a quadruple-ζ basis set [91]. The presented excellent performance of the local correlation approaches is partly due to the moderate system size and IE strength of the S66 dimers. However, approaching the rigorous few percent relative accuracy in the CCSD(T) IEs is already important for the most complicated but medium-sized S66 dimers. Such accuracy enables rigorous and informative comparisons with alternative methods capable of high-accuracy, including local correlation approaches, DMC, SAPT, RPA, and maybe even the leading DFT methods [25]. Moreover, as the system size grows to a very large scale, the IE has a faster than linear scaling component in terms of the size of the interacting surface and has been documented to reach up to 100 kcal/mol CCSD and 10-20 kcal/mol (T) contributions for a 132-atom fullerene-in-a-nanoring [27] and a 1023-atom protein-ligand complex [91] studied recently. In such cases, a relative error of even better than a few percent is required for well-converged IEs. Therefore, when more challenging and more extended molecular complexes are studied, we recommend cautious convergence tests and extrapolation toward the local approximation free CCSD(T) limit to estimate [91] and to avoid potentially high local correlation errors.