Appointing silver and bronze standards for noncovalent interactions: A comparison of spin-component-scaled (SCS), explicitly correlated (F12), and specialized wavefunction approaches

Noncovalent interactions (NCI) are fundamental to self-assembly, organic crystal structure, protein folding, drug binding, and many other topics of relevance for physics, chemistry, and biology.^1–3 Unfortunately, NCI can be difficult to model accurately. Quantum mechanical methods such as Hartree–Fock (HF) or standard functionals from density functional theory (DFT) can fail completely for weakly bound systems due to inadequate treatment of the electron correlation effects that give rise to long-range London dispersion interactions. Even seemingly sophisticated post-HF wavefunction methods can provide a poor description of NCI. Although second-order Møller–Plesset perturbation theory (MP2) provides a rather good description of hydrogen-bonding interactions in water⁴ or weak dispersion interactions in methane dimer,⁵ it greatly overbinds more polarizable systems such as benzene dimer.⁶ Disappointingly, the even more sophisticated coupled-cluster (CC) singles and doubles (CCSD) method⁷ tends to underbind polarizable complexes, sometimes by about as much as MP2 overbinds them.^8–10 Within the standard hierarchy of wavefunction-based quantum chemical methods, the gold standard, coupled-cluster with perturbative triple excitations [CCSD(T)],¹¹ provides consistently reliable results.^9,12 However, due to the substantial computational expense of standard CCSD(T) (which scales as o³v⁴, where o and v are the number of occupied and virtual orbitals, respectively), there is a strong motivation to find more approximate methods that can still model NCI accurately and reliably.

Numerous studies have compared various approximate quantum chemical methods for weakly bound dimers, but they often employ only a handful of test cases, or perhaps only a single standard test set, such as the S22 database of Hobza and co-workers,¹³ a collection of 22 van der Waals dimers ranging in size from (H₂O)₂ to the adenine·thymine complex. A larger number of test cases, sampling a wider variety of interaction motifs, is desirable for reliably assessing the performance of various approximate methods. Moreover, some approaches developed over the past few years are accurate enough that assessment of their performance requires the generation of even more accurate benchmarks. We have found^14,15 that a focal-point scheme^16,17 combining extrapolations to the MP2 complete-basis-set (CBS) limit using the aug-cc-pVTZ and aug-cc-pVQZ basis sets,^18,19 along with a higher-order correlation correction evaluated as the difference between CCSD(T) and MP2 in a triple-ζ basis set augmented by diffuse functions, provides quite reliable benchmarks closely approaching the CCSD(T)/CBS limit. We denote this scheme as MP2/aTQZ +

$δ MP 2 CCSD (T)$/aTZ or, more succinctly, CCSD(T)/[aTQZ; δ:aTZ]. Applying this or similar approaches, we have developed three benchmark databases for NCI: (1) NBC10,¹⁴ which contains ten non-bonded potential curves for systems like methane dimer, pyridine dimer, and various configurations of the benzene (Bz) dimer; (2) HBC6,^14,20 which provides fully relaxed dissociation curves for six doubly hydrogen-bonded complexes consisting of formic acid (FaOO), formamide (FaON), and formamidine (FaNN) monomers; and (3) HSG,^14,21 which considers intermolecular contacts for fragments taken from the crystal structure of indinavir bound to HIV-II protease. In addition, even higher quality interaction energies (IE) have recently been obtained^14,22,23 to improve the available reference values for the S22 database.¹³ The four databases considered in this work comprise 345 CCSD(T)/CBS energies. Other notable databases of similar benchmark IE quality are S66 (revised)²⁴ and A24 (Ref. 12) by Hobza and co-workers.¹³ 

For this study, we employ the revision A benchmarks for NBC10, HBC6, and HSG databases and the revision B benchmark for S22, all tabulated in Ref. 14. The bimolecular complexes in these databases are displayed in Fig. 1, which also introduces a color-coding scheme based on symmetry-adapted perturbation theory (SAPT)²⁵ analysis in which red designates electrostatically dominated complexes, typically those with hydrogen-bonding (HB), blue designates dispersion-dominated (DD) complexes, and green represents in-between or “mixed-influence” (MX) cases. Table I and the ternary diagram in Fig. 1, which places a colored dot for each system based on the relative importance of electrostatics, dispersion, and induction toward the IE, indicates that our four databases provide reasonably good coverage of possible NCI motifs.

Databases are used in this work to assess a large number of wavefunction methods, ranging in complexity from MP2 to CCSD(T), with and without incorporation of explicit correlation (F12).²⁶ Of special interest are the numerous modifications of standard methods that have been proposed to improve treatment of NCI. These include spin-component-scaling (SCS) of MP2 [SCS-MP2, SCS(N)-MP2, and SCS(MI)-MP2 (plain, nucleobase, and molecular interactions flavors, respectively)].^27–29 and CCSD [SCS-CCSD and SCS(MI)-CCSD].^8,30 We also examine the average of MP2 and MP3, dubbed MP2.5,³¹ the MP2 coupled method (MP2C),^32,33 and dispersion-weighted versions of MP2 (DW-MP2)³⁴ and CCSD(T)-F12 [DW-CCSD(T)-F12].³⁵ These methods are combined with various standard correlation-consistent basis sets^18,19 from aug-cc-pVDZ to aug-cc-pV5Z where feasible. We compare the accuracy of these approaches as well as their computational cost.

A similar study was performed by Hobza and co-workers for the S66 test set when it was introduced,¹⁰ and the performance of some of the methods surveyed was re-analyzed²⁴ when the benchmark energies in S66 were improved from MP2/aTQZ +

$δ MP 2 CCSD (T)$/aDZ to MP2/aTQZ +

$\delta _{\rm{ MP2}}^{\rm{CCSD(T)}}$

$δ MP 2 CCSD (T)$/aDTZ. Those studies focused on MP2, MP2C, SCS-MP2, SCS(MI)-MP2, DW-MP2, MP3, MP2.5, CCSD, SCS-CCSD, and SCS(MI)-CCSD in the CBS limit, and they also reported error statistics for the S66x8 and S66a8 databases, which distort the S66 molecules by changing the intermolecular distance or the intermolecular orientation, respectively (although energies for these extended databases were obtained with the somewhat less accurate

$\delta _{\rm{ MP2}}^{\rm{CCSD(T)}}$

$δ MP 2 CCSD (T)$/aDZ approach like the original S66 database). The present study is similar in spirit and differs primarily by using an alternative, larger test set at the highest quality (345 energies vs. 66); by extending the analysis to explicitly correlated (F12) methods, including DW-CCSD(T)-F12; and by including in the analysis various basis sets, pairs of basis sets for CBS extrapolation, and focal-point schemes, affording a wider array of model chemistries in the search for good Pauling points.³⁶ Additionally, we recently used these same databases to evaluate several popular density functional approximations³⁷ and various truncations of symmetry-adapted perturbation theory.³⁸ Hence, combined with our previous work, the present study allows for a head-to-head comparison of various wavefunction methods, dispersion-corrected DFT methods, and SAPT approaches for NCI. (Our previous DFT study³⁷ was published prior to our most recent revisions of the benchmark interaction energies,¹⁴ and hence the supplementary material³⁹ contains updated error statistics for the previously studied DFT methods.)

Further studies of similar character to this include the halogen-containing and halogen-bonded surveys of Hobza⁴⁰ and Martin and Kozuch,⁴¹ respectively. Additional review of S66 by Riley et al.⁴² for MP2, SCS-MP2, and SCS-S66-MP2 [S66-trained SCS(MI)-MP2] considered local and explicitly correlated variations on MP2. While the present work also deals with MP2-F12- and CC-F12-based methods and incorporates density-fitting where available, many further promising algorithmic approaches to reducing the computational cost of CCSD(T) such as local,^43,44 fragment,^45,46 limited virtual orbital space,^47,48 domain-based pair-natural orbitals,⁴⁹ and tensor hypercontraction⁵⁰ are not touched upon.

For the purposes of evaluating or parameterizing next-generation force-fields, even larger collections of high-quality benchmark values for NCI would be desirable. For a database with thousands of reference energies, or for test cases with substantially larger numbers of atoms, our currently recommended gold standard prescription for estimating the CCSD(T)/CBS limit may be too computationally expensive. Hence, a goal of the present work is to identify substitutes for the gold standard that remain reliable enough to be used for meaningful benchmarking purposes. After reviewing the accuracy of 394 model chemistries in Secs. III A–III D, we examine their relative efficiency in Sec. III E and define “silver,” “bronze,” and “pewter” standards for NCI.

Spin-component-scaling has been popular for noncovalent interactions. Originally proposed by Grimme in 2003 as a simple procedure for improving the average quality of MP2,²⁷ SCS separates the same-spin (SS) and opposite-spin (OS) components of the correlation energy (E_corr) and multiplies each by a separate scaling factor. Parameters were obtained by fitting to 51 high-quality reaction energies, but the same parameters greatly improved the IE for two configurations of Bz₂,²⁷ correcting the well-known overbinding of π-π contacts by MP2.⁶ These promising results led to initial excitement about SCS-MP2 as a general method for NCI, but this excitement was dampened after finding that SCS overcorrected and spoiled the good performance of MP2 for systems like (CH₄)₂^51,52 and led to somewhat larger errors for hydrogen-bonded complexes.⁵¹ 

In 2007, Hill and Platts²⁸ re-parameterized SCS-MP2 by fitting to the interaction energies of ten stacked nucleobases, yielding the SCS(N)-MP2 approach with significantly reduced errors for the S22 database. Around the same time, Distasio and Head-Gordon²⁹ also re-examined SCS for applications to NCI, and they formulated the “molecular interactions” MP2 method, or SCS(MI)-MP2, by fitting to S22. Unlike previous efforts, Distasio and Head-Gordon²⁹ reduced basis set sensitivity by obtaining separate scaling parameters for each basis set and extrapolation, and moreover they employed non-augmented correlation-consistent basis sets, cc-pVXZ (XZ).¹⁸ Because the parameters are tuned to the basis set, this permits dispensing with the diffuse functions that would ordinarily be important in NCI computations. Marchetti and Werner in 2009 sought to balance the advantages of MP2 for electrostatics-dominated systems and SCS-MP2 for dispersion-dominated systems by switching between the methods according to a “dispersion-weighting,” yielding DW-MP2.³⁴ These authors formulated DW-MP2 using explicitly correlated (MP2-F12) calculations, but the present work uses the label DW-MP2-F12 in that context. Switching in Eq. (1) is accomplished by means of a variable ω defined in Eq. (2) that depends on the ratio of HF and MP2 IE and on two parameters fit to the S22 database,

In this work, SCS parameters are collected into Table II, and density-fitted (DF) MP2 quantities are used when conventional values from a higher level calculation, e.g., MP3, are not available.

Using a Casimir–Polder transformation, the MP2 dispersion energy (⁠

$E disp (20)$⁠) can be rewritten as an integral over frequency-dependent response functions for the component monomers determined at the uncoupled Hartree–Fock level. The use of uncoupled response functions can lead to significant overestimations of the overall dispersion energy, particularly in π–π stacked complexes.^33,53 Hence, Hesselmann and co-workers excised the troublesome uncoupled dispersion energy component from MP2 theory (by computing it independently from SAPT) and replaced it with a more accurate TDDFT quantity to form MP2 coupled, or MP2C,^32,33

For S22, MP2C very encouragingly reduces mean absolute error (MAE) to 0.17 kcal/mol from 0.78 for MP2.³² In this work

$E disp (20) SAPT 0$ is computed using Psi4 and

$E_{\rm{disp(20)}}^{\rm{TDDFT}}$

$E disp (20) TDDFT$ is computed from the exchange-only ALDA functional using local HF orbitals with MOLPRO.⁵⁴ 

Observing the happy bracketing of stacked π-π errors by MP2 and MP3 levels of theory, Pitoňák et al.³¹ formalized the prescription in MP2.5,

By varying the relative MP2 and MP3 weightings over the S22 database and applying a mixture of basis set treatments, the authors found optimal values of 0.4–0.6 and settled upon the outright average, which yields RMSE 0.29 kcal/mol, an improvement over its ingredients and SCS-MP2 variants at a computational cost below that of CCSD.

Following the introduction of SCS-MP2, Takatani, Hohenstein, and Sherrill⁸ applied a similar independent scaling of same-spin and opposite-spin components of the doubles correlation energy in the CCSD method to define SCS-CCSD (the singles contribution is unscaled). Optimization of the scaling parameters on a set of 48 reaction energies (the SCS-MP2 set less three) led to a near halving of MAE for these systems to 1.1 kcal/mol (compare 1.8 kcal/mol for SCS-MP2). For NCI, SCS-CCSD showed low non-parallelity errors (failure to align the shapes of dissociation curves, particularly the minimum-energy intermonomer distance, between reference and examined model chemistries) for sandwich Bz₂ and (CH₄)₂, two systems that often resist simultaneous good treatment. Shortly thereafter, Hobza and co-workers re-parameterized SCS-CCSD based on the S22 database to yield SCS-CCSD for molecular interactions, SCS(MI)-CCSD.³⁰ SCS(MI)-CCSD cuts MAE to a quarter of SCS-CCSD values both for the database itself and for displacements along dissociation curves.

In addition to the standard supermolecular approaches described above, this work examines several of their explicitly correlated counterparts. Explicitly correlated methods introduce to the theoretical ansatz terms (denoted generically by F12) that depend explicitly on the distance between pairs of electrons and thereby accelerate convergence toward the complete basis set limit.^26,55,56 MP2-F12 computations were performed as implemented in MOLPRO,⁵⁴ using the default geminal β = 1.0 parameter and the 3C(FIX) ansatz which is orbital invariant, size consistent, and free of geminal basis set superposition error.^57–59 Unless figures were available as a byproduct of CCSD(T)-F12, MP2-F12 computations applied DF at all stages. MP2-F12 and CC-F12 computations performed in MOLPRO correct for basis set incompleteness in the Hartree–Fock energy by adding to it a correction for single excitations into the complementary auxiliary basis set⁶⁰ (CABS) space.⁶¹ It should be noted that the dispersion weighting approach was originally developed for explicitly correlated calculations so DW-MP2 in Ref. 34 corresponds to DW-MP2-F12 here. To form MP2-level methods DW-MP2-F12 and MP2C-F12, quantities

$E IE HF$⁠,

$E_{\rm{IE}}^{\rm{MP2}}$

$E IE MP 2$⁠, and

$E_{\rm{IE}}^{\rm{SCS\hbox{-}MP2}}$

$E IE SCS - MP 2$ in Eqs. (1)–(3) are simply replaced by

$E_{\rm{IE}}^{\rm{HF\hbox{-}CABS}}$

$E IE HF - CABS$⁠,

$E_{\rm{IE}}^{\rm{MP2-F12}}$

$E IE MP 2−F12$⁠, and

$E_{\rm{IE}}^{\rm{SCS\hbox{-}MP2\hbox{-}F12}}$

$E IE SCS - MP 2-F12$⁠, respectively.

Explicitly correlated coupled-cluster methods (CCSD-F12, SCS-CCSD-F12, SCS(MI)-CCSD-F12, and CCSD(T)-F12) have been accessed for the most commonly used ansatzë, F12a and F12b.^57,59,61,62 The triples correction is adapted for F12 computations by the approach of Werner and co-workers^34,61 through scaling by the ratio of MP2-F12 and MP2 energies:

The two asterisks in (T**) indicate that the dimer scale factor was used for each component of the interaction energy to preserve size-consistency.³⁴ This is the scheme used throughout this work instead of (T*) where dimer and monomers are scaled independently.

Recently, Marshall and Sherrill³⁵ extended dispersion weighting to coupled cluster theory, mixing CCSD(T)-F12a (performs well for hydrogen-bonding systems) and CCSD(T)-F12b (superior for dispersion-dominated complexes) in the same way that MP2 and SCS-MP2 are mixed in DW-MP2.

Equations (6) and (1) are analogous, and both share Eq. (2) (after appropriate F12 substitutions). By fitting the two parameters that determine the mixing factor to S22B,¹⁴ DW-CCSD(T**)-F12 is able to halve the MAE associated with either of its component methods.

Wavefunction methods were evaluated with the double- (aDZ), triple- (aTZ), quadruple- (aQZ), and quintuple-ζ (a5Z) basis sets of Dunning augmented by diffuse functions on all atoms^18,19 (aug-cc-pVXZ abbreviated aXZ throughout). Generally, all ζ-levels are reported for MP2-based techniques, while MP2-F12 runs through Qζ, conventional CC is confined to Dζ and Tζ, and CC-F12 uses only Dζ. Of special note, results for several CC methods [i.e., CCSD, SCS-CCSD, SCS(MI)-CCSD, and CCSD(T)] are derived from quantities computed in obtaining the original CCSD(T) benchmarking energies.¹⁴ These quantities involve mild basis set variations^63,64 but are presented under the unified labels aDZ or aTZ. All IE were computed using the counterpoise (CP) correction scheme of Boys and Bernardi,⁶⁵ unless explicitly stated otherwise. Systematic truncations of the diffuse space for Dunning basis sets (dubbed calendar⁶⁶ sets) have also been examined for MP2- and MP2-F12-based methods. Consult Table I and Sec. IIC of Ref. 38 for description and computational details.⁶⁷ 

Several of the wavefunction methods considered employ density-fitting and thus require specification of auxiliary basis sets. Generally, the defaults for Dunning basis sets suggested by MOLPRO and Psi4 were employed. For MP2, this is aug-cc-pVXZ-JKFIT⁶⁸ for Coulomb/exchange integrals in the SCF and aug-cc-pVXZ-RI⁶⁹ for fitting in the MP2 portion. Methods MP2-F12 and CCSD(T)-F12 (CC reference is not DF) additionally use cc-pVXZ-JKFIT for Fock and exchange matrices and for forming the many-electron integrals in the CABS approach; remaining integral quantities are fit with aug-cc-pVXZ-RI. MP2C uses several auxiliary basis sets with an aXZ orbital basis: aug-cc-pVXZ-JKFIT (DF roles in IE and monomer DF-HF Coulomb/exchange and monomer local HF), aug-cc-pVXZ-RI (MP2

$E disp (20)$ dispersion), and cc-pVXZ-JKFIT (monomer Coulomb/exchange). SAPT computations in Psi4 utilize aug-cc-pVXZ-JKFIT for the underlying supermolecular DF-HF and aug-cc-pVXZ-RI for density-fitting within SAPT terms.

Since correlation energy converges slower with respect to basis set size than the HF reference energy (⁠

$E total HF$⁠) underlying any correlated method, basis set extrapolation procedures have been developed⁷⁰ for wavefunction methods to treat the two independently. When the correlation energy is modeled by

$E_{\rm{corr}}^{\rm{X}\zeta } = E_{\rm{corr}}^{\rm{CBS}} + A/\ell ^{\rm{pow}}_{\rm{max}}$

$E corr Xζ=E corr CBS +A/ℓ max pow$⁠, solution for ℓ_max = 3, 4, or 5 (a DT, TQ, or Q5 extrapolation, respectively) combined with the single largest basis HF energy involved leads to the formula (TQ demonstrated)

Where adjacent ζ-level basis sets are available, this work uses two-point Helgaker⁷¹ extrapolations where pow = 3 (unless otherwise stated) or (when stated explicitly) Hill⁵⁶ extrapolations where pow varies by basis family, ζ-level, and method. Extrapolated quantities will be referred to as level-1 basis treatments and solitary basis sets as level-0.

Although correlation energy generally converges slowly with respect to basis ζ-level, it has been observed that, particularly for CCSD(T), the quantity

$δ MP 2 post - MP 2=E corr post - MP 2−E corr MP 2$ is comparatively static with basis set size (though often not monotonically convergent until aQZ^14,15). This enables level-2 basis set refinements to be defined, wherein a MP2 energy in a large or extrapolated basis is subjoined by a δ-correction evaluated in a smaller, more accessible basis (or basis set extrapolation). Just as the level-1 refinement is valid for post-HF approaches that can separate a correlation energy from the reference energy, level-2 basis treatments can be practiced on any post-MP2 technique. In this work, MP2C, MP3, MP2.5, CCSD, SCS-CCSD, SCS(MI)-CCSD, CCSD(T), MP2C-F12, and the CC-F12 variants are eligible. The following equation demonstrates the assembly of a level-2 refinement for the highest basis treatment here considered, [Q5ζ; δ:Tζ],

The notation [Aζ; δ:Bζ] is meant to succinctly denote a post-MP2 δ-correction at basis B (generally D, DT, or T), appended to a MP2 energy or extrapolation at basis A (generally DT, T, TQ, or Q5), appended to the largest single-basis reference energy encompassed by A. Equation (10) shows a rearrangement of Eq. (8) equating a MP2 energy plus small-basis “coupled-cluster correction” with a small-basis post-MP2 energy plus “basis-set correction.”

SAPT0 values for MP2C have been accessed through the quantum chemistry code Psi4.⁷² The remaining wavefunction calculations were performed within MOLPRO⁵⁴ 2009.1 and 2010.1. Methods with empirical adjustments employ the parameters in Table II. Timing results were generated on a dedicated workstation featuring a Xeon Core i7-3930K processor (6 cores, overclocked to 3.9 GHz) and 64 GB memory. The system used for comparing timings is S22-7, a 30-atom Watson–Crick-bonded adenine·thymine complex of C₁ symmetry. The dimer for MP3/aTZ required a serial run;⁷³ all others were fully threaded across six processors through Intel MKL 10.3. Timings results for MP2 (DF applied to both SCF and MP2) were extracted from Psi4, being the faster code, in preference to MOLPRO. Regrettably, comparable CC/Tζ numbers are not available since that calculation²³ employed massively parallel NWCHEM⁷⁴ on a Cray, rather than MOLPRO on a workstation as here. In contrast to the efficiency study in Ref. 37 (which also employed far inferior hardware), tabulated wall times are for all components of a supermolecular interaction energy (generally, sum of bimolecular complex and two CP-corrected monomers) and thus are directly comparable to those in Ref. 38.

The performance of a method and basis set (together, a theoretical model chemistry⁷⁵) for a database is evaluated by comparing IE (always in kcal/mol) for each individual system to that of its benchmarked value and is summarized by the quantity mean absolute error. MAE is used in this work in preference to mean absolute percent error (MA%E) though each has advantages and disadvantages for comparisons across NCI binding motifs or dissociation curves that are discussed fully in Sec. IIC of Ref. 37. Other common measures of error are compiled in the supplementary material,³⁹ including definitions, relative error counterparts to Tables III–VII, and statistics for each model chemistry and database combination. In this work, figures and statistics in the text use the “near-equilibrium” subset, including all S22 and HSG and the five points nearest equilibrium from the dissociation curves of NBC10 and HBC6, while tables use the entire 345 systems from four databases.

While summary statistics are valuable, a model chemistry can be more completely assessed by a visual representation of errors from all database members. Results for several of the best methods examined in this work are illustrated by strip charts included in Tables IV–VII, the conventions for which are outlined here. Each horizontal strip represents the results for the four databases (S22, NBC10, HBC6, and HSG) with a given model chemistry. Thin vertical lines plot the error in interaction energy (kcal/mol) for each member of the database in either the underbound (right of the zero-error line) or overbound (left of the zero-error line) sector of the chart. Individual quantities lying beyond the range of the graph are omitted without annotation. The vertical lines are colored to reflect the NCI bonding character as revealed by SAPT (generally, red: hydrogen-bonded; blue: dispersion-dominated; green: mixed). Color maps the relative dominance of dispersion vs. electrostatics contributions and is shown for each system by the rectangular outlines in Fig. 1. A black rectangular marker indicates the MAE over all the databases, with its position being fixed in the “overbound” sector of the graph for ease of comparison to the zero line, regardless of the preponderance of individual subset markings. For databases composed of potential energy curves (i.e., NBC10 and HBC6), only five points centered on the equilibrium geometry of each curve are plotted and used to compute the MAE marker position; this strategy avoids the visual confusion that would arise from the many small-error points at the tail of a curve that are not directly comparable to the minimum-energy values in other databases.

The most condensed summaries of wavefunction-based model chemistries are the MAE statistics averaged over four databases in Table III. Breakdowns by individual database and binding motif are presented in Tables IV–VII, arranged by common basis set. Illustrations of overall and binding motif error arranged by common method are collected by computational family in Figs. 2 and 3: (a) MP2-, (b) CC-, (c) MP2-F12-, and (d) CC-F12-based methods. Timings are presented in Fig. 4 and are tabulated for model chemistries in Tables IV–VII. Comparison of the very best wavefunction model chemistries examined to their competitors in SAPT and DFT is given in Fig. 5. The authors have striven to make the figures and tables informative and self-contained. Readers can examine them, then re-join the text at Sec. III E, as an abridged route through this work. The detailed computational figures underlying this study can be accessed externally, either through static tables³⁹ or an interactive website.⁷⁶ 

From HF theory which predicts strongly underbound IE with very large spread (individual errors 0.7–15.1 kcal/mol), the introduction of electron pair correlation in MP2 commences effective treatment of hydrogen-bonding NCI. As seen in Fig. 2(a.i), HB (red bars) converges strongly with basis set; aDTZ and aQZ reach tolerable error statistics (signed error mean 0.31 kcal/mol and standard deviation 0.25 kcal/mol for near-equilibrium systems; hereafter {+0.31; 0.25}) and aTQZ improves further with respect to centering about the zero-error line, {+0.04; 0.24}. In contrast, mixed-influence and dispersion-dominated systems exhibit two patterns of behavior (neither is ζ-level sensitive): (i) modest error at ⩾aTZ {−0.23; 0.29} for the less-polarizable cases, especially HSG (all) and Bz·CH₄ and (CH₄)₂ curves in NBC10 and (ii) highly overbound and large error span {≈−2; ≈0.6} and short-shifted equilibrium inter-monomer distances (R_eq) for the more-polarizable π-stacking systems abundant in S22 and NBC10.⁷⁷ The disparate performance of MP2 for MX/DD is well-documented^5,6,31 and follows from implicit use of uncoupled (HF) monomer polarizabilities. Unless pre-validated for the system and basis set in question, MP2 cannot be recommended for NCI.

Tuned scaling of same- and opposite-spin contributions in SCS-MP2 proves (Fig. 2(a.ii)) to be a trade-off between HB and certain MX/DD systems. Previous work⁵¹ showed this method to disfavor HB and alkane systems. Here, HB is wildly underbound and spread in error, converging between aDZ {+3.42; 1.22} and aTQZ {+1.60; 0.53} slowly toward an underbound mean. Also disrupted is MP2's success for less-polarizable MX/DD; though slightly reduced in spread, their mean is shifted to {≈ +0.44; ≈0.16} for ⩾aDTZ. Treatment of the more-polarizable subset is vastly improved (excepting S22-13 and 15), reaching optimal {+0.01; 0.05} at aTZ then drifting slowly overbound. NBC10 shows error bifurcation between sandwich and parallel-displaced configurations (overbound) and T-shaped and less-polarizable members (underbound), yet each dissociation curve lacks non-parallelity errors. IE errors for MX/DD in SCS-MP2 are fairly static with respect to ζ-level, allowing more-polarizable cases (including all NBC10 {+0.21; 0.19}) to be computed accurately at ⩾aTZ.

SCS(N)-MP2 improves the apparent strong performance of SCS-MP2 for MX/DD, as illustrated in Fig. 2(a.iii). For less-polarizable MX/DD, all basis treatments ({+0.12; 0.20} at aTZ) are estimated better than SCS-MP2; compared to MP2, error distributions are shifted underbound in SCS(N)-MP2 so that the latter is inferior for aDZ but is better centered for larger basis sets. Meanwhile, more-polarizable MX/DD are estimated far better than MP2 and slightly worse (except aDZ at {+0.03; 0.21}) than SCS-MP2, yielding {[−0.54, −0.37]; ≈0.20} at ⩾aTZ. These modestly positive figures are countered by non-parallelity errors in NBC10. HB is satisfactory at aDZ {+0.47; 0.35}, far better than MP2 or SCS-MP2. At larger basis sets, HB errors spread and means diverge so that for ⩾aDTZ, MP2 is superior. As the method was not fit to HB cases, such off-kilter results are not surprising. Overall, SCS(N)-MP2/aDZ is a good prospect for a cheap and accurate IE for a general NCI complex, with error profile {+0.28; 0.26} and most individual errors [−0.3, +1.0] kcal/mol, yet neither its energy curves nor its basis trends are sound.

By defining different parameters for each basis treatment, SCS(MI)-MP2 largely suppresses the rival trends among NCI binding motifs that trouble previous MP2-based methods and achieves comparable energy spreads for HB and MX/DD classes and all basis sets. Though parameterized for non-diffuse-augmented basis sets, patterns for aXZ and XZ are very similar (excepting TZ is strongly preferred over aTZ), and the latter is discussed here; both are shown in Fig. 2(a.iv). HB complexes are well-centered and exhibit moderate error spread for TZ ({+0.08; 0.46} with FaOO·FaOO and FaOO·FaNN curves being the worst offenders, >1 kcal/mol overbound), contracting slightly for larger sets. Hobza and co-workers¹⁰ also noted SCS(MI)-MP2's good HB performance over MP2 at ζ-levels where the latter is far from converged to the basis set limit. MX/DD complexes also show little basis-treatment variation, and most individual errors lie within −0.5 and +0.7 kcal/mol, {+0.10; 0.27} at TZ. As SCS(MI)-MP2 also is free of NBC10 and HBC6 non-parallelity errors, it can be generally recommended for NCI at all defined basis set treatments.

As illustrated in Fig. 2(a.v), DW-MP2 circumvents contrary HB and MX/DD ζ-trends by invoking on a per-system basis the method (MP2 in Fig. 2(a.i) or SCS-MP2 in (a.ii)) likelier to perform well. For HB, DW-MP2 achieves error distributions very like to MP2 ({+0.41; 0.26} at aDTZ) only shifted ≈0.1 underbound due to ω being not quite 1 in Eq. (1). The least- (e.g., (CH₄)₂) and most-polarizable (e.g., sandwich Bz₂) MX/DD systems face similar compromises since the former have low ω and so receive the less accurate SCS-MP2 and the latter have ω not quite 0 and so do not receive the full benefit of SCS-MP2. Remaining MX benefit considerably from method interpolation so that MX/DD overall with DW-MP2/aDTZ is {+0.08; 0.24}. The mechanism of dispersion weighting can lead to discontinuity in dissociation curves. Although statistics for NBC10 are exemplary (individual errors ±0.4 kcal/mol for ⩾aDTZ), DW-MP2 shows ≈0.1 Å lengthening of R_eq and accompanying non-parallelity errors present in neither of the parent approaches and correlated with ω shifts of [0.2, 0.5] around the R_eq geometry. Judging by error distribution, DW-MP2 can be recommended for general NCI at ⩾aDTZ, but one should be mindful that HB fares best at the costliest aTQZ treatment and MX/DD has underlying non-parallelity issues.

The MP2C method, illustrated in Fig. 2(a.vi), exhibits considerable (and non-NCI-parameterized!) improvement over MP2 for mixed and dispersion systems. Over level-0 and level-1 basis treatments, more-polarizable MX/DD (excepting outliers S22-13 and 15 which behave more like HB) exhibit exemplary error convergence and minimal spread, from overbound at aDZ {+0.71; 0.19}, to optimal at aDTZ and aQZ ({≈−0.01; 0.04}) to slightly overbound at aTQZ. Less-polarizable MX/DD are also much improved over MP2, being optimal at aTZ {+0.03; 0.09}. Expectedly, if regrettably, HB in MP2C is not improved over MP2, with the dispersion correction generally making MAE slightly better and error spread slightly worse at each basis or extrapolation. However, level-2 basis treatments recover the HB results of high-ζ MP2 at low-ζ MP2C such that [aTQZ; δ:aDZ], [aTQZ; δ:aDTZ], and [aTQZ; δ:aTZ] give error distributions {[−0.08, +0.16]; 0.24} like to the underlying MP2/aTQZ {+0.04; 0.24}. MX/DD systems also benefit by level-2 basis treatments, reaching {[−0.10, +0.00]; 0.03} for more- and {[−0.13, −0.07]; 0.09} for less-polarizable subsets (this and further MP2C discussion pertains to the three model chemistries listed for HB). MP2C exhibits many features of a promising method, including absence of non-parallelity errors over NBC10 and HBC6 dissociation curves (with the exception of aDZ, for which R_eq is shifted to longer values compared to the reference); its main lack is a large individual-error span (1.4 kcal/mol), especially for HB. Notwithstanding that one weakness, MP2C is recommended for general NCI, particularly at the most computationally efficient [aTQZ; δ:aDZ] basis treatment {+0.04; 0.22} with MAE 0.17 kcal/mol.

Further approaches for improving or expediting MP2-based computations through basis truncation or extrapolation have also been considered. Figure 3(a) is an expanded version of Fig. 2(a) that illustrates the convergence patterns of calendar basis sets from which diffuse functions have been systematically truncated. Considering only binding motif subsets where fully augmented basis sets yield MAE < 0.5 kcal/mol, calendar basis sets at light truncation largely recover the aXZ value. Across MP2, SCS-MP2, SCS(N)-MP2, and DW-MP2, MAE differs by +0.02–+0.10 for haDZ, by −0.04–+0.05 for haTZ, and by −0.02–+0.05 kcal/mol for haQZ. Truncating further, jaTZ differ by −0.16–+0.22 kcal/mol, and jaQZ and maQZ are also acceptable. SCS(MI)-MP2 in Fig. 3(a.iv) is largely steady across diffuse levels ⩾DTZ, varying by no more than ±0.04 kcal/mol for HB and less-polarizable MX/DD but by >0.5 kcal/mol for the more-polarizable subset. MP2C in Fig. 3(a.vi) shows HB convergence behavior similar to MP2 for level-0 and -1 basis sets. For the more practical level-2 [aTQZ; δ:aDZ] basis treatment, δ-corrections with haDZ and jaDZ provide nearly as good results as aDZ (subset MAE difference −0.03–+0.05 kcal/mol). Also investigated briefly were alternate extrapolation parameters for MP2 from Table VIII of Hill et al.⁵⁶ Although for DTζ these moderately improve HB statistics for MP2 and SCS-MP2 (and hence DW-MP2) at all calendar levels, there is a cost to MX/DD subsets. At TQζ, both the benefit and harm are damped, and by Q5ζ, the Helgaker and Hill schemes are equivalent.

Figure 2(b) illustrates the progression of (non-explicitly-correlated) post-MP2 wavefunction approaches from MP3 to CCSD(T), as well as the basis ζ-level convergence within each. It is apparent among the conventional methods in panels (b.i) MP3, (b.iii) CCSD, and (b.vi) CCSD(T) that level-0 basis treatments are performing no better [excepting CCSD(T)/aTZ] than their MP2 counterparts in Fig. 2(a.i), giving MAE of 1.6, 2.0, 1.3 and 0.8, 1.3, 0.5 kcal/mol for aDZ and aTZ, respectively, despite the incorporation of higher order correlation. Unsurprisingly, HB is the worst subset (MAE > 2.0 at aDZ), given that the dominating electrostatic effects are primarily addressed at MP2 which itself exhibits large errors at aDZ and aTZ then converges swiftly and monotonically thereafter. For all NCI classes, these post-MP2 methods are pervasively underbound (individual errors 0.0–6.5 kcal/mol) and accordingly predict overly long equilibrium separations for complexes in NBC10 and HBC6. A basis set correction (Eq. (10)) based upon sufficiently converged (that is, beyond level-0 aTZ) MP2 ably counteracts these deficiencies (MP2 itself overbinds IE and over-shortens R_eq). The remaining five level-1 and -2 basis treatments (appearing in Fig. 2(b)) are analyzed for both conventional methods and those designed for NCI: MP2.5, SCS-CCSD, SCS(MI)-CCSD.

Hydrogen-bonded complexes in MP3 have similar basis-set convergence patterns to MP2, reaching passable at best {[+0.08, +0.28]; 0.27}. In contrast, MX and DD systems shift from overbound with large error span (individual errors [≈0, 3] kcal/mol) for MP2 to underbound with large error span for MP3 (particularly parallel-displaced Bz₂ curves). All level-1 and -2 treatments perform similarly, yielding {+0.29; 0.20} for less-polarizable and {+1.46; 0.38} for more-polarizable subsets ([aTQZ; δ:aTZ]; S22-13 and 15 excluded to compare with MP2.5, for which they are outliers). The remarkable error balance between MP2 and MP3 for dispersion-dominated complexes inspired the creation of their average, MP2.5, shown in Fig. 2(b.ii). MX/DD are indeed substantially corrected, showing better centering for the less-polarizable {−0.04; 0.09} and far less spread for the more-polarizable {−0.46; 0.08}, albeit firmly overbound. HB also improves, especially in error span, with the best now yielding {[+0.06, +0.16]; 0.13}. Due to the spread in MX/DD errors, MP3 cannot be recommended in general for NCI, but MP2.5 offers good (and nearly identical) error profiles for all level-2 treatments built upon MP2/aTQZ. By both computational scaling and statistics, [aTQZ; δ:aDZ] is the choice for MP2.5, with overall {−0.04; 0.18} and MAE 0.18 kcal/mol; Hobza and co-workers¹⁰ concur in recommending this same basis treatment.

CCSD amplifies the defects of MP3, with all bonding motifs showing wide spreads of underbinding errors. No basis treatment alleviates the acute underbinding, and the “best” model chemistries [aTQZ; δ:aDZ] and [aTQZ; δ:aTZ] have overall error ≈{+0.8; 0.4}. SCS-CCSD displays better basis set convergence than CCSD for HB along aDZ to aTZ to aDTZ, but even the best [aTQZ; δ:aDTZ] at {+0.33; 0.16} is considerably underbound and spread for so costly a computation. In contrast, SCS-CCSD is excellent for MX/DD, especially for less-polarizable complexes, with all subset MAE < 0.25 kcal/mol for all databases. Although SCS-CCSD/[aTQZ; δ:aDZ] performs so well for MX/DD at {+0.03; 0.09}, its treatment of HB leads to overall {+0.40; 0.29} and MAE 0.42 kcal/mol.

By (roughly) exchanging the p_OS and p_SS parameters from SCS-CCSD, SCS(MI)-CCSD retains the former's good performance for MX/DD while largely ameliorating errors for HB. For the four level-2 basis treatments considered, all MX and DD subset MAE are < 0.15 kcal/mol and are centered, condensed, and show good non-parallelity. Less-polarizable systems are especially good at {[−0.04, −0.02]; [0.03, 0.05]}, while for more-polarizable MX/DD, double-ζ treatments are slightly better ({−0.02; 0.03} at [aTQZ; δ:aDZ]) than triple-ζ ({+0.16; 0.04} at [aTQZ; δ:aTZ]). HB in SCS(MI)-CCSD is more rapidly convergent with respect to basis set than SCS-CCSD, and errors are subdued by [aTQZ; δ:aDZ] at {+0.25; 0.11}, though still underbound, and outstanding by [aTQZ; δ:aTZ] at {+0.01; 0.05}. Hobza and co-workers¹⁰ also observed the balanced treatment of NCI motifs by SCS(MI)-CCSD. Although the importance of connected triples excitations for dispersion is clearly shown by the improvement of CCSD(T) over CCSD (MAE 0.15 vs. 0.99 kcal/mol at aDTZ), SCS(MI)-CCSD contrives to imitate the basis set convergence patterns of CCSD(T) very closely for all bonding motifs (compare Figs. 2(b.v) and 2(b.vi)). Insofar as the chemical application can tolerate the modest HB error of CCSD(T)/[aTQZ; δ:aDZ], SCS(MI)-CCSD at the same basis treatment (overall MAE 0.15 kcal/mol) can be recommended as a cheaper method (avoiding triples computation and triple-ζ CC). At [aTQZ; δ:aTZ], though the difference between SCS(MI)-CCSD and CCSD(T) (essentially reference quality) is greater (aforementioned underbinding of more-polarizable MX/DD), SCS(MI)-CCSD is recommended without qualification with overall MAE 0.06 kcal/mol.

Explicitly correlated MP2 at all basis ζ-levels is very similar to high-ζ conventional MP2 (compare Figs. 2(a.i) and 2(c.i)), with individual database member errors differing by only 0.05 kcal/mol on average between MP2-F12/aDZ and MP2/aTQZ or aQ5Z. For MX/DD, all five MP2-F12 basis treatments and MP2/aTQZ share error profiles for less-polarizable of {[−0.38, −0.32]; [0.31, 0.34]} and more-polarizable of {[−2.26, −2.21]; [0.62, 0.64]}. Though errors of overbinding and spread are more extreme for the latter subset, nearly all NBC10 curves show unduly short R_eq. Hydrogen-bonding IE, too, are constant for aDZ through aTQZ (in contrast to conventional MP2) at <2% MA%E. Although HB errors are well-centered and mostly within the ±0.3 kcal/mol range, some FaNN systems are considerably (approaching 1 kcal/mol) overbound. Given the remarkable uniformity of error distributions in MP2-F12 model chemistries, selecting a basis beyond aDZ appears unwarranted. Like MP2, the explicitly-correlated method should not be used for mixed-influence or dispersion-dominated complexes due to significant overbinding.

Application of spin-component scaling to MP2-F12 at any basis ζ-level unsurprisingly reproduces the high-ζ conventional analog, as exhibited in Fig. 2(c.ii–v) (excepting (c.iv), vide infra). Accordingly, SCS-MP2-F12 inherits most of the undesirable properties of SCS-MP2, and though greatly improved over MP2-F12 for more-polarizable MX/DD cases at the cost of significantly larger HB errors, it is not recommended for general NCI (MAE 0.92 kcal/mol). SCS(N)-MP2-F12 fares no better since, by emulating SCS(N)-MP2/aTQZ, it shows good results for less-polarizable MX/DD and bad outliers for HB (approaching −3 kcal/mol). Thus, SCS(N)-MP2-F12/aDZ (MAE 0.53 kcal/mol) is inferior to its non-explicitly correlated counterpart, SCS(N)-MP2/aDZ (MAE 0.33 kcal/mol). SCS(MI)-MP2-F12 employs different scaling parameters at each ζ-level so Fig. 2(c.iv) alone of the MP2-F12-based approaches shows pronounced basis-set trends. Since the parameters were designed for their non-explicitly correlated analogs, no single SCS(MI)-MP2-F12 model chemistry (MAE 0.41 kcal/mol for aTQZ) performs well for all binding motifs and databases, nor any better than than SCS(MI)-MP2/TZ or QZ (MAE 0.32 kcal/mol). Finally, DW-MP2-F12 combines good performances for HB by MP2-F12 and for MX/DD by SCS-MP2-F12 with convergence at aDZ (the parent DW-MP2 requires at least aDTZ) to yield quite small errors ({+0.07; 0.28}; MAE 0.28 kcal/mol) relative to computational complexity. However, like its non-F12 analog, DW-MP2-F12 sometimes exhibits R_eq shortening and non-parallelity errors associated with the dispersion-weighting mechanism, calling into question its suitability for potential energy curves.

The tantalizing potential to combine within MP2C-F12 the performance of MP2C for MX/DD subsets and the ζ-convergence-acceleration properties of MP2-F12 for the HB subset is realized in Fig. 2(c.vi). Among hydrogen-bonded systems, both conventional and coupled MP2-F12 at all basis treatments exhibit the same ≈1.3 kcal/mol range of individual errors that, though broad, are low on a relative scale: MA%E ≈1.3%. MP2-F12 (and accordingly MP2C-F12) at aDZ treats HB systems as well as any (rigorously ab initio) model chemistry short of CCSD(T). Less-polarizable MX/DD complexes show excellent centering, little error spread, and no non-parallelity problems. All basis treatments are encompassed by {[−0.14, −0.02]; [0.07, 0.09]}, with larger basis sets tending slightly toward the overbound so that, happily, best is aDZ at {−0.02; 0.07} and MAE 0.06 kcal/mol. Once again excepting outliers S22-13 and 15, MP2C-F12 performance for more-polarizable MX/DD is outstanding. All basis treatments are within {[−0.14, +0.05]; [0.03, 0.06]}, with [aDTZ; δ:aDZ] and [aTQZ; δ:aDZ] being especially good and aDZ at {+0.05; 0.06} and aTZ being scarcely worse. Overall, MX/DD (including S22-13 and 15) are handled slightly better by model chemistries built around MP2C computations at aDZ than at aTZ. After seeking to maximize performance at MP2C's DD strength (and seeking to minimize the number of computations needed to assemble an interaction energy, MP2C's weakness), MP2C-F12/aDZ is recommended for general NCI for all but the most challenging electrostatics-dominated systems (e.g., double hydrogen bonds and ionic complexes), as it supports MAE 0.27 kcal/mol and MA%E 1.7% for HB and 0.08 kcal/mol and 3.9% for MX/DD.

Also surveyed briefly are suggestions from the literature for improving MP2-F12/aXZ quantities. Figure 3(c) is an expansion of Fig. 2(c) showing convergence with respect to both ζ-level and degree of truncation of the diffuse space through calendar basis sets. Within a ζ-level, depleting diffuse functions generally shifts IE errors toward underbinding. These shifts are sizeable for the smallest basis sets DZ and jaDZ, but haDZ, jaTZ, and haTZ differ from their fully augmented counterparts by MAE ⩽ 0.03 kcal/mol for all subsets in Fig. 3(c.i–v) (exception being more-polarizable MX/DD at triple-ζ where larger MAE differences favor the truncated basis sets). The cc-pVDZ-F12 and cc-pVTZ-F12 basis sets developed⁷⁸ for explicitly correlated computations also were tested (on S22) and found to shift errors by ≲0.1 kcal/mol toward underbound compared to their Dunning counterparts, which, given the tendency of MP2-F12, SCS(N)-MP2-F12, SCS(MI)-MP2-F12, and DW-MP2-F12 to be overbound, confers a MAE advantage for these methods of 0.03–0.08 over aDZ and 0.01–0.04 kcal/mol over aTZ. Finally, basis set extrapolation exponents of Hill et. al tailored to the MP2-F12(3C-FIX) ansatz (Table IX of Ref. 56) were applied to S22 but found to vary inappreciably (MAE subset differences ⩽0.03 kcal/mol) from the traditional Helgaker extrapolation for aTQZ, aDTZ, and haDTZ, whereas jaDTZ and DTZ slump by up to 0.15 kcal/mol. None of calendar basis, F12-basis, or Hill extrapolation boost the performance of MP2C-F12/aDZ; indeed haDZ or cc-pVDZ-F12 degrades its S22 overall MAE by 0.04 or 0.09 kcal/mol, respectively. Figure 3(c.vi) does show that calendar basis sets in level-2 treatments can maintain the performance quality of fully augmented counterparts.

Like its conventional counterpart, CCSD-F12a/b is consistently underbound at the double-ζ level, averaging [0.43, 1.73] kcal/mol for all subsets and basis treatments. Indeed, from inspection of Figs. 2(b.iii) and 2(d.i), it is unsurprising that error distribution of the best conventional model chemistry, CCSD/[aTQZ; δ:aDZ] at {+0.78; 0.38}, is comparable to the best explicitly correlated (level-0) scheme, CCSD-F12a/aDZ at {+0.85; 0.40}. For all binding motifs, CCSD-F12b has average errors ≈[0.1, 0.3] kcal/mol larger than CCSD-F12a, consistent with reports⁶¹ that the F12 “a” variant is superior for small basis sets. Focal-point (level-2) basis treatments show no improvement beyond CCSD-F12/aDZ.

SCS-CCSD-F12a in Fig. 2(d.ii lower) also mirrors its conventional counterpart in that spin-component-scaling has a far greater effect upon mixed and dispersion-dominated systems than on HB, reducing errors by ≈0.7 kcal/mol to {+0.01; 0.07} for MX/DD, there being no difference in treatment of more- vs. less-polarizable systems, and mending the over-long NBC10 R_eq given by CCSD-F12. HB errors decrease by ≈0.5 kcal/mol to {+0.43; 0.20} at aDZ— still underbound and with non-trivial error spread. SCS-CCSD-F12b follows the same bonding-motif patterns as “a,” but subset MAE are all larger (Fig. 2(d.ii upper)), and more-polarizable MX/DD are underbound, {+0.38; 0.13}. Again, neither variant derives appreciable benefit from level-2 basis treatments.

SCS(MI)-CCSD-F12 in Fig. 2(d.iii) is the first CC-F12 method to yield excellent results for both HB and MX/DD complexes. For the “a” variant, all individual errors are within ±0.2 kcal/mol, more- and less-polarizable MX/DD subsets are treated evenly, and curves are free of non-parallelity errors. The three binding motifs each shift slightly (≈0.05 kcal/mol) toward stronger interactions between aDZ ({+0.03; 0.05} for HB and {−0.02; 0.03} for MX/DD) and the level-2 basis treatments. For both HB and more-polarizable MX/DD, SCS(MI)-CCSD-F12a achieves MA%E < 1%, while less-polarizable is 3.9%. Compared to SCS-, SCS(MI)-CCSD-F12b has a similar error profile for MX/DD {+0.14; 0.12} and far improved for HB {+0.22; 0.08}. While free of non-parallelity errors, “b” is generally underbound, especially for HBC6 and sandwich and parallel-displaced NBC10 structures. Less-polarizable MX/DD systems {+0.07; 0.06} are more accurate than more-polarizable {+0.34; 0.07}. Although subsets show somewhat reduced error with level-2 basis treatments, aDZ basis provides the best error and efficiency balance, just as for the F12a variant.

The incorporation of triples excitations into CCSD-F12 has a salutary effect on all classes of interactions, reducing MAE from 0.94, 0.49, and 1.45 kcal/mol for HB, less-, and more-polarizable MX/DD subsets in CCSD-F12a (Fig. 2(d.i)) to 0.07, 0.05, and 0.31 kcal/mol in CCSD(T**)-F12a (Fig. 2(d.iv)). Hydrogen-bonding systems fare better with the “a” variant, having an error distribution of {+0.03; 0.06} at aDZ (similar magnitude only overbound for the level-2 basis treatments). More-polarizable MX/DD systems exhibit the opposite trend with respect to F12 variant. Whereas F12a is consistently overbound, {−0.31; 0.08} at aDZ, with some S22 individual errors reaching −0.6 kcal/mol, F12b performs better than any previously discussed CC-F12 approach, with {−0.03; 0.04} at aDZ. Less-polarizable MX/DD complexes are treated equally well at aDZ by F12a {−0.04; 0.04} and by F12b {+0.05; 0.04}. However, since level-2 model chemistries shift individual errors toward overbinding, these systems are harmed by F12a, yet computed exceedingly well by F12b. Although still showing appreciable underbinding for HB (up to +0.55 kcal/mol), CCSD(T**)-F12b/[aDTZ; δ:aDZ] is the best model chemistry of this technique (overall MAE 0.09 kcal/mol, MA%E 1.3%), exhibiting (in common with the entirety of the method) no non-parallelity difficulties and being actually very similar in error distribution to CCSD(T)/[aTQZ; δ:aDZ].

The DW-CCSD(T)-F12 method takes advantage of CCSD(T)-F12a/b bracketing IE errors between centered (F12a) and underbound (F12b) for HB and centered (F12b) and overbound (F12a) for MX/DD. Its success is apparent in Fig. 2(d), where the best performances of the two variants in Fig. 2(d.iv) are united in Fig. 2(d.v). For hydrogen-bonding, DW-CCSD(T**)-F12 matches the excellent results from parent CCSD(T**)-F12a (itself the best result for this binding motif short of triple-ζ CCSD(T)), showing {+0.04; 0.06} for aDZ and all database subset MAE < 0.10 kcal/mol (<0.6% MA%E). Mixed and dispersion-dominated systems are treated with comparable accuracy by DW-CCSD(T**)-F12/aDZ at {+0.01; 0.05} with all database subset MAE < 0.07 kcal/mol (<3.5% MA%E). Across all near-equilibrium systems, individual errors range from [−0.12, +0.30] kcal/mol. Applying level-2 basis treatments tightens the error spread across most subsets, but the changes are ⩽0.03 kcal/mol MAE and <1% MA%E and so are not sufficient to warrant the extra assembly required. Dissociation curves for NBC10 and HBC6 show occasional discontinuities in the dispersion weighting parameter at the tails or walls, but these are sufficiently minor that they are not apparent in actual energy traces. DW-CCSD(T**)-F12/aDZ is commensurate with the other standout among CC-F12 methods, SCS(MI)-CCSD-F12a/aDZ, having binding motif subsets which vary by MAE ⩽ 0.05 kcal/mol and <2% MA%E. Both are recommended for general NCI applications.

To gauge computational efficiency for the wavefunction methods considered here, timings for a single system of appreciable size (adenine·thymine) are represented in Fig. 4 and Tables IV–VII. It is readily apparent that costly computations are by no means guarantors of accurate results. Instead, model chemistries that present a good balance of accuracy and expense (Pauling points) are sought.

With the highly efficient and scalable DF-MP2 codes existing today, most any basis set is accessible for system sizes under review in this work [see the flat MP2 curves in Fig. 4 (pink traces) showing that DZ to aQZ requires 1–70 min]. Thus, MP2-based model chemistries can be selected primarily by error profiles, rather than computational cost. Though increasing basis set size generally returns better error statistics within method for the inexpensive MP2-based techniques in Fig. 2(a.i–v), acceptable errors profiles over all NCI motifs are best represented by: SCS(N)-MP2/aDZ (overall MAE: 0.33 kcal/mol; overall MA%E: 20.2%; time for A·T: 3 min), SCS(MI)-MP2/QZ (0.32; 16.5%; 21 min), DW-MP2/aDTZ (0.33; 13.6%; 16 min), and SCS(MI)-MP2/TZ (0.36; 17.3%; 4 min). Among these, a winner (albeit not clear-cut) can be selected for a pewter standard (connoting everyday or non-benchmark use) of SCS(MI)-MP2/QZ, which has the advantages of having been constructed taking all classes of NCI into account, requiring only two straightforwardly defined parameters, and supporting no non-parallelity issues. The cost of explicitly correlating MP2 (maroon traces in Fig. 4) is 20–90 × between (respectively) aQZ and aDZ, so that the best MP2-F12-based model chemistries, SCS(MI)-MP2-F12/DTZ (0.39; 15.7%; 9.1 h) and DW-MP2-F12/aDZ (0.28; 13.9%; 4.1 h), are not sufficiently improved over their MP2 counterparts (albeit without refitting SCS parameters) to justify the extra expense.

Rigorously ab initio wavefunction methods less complete than CCSD(T) (i.e., MP2, MP3, and CCSD) are generally poor for NCI; unsurprising since high-quality treatments of both electrostatics and correlation are required. Moreover, even the cleverest tweaks, scalings, averagings, and contortions to low-scaling MP2 only deliver methods comparable to B3LYP-D in quality (see Table III). Thus, though this work supplies recommendations for everyday use, or pewter standard (for readers of Ref. 38, this role is comparable to sSAPT0/jaDZ, zeroth-order exchange-scaled SAPT with a high-angular-momentum-truncated aDZ basis), the important role of wavefunction theory for NCI is to provide benchmark-quality interaction energies, a role unassailable by present DFT approaches.

An immediate candidate for benchmark standard is MP2C in Figs. 2(a.vi) and 2(c.vi) (salmon traces in Fig. 4). Although requiring considerable time for large basis sets, aDZ to aQZ entailing 1.5–46 h, MP2C converges quickly for MX/DD systems with respect to basis ζ-level and so performs very well at aDZ with a basis set correction provided by MP2/aTQZ. By building MP2C upon MP2-F12, one can dispense with the focal-point scheme and use aug-cc-pVDZ outright. Selecting between MP2C/[aTQZ; δ:aDZ] (0.17; 4.2%; 2.9 h) and MP2C-F12/aDZ (0.16; 3.1%; 5.5 h) error-wise is indistinguishable, timings-wise is comparable, and may come down to manual-labor-wise: eight input files versus five. Both are recommended as our bronze standard for benchmarking NCI. Its weakness is in spread of errors for electrostatics dominated systems; work in this group has shown that for ionic NCI complexes (monomer charges +/− or −/−), MP2C-F12/aDZ individual errors can approach 1 kcal/mol.

Computational cost increases substantially for methods of greater complexity than MP2; although MP3 (green traces in Fig. 4), CCSD (blue), and CCSD(T) (purple) all fit within the week time frame at aDZ, only MP3 is expected to do so at aTZ. The proximity of light and dark pairs of blue and purple markers at aDZ shows that F12 only incurs a modest additional cost for CCSD (5.5 h) and CCSD(T) (7.3 h). Among the model chemistries showing the lowest errors in Figs. 2(b) and 2(d) are MP2.5/[aTQZ; δ:aDZ] (0.18; 7.1%; 3.8 h), SCS(MI)-CCSD/[aTQZ; δ:aDZ] (0.15; 2.7%; 35.8 h), SCS(MI)-CCSD-F12a/aDZ (0.05; 2.3%; 39.9 h), SCS(MI)-CCSD/[aTQZ; δ:aTZ] (0.06; 2.9%; comparable timing not available), and DW-CCSD(T**)-F12/aDZ (0.06; 1.8%; 97.3 h). Of these, MP2.5 offers no improvement over bronze (and shows slight non-parallelity errors in NBC10), double-ζ SCS(MI)-CCSD, too, is no better than bronze while being appreciably costlier, and triple-ζ SCS(MI)-CCSD, though accurate, is off-the-scale in cost. The remaining methods, SCS(MI)-CCSD-F12a and DW-CCSD(T**)-F12, are highly accurate and differ most noticeably in cost. The authors regard mere doubling (at A·T system size) of computational effort as a worthy exchange for the incorporation of more physics [(T) term] in benchmarking and appoint DW-CCSD(T**)-F12/aDZ as our silver standard for benchmarking NCI, with SCS(MI)-CCSD-F12a being a runner-up that might become preferable for computational efficiency in larger systems.

Future investigations on NCI will involve systems larger in size or larger in number than those accessible by gold standard (⩾CCSD(T)/[aTQZ; δ:haTZ]) levels of theory. While DFT is an invaluable tool, its growing collection of Jacob's ladders^79,80 cannot yet reach systematic convergence, and benchmarks of reliable or excellent quality continue to make recourse to wavefunction theory. Continuing our efforts to provide directly comparable assessments for computational approaches to NCI, this work adds 394 wavefunction model chemistries to existing DFT³⁷ and SAPT³⁸ publications. In particular, here are examined MP2 and associated spin-component-scaling, dispersion-weighted, coupled, and explicitly correlated variants, CCSD and associated spin-component-scaling and explicitly correlated variants, and CCSD(T) and associated dispersion-weighted and explicitly correlated variants. All methods are evaluated with basis sets augmented by diffuse functions and with Helgaker extrapolations thereof; post-MP2 methods are assessed with focal-point schemes built upon MP2[-F12]. Additionally, MP2- and MP2-F12-based techniques are investigated with diffuse-function-truncated (calendar) basis sets. Altogether, several items of interest and Pauling points for routine use and benchmarking are extracted.

• For levels of theory on the order of MP2, accuracy is hit-or-miss between tolerable and poor. (Overall MAE in kcal/mol and any liabilities in parentheses.)

  • For hydrogen-bonded systems, judging by MAE values, hits include MP2/⩾aDTZ (0.21–0.37), SCS(N)-MP2/aDZ (0.51; maxE >1 underbound), SCS(MI)-MP2/⩾TZ (0.33–0.39; maxE >1 overbound), and DW-MP2/⩾aDTZ (0.23–0.42).

  • For mixed-influence and dispersion-dominated systems, judging by MAE values, hits (0.27–0.32) include SCS(N)-MP2/⩾aDZ (long R_eq for aDZ) SCS(MI)-MP2/⩾TZ, and DW-MP2/⩾aTZ (long R_eq, maxE >2 overbound for close contacts).

  • For overall NCI, judging by MAE values, hits (0.27–0.36) include SCS(N)-MP2/aDZ (long R_eq, maxE >1 underbound), SCS(MI)-MP2/⩾TZ (maxE >1 overbound), DW-MP2/⩾aDTZ (long R_eq, maxE >2 overbound for close contacts).

  • Clearly, no single MP2-based method is consistently sound. SCS(MI)-MP2/QZ is recommended as a wavefunction technique for everyday applications since it, among the model chemistries with tolerable error quantities, is expected to disappoint only in accuracy rather than by more invidious means (e.g., non-parallelity errors or classes of NCI outliers). In statistics, TT (overall): MAE 0.32 kcal/mol (MA%E 16.5%); HB (hydrogen-bonding): 0.37 (2.7%); MX/DD (mixed-influence and dispersion-dominated): 0.30 (24.5%); A·T (timing for adenine·thymine): 0.3 h.

• Levels of theory extending MP2 can provide considerably more accuracy than MP2 variants at modestly higher cost but also considerably worse accuracy or enormously higher cost.

  • The averaging of MP2 and MP3 provides even performance across binding motifs and is quite affordable with application of a post-MP2 δ-correction, as in MP2.5/[aTQZ; δ:aDZ]. In statistics, TT: 0.18 (7.1%); HB: 0.18 (1.2%); MX/DD: 0.21 (10.3%); A·T: 3.8 h.

  • Comparable to MP2.5 in overall error, MP2C performs especially well for dispersion-dominated cases. It is weaker for electrostatically dominated systems, especially charged ones. Atop cheap conventional MP2, a δ-correction elicits excellent performance with MP2C/[aTQZ; δ:aDZ]. In statistics, TT: 0.17 (4.2%); HB: 0.25 (1.6%); MX/DD: 0.10 (5.7%); A·T: 2.9 h. Atop explicitly correlated MP2, a small basis can be used outright with MP2C-F12/aDZ. In statistics, TT: 0.16 (3.1%); HB: 0.27 (1.7%); MX/DD: 0.08 (3.9%); A·T: 5.5 h. Both are recommended as a “bronze standard” for benchmarking, being perhaps only an order of magnitude more expensive than everyday DFT and wavefunction approaches for systems of the size considered here.

  • As MP2.5 scales as non-iterative O(N⁶) and MP2C is presently dominated^33,81 by the iterative O(N⁴) with high prefactor local HF and TDDFT portion, the balance in computational cost seen for A·T may tip in favor of the latter for larger systems.

• CC levels of theory are unabashedly poor for all NCI motifs until triple excitations are included at high-ζ basis sets. The “gold standard” CCSD(T)/[aTQZ; δ:aTZ] used here for reference interaction energies achieves accuracy through including triples and employing CC δ-corrections to mimic large basis set results.

  • By applying explicit correlation to elicit accurate quantities at low-ζ and by mixing ansatzë with different convergence properties (F12a for HB and F12b for MX/DD), DW-CCSD(T**)-F12/aDZ is recommended as a benchmarking “silver standard” for delivering near golden accuracy at considerably lower cost. In statistics, TT: 0.06 (1.8%); HB: 0.07 (0.5%); MX/DD: 0.04 (2.4%); A·T: 97.3 h.

  • The important effect of triples excitations can be imitated, at least for systems of moderate size as surveyed in this work, by a parameterized scaling of singlet and triplet components of the CCSD correlation energy. Though bypassing the triples cost, the best level of theory with conventional CC, SCS(MI)-CCSD/[aTQZ; δ:aTZ], yet requires a large basis set to obtain silver quality results. In statistics, TT: 0.06 (2.9%); HB: 0.05 (0.3%); MX/DD: 0.07 (4.3%); A·T: comparable timing not available. With explicit correlation, the same accuracy is achieved more efficiently through SCS(MI)-CCSD-F12a/aDZ. In statistics, TT: 0.05 (2.3%); HB: 0.08 (0.5%); MX/DD: 0.04 (3.4%); A·T: 39.9 h.

The authors are grateful to Robert M. Parrish who accelerated the DF-MP2 code in PSI4 to such a degree over course of this project that our recommendations had to adapt to suit. This work was performed under the auspices of grants provided by the United States National Science Foundation (NSF) (Grant No. CHE-1300497). The Center for Computational Molecular Science and Technology is funded through a NSF CRIF award (Grant No. CHE-0946869) and by Georgia Institute of Technology.