Understanding Discrepancies of Wavefunction Theories
for Large Molecules

Abstract

Quantum mechanical many-electron calculations can predict properties of atoms, molecules and even complex materials. The employed computational methods play a quintessential role in many scientifically and technologically relevant research fields. However, a question of paramount importance is whether approximations aimed at reducing the computational complexity for solving the many-electron Schrödinger equation, are accurate enough. Here, we investigate recently reported discrepancies of noncovalent interaction energies for large molecules predicted by two of the most widely-trusted many-electron theories: diffusion quantum Monte Carlo and coupled-cluster theory. We are able to unequivocally pin down the source of the puzzling discrepancies and present modifications to widely-used coupled-cluster methods needed for more accurate noncovalent interaction energies of large molecules on the hundred-atom scale. This enhances the reliability of predictions from quantum mechanical many-electron theories across a wide range of critical applications, including drug design, catalysis, and the innovation of new functional materials, such as those for renewable energy technologies.

1 Introduction

The prediction of electronic transition energies for a single hydrogen atom in 1926 marks the beginning of an incredible successful era of quantum mechanics [1, 2]. Shortly after this breakthrough, Paul Dirac famously noted that the underlying physical laws necessary for much of physics and all of chemistry are now completely known. However, he also pointed out that the exact application of these laws leads to equations that are much too complicated to solve [3]. This paradigm of theoretical chemistry and physics is prevalent until today. In particular, the exponential growth of the computational complexity of the many-electron problem with system size still makes an exact solution of the electronic Schrödinger equation for more than a few atoms impossible. As a consequence, a hierarchy of increasingly accurate methods that is capable of producing reference results at the expense of tractable yet high computational cost has emerged. These reference results are pivotal in order to develop, assess, and further improve computationally more efficient but in general less accurate approximations. In this context, a prime example was the numerical prediction of highly accurate ground state energies for the uniform electron gas using diffusion Monte Carlo methods, leveraging the development of approximate exchange-correlation functionals that ultimately led to the breakthrough of density functional theory in computational materials science during the last decades [4, 5, 6].

At present, quantum mechanical many-electron calculations of systems containing more than 100 atoms have become possible thanks to methodological developments and considerable growth in computing power. These methodological improvements are often based on taking advantage of the relative short-rangedness of many-electron correlation effects  [7, 8, 9]. In this manner, the scaling of the computational complexity with respect to system size can be lowered. However, recently several works [10, 11, 12, 13] showed that there exist alarming discrepancies between predicted interaction energies for large molecules when using two of the most widely-trusted highly accurate many-electron theories: DMC and CCSD(T), which stand for diffusion Monte Carlo and coupled-cluster theory using single, double, and perturbative triple particle-hole excitation operators, respectively. These observations are a source of great concern in the electronic structure theory community because, in the case of noncovalent interactions between molecular complexes, both CCSD(T) and DMC are considered highly reliable benchmark methods [14, 15]. Furthermore, the observed discrepancies are large enough to cause qualitative differences in calculated properties of materials, which can have scientific, technological, and even clinical implications. For example, accurate crystal structure predictions are crucial in drug design to differentiate between harmful and effective polymorphs [16, 17, 18]. Similarly, reliable reference methods are essential for discovering and designing new functional materials for applications such as renewable energy storage and conversion, including catalysis, or solar cells [19, 20, 21]. As machine learning increasingly pervades all areas of computational first-principles physics, the accuracy of these reference methods, which provide the training data, becomes even more critical [22, 23, 24].

In the following, we analyze a set of large molecular systems where large discrepancies between approximated versions of DMC and CCSD(T) were observed [10, 11]. Importantly, a direct experimental measurement of the computed interaction energies of these systems is complicated and prone to significant uncertainties. Therefore, it is an open challenge to identify the origin of the observed deviations for the employed highly accurate yet approximate theoretical approaches. A question of paramount importance is, are errors introduced by the fixed-node approximation in DMC and the truncated cluster ansatz of CCSD(T) causing these discrepancies? Or, are approximations used to reduce the scaling of the computational complexity not accurate enough to make reliable predictions?

