Benchmarking DFT methods with small basis sets for the calculation of halogen-bond strengths

In recent years, halogen bonding has become an important design tool in crystal engineering, supramolecular chemistry and biosciences. The fundamentals of halogen bonding have been studied extensively with high-accuracy computational methods. Due to its non-covalency, the use of triple-zeta (or larger) basis sets is often recommended when studying halogen bonding. However, in the large systems often encountered in supramolecular chemistry and biosciences, large basis sets can make the calculations far too slow. Therefore, small basis sets, which would combine high computational speed and high accuracy, are in great demand. This study focuses on comparing how well density functional theory (DFT) methods employing small, double-zeta basis sets can estimate halogen-bond strengths. Several methods with triple-zeta basis sets are included for comparison. Altogether, 46 DFT methods were tested using two data sets of 18 and 33 halogen-bonded complexes for which the complexation energies have been previously calculated with the high-accuracy CCSD(T)/CBS method. The DGDZVP basis set performed far better than other double-zeta basis sets, and it even outperformed the triple-zeta basis sets. Due to its small size, it is well-suited to studying halogen bonding in large systems.

section_content: A halogen bond is a noncovalent interaction between an electrophilic region on a halogen atom (σ-hole [1] ) and a nucleophilic region on another, or the same, molecule [2]. As is shown in recent reviews [3] [4] [5] [6] [7], it is an increasingly important tool in many fields of research, such as supramolecular chemistry, crystal engineering and, more recently, biosciences and medicinal chemistry. Much of the early work on halogen bonding was done on solid phase. However, in recent years, halogen bonding in solution has gained significant interest [8] [9] [10] [11] [12] [13] [14] [15] [16]. The theoretical aspects of halogen bonding have been studied extensively with the help of high-level and highaccuracy post-Hartree-Fock (HF) methods, such as those based on the coupled cluster (CC) technique [17] [18] [19] [20] [21] [22] [23] [24] [25] [26] [27]. Although halogen bonding is best described using high-level methods, it is recognized that such methods greatly limit the size of the model systems, posing severe limitations, e. g., in biosciences [4]. Therefore, studies with high-accuracy methods have focused mainly on small complexes of only a few atoms. Molecular mechanics methods, based on force fields, are much more suitable for large systems, and are routinely used to study e. g., protein-ligand interactions. However, partly due to its anisotropicity, halogen bonding has not been well-described with traditional force fields [28], which have only recently been modified to better describe halogen bonding [28] [29] [30] [31]. 

In terms of speed and accuracy, methods based on density functional theory (DFT) fall in between post-HF methods and molecular mechanics. Although much slower than the forcefield-based methods, they scale much better than the high-level Electronic supplementary material The online version of this article (doi:10. 1007/s00894-017-3212-4) contains supplementary material, which is available to authorized users. 

post-HF methods and are therefore better suited for studying large systems. As a downside, the accuracy of the DFT methods cannot be increased systematically the same way as in CC methods, where the accuracy can be improved by accounting more precisely for the electron correlation. Also the exact form of the exchange-correlation functional is not known. Several approximations have been made, resulting in a number of different functionals. Although some hierarchy from less accurate to more accurate can be built, there are no clear rules on which functional is the most suitable for a given problem. To find a suitable functional, it is useful to compare the results calculated with different functionals to those calculated with high-level post-HF methods such as CCSD(T)/CBS which, as stated by Rezác et al. [32], Bgives results very close to the highest reasonable computational methods available at the moment^. 

Besides the functionals, another important factor in determining the speed and accuracy of the calculations is the quality of the basis set. Although some benchmarking of DFT has been made for halogen bonding, the studies have focused mainly on comparing different functionals and less on comparing the basis sets [21, 33]. Since large basis sets, in general, yield more accurate results, triple-zeta basis sets or larger have been employed routinely when comparing different functionals. However, in order to study halogen bonding in large systems, such basis sets can be too large. Smaller basis sets, such as those of double-zeta quality, make calculations faster but have been considered to be too small to accurately describe noncovalent interactions, including halogen bonding. However, with large systems they may be the only DFT-based option, and a fast but accurate DFT method would be of great help in studying halogen bonding in supramolecular chemistry and biosciences. 

In this study, we used two data sets of halogen-bonded complexes (Fig. 1 ) to compare DFT methods with small, double-zeta basis sets for their accuracy to estimate halogen-bond strengths. Larger, triple-zeta basis sets were also included for comparison. The reference values were taken from the literature and have been calculated with CCSD(T)/CBS [21, 22]. The smaller data set XB18 [21] (Fig. 1a ) contained complexes of nine halogen diatomics with hydrogen cyanide or formaldehyde, and was used to compare 46 DFT methods for their accuracy in reproducing the reference dissociation energy values. The small sizes of the complexes made it practical to evaluate a large number of methods. Also, the range of dissociation energies, from 1. 41 to 10. 09 kcal mol -1, covered interactions from weak to strong halogen bonds. A larger data set of 33 complexes [22] (Fig. 1b ) between iodobenzene derivatives and pyridine (Iodobenzene data set) was then used to further evaluate several methods. The latter set contained the types of complexes more relevant to supramolecular chemistry, crystal engineering and biosciences. Additionally, it contained molecules that are structurally very similar to each other, and is therefore a good set for testing how well the DFT methods can estimate changes in halogen-bond strength when the bond donor is only slightly modified, for example by adding a chloro-substituent on the aromatic ring. Such small changes are very common in drug design, and could also be useful in fine-tuning halogen bonding in supramolecular chemistry or crystal engineering. The dissociation energies in the set ranged from 3. 43 kcal mol -1 to 6. 65 kcal mol -1, which covers typical halogen-bond strengths quite well.
section_name: Introduction

section_content: All calculations were performed with Gaussian 09 Revision D. 01 (Gaussian Inc., Wallingford CT). Ultrafine grid was used for all calculations. Default values were used for other parameters. The basis sets not available in Gaussian were obtained from the EMSL Basis Set Exchange [34, 35] (http://bse. pnl. gov/). All geometries were optimized on a potential energy surface corrected for basis set superposition error (BSSE) using the counterpoise method by Boys and Bernardi [36] as implemented in Gaussian 09. Frequency calculations were performed for all complexes after geometry optimization to confirm that the obtained geometries are true minima (i. e., no imaginary frequencies). Pseudopotentials were used for iodine (cc-basis sets), bromine and iodine (def2basis sets) or chlorine, bromine and iodine (LanL2DZ and LanL2DZdp basis sets). The reference values for the halogen-bonded complexes were taken from the literature, and have been calculated with CCSD(T)/CBS [21, 22]. The more technical side of benchmarking noncovalent Fig. 1 Structures of the complexes in the XB18 [21] and iodobenzene [22] data sets interactions, including halogen bonding, has been reviewed previously [23, 37, 38] and is outside the scope of this paper. 

Our main focus was to test how well double-zeta basis sets performed in describing halogen bonding energies. Additionally, to study how the best double-zeta basis sets compared to larger basis sets, several triple-zeta basis sets, with and without diffuse functions, were included. Altogether, 16 basis sets were studied (see Table 1 ): def2-series: [39] [40] [41] [42] def2-SVP, def2-SVPD, def2-TZVP, def2-TZVPD; cc-series: [42] [43] [44] [45] [46] [47] [48] cc-pVDZ-PP, aug-cc-pVDZ-PP, cc-pVTZ-PP, aug-cc-pVTZ-PP; others with pseudopotentials: LanL2DZ [49], LanL2DZdp [50] ; allelectron basis sets: DGDZVP [51, 52] and those by Jorge and coworkers: [53] [54] [55] [56] [57] [58] [59] [60] [61] [62] DZP, ADZP, TZP and ATZP. We focused mainly on the functional M06-2X [63], which has been shown previously to perform very well in estimating halogen bond strengths when combined with a large basis set aug-cc-pVTZ [21, 64]. The doublezeta basis sets were also tested with two other functionals, PBE0 [65] and B3LYP [66]. The former was chosen because it is a popular functional and has often been used to study halogen bonding, especially to rationalize experimental results [67] [68] [69] [70]. Although the latter is known to be quite poor in describing halogen bonding interactions [4, 33], it has been, and still is, often used in various other DFT studies and was therefore chosen for comparison. Also, attempts have been made to improve its accuracy, for example by including empirical dispersion [71] [72] [73] (B3LYP-D3), long-range correction (CAM-B3LYP [74] ) or both (CAM-B3LYP-D3). These variants and some other functionals with global, short-range and/or long-range HF exchange, and in some cases with dispersion correction, were also tested: ωB97 [75], ωB97X [75], ωB97XD [76], B97-D3 [73, 77], M06 [63], M06-L [78], M06-HF [79], M11 [80] and MN12-SX [81]. 

To analyze the data, the calculated energies and the reference energies were compared by several measures. Rootmean-square of deviation (RMSD) depicts the average deviation, mean absolute error (MAE) reveals if the method is systematically over-or underestimating dissociation energies, and maximum deviation is the absolute value of the biggest error among the complexes. The R 2 -value, when plotting a linear correlation between the calculated values and the reference values, is another measure of accuracy, and, for a good correlation, the value should be close to 1. It is also useful to test if there is a linear correlation between the errors and the reference energy values. For a good method, when plotting the errors as a function of the reference energies, the slope should have a value close to 0, meaning that the error is not dependent on the dissociation energy value. A value >0 indicates that the method systematically underestimates low energies and overestimates the high energies, or that it overestimates (underestimates) low energies less (more) than higher energy values. A value <0 indicates the opposite.
section_name: Methods

section_content: XB18 data set
section_name: Results and discussion

section_content: The statistical data for the M06-2X functional, when applied to the XB18 data set (see Fig. 1 ) using different basis sets, is shown in Table 1. The five double-zeta basis sets without diffuse functions showed mixed results. The most important observation is that the all-electron basis set DGDZVP is by far the most accurate. The correlation with the reference values is excellent (R 2 = 0. 997), the mean absolute error (MAE) is practically zero (-0. 001 kcal mol -1 ), and the maximum deviation is only 0. 3 kcal mol -1. The RMSD-value is the lowest among all the basis sets tested, including the triple-zeta basis sets. 

On the contrary, the other all-electron basis set DZP performs poorly, and consistently underestimates the dissocia t i o n e n e r g i e s. O f t h e b a s i s s e t s e m p l o y i n g pseudopotentials, the often-used def2-SVP and cc-pVDZ-PP perform satisfactorily, but the correlation is worse and the maximum deviation higher than for DGDZVP. LanL2DZ, the smallest of the basis sets, shows the poorest results, with R 2 -value of only 0. 945, maximum deviation of 2. 36 kcal mol -1, and significantly higher RMSD-value than for the other basis sets. Apart from LanL2DZ, the slope values for these basis sets are close to zero, i. e., they do not significantly underestimate low dissociation energies or overestimate high energies. 

Overall, the order of accuracy for these five basis sets is DGDZVP > def2-SVP ≈ cc-pVDZ-PP > DZP > LanL2DZ. 

Since there are no diffuse functions available for the DGDZVP basis set, only the other four basis sets were tested with additional diffuse functions. Although the diffuse functions should, in principle, improve the results, for some basis sets, the opposite is actually observed. Although def2-SVPD and aug-cc-pVDZ-PP do show better correlation with the reference values (i. e., higher R 2 ), the RMSD, MAE, maximum deviations, and slope values are actually higher. On the contrary, augmenting LanL2DZ with diffuse functions improves the results greatly, and the resulting basis set LanL2DZdp shows the second-best results from all the double-zeta basis sets studied, being outperformed only by DGDZVP. Augmenting DZP with diffuse functions improves the results slightly, although the R 2 -value remains low. 

As expected, the triple-zeta basis sets are, in general, much more accurate than the double-zeta basis sets, with the important exception of DGDZVP. The cc-pVTZ-PP basis set shows significantly better results than def2-TZVP, and is very similar in accuracy to DGDZVP. Unlike with the double-zeta basis set DZP, adding diffuse functions to TZP decreases the accuracy significantly, and ATZP gives exceptionally bad results, worse than any other basis set tested. Similarly to double-zeta basis sets, augmenting triple-zeta basis sets with diffuse functions increases the slope values. The effect is more pronounced for cc-pVTZ-PP than for def2-TZVP. 

Overall, DGDZVP is much more accurate than the other double-zeta basis sets and, despite its smaller size, is comparable to, and even better than, the much more computationally demanding triple-zeta basis sets.
section_name: M06-2X

section_content: The statistical data for the PBE0 and B3LYP functionals, when applied to the XB18 data set (see Fig. 1 ) using different basis sets, is shown in Table 2. DGDZVP again gives the best results, and with the other basis sets the addition of the diffuse functions reduces accuracy. The dissociation energies are approximately 0. 10 -0. 25 kcal mol -1 lower than with M06-2X, depending on the basis set. Although PBE0 is fairly accurate, it tends to underestimate low dissociation energies and slightly overestimate higher energies. Apart from ADZP, this is especially apparent with basis sets augmented with diffuse functions (slope = 0. 22 -0. 31). 

As expected, B3LYP shows poor results, and the dissociation energies are underestimated with all the basis sets. However, with B3LYP, the diffuse functions clearly improve the results and def2-SVPD and aug-cc-pVDZ-PP are slightly more accurate than the other basis sets. On the other hand, similarly to PBE0, the diffuse functions cause a significant increase in the slope values. Unlike with M06-2X and PBE0, DGDZVP did not perform especially well with B3LYP.
section_name: PBE0 and B3LYP

section_content: The tendency of B3LYP to underestimate noncovalent interaction energies is well known, and attempts have been made to improve its accuracy. Several modifications were tested with the DGDZVP basis set, and the results are presented in Table 3. Adding dispersion correction (-D3) to B3LYP improves the performance significantly. The energies on average are increased by 1. 1 kcal/mol and the MAE is only 0. 14 kcal/mol. The RMSD and maximum deviations are also significantly lower. Another approach to improve B3LYP is to include a long-range correction, for example by the coulomb-attenuating method. The resulting functional CAM-B3LYP shows improvement over B3LYP but it still underestimates the dissociation energies. Adding dispersion correction (-D3) to CAM-B3LYP further improves the accuracy slightly. 

Several other functionals with the DGDZVP basis set were also tested. The long-range-corrected functional ωB97, which does not contain exact short-range exchange, is between PBE0 and B3LYP in accuracy. The variant ωB97X, which contains a small fraction of exact shortrange exchange performs slightly better. Adding Grimme's D2 dispersion correction to ωB97X (ωB97X-D), does not change the overall accuracy, but, surprisingly, decreases the dissociation energies and the mean absolute error changes from 0. 44 to -0. 44 kcal mol -1. These functionals do not show bias to overestimate higher dissociation energies, and the slope values are all close to zero.
section_name: Dispersion and long-range corrections

section_content: The other M06 functionals of the Minnesota series, namely M06, M06-L and M06-HF, do not perform much better than B3LYP (Table 3 ). M06-HF, which contains 100% exact exchange, is slightly less accurate than the other two. Similarly to what has been previously reported with a large basis set [21], M06 and M06-L overestimate the dissociation energies, whereas M06-HF underestimates them. The newer variant M11, which also contains 100% exact exchange, performs similarly to M06-HF. MN12-SX was similar in accuracy to M06, but it underestimates the dissociation energies instead of overestimating them, and the underestimation is systematically bigger for higher dissociation energies, as shown by the negative slope value.
section_name: Minnesota functionals

section_content: As DGDZVP was shown to perform very well with the XB18 data set, it was used to study also the larger iodobenzene data set. Three basis sets of the cc-pVXZ series were also included for comparison. The results are summarized in Table 4. In contrast to the XB18 data set, PBE0 showed better results than M06-2X with all the basis sets. B3LYP again underestimated the dissociation energies significantly (MAE = -1. 24 kcal mol -1 ). As seen with B3LYP versus B3LYP-D3, the dispersion correction again improves the results, although less than with the XB18 data set, and B3LYP-D3 clearly overestimates dissociation energies. 

The cc-pVTZ-PP basis set, like the XB18 data set, shows results very similar to DGDZVP. With M06-2X, the doublezeta basis sets cc-pVDZ-PP and aug-cc-pVDZ-PP show clearly worse results than DGDZVP or cc-pVTZ-PP, but with PBE0 the difference is smaller. 

With PBE0, all basis sets showed better results than M06-2X/DGDZVP, but, as previously mentioned, PBE0 tends to underestimate low energies and overestimate high energies. The tendency to overestimate higher energies can be seen clearly with all basis sets, and the slope values are consistently higher for PBE0 than for M06-2X (Table 4 ). The energies in this data set vary from 3. 43 kcal mol -1 to 6. 65 kcal mol -1, which seems to be a good range for PBE0. Presumably due to the lack of very low or very high dissociation energies in this data set, PBE0 performs well. 

Unlike the other basis sets, apart from LanL2DZ, DGDZVP does not contain polarization functions on hydrogens. It was recently shown [82] that they have very little effect on the calculated electrostatic potential maxima or minima on molecular surfaces, and are therefore not expected to have much effect on the halogen-bond dissociation energies. They do, however, increase the number of basis functions and therefore the computation time, especially when dealing with large molecules containing many hydrogens.
section_name: Iodobenzene data set

section_content: Our aim was to benchmark DFT methods with small basis sets for their ability to reproduce the halogen-bond strengths previously calculated with the high-accuracy computational method CCSD(T)/CBS. Initially, 46 DFT methods (15 functionals and 16 basis sets) were evaluated using a data set of 18 small halogen-bonded complexes (XB18 data set [21] ); 12 of the methods were also tested using a larger data set of 33 complexes between iodobenzene derivatives and pyridine (iodobenzene data set [22] ). Although the main focus was in evaluating double-zeta basis sets, several triple-zeta basis sets were included for comparison. 

The dissociation energies of the XB18 data set were best reproduced with the DGDZVP basis set, which was more accurate than even the much larger triple-zeta basis sets, although cc-pVTZ-PP was quite close in accuracy. M06-2X was clearly the best functional for this data set, and the combination M06-2X/DGDZVP reproduced the reference dissociation energies very well. Another commonly used functional, PBE0, gave quite good results although it had a tendency to underestimate low energies and overestimate high energies. In the case of the iodobenzene data set, PBE0 performed better than M06-2X with all the basis sets, although it again underestimated low energies and overestimated high energies. As with the XB18 data set, DGDZVP was the most accurate, followed closely by cc-pVTZ-PP. As expected, B3LYP gave poor results for both data sets, while adding empirical dispersion correction improved the results significantly. 

In conclusion, for both data sets the double-zeta basis set DGDZVP performed better than the other double-zeta basis sets, and at least as well as triple-zeta basis sets. Augmenting basis sets with diffuse functions did not always improve the results, and, in some cases, made them much worse. Due to its small size, DGDZVP is well-suited for calculating halogenbond strengths on large complexes, or on a large set of smaller complexes, and is therefore an excellent choice when studying halogen-bonding in the design of supramolecular functional materials or its role in biosciences.
section_name: Conclusions

section_content: Acknowledgements A. P. greatly acknowledges the Academy of Finland Research Fellowship Program (Decision Nos. 277091 and 284553 ) and the Emil Aaltonen Foundation for financially supporting this work.