Assessing the Applicability of the Geometric Counterpoise Correction in B2PLYP/Double-ζ Calculations for Thermochemistry, Kinetics, and Noncovalent Interactions*

Keywords: density functional theory, double-hybrid density functionals, noncovalent interactions, thermochemistry, London dispersion, atomic-orbital basis sets, computational chemistry, theoretical chemistry.

When using finite atomic orbital (AO) basis sets in an electronic structure calculation, the computed properties are subject to the basis set error.^[1–4] The basis set error is classified into two types: the basis set superposition error (BSSE) and the basis set incompleteness error (BSIE).^[1] Our present study focuses solely on the BSSE, for which several corrections have been proposed.^[1,5–12] BSSE arises from an unbalanced basis set expansion of supramolecular complexes and their constituent monomers. In a dimer complex with subsystems A and B, the unoccupied basis functions of subsystem A can be used by subsystem B (and vice-versa) to lower its energy, which leads to the over-binding of complexes (Fig. 1). That artificial stabilisation can be removed by the counterpoise correction (CP) scheme, which was first introduced by Boys and Bernardi (the BB-CP scheme).^[13] In this approach, the energy of each fragment is calculated with the basis functions of the full complex, but without electrons and nuclei of the other fragments; orbitals located at the positions of the other fragments are therefore also known as ‘ghost orbitals’. For example, the BB-CP correction for a dimer can be presented as:

Fig. 1.  Schematic representation of the basis set superposition error (BSSE).

Here, a and b are the basis sets of subsystems A and B in their frozen complex geometries, with energies E(A) and E(B), respectively; the dimer’s basis set is represented by ab. The CP-corrected binding energy (BE) can then be obtained as:

Similarly, the BB-CP correction can be applied to multimer systems beyond dimers.

The BB-CP scheme has two major shortcomings. The standard scheme can only calculate the intermolecular BSSE for noncovalently bound complexes. This is because it requires partitioning a system into all possible fragments. For some simple systems, this partitioning is straightforward, but for systems involving intramolecular noncovalent interactions (NCIs), this is not the case. The second problem is the computational cost, which increases with the number of fragments due to a larger number of separate calculations.

In 2006, Galao and Alvarez-Idaboy addressed intramolecular BSSE by introducing an ‘atom-by-atom’ scheme.^[14] Jensen later generalised this into the ACP-n approach (atomic counterpoise method) where the total BSSE is calculated as the sum of atomic contributions.^[15] A few empirical models for BSSE calculations—models that do not require additional quantum-chemical calculations, but rely on fitted statistical models—were introduced in Refs [16 and 17]. DiLabio and co-workers’ atom-centred potentials^[18–20] compensate for basis set limitations in treating NCIs with dispersion-corrected Hartree-Fock (HF) and Density Functional Theory (DFT) methods, the latter of which we here refer to as density functional approximations (DFAs).

In 2012, Kruse and Grimme introduced a semi-empirical correction scheme called ‘geometrical counterpoise correction (gCP)’.^[21] The gCP scheme is based on an atom pairwise approach. This scheme only requires geometric information as an input and no wave function data is needed. gCP has gained considerable attention in the computational chemistry community because of its efficiency, ease of use, and transferability. The gCP correction has to be parametrised for each basis set. For a given basis set, there is one set of parameters for HF and another set that applies to all GGA, meta-GGA, and hybrid DFT methods. In addition, the correction has also been used in specialised cases, such as HF-3c,^[22] HSE-3c,^[23] PBEh-3c,^[24] B97-3c,^[25] HFsol-3c,^[26] PBEsol0-3c,^[26] HSEsol-3c,^[26] and r²SCAN-3c.^[27] Those are DFT methods with specialised basis sets and other corrections to allow fast and accurate calculations. Regardless of the underlying DFT methods, the gCP correction has to be combined with a London-dispersion correction. Usually, this has been the DFT-D3^[28,29] or DFT-D4^[30,31] corrections.

In 2017, Head-Gordon and co-workers introduced the DFT-C^[32] correction for BSSE, which is different from Kruse and Grimme’s gCP approach in a few critical areas; see Ref. [32] DFT-C for more details on those differences. Although DFT-C showed improvement over gCP, this approach is very specialised and implemented in fewer software packages, therefore, beyond the scope of this study.

DFT currently offers the highest cost-accuracy ratio for thermochemistry, kinetics, and NCIs. Among the various DFAs, double hybrid density functionals (DHDFs) stand out for their accuracy, reliability, and overall robustness.^[33–40] DHDFs are a combination of conventional DFT and second-order Møller-Plesset perturbation (MP2) theory. It is known that conventional wave function electron correlation methods converge more slowly to the complete basis set (CBS) limit than HF or conventional DFT. Basis set convergence and BSSE for lower-rung DFT methods and double hybrids have been extensively studied for the S66x8^[41] benchmark set for noncovalent interactions in Ref. [42]. While the basis set convergence of a DHDF calculation is faster than for wave function electron correlation methods,^[43] there is a slower basis set convergence compared to conventional DFT due to the MP2 part. In electronic structure calculations, using a smaller basis set (e.g. a double-ζ AO basis) may become necessary for the treatment of large systems due to limited computing resources. Explicitly correlated approaches, such as R12 and F12 variants,^[44–48] allow using smaller AO basis sets with faster convergence to the CBS limit, and this has also been demonstrated for DHDFs.^[43] The other alternative to allow for faster DHDF calculations is a computationally efficient split-valence basis set named ‘DH-SVPD’,^[49–51] which has been especially designed for DHDFs and the treatment of weakly bound systems.

We are motivated by the question if the gCP scheme is beneficial for double-hybrid DFT calculations, as they have so far been excluded from the gCP scheme, most likely due to MP2 having a different basis set convergence behaviour. Herein, we present a brief proof-of-concept study to see if any developments in such a direction are justified. For that reason we choose the B2PLYP^[33] functional to test the impact of gCP. In the first step of our study, we analyse the major shortcomings of using B2PLYP with a small basis set (def2-SVP^[52]), namely, BSSE and missing London dispersion. We will then present the B2PLYP-gCP-D3(BJ)^[29]/def2-SVP scheme and compare it with B2PLYP-D3(BJ)/def2-SVP, B2PLYP/def2-SVP, and B2PLYP-D3(BJ)/quadruple-ζ results. The individual contributions from gCP and DFT-D3(BJ) will be analysed in all cases. This analysis will be carried out on the S66x8 set, followed by the GMTKN55^[35] database for general main-group thermochemistry, kinetics, and NCIs.

This manuscript illustrates the benefits of using DHDF-gCP-D3(BJ)/double-ζ over DHDF/double-ζ and DHDF-D3(BJ)/double-ζ for electronic structure calculations. Our conclusions will reveal useful insights for new approaches, and enable potential future developments beyond B2PLYP.

Our starting point for determining improved electronic energies is a standard DHDF calculation (B2PLYP in this study) with a small basis set. Two terms are added to the B2PLYP energy (E^{B2PLY P}) to account for BSSE and missing London dispersion forces. The total energy for the resulting B2PLYP-gCP-D3(BJ) method can be written as:

where and are the energy contributions from the gCP and DFT-D3(BJ) corrections, respectively, which will be described in more detail in the following two sections.

All electronic structure calculations were carried out with ORCA 4.0.0 and 4.0.1,^[53,54] while the DFT-D3(BJ) and gCP energy corrections were obtained with the standalone programs by Grimme and co-workers.^[55,56] The generally observed trends for B2PLYP were compared to BLYP^[57–59] and B3LYP,^[60,61] which are GGA and hybrid DFAs, respectively, with the same underlying DFT exchange and correlation components as B2PLYP. All BLYP calculations were treated with the resolution of the identity (RI-J) approximation for Coulomb integrals and appropriate auxiliary basis sets,^[62,63] while hybrid, double-hybrid, and MP2-type calculations were done with the combined resolution-of-the-identity for Coulomb and the chain-of-sphere approximations for exchange integrals (RIJCOSX).^[64] We also used the frozen core approximation for the second-order perturbative steps of the DHDF and MP2-type calculations to prevent core-core BSSE,^[65] as well as the RI approximation with appropriate auxiliary basis sets to speed up those steps.^[63] All calculations were carried out with ORCA’s quadrature grid ‘3’, followed by a non-iterative step with the larger grid ‘4’. The SCF convergence criterion was set to 10⁻⁷ E_h. The non-hydrogen atoms in the GMTKN55 sets G21EA,^[66,67] AHB21,^[68] and IL16^[68] were augmented with diffuse s and p functions from the Dunning class of basis sets,^[69] while hydrogens were subjected to diffuse s functions. Oxygens in WATER27^[35,70] were also augmented with diffuse s and p functions. While the gCP correction in our case has not been fitted for diffuse functions, it was necessary to add them in the case of the aforementioned four benchmark sets to make our results commensurate with our previous GMTKN55 studies and to properly treat the negatively charged species^{[35,38,71–73]} in those sets. Similar to our previous GMTKN55 papers, we use the prefix aug’ in front of a basis set name. For instance, “aug’-def2-SVP” indicates that diffuse functions were used for a specific benchmark set. If we have a more general discussion “(aug’)-def2-SVP” indicates that some results are based on basis sets with diffuse and others without diffuse functions.

The basis functions for core electrons of heavy elements were replaced with def2-ECP^[52] effective core potentials in the GMTKN55 sets HAL59,^[74,75] HEAVY28,^[28] and HEAVYSB11.^[28]

The gCP correction is defined as:^[21]

where represents the difference in atomic energies between a large basis set and the target basis set for free atom A, S_AB is the Slater-type overlap integral assessed over s-type orbitals, represents the number of virtual orbitals on atom B, and σ, α and β are scaling parameters. Slater exponents in S_AB are weighted by a scaling parameter dubbed η. Those four scaling parameters need to be fitted for a given level of theory; in our case this meant that we fitted them for the B2PLYP/def2-SVP level of theory. The fitting procedure for the gCP correction uses the S66x8 benchmark set, which is an extension of the S66 dataset.^[41] S66 contains 66 noncovalently bound dimers and is separated into systems that are dominated by electrostatics (hydrogen bonds), dispersion, or a mixture of both. To allow for an analysis of interaction energies in non-equilibrium situations, the S66x8 extension considers eight different intermolecular distances for each dimer, namely at 0.90r_e, 0.95r_e, r_e, 1.05r_e, 1.10r_e, 1.25r_e, 1.5r_e, and 2r_e, where r_e is the dimer-specific equilibrium distance. This means that S66x8 contains 528 dimers whose geometries we took from Ref. [41]. For our target level of theory (B2PLYP/def2-SVP), the standard BB-CP corrections were calculated for all 528 dimers. The four gCP scaling parameters were determined in a least-squares fit against those BB-CP values. For the dimers with the shortest intermolecular distances (i.e. 0.90r_e), the weights of the deviations from that correction were reduced to 0.5 during our fit, as recommended in the original gCP publication.^[21]

In passing, we would also like to mention that we additionally analysed a damped version of gCP, which was first used in the context of PBEh-3c.^[24] The damped-gCP variant uses the basic form of gCP, but the short-range part of the correction is damped according to:

where are radii introduced in the DFT-D3(0) correction with zero-damping.^[28] During initial tests, we saw the statistical superiority of gCP over damped-gCP (see Table S4 in the Supplementary Material). Therefore, we decided to keep the basic form of gCP and the damped version will not be discussed any further.

The two-body London dispersion energy according to the DFT-D3(BJ) scheme is defined as:^[29]

with chemical-environment-dependent^[28] dispersion coefficients C_n (n = 6, 8), internuclear distances R_AB, order-dependent scaling factors s_n, cut-off radii , and empirical parameters a₁ and a₂ that control the magnitude of added dispersion at short and medium interatomic distances. For more details on this correction, see Refs [29, 76, and 77]. As we will outline in the results section, we observed benefits of refitting DFT-D3(BJ) for our purposes over using the standard parametrisation, which is based on B2PLYP/def2-QZVP. In this section we cover briefly how we performed such refits.

The s₆ parameter is set to unity for most functionals to ensure the correct asymptotic behaviour of the dispersion energy. The exception are double hybrids, which have an additional long-range dispersion contribution from their MP2 part. Herein, we determined s₆ for the B2PLYP/def2-SVP level of theory according to the procedure presented in Ref. [78]. For that purpose, we considered the Ne₂, Ar₂, and Kr₂ dimers at interatomic separations of 7 Å, 10 Å, and 10 Å, respectively. The MP2 contributions to the interaction energies of the three dimers were calculated and compared to CCSD(T)^[79]/aug-cc-pVTZ^[69,80] dispersion energies. This allowed us to calculate the average amount of long-range dispersion recovered by B2PLYP/def2-SVP. This average was then subtracted from unity to get the value s₆. Our resulting s₆ is 0.830.

Having determined s₆, the remaining three parameters s₈, a₁, and a₂ in Eqn 6 were determined by a least-squares fit against the reference interaction energies of the combined S66x8, S22x5,^[81] and NCIBLIND^[82] benchmark sets. We carried out such fits for different BSSE-corrected B2PLYP/def2-SVP interaction energies, as outlined in the results section.

To further assess the performance of B2PLYP/def2-SVP with the DFT-D3(BJ) and gCP corrections, we used the GMTKN55 database, which comprises 2462 total single-point energies distributed over 55 datasets. GMTKN55 is classified into five categories. The first category, ‘basic properties and reactions of small systems’, consists of 18 sets and covers reaction energies for small systems, total atomisation energies, ionisation potentials, electron affinities, and self-interaction-error-related problems. The second category offers nine test sets for isomerisation and reaction energies of larger systems. Barrier height related problems are dealt with in the third category (seven test sets in total). The fourth category comprises 12 benchmark sets dealing with intermolecular noncovalent interactions, while intramolecular interactions are described by the nine benchmark sets in the fifth category. A summary of the GMTKN55 database and the names of all individual 55 benchmark sets can be found in Table 1.

Table 1.  Overview of the GMTKN55 database and its five categories: basic properties and reactions of small systems (A), reaction energies of larger systems and isomerisation (B), barrier heights (C), intermolecular noncovalent interactions (D), and intramolecular noncovalent interactions (E) (for more details, see Ref. [35])

First, we analyse the behaviour of BSSE for the S66x8 benchmark set. For def2-SVP, the average BB-CP corrections are reported for BLYP, B3LYP, and B2PLYP in Figs 2 and S1; the individual BB-CP values are shown in Table S1. The figures show the CP averages (CP_avg) for the hydrogen-bonded, dispersion-dominated, and mixed complexes, as well as the CP_avg values for the eight different intermolecular distances. Hydrogen-bonded complexes show the largest CP_avg (BLYP = 3.10 kcal mol⁻¹, B3LYP = 2.65 kcal mol⁻¹, and B2PLYP = 2.76 kcal mol⁻¹), followed by mixed (BLYP = 1.36 kcal mol⁻¹, B3LYP = 1.16 kcal mol⁻¹, and B2PLYP = 1.34 kcal mol⁻¹), and dispersion-dominated complexes (BLYP = 1.22 kcal mol⁻¹, B3LYP = 1.05 kcal mol⁻¹, and B2PLYP = 1.31 kcal mol⁻¹) (see Fig. 2 for details). The complexes with the shortest intermolecular distances show the largest CP corrections; e.g. the CP_avg values for B2PLYP are 3.19, 2.79, 2.45, 2.16, 1.91, 1.30, 0.69, and 0.11 kcal mol⁻¹ for the dimers for the eight distances 0.90r_e, 0.95r_e, r_e, 1.05r_e, 1.10r_e, 1.25r_e, 1.5r_e, and 2r_e, respectively (Fig. S1).

Fig. 2.  Average BB-CP corrections for the S66x8 benchmark set and its categories. The def2-SVP basis set was used for all calculations.

Upon closer inspection of our results, we found that GGAs have a larger BSSE than HF theory. This also explain why a hybrid, such as B3LYP has smaller BSSE than BLYP. We note that the same can be seen in the discussion of the original gCP paper in Ref. [21].

Furthermore, B2PLYP shows a tendency to yield larger CP corrections than the B3LYP hybrid DFA. For instance, the CP_avg values for the entire S66x8 set are 1.82, 1.64, and 1.91 kcal mol⁻¹ for B2PLYP, B3LYP, and BLYP, respectively. This increase in CP_avg for B2PLYP is most likely due to its MP2 component. We validated this hypothesis by comparing our results for B2PLYP to those of other DHDFs with different degrees of MP2 contribution, as well as the MP2 method itself. The average CP correction for PBE0-2^[134] (79 % Fock exchange and 50 % MP2 correlation) and PBE0-DH^[135] (50 % Fock exchange and 12.5 % MP2 correlation) are 1.75 and 1.40 kcal mol⁻¹. The CP_avg values for MP2 and spin-component-scaled SCS-MP2^[136] are 2.31 and 2.23 kcal mol⁻¹, respectively. This shows that the amount of MP2 correlation is the major influencing factor for double hybrids. Our results in Figs 2 and S1 reveal that the BSSE of a double hybrid differs from that of the hybrids and pure DFAs. This justifies why the globally adjusted gCP parameters for DFT from Ref. [21] should not be used for double hybrids, which is why we refit them here for B2PLYP.

Although the BB-CP and gCP corrections were developed as independent procedures, London dispersion corrections play an important role in the DFT-based treatment of geometries, thermochemistry, kinetics, and noncovalent interactions.^{[29,35,71,72,133,137–146]} Herein, we will study in detail if a refit of the DFT-D3(BJ) parameters has any positive effects. We find a significant improvement in performance when the DFT-D3(BJ) parameters s₈, a₁, and a₂ in Eqn 6 are re-determined in the presence of the BB-CP correction with the procedure discussed earlier. This reparametrised version of B2PLYP-CP-D3(BJ) performs extremely well, with a decrease in the MAD from 0.73 to 0.26 kcal mol⁻¹ for the S66x8 benchmark set when compared to a BB-CP based version that uses the original DFT-D3(BJ) parameters. Consequently, we decided to refit the DFT-D3(BJ) parameters for the gCP scheme as well.

When the damping parameters for the DFT-D3(BJ) correction are optimised in the presence of the gCP correction, we see a significant improvement in performance and the mean absolute deviation (MAD) is reduced by 0.39 kcal mol⁻¹ compared to using the standard DFT-D3(BJ) parametrisation for B2PLYP (see Table 2 for details). Following that, we investigate whether we can apply the global gCP parameters implemented for DFT (Ref. [21]) to B2PLYP. As shown in Table 2, the accuracy obtained by re-fitting both the DFT-D3(BJ) and gCP parameters for B2PLYP is not negligible, which is why this is our chosen strategy for this study.

Table 2.  Mean absolute deviations (MADs) and mean deviations (MDs) in parentheses (in kcal mol⁻¹) for the S66x8 benchmark set and its categories
‘D3(BJ)’ and ‘gCP’ denote the addition of the gCP and dispersion corrections. The def2-SVP basis set was used for all calculations

The various fitted parameters can be found in Tables S2 and S3 in the Supplementary Material. To avoid confusion, we would like to clarify that ‘B2PLYP-gCP-D3(BJ)’ in this manuscript always refers to the refitted gCP and DFT-D3(BJ) parameters, whereas ‘B2PLYP-D3(BJ)’ uses the original DFT-D3(BJ) damping parameters for B2PLYP regardless of the chosen basis set.

In Fig. 3, we illustrate the effect of the gCP and dispersion corrections against the S66x8 reference interaction energies, which are based on high-level data at the CBS limit. Note that we use the reference values from the original publication^[41] to make our analysis commensurate with previous studies,^[21,147] particularly our GMTKN55 ones of the related S66.^{[35,38,71–73]} For slightly changed reference values, which would not change our main conclusions, see Refs [42 and 148]. Dispersion and gCP contributions to B2PLYP-gCP-D3(BJ) are reported separately to demonstrate the energetic effects of the corrections. A positive mean deviation (MD) indicates an overestimation of the interaction energies. Starting from plain B2PLYP/def2-SVP [MD = −0.15 kcal mol⁻¹ and MAD = 1.46 kcal mol⁻¹], the inclusion of the dispersion correction is expected to make the dimers of S66x8 more stable, which is indeed the case with MD/MAD = 1.14/1.22 kcal mol⁻¹ for B2PLYP-D3(BJ)/def2-SVP. These values are contaminated by BSSE. Additionally adding the gCP correction to the dispersion-corrected result induces the expected destabilising effect and brings the interaction energies closer to the reference values (MD/MAD = 0.11/0.56 kcal mol⁻¹). The hydrogen-bonded systems in the S66x8 benchmark set show overestimation of the interaction energies for plain B2PLYP (MD = 1.74 kcal mol⁻¹). Adding the dispersion correction [B2PLYP-D3(BJ)] leads to further over-stabilisation (MD = 2.47 kcal mol⁻¹). This in turn is compensated for by the gCP correction and results in a final MD of 0.62 kcal mol⁻¹ for B2PLYP-gCP-D3(BJ). The situation is different for the dispersion-driven systems of the S66x8 benchmark set, where B2PLYP yields interaction energies that are too low (MD = −1.72 kcal mol⁻¹). Correcting the energies for dispersion and BSSE brings interaction energies closer to the reference interaction energies (MD = 0.39 kcal mol⁻¹). A look at the third category of S66x8 shows that the gCP-D3(BJ)-corrected B2PLYP yields an MD that is better than the MD of plain B2PLYP (0.10 versus −0.51 kcal mol⁻¹).

In conclusion, most of the results shown in this manuscript can be explained by the advantageous interplay of the opposing effects of the dispersion correction (stabilisation) and the addition of gCP correction (destabilisation). Our results in this section therefore validate small basis set calculations for NCIs. Although our herein presented approach is not competitive with the standard B2PLYP-D3(BJ) near the basis set limit [MAD(B2PLYP-D3(BJ)/def2-QZVP^[52] = 0.16 kcal mol⁻¹],^[147] it certainly performs better than B2PLYP/def2-SVP and B2PLYP-D3(BJ)/def2-SVP without requiring additional computational resources.

Next, we look at the GMTKN55 database and its categories. The benchmarking is done against the original GMTKN55 reference values which mostly consist of high-level data at the complete basis set limit. We focus mostly on the NCI benchmark sets but also discuss the other categories briefly. Individual statistics for all 55 sets are shown in Tables S5–S7 in the Supplementary Material.

Fig. 4 shows the MDs of B2PLYP, B2PLYP-D3(BJ), and B2PLYP-gCP-D3(BJ) for the noncovalent interaction datasets of GMTKN55. Intermolecular NCIs are described by 12 benchmark sets (Table 1) and we discuss a few, select examples. The interaction energy trends for the n-alkane dimer set ADIM6^[28,35] are similar to those of S66x8, meaning that interaction energies are insufficiently described by plain B2PLYP (MD = −2.88 kcal mol⁻¹). The improvement by using gCP in conjunction with the reparametrised DFT-D3(BJ) is even better than for S66x8 and lowers the magnitude of the MD by 1.90 kcal mol⁻¹. The interaction energy trends for the bound dimers in S22 and S66 are similar to those of ADIM6 and S66x8, meaning that the introduction of dispersion into the treatment leads to an over-stabilisation of dimers and yields MDs of 1.81 kcal mol⁻¹ and 1.69 kcal mol⁻¹, respectively. This in turn is mostly compensated by additionally adding the gCP correction to yield nearly perfect MDs of −0.01 and 0.21 kcal mol⁻¹.

For RG18^[35] (interaction energies in rare gas complexes) the gCP and dispersion corrections show the same trends as before, but both corrections have different impacts on the interaction energies and the overall destabilising effect of gCP is more prominent [MD(B2PLYP) = −0.09 kcal mol⁻¹, MD[B2PLYP-D3(BJ)] = 0.42 kcal mol⁻¹, and MD[B2PLYP-gCP-D3(BJ)] = −0.83 kcal mol⁻¹]. The MD value of CARBHB12^[35] (hydrogen-bonded complexes of carbene analogues with small molecules) shows that dispersion- and BSSE-uncorrected B2PLYP slightly overestimates the interaction energies (MD = 0.93 kcal mol⁻¹). These interactions are overestimated even further when the dispersion correction is included with an MD value of 1.58 kcal mol⁻¹, which is reduced to 1.16 kcal mol⁻¹ when gCP and DFT-D3(BJ) are both applied. The halogen bonding interaction set HAL59^[35,74,75] also yields good results when the gCP and DFT-D3(BJ) corrections are applied, with the MD of 1.51 kcal mol⁻¹ for B2PLYP-D3(BJ) improving to 0.82 kcal mol⁻¹ for B2PLYP-gCP-D3(BJ).

We remind the reader that def2-SVP was enhanced by diffuse functions to take anionic systems into account in benchmark sets such as WATER27^[35,70] [binding energies in (H₂O)_n, H⁺(H₂O)_n, and OH^–(H₂O)_n] or IL16^[35,68] [interaction energies in anion-cation dimers] (see Computational Details section). Those diffuse functions were added to make our results consistent with previous GMTKN55 studies and to allow for an analysis of the entire database. We notice a strong underestimation of the interaction energies in WATER27 for B2PLYP-gCP-D3(BJ) (MD = −13.13 kcal mol⁻¹) and IL 16 (MD = 5.13 kcal mol⁻¹). In future work, reparametrising gCP for such a slightly changed version of def2-SVP might lead to improvements, but problems with the basis set cannot be ruled out either at this stage. That being said, the interaction energies of the cation-neutral dimers in CHB6^[35,68] are overestimated by 2.64 kcal mol⁻¹ for B2PLYP-gCP-D3(BJ) relative to B2PLYP. CHB6 does not require any diffuse functions and the fact that it is an outlier might also hint at basis set insufficiencies, which might require further study in the future.

When analysing intramolecular NCIs, we noticed that one of the largest deviations we observed for plain B2PLYP were for BUT14DIOL^[35,132] (relative energies in butane-1,4-diol conformers) and SCONF^[35,67,130] (sugar conformers) with MD values of 2.28 kcal mol⁻¹ and 1.22 kcal mol⁻¹, respectively. Both sets strongly benefit from the gCP-D3(BJ) combination (MDs = 0.89/0.02 kcal mol⁻¹), which is not surprising as conformer stabilities in those sets are driven by intramolecular hydrogen bonding; we have seen in our discussion of S66x8 how hydrogen-bonded systems benefit from both corrections (see Fig. 3). For the intramolecular NCI category of GMTKN55, conformational energies are sometimes difficult to determine, but generally acceptable at the B2PLYP/large-basis level. B2PLYP-D3(BJ) with a small AO basis fails to establish the correct order. The benchmark reference of one of the systems in the SCONF test set, for example, predicts the two most stable conformers to be 5.54 kcal mol⁻¹ apart from each other. B2PLYP-D3(BJ) predicts a difference of −1.25 kcal mol⁻¹ between the two, where the minus sign indicates a reversed energetic order of the two conformers. The order is correct for B2PLYP-gCP-D3(BJ), with an energy difference of 3.96 kcal mol⁻¹.

Fig. 3.  Mean and mean absolute deviations for the S66x8 benchmark set and its categories. The def2-SVP basis set was used for all calculations. A positive mean deviation indicates overestimation of the interaction energies.

Fig. 4.  Mean deviations (MDs) for the noncovalent interaction subsets of GMTKN55. The (aug’)-def2-SVP basis set was used for all calculations. A positive MD indicates overestimation of the interaction energies, except for AHB6, CHB21, and IL16.

Having seen the impact of the refitted gCP and DFT-D3(BJ) corrections on the NCI test sets, we now turn our attention to the non-NCI subsets of GMTKN55 (see Fig. 5). We discuss the reaction barrier heights first (seven sets, see Table 1). In the PX13^[116] test set (proton-exchange barriers in small-molecule clusters), B2PLYP and B2PLYP-D3(BJ) yield barriers that are on average too low (MDs of −6.19 and −6.51 kcal mol⁻¹, respectively). Correcting for BSSE in B2PLYP-D3(BJ) raises the barriers significantly (MD = −0.81 kcal mol⁻¹). For BHDIV10^[35]—a test set for barrier heights of ‘diverse reactions’—the addition of the dispersion correction leads to a general underestimation of the barriers (MDs of −0.39 to −0.94 kcal mol⁻¹, respectively), which is only partially compensated for by the gCP correction (MD = −0.89 kcal mol⁻¹).

Fig. 5.  Mean deviations (MDs) for the non-NCI test sets of the GMTKN55 database. The (aug’)-def2-SVP basis set was used for all calculations.

For BHROT27^[35]—barrier heights for rotation around single bonds—there is little statistical difference between the tested methods. For BHPERI^{[67,111–113]}—barrier heights of pericyclic reactions—B2PLYP shows better results than B2PLYP-D3(BJ) (MDs of −0.17 and −2.32 kcal mol⁻¹, respectively) and this time the addition of gCP does not improve the results (MD = −2.46 kcal mol⁻¹), which indicates that this set is more challenging then the previously discussed ones.

For several benchmark sets belonging to the remaining categories of GMTKN55, such as total atomisation energies for small systems (W4-11^[83]), dissociation and other reactions of alkaline compounds (ALK8^[35]), reaction energies of selected G2/97^[149] systems (G2RC^[35,67]), or oligomerisation and H₂ fragmentations of NH₃/BH₃ systems (NBPRC^{[35,67,78,89]}), B2PLYP or B2PLYP-D3(BJ) are statistically superior to B2PLYP-gCP-D3(BJ). A likely reason could be a seemingly better result due to error compensation between remaining errors of the method and the basis set incompleteness error. Alternatively, it can also mean that the reparametrised ‘gCP-D3(BJ)’ corrections might not be suitable for the properties covered by those benchmark sets. Additionally, atomisation energies in W4-11 are poorly described by a method relying on a small-basis MP2 contribution, such that it is not expected that either of the two corrections can improve the results. Having established the impact of the combined gCP and DFT-D3(BJ) corrections on B2PLYP/(aug’-)def2-SVP numbers, the remaining question is how close those numbers get to B2PLYP-D3(BJ) near the CBS limit, which for GMTKN55 is B2PLYP-D3(BJ)/(aug’)-def2-QZVP. In the following, we will collectively use the acronym ‘QZ’ for the larger basis set and ‘DZ’ for the smaller. In order to answer that remaining question, we calculated the differences in MADs between DZ-based levels and published^[35] B2PLYP-D3(BJ)/QZ results; this difference is called ‘ΔMAD’. The resulting values are shown for each benchmark set in Figs 6 and 7.

Fig. 6.  Mean absolute deviation difference (ΔMAD) between DZ-based results and B2PLYP-D3(BJ)/QZ for the NCI test sets in GMTKN55.

Fig. 7.  Mean absolute deviation difference (ΔMAD) between DZ-based B2PLYP results and B2PLYP-D3(BJ)/QZ for the non-NCI test sets in GMTKN55.

A positive ΔMAD means that the investigated method gives a larger average error than B2PLYP-D3(BJ)/QZ. A negative value means that we observe fortuitous compensation between shortcomings of the tested method and the basis-set incompleteness error, which make the MADs at the DZ level smaller than at the QZ level.

As can be seen in Figs 6 and 7, most of the NCI benchmark sets of GMTKN55, as well as a significant number of non-NCI test sets, benefit from the inclusion of the gCP correction. For example, the ΔMADs for S22 and S66 improve by 1.26 kcal mol⁻¹ and 1.06 kcal mol⁻¹, respectively. As expected, however, most ΔMAD values are positive.

The first GMTKN55 study introduced so-called weighted total mean absolute deviation (WTMAD) schemes. Details on the two schemes called ‘WTMAD-1’ and ‘WTMAD-2’ can be found in Ref. [35]. For our purpose, it is only important to know that the idea behind a WTMAD is to combine the MADs of the 55 individual benchmark sets to into one number, which can be used to conveniently assess functional robustness and allows to establish a functional ranking. The interested reader can find WTMAD-2 values in Table S8 in the Supplementary Material. In Ref. [150], an article showing why the popular B3LYP/6-31G*^[151] model chemistry should not be used, the authors introduced the WTΔMAD scheme for the GMTKN30^[78] predecessor of GMTKN55. According to that idea the individual MADs are replaced with ΔMADs, but otherwise the WTMAD formalism stays the same. We did the same here for our GMTKN55 results based on GMTKN55’s WTMAD-2^[35] scheme. Such WTΔMAD-2 values are listed in Table 3 for GMTKN55 and its individual categories. This allows us to compare the general applicability of B2PLYP-gCP-D3(BJ)/DZ relative to B2PLYP-D3(BJ)/QZ.

Table 3.  WTΔMAD-2 values (in kcal mol⁻¹) for the GMTKN55 database and its categories: basic properties and reactions of small systems (A), reaction energies of larger systems and isomerisation (B), barrier height (C), intermolecular noncovalent interactions (D), intramolecular noncovalent interactions (E), all noncovalent interactions (F) and the entire database (G)

For the entire database, we obtain WTΔMAD-2 values of 8.28 kcal mol⁻¹ for B2PLYP, 8.60 kcal mol⁻¹ for B2PLYP-D3(BJ), and 6.17 kcal mol⁻¹ for B2PLYP-gCP-D3(BJ), respectively, which shows an overall improvement for the latter. This is mostly due to the barrier height and NCI categories (see Table 3). We particularly notice a significant improvement when gCP is applied for intramolecular NCIs: B2PLYP-gCP-D3(BJ) has an WTΔMAD-2 of 7.59 kcal mol⁻¹ compared to 16.33 and 21.87 kcal mol⁻¹ for plain B2PLYP and B2PLYP-D3(BJ), respectively.

Overall, the accuracy of our developed scheme is limited due to the basis set incompleteness error associated with the small basis. That being said, gCP has been developed as a general-purpose tool and can be very easily adapted to different basis sets to limit that problem. Nevertheless, we have shown that the gCP scheme is in principle adaptable to double-hybrid DFT methods.

In this study, we investigated the applicability of Grimme’s geometric counterpoise correction (gCP) for the basis set superposition error (BSSE) in double-hybrid DFT calculations. This scheme only requires geometrical information as an input and no orbital/density information is needed. Therefore, this correction is almost free compared to the cost of the underlying electronic structure calculation. The choice of the functionals and basis sets was not the scope of this manuscript, and we exemplified our analysis with B2PLYP/def2-SVP. Our conclusions should be considered as a proof-of-concept study.

Our work revealed the error compensation between missing London dispersion and BSSE in B2PLYP/def2-SVP results for thermochemistry, kinetics, and noncovalent interactions. We suggest B2PLYP-gCP-D3(BJ) as the more reliable alternative, with both corrections explicitly parameterised for the chosen level of theory. In fact, reparametrising DFT-D3(BJ) for B2PLYP-gCP/def2-SVP was beneficial. Benchmark results on the S66x8 benchmark set and GMTKN55 database showed statistical improvements of B2PLYP-gCP-D3(BJ) over plain B2PLYP and B2PLYP-D3(BJ). B2PLYP just with dispersion correction [B2PLYP-D3(BJ)] showed a significant overestimation of interaction energies, barrier heights with larger deviations from reference values, and wrong relative stabilities in conformers. Those problems could be removed or their severity at least reduced by adding gCP. The error cancellation in plain B2PLYP is unsystematic and depends on the chemical nature of the test set. The improvement by gCP over B2PLYP-D3(BJ) is largest for noncovalent interactions (in particular intramolecular NCIs).

In summary, the addition of a specifically optimised gCP-D3(BJ) correction removes two major deficiencies (missing London dispersion and BSSE) and consequently shows a more robust and reliable picture of B2PLYP with a small basis set. Our observation of the opposing effects of the dispersion correction (stabilisation) and gCP correction (destabilisation) is commensurate with a similar picture drawn for the B3LYP/6-31G* level of theory in Ref. [150].

Despite the promising trends presented herein, a comparison of B2PLYP-gCP-D3(BJ)/def2-SVP and B2PLYP-D3(BJ)/def2-QZVP shows that the approach's accuracy is limited, mostly due to remaining basis set incompleteness errors, which again is something that can also be seen in Ref. [150]. However, if the application of larger basis sets is not feasible, we strongly recommend using the gCP-D3(BJ) variant. This study encourages the possibility of extending this method to other double hybrids and basis sets.

Additional information on BSSE in the S66x8 set, the optimised DFT-D3(BJ) and gCP parameters, an analysis of the damped-gCP variant, and all statistical results for GMTKN55 are available on the Journal’s website.

The data that support this study are available in the article and accompanying online supplementary material. More detailed information on specific energies can be obtained from the corresponding author upon request.

The authors declare no conflicts of interest.

This research did not receive any specific funding.