2 Results

We present an approach to CCSD(T) calculations of large molecules which builds upon our previous work. This approach enables us to unambiguously test if the employed approximations for DMC and CCSD(T) cause the puzzling discrepancies between their predictions. In particular, our methodology exhibits three striking advantages. Firstly, due to its efficient and massive computational parallelization, we omit any local correlation approximation, as was employed for the CCSD(T) calculations in Refs. [11, 12, 13]. Secondly, we use a plane wave basis set to enable an unbiased assessment of the quality of previously employed tabulated atom-centered basis functions. Thirdly, we are able to study the influence of higher-order contributions to the many-electron perturbation expansion beyond CCSD(T) theory for large molecular complexes.

2.1 Assessing the representation of the wavefunction

In practice, every calculation employing a wavefunction ansatz expands the true wavefunction in a finite basis. Tabulated atom-centered localized Gaussians are the most widely-used basis set for the expansion of many-electron wavefunctions in quantum chemical calculations of large molecules. They achieve a transferable and well-balanced description of electronic correlation effects for a wide range of systems. However, despite their computational advantages, they suffer from one major drawback. As they approach completeness by including more and more diffuse functions, needed to capture polarization effects, their near-linear dependence increases, especially for densely packed structures. Although there exist techniques to ameliorate this problem [25, 26, 27], it remains difficult to employ large basis sets for large systems as done in this work [26].

In contrast to localized basis sets, plane waves form a systematically improvable orthonormal basis and are typically employed in combination with finite simulation cells and periodic boundary conditions. Therefore, computed interaction energies need to be converged with respect to both cell size and plane wave basis set size. Furthermore, compared to the number of Gaussians, the required number of plane wave basis functions is significantly larger, making the latter approach computationally more expensive. To remedy the computational bottlenecks associated with plane wave basis calculations of molecules, we have developed efficient techniques that transform the large plane wave basis set to a system-specific and significantly more compact as well as systematically improvable natural orbital basis [28, 29]. Although the computational cost to determine these transformations is significant, it remains small compared to the computational cost of CCSD(T) calculations for large systems. A detailed describtion of our plane wave approach can be found in the supplementary information.

We now seek to demonstrate that our plane wave basis approach works reliable for the study of noncovalent interactions between molecules and combines the best of two worlds: compactness and systematic improvability without linear dependencies. To this end we turn to the discussion of the computed interaction energy of the parallel displaced benzene dimer on the level of CCSD(T) theory and compare to results from basis set converged Gaussian calculations (see supplementary information).

(a)  

(b)  

Figure 1: Convergence behavior of the interaction energy for the benzene parallel displaced dimer system. $\Delta E$ is the difference between our plane wave based approach and the reference results from Gaussian basis calculations extrapolated to the CBS limit (see supplementary information). a, shows the convergence with respect to the number of natural orbitals for a fixed volume of about $4018\,\text{\AA}^{3}$. b, displays the convergence with respect to the volume of the simulation cell with a fixed number of natural orbitals of $N_{\mathrm{v}}/N_{\mathrm{o}}=15$.

The total CCSD(T) energy is composed of three terms, the HF total energy, the CCSD correlation energy and the perturbative triples contribution which we denote as (T). Fig. 1(a) depicts the convergence of the CCSD and (T) correlation energy contributions to the computed CCSD(T) interaction energy for a fixed box size with respect to the number of basis functions (natural orbitals) per occupied state ($N_{\mathrm{v}}/N_{\mathrm{o}}$). We include a recently introduced correction to accelerate the convergence of correlation energies to the complete basis set limit (CBS) [30]. Our findings show that a basis set size of $N_{\mathrm{v}}/N_{\mathrm{o}}=15$ suffices to achieve convergence to within a fraction of $0.1\,\text{kcal/mol}$. We employ this basis set to compute the CCSD(T) interaction energies and its Hartree–Fock (HF), CCSD and (T) correlation energy contributions for different simulation cell sizes. Fig. 1(b) shows that these contributions converge rapidly. Our fully converged estimate of the CCSD(T) interaction energy for the parallel displaced benzene dimer is $-2.62\,\text{kcal/mol}$, which is in excellent agreement with results obtained using Gaussian basis sets of $-2.70\,\text{kcal/mol}$.

Table 1: Interaction energy in kcal/mol of C2C2PD obtained at different levels of theory including MP2, CCSD, CCSD(T), CCSD(cT) and DMC. The stochastic uncertainty of the DMC estimates identifies a 95% confidence interval and is given in parenthesis notation. The uncertainty of this work’s are dominated by the remaining basis set error and the uncertainty of the box size extrapolation. The uncertainty of the referenced LNO and DLPNO results are estimated from the tightness parameter controlling the local approximation.

Theory	Interaction energy	Ref.
MP2	-38.5 $\pm$ 0.5	this work
CCSD	-13.4 $\pm$ 0.5	this work
CCSD(T)	-21.1 $\pm$ 0.5	this work
LNO-CCSD(T)	-20.6 $\pm$ 0.6	Ref. [11]
DLPNO-CCSD(T0)	-20.9 $\pm$ 0.4	Ref. [13]
DMC	-18.1(8)	Ref. [11]
DMC	-17.5(14)	Ref. [10]
CCSD(cT)	-19.3 $\pm$ 0.5	this work

Having established a reliable and systematically improvable framework to achieve CBS limit estimates of CCSD(T) interaction energies, we now turn to the discussion of results for the parallel displaced coronene dimer (C2C2PD) also studied in Refs. [11, 10, 13]. In this case, significant discrepancies between CCSD(T) and DMC of about $1\,\text{kcal/mol}$ after subtracting error bars were reported. The convergence with respect to box size and number of natural orbitals is presented in detail in the supplementary information. Table 1 summarizes the computed interaction energies on different levels of theory. We find that our canonical CCSD(T) estimate agrees to within $0.5\,\text{kcal/mol}$ with domain based local pair-natural orbital (DLPNO-CCSD(T₀)) and local natural orbital (LNO-CCSD(T)) findings from Refs. [11, 13]. Therefore, basis set incompleteness errors cannot be the cause for the discrepancy between DMC and CCSD(T) in this case. Likewise, local approximations used for LNO and DLPNO approaches can also be ruled out. However, it is noteworthy that the $-21\,\text{kcal/mol}$ CCSD(T) interaction energy contains a large (T) contribution of about $-8\,\text{kcal/mol}$, indicating that the correct treatment of triple particle-hole excitation effects for the electronic correlation plays a crucial role.

2.2 All that glitters is not gold: overcorrelation of CCSD(T)

The most important question of the present work is: what is the origin of the discrepancy between CCSD(T) and DMC molecular interaction energies? Having ruled out errors from local approximations and incomplete basis sets for the parallel displaced coronene dimer (C2C2PD), we now seek to assess the (T) approximation, which makes a significant contribution to the interaction energy of C2C2PD. In passing we anticipate that the (T) contribution to the interaction energy is also relatively large for all other systems with a significant discrepancy reported in Ref. [11] (see supplementary information).

The (T) approximation was introduced in the seminal work by Raghavachari et. al. [31]. Since then, it has become one of the most widely-used benchmark methods –sometimes referred to as the ’gold standard’ of molecular quantum chemistry– for systems that can be well approximated using single reference based approaches. However, we argue that the partly significantly too strong interaction energies in CCSD(T) theory are caused by the employed truncation of the perturbation series in the approximation of the triple particle-hole excitation operator. As already demonstrated by Nguyen et. al., even weak intermolecular interactions cannot be approximated in a systematically improvable manner by truncated high-order perturbation theories for systems with large polarizabilities [32]. This is also reflected in the too strong interaction energies on the level of MP2 theory, as can be observed for C2C2PD in Table 1. In the extreme case of an infinite polarizability, as it occurs in metallic systems, MP2 and CCSD(T) even yield divergent correlation energies in the thermodynamic limit, which is referred to as infrared catastrophe [33, 34]. In contrast, perturbation theories based on the resummation of certain terms to infinite order can yield interaction energies with an accuracy that is less dependent on the polarizability. Prominent examples for such approaches include the CCSD theory as well as the random-phase approximation. We have recently presented a method, denoted as CCSD(cT), that averts the infrared catastrophe of CCSD(T) by including selected higher-order terms in the triples amplitude approximation without significantly increasing the computational complexity [34].

Understanding the discrepancy

For the present work it is important to note that the main difference between CCSD(cT) and CCSD(T) theory originates from the employed approximation to the triple particle-hole excitation amplitudes. The triple amplitudes of the (cT) approximation are given in diagrammatic and algebraic form by [34]

(1)

where $\hat{V}$ and $\hat{T}_{2}$ stand for the Coulomb interaction and the double particle-hole excitation operator, respectively. For brevity, the contributions from the single excitation operator are not included and only one additional ‘direct’ diagram is depicted. In here, $\Delta^{ijk}_{abc}=\varepsilon_{i}+\varepsilon_{j}+\varepsilon_{k}-\varepsilon_{a}-\varepsilon_{b}-\varepsilon_{c}$, with $\varepsilon$’s being one-electron HF energies. The bra- and ket-states correspond to a triple excited and reference state, respectively. The (T) approximation disregards the term $[[\hat{V},\hat{T}_{2}],\hat{T}_{2}]$, which is included in (cT) and also occurs in full CCSDT theory. This term effectively screens the bare Coulomb interaction of the $[\hat{V},\hat{T}_{2}]$ term and has an opposite sign, making it crucial for systems with large polarizability. As discussed by Nguyen et. al. [32], truncated finite-order perturbation theories can exhibit large errors in noncovalent interactions due to incomplete screening. However, for small and weakly polarizible systems the $[[\hat{V},\hat{T}_{2}],\hat{T}_{2}]$ contribution is small, making the (T) and (cT) approximation agree, as it was already shown for a small set of molecules [34].

We now demonstrate that using CCSD(cT) instead of CCSD(T) theory restores excellent agreement for noncovalent interaction energies with DMC findings. First, we consider again the coronene dimer. Table 1 shows that the binding energy for the coronene dimer calculated on the level of CCSD(cT) theory is by almost $2\,\text{kcal/mol}$ closer to the DMC estimate compared to CCSD(T) theory, achieving chemical accuracy ($1\,\text{kcal/mol}$) in comparison to DMC after subtracting error bars. Next, we investigate the accuracy of CCSD(T) and CCSD(cT) compared to DMC for noncovalent interactions in smaller molecules. To this end, we study a set of dimers containing up to 24 atoms that were also investigated in Ref. [11]. This gives us another opportunity to assess the effect of local approximations at the level of CCSD(T) theory.

Figure 2: Deviations of coupled cluster results from DMC results for a set of noncovalently bound dimers with up to 24 atoms. The CCSD(T) and CCSD(cT) values are CBS estimates obtained from basis set extrapolation using aug-cc-pVTZ and aug-cc-pVQZ basis sets [35, 36] (details in the supplementary information). The LNO-CCSD(T) and DMC data are from Ref. [11]. The uncertainty measures are as described in Table 1. The uncertainty of DMC is shown by the blue area.

Fig. 2 depicts the deviations of all computed interaction energies from DMC reference values taken from Ref. [11]. It should be noted that DMC references and differences to LNO-CCSD(T) interaction energies are shown with error bars derived from statistical noise and the estimated local approximation error, respectively [11]. Using our massive computational parallelization approach, we are able to add canonical CCSD(T) interaction energies extrapolated to the CBS limit to the comparison to DMC. We stress that for these relatively small molecules, the employed atom-centered Gaussian basis sets together with a standard extrapolation procedure achieve excellent agreement with CBS limit estimates as verified on the level of MP2 theory and summarized in the supplementary information. Importantly, our canonical CCSD(T) results are in good agreement with LNO-CCSD(T) findings to within its error bars. The only minor exception is observed for the parallel displaced uracil dimer, where canonical CCSD(T) predicts a slightly stronger interaction. A comparison to DMC reveals that CCSD(T) theory predicts on average about $0.3\,\text{kcal/mol}$ stronger interaction energies. Based on LNO-CCSD(T) and DMC data alone such a statement cannot be made due to the relatively large and mostly overlapping error bars. However, our well converged canonical CCSD(T) findings allow drawing such conclusions. Only for the T-shaped pyridine and benzene dimers, DMC and CCSD(T) binding energies agree to within the DMC errors. Note that these systems have a smaller (T) contribution to the intereaction energy, compared to the parallel displaced systems. All other systems exhibit small but significant discrepancies between CCSD(T) and DMC results, which is consistent with the even larger discrepancies reported for the larger molecules in Ref. [11]. Similar to our findings for the coronene dimer reported in Table 1, Fig. 2 shows that CCSD(cT) interaction energies agree significantly better with DMC values than their CCSD(T) counterparts.

Figure 3: Comparison between the full triples and the perturbative triples approaches, (cT) and (T) for a set of molecules contained in the S22 data set [37]. a, Comparing the total triples correlation energy $E$ of the molecular systems. b, Comparing the triples contribution to the interaction energy $E_{\text{int}}$ of the dimers. The dashed line guides the eye to estimate the deviations from the full T result.

Given the good agreement between DMC and CCSD(cT) for the systems studied above, an important question to ask is if CCSD(cT) is really more accurate than CCSD(T) for noncovalent interaction energies? To answer this question, following the paradigm of theoretical quantum chemistry, we turn to more accurate reference methods. CC methods form a rapidly convergent series of approximations to the true ground state wavefunction and energy, especially for weakly interacting systems [38]. This implies that the following hierarchy applies to absolute energy differences: CCSD-HF $>$ CCSDT-CCSD $>$ CCSDTQ-CCSDT. Therefore, we first seek to assess if CCSD(cT) is a better approximation to CCSDT than CCSD(T) theory. Our massive computational parallelization approach allows us to compute CCSDT correlation energies for a set of molecules contained in the S22 set employing cc-pVDZ basis sets [35], which suffices for the following qualitative discussion. The respective CCSDT and CCSD(T) correlation energies are depicted in Fig. 3a. Our findings show that CCSD(T) and CCSDT energies are in good agreement even in the case of large molecular complexes. However, in contrast to total correlation energies, a comparison between CCSDT and CCSD(T) interaction energies of the same molecular complexes draws a completely different picture. Fig. 3b shows that the larger the triples contribution to the intermolecular interaction energy gets in magnitude, the larger CCSD(T) overestimates compared to CCSDT. We suggest the following conclusion: CCSD(T) and CCSDT total correlation energies agree fortuitously well, which fails in the case of physically more relevant strong interaction energies. The opposite applies to CCSD(cT), which agrees less good with CCSDT for absolute correlation energies, but very well in the case of interaction energies depicted in Figs. 3a & 3b, respectively.

2.3 Estimating the overcorrelation of (T) for weak interactions

In summary, Table 1 and Figs.2 & 3 demonstrate that CCSD(cT) theory achieves excellent agreement for noncovalent interaction energies between molecular complexes compared to DMC and CCSDT theory. However, we stress that the CC series of methods (CCSD, CCSDT and CCSDTQ) is observed to yield monotonic and rapidly converging interaction energies for small and weakly bound complexes [38]. Based on this knowledge, we emphasize that the Q contribution to the interaction energies can be expected to be smaller than its T counterpart, but could possibly yield a significant contribution. Indeed, this is part of the reason for the success of the CCSD(T) approximation for very small molecules, where CCSD(T) is often fortuitously closer to CCSDTQ than CCSDT [38]. Here, we argue that this error cancellation no longer functions in the case of large molecular complexes involving strongly polarizable systems such as C2C2PD, C3GC and C60@[6]CPPA. A similar problem is known to occur in MP2 theory, where the truncation of the perturbation series also leads to significantly too strong interaction energies for systems with large polarizability, although MP2 yields relatively accurate interaction energies for systems with an intermediate polarizability [32].

Figure 4: Correlation between the ratio of (T) and (cT) with the ratio of the MP2 and CCSD correlation energy contributions to interaction energies of a set of dispersion-dominated complexes from the S22, L7 and S66 benchmark datasets [37, 39, 40]. Selected cases are labeled and visualized: Methane dimer, GCGC (guanine-cytosine tetramer), BeBePD (benzene dimer parallel displaced), C2C2PD (coronene dimer parallel displaced). All the data is available in the supplementary information.

To quantify and support the statements above, Fig. 4 illustrates that there exists a correlation between the ratio of (T) and (cT) with the ratio of the MP2 and CCSD correlation energy contributions to the interaction energies of all studied molecules in this work with dispersion-dominated interactions. This demonstrates that (T) exhibits a tendency to overestimate the absolute binding energy in a similar manner as MP2 theory for more polarizable systems. Although an overestimation of the (T) binding energy contribution compared to its (cT) counterpart by about 10% might yield a fortuitously better agreement between CCSD(T) and CCSDTQ, we argue that 20–30% overestimation is expected to yield significantly too strong interaction energies. For example, the values of (T)/(cT) for the Benzene-Benzene PD and coronene dimer are 1.2 and 1.3, respectively.

Having demonstrated and explained the reasons for the overestimation of absolute interaction energies on the level of CCSD(T) theory for small molecules with up to 24 atoms and the C2C2PD system, we now want to turn to the discussion of the remaining large molecular complexes where significant absolute discrepancies between DMC and CCSD(T) have been observed. These systems include C3GC from the L7 data set and the C60@[6]CPPA buckyball-ring. Here, substantial differences in the binding energies of $2.2\,\text{kcal/mol}$ and $7.6\,\text{kcal/mol}$ were reported, respectively after subtracting error bars. Although CCSD(cT) calculations for systems of that size are currently not feasible using our approach, we now introduce a simplified model that allows us to estimate the change in interaction energies from CCSD(T) to CCSD(cT) in an approximate manner. Given the linear trend between the different correlation energy contributions to the interaction energy depicted in Fig. 4, it is possible to estimate the (cT) contribution for systems where only MP2, CCSD and (T) are known. These numbers can be calculated using a computationally efficient LNO-CCSD(T) implementation [41]. Results computed in this manner are denoted as CCSD(cT)-fit. Details about this procedure and the corresponding error estimates are provided in the supplementary information.

Table 2: Comparison of the interaction energy for large molecular complexes in kcal/mol as calculated by different levels of theory. Showcasing partially large discrepancies between CCSD(T) and DMC on the one hand, and an excellent agreement between CCSD(cT) and DMC results for complexes up to the 100-atom scale on the other hand. CCSD(T) and CCSD(cT) results are obtained using our plane wave approach. The calculation and the uncertainty of CCSD(cT)-fit is explained in the supplementary information.

System	CCSD(T)	LNO-CCSD(T) [11]	CCSD(cT)	CCSD(cT)-fit	DMC [11]	DMC [10]
GGG	-1.5 $\pm$ 0.5	-2.1 $\pm$ 0.2	-1.2 $\pm$ 0.5	-1.8 $\pm$ 0.2	-1.5(6)	-2.0(8)
GCGC	-13.1 $\pm$ 0.5	-13.6 $\pm$ 0.4	-12.5 $\pm$ 0.5	-12.8 $\pm$ 0.5	-12.4(8)	-10.6(12)
C2C2PD	-21.1 $\pm$ 0.5	-20.6 $\pm$ 0.6	-19.3 $\pm$ 0.5	-18.9 $\pm$ 0.7	-18.1(8)	-17.5(14)
C3A		-16.5 $\pm$ 0.8		-15.3 $\pm$ 0.9	-15.0(10)	-16.6(18)
PHE		-25.4 $\pm$ 0.2		-25.0 $\pm$ 0.2	-26.5(13)	-24.9(12)
C3GC		-28.7 $\pm$ 1.0		-26.7 $\pm$ 1.1	-24.2(13)	-25.1(18)
$\text{C}_{60}$@[6]CPPA		-41.7 $\pm$ 1.7		-35.6 $\pm$ 2.0	-31.1(14)

Table 2 gives our estimated CCSD(cT) interaction energies in comparison to CCSD(T) and DMC findings for seven large molecular complexes. A comparison between CCSD(cT)-fit and the explicitly calculated CCSD(cT) results for GGG, GCGC and C2C2PD shows that the linear regression model is sufficiently reliable for the systems studied in this work. For comparison Table 2 also summarizes the DMC interaction energies from Refs. [10] and [11], which agree to within at least $1\,\text{kcal/mol}$ for GGG, C2C2PD, PHE and C3GC. For the remaining systems the DMC estimates show a larger discrepancy and for C60@[6]CPPA only one DMC estimate is available. Although the DMC binding energies have overlapping error bars, the remaining uncertainties are relatively large, illustrating that obtaining highly accurate interaction energies for these large molecules is also challenging for DMC.

As already discussed in Ref. [11], CCSD(T) interaction energies listed Table 2 exhibit large discrepancies compared to DMC for C2C2PD, C3GC and C60@[6]CPPA. In contrast, CCSD(cT)-fit resolves these discrepancies for all systems on the hundred-atom scale, achieving excellent agreement with DMC estimates of Hamdani et. al. [11] to within chemical accuracy ($1\,\text{kcal/mol}$) after subtracting the error bars. Even for C60@[6]CPPA, which contains 132 atoms, a discrepancy of only $1.1\,\text{kcal/mol}$ remains, although the error bar of CCSD(cT)-fit is relatively large in this case. We argue that the remaining discrepancies are potentially caused by uncertainties in DMC, CCSD(cT)-fit and the underlying LNO-CCSD(T) calculations. It should be noted that the error bars of LNO-CCSD(T) interaction energies are in some cases underestimated, as exemplified for the Uracil-Uracil PD dimer by the comparison between canonical CCSD(T) and LNO-CCSD(T) interaction energies shown in Fig. 2. Furthermore, the DMC interaction energy of C60@[6]CPPA has not yet been verfied independently using a different DMC implementation as it was done for all other systems listed in Table 2. We also stress that in some cases the differences between the DMC estimates are larger than their respective error bars.

3 Conclusion

Our work unequivocally demonstrates that, due to the employed truncation of the many-body perturbation series expansion, one of the most widely-used and accurate quantum chemistry approaches – CCSD(T) theory – in certain cases binds noncovalently interacting large molecular complexes too strongly. Our findings show that a simple yet efficient modification denoted as CCSD(cT) remedies these shortcomings. This paves the way for highly reliable benchmark calculations of large molecular complexes on the hundred-atom scale that play a crucial role in scientific and technological problems, for example, drug design and surface science. We stress that the more accurate CCSD(cT) approximation can directly be transferred to computationally efficient low-scaling and local correlation approaches, which will substantially advance the applications of theoretical chemistry as well as physics in all areas of computational materials science where highly accurate benchmark results are urgently needed. We are witnessing an unremitting expansion of the frontiers of accurate electronic structure theories to ever larger systems which when combined with machine-learning techniques, has the potential to transform the paradigm of modern computational materials science.

• Funding T.S. acknowledges support from the Austrian Science Fund (FWF) [DOI: 10.55776/ESP335]. The computational results presented have been largely achieved using the Vienna Scientific Cluster (VSC). Andreas Irmler and Alejandro Gallo acknowledge support from the European Union’s Horizon 2020 research and innovation program under Grant Agreement No. 951786 (The NOMAD CoE). Support and funding from the European Research Council (ERC) (Grant Agreement No. 101087184) is gratefully acknowledged.

• Conflict of interest The authors declare no competing interests.

• Ethics approval Not applicable

• Consent to participate Not applicable

• Consent for publication Not applicable

• Availability of data and materials The data that support the findings of this study are available from the corresponding authors upon reasonable request.

• Code availability The VASP code is a copyrighted software and can be obtained from its official website. The CC4S code is free and open-source and can be downloaded from its website.

• Authors’ contributions A.I. and T.S. performed the calculations. A.I. and A.G. implemented the methods in Cc4s and the employed quantum chemistry packages. A.G. advised and implemented related interface methods in the VASP code. All authors wrote the paper.