Density Functional Theory for Microwave Spectroscopy of Non

Article pubs.acs.org/JPCA

Cite This: J. Phys. Chem. A 2018, 122, 4894−4901

Density Functional Theory for Microwave Spectroscopy of Noncovalent Complexes: A Benchmark Study P. Kraus* and I. Frank Institut für Physikalische Chemie und Elektrochemie, Leibniz Universität Hannover, Callinstraße 3A, 30167 Hannover, Germany

Downloaded via KAOHSIUNG MEDICAL UNIV on July 25, 2018 at 04:24:54 (UTC). See https://pubs.acs.org/sharingguidelines for options on how to legitimately share published articles.

S Supporting Information *

ABSTRACT: In this work, we compare the results obtained with 89 computational methods for predicting noncovalent bond lengths in weakly bound complexes. Evaluations for the performance in noncovalent interaction energies and covalent bond lengths obtained from five other data sets are included. The overall best performing density functional is the ωB97M-V method, achieving balanced results across all three categories. For noncovalent geometries, the best methods include B97M-V, B3LYP-D3(BJ) and DSD-PBEPBE-D3(BJ). The effects of systematic improvement of the density functional approximation and of dispersion corrections are also discussed.

1. INTRODUCTION For noncovalently bound complexes, a large portion of recent efforts in research, development, and benchmarking of computational methods focuses on obtaining accurate interaction energies.1−4 A great amount of progress has been achieved in this field over the recent years. At the highest levels of wave function theory (WFT), it has been shown that interaction energies computed with coupled cluster calculations including singles, doubles, and pertubative triples used with a basis set extrapolated to the complete basis set limit (CCSD(T)/CBS) are in excellent agreement with calculations also including core−core and core−valence correlation, relativistic effects and quadruple excitations.3 On the other hand, the advent of easy-to-use empirical dispersion corrections,5 now including three-body interactions and effects of local chemical environment,6 dynamic polarizabilities of atoms,7 or nonlocal correlation,8 allows density functional theory (DFT) methods to achieve a comparatively good performance in the description of interaction energies, provided large enough basis sets are applied,4 even when comparatively dated density functionals are used. Additionally, with double-hybrid density functional methods, which include a part of MP2 correlation in addition to the part of Hartree−Fock (HF) exchange used in single-hybrid DFT methods, accuracy comparable to CCSD(T)-based benchmarks is possible.9,10 However, accurate description of energies is only one aspect that makes a computational method generally applicable. For certain structure-sensitive applications, including microwave spectroscopy, an accurate prediction of geometries is equally, if not more important. The data acquisition in such experiments generally starts with a computation of a predicted structure, to obtain a first guess of the rotational constants. If the predicted rotational constants differ significantly (>10 MHz) from the true values, the potential to waste valuable instrument time is high, especially if a broadband instrument is unavailable. In addition to the data collection issues, accurate structural © 2018 American Chemical Society

predictions and ab initio based spectra can decrease the difficulty in interpretation of experimental data. Furthermore, there is currently little guidance for an unbiased choice of computational methods, as accuracy for noncovalent structures does not immediately follow from an accurate prediction of energies. Therefore, in the current work, we attempt to answer the question: Which affordable computational method is currently the best for predicting geometries of noncovalently bound complexes? To help us answer this question, we have used the recently developed data set of semiexperimental noncovalent bond lengths in dimers and trimers (NCDT16)11 and 89 computational methods, including the HF, MP2, CCSD and CCSD(T) wave function-based methods. On the density functional front, 16 generalized density gradient approximation (GGA) methods, 11 meta-GGA functionals depending also on the kinetic energy density, 28 hybrid GGA functionals incorporating HF-exchange, 20 hybrid meta-GGA functionals, and 8 double-hybrid functionals also incorporating MP2 correlation are compared. Additionally, two composite methods are included: an “affordable” PBEh-3c method intended to predict good geometries12 and a “near gold-standard” MP2/ CBS + δCCSD(T)/3-ζ extrapolated-basis method (shortened to MP2 + δ(T) in the following). The NCDT16 data set comprises 45 noncovalent bond lengths between the heavy atoms of 16 dimers and trimers, obtained using experimental rotational constants corrected for nonequilibrium effects at the double-hybrid level of DFT. While angles are not included in the database, angular information is included by considering more than one bond length for all nonsymmetric complexes.11 The current data set is somewhat skewed toward dispersion interactions, with 9 out of 16 complexes having interaction Received: April 9, 2018 Revised: May 7, 2018 Published: May 11, 2018 4894

DOI: 10.1021/acs.jpca.8b03345 J. Phys. Chem. A 2018, 122, 4894−4901

Article

The Journal of Physical Chemistry A

Table 1. Overview of the Computational Methods Used in This Work, Grouped According to Their Classification (“Jacob’s Ladder”): GGAs (F(∇ρ[r])), Meta-GGAs (F(∇ρ[r], τ[r])), Single Hybrids, Single Meta-GGA Hybrids, Double Hybrids, and Coupled Cluster Methods GGA

meta-GGA

XHF + GGA

XHF + meta-GGA

XHF + CMP2 + GGA

CC

BLYP BLYP-D2 BLYP-D3(BJ) BLYP-D3m(BJ) PBE PBE-D2 PBE-D3(BJ) PBE-D3m(BJ) HCTH/120 HCTH/120-D3(BJ) VV10 B97-D2 B97-D3(BJ) B97-D3m(BJ) N12 N12-D3(BJ)

M06-L M06-L-D3 M11-L M11-L-D3(BJ) MN12-L MN12-L-D3(BJ) MN15-L MN15-L-D3 TPSS TPSS-D3(BJ) B97M-V

HF B3LYP B3LYP-D2 B3LYP-D3 B3LYP-D3(BJ) B3LYP-D3m(BJ) B3PW91 B3PW91-D3(BJ) ωPBE ωPBE-D3(BJ) ωPBE-D3m(BJ) PBE0 PBE0-D2 PBE0-D3 PBE0-D3(BJ) PBE0-D3m(BJ) ωPBE0 LC-VV10 B97 ωB97X ωB97X-D ωB97X-V BHandH BHandHLYP HSE06 HSE06-D3(BJ) SOGGA11-X SOGGA11-X-D3(BJ) PBEh-3c

M05 M05-D3 M05-2X M05-2X-D3 M06 M06-D3 M06-2X M06-2X-D3 M06-HF M06-HF-D3 M08-HX M08-HX-D3 M08-SO M11 M11-D3(BJ) MN15 MN15-D3(BJ) TPSSH TPSSH-D3(BJ) ωB97M-V

MP2 B2PLYP B2PLYP-D3(BJ) B2PLYP-D3m(BJ) PBE0-2 DSD-BLYP-D3(BJ) DSD-PBEP86-D3(BJ) DSD-PBEPBE-D3(BJ)

CCSD CCSD(T) MP2 + δ(T)

All electronic structure calculations were carried out with the density-fitted variants of WFT (including MP2, CCSD or CCSD(T)) and DFT methods. An “ultra fine” (99, 590) integration grid was used; for cases where SCF convergence was difficult, an even larger (250, 974) grid was applied with a density damping of 20%. The convergence thresholds in geometry optimizations were tightened to ΔE < 1 μEh, Fmax < 60 μEh/a0, FRMS < 40 μEh/a0, ΔRmax < 1 ma0, and ΔRRMS < 1 ma0. In some complexes, such as for the Ne−OCS dimer, results vary significantly when the default (looser) convergence thresholds are used. For all methods with the exception of coupled cluster and composite calculations, the def2-TZVPD basis set was used17,18 due to its availability for heavier elements. For the CCSD and CCSD(T) calculations, the maycc-pVTZ basis set was used19 to reduce computational cost. The PBEh-3c composite method has been applied with its modified double-ζ basis set.12 Finally, the extrapolated composite method denoted “MP2 + δ(T)” comprises a frozen-core density-fitted MP2 calculation extrapolated toward the basis set limit20 by using triple- and quadruple-ζ correlation-consistent basis sets (i.e. MP2/aug-cc-pV[TQ]Z),21,22 with further correlation effects (“δ(T)”) approximated from the difference between frozen-core density-fitted CCSD(T) and MP2 energies at the may-cc-pVTZ level. The use of quadruple-ζ basis in the MP2 extrapolation is critical, as the extrapolation from double- and triple-ζ basis sets (ie. MP2/augcc-pV[DT]Z) performs considerably worse than the MP2/ def2-TZVPD method. An overview of the computational

energies below 4 kJ/mol. However, three complexes with dominant π−π interactions and four H-bonded systems ensure that the interaction energy ranges of 5−10 kJ/mol and >15 kJ/ mol, respectively, are also considered. Interaction energies for clusters in the NCDT16 data set, obtained with the bestperforming methods, are presented at the semiexperimental (rSE e ) geometries and compared to Eint results at the optimized geometries. To make our analysis more robust, we have carried out further calculations on five other previously published benchmark data sets: a set of semiexperimental bond lengths in small covalently bound compounds (CCse22),13 its ab initio bond length counterpart (CS20),14 two data sets of highly accurate ab initio interaction energies at the prescribed geometries (NCCE3114 and A243), and a database of 21 equilibrium CCSD(T)/CBS structures (A21G)3 to serve as a counterpart to our NCDT16 set.

2. COMPUTATIONAL METHODS All ab initio calculations were performed using a development version (ver. > 1.2a1.dev781) of the electronic structure code ψ4,15 coupled to the DFT exchange-correlation library LIBXC v. 3.0.0.16 Our modifications to the ψ4 codebase consisted of adding numerous density functionals which are a part of LIBXC v. 3.0.0; they are now included in the development version. The use of such open-source programs is in our view critical for reproducibility and auditing of computational data. As a result, the driver routines used in the current work are available in the Supporting Information code archive. 4895


Article


Figure 1. MAE of various computational methods in selected noncovalent bond lengths averaged over the 16 complexes in the NCDT16 data set. The color denotes the class of the method (dark gray, GGA; light gray, meta-GGA; blue, global hybrid; light blue, range-separated hybrid; purple, meta-GGA hybrid; pink, range-separated meta-GGA hybrid; orange, double hybrid; yellow, coupled cluster). Accuracy decreases from left to right. Error bar shows one standard error (σ(|ΔreSE|)/ 16 ).

Figure 2. MAE of the ten best computational methods in subsets of the NCDT16 data set: (A) dispersion dominated complexes including Ne, (B) other dispersion dominated complexes, (C) complexes with significant π−π interactions, and (D) H-bonded complexes. Colors are as in Figure 1.

than CCSD(T), as well as the extrapolated MP2 + δ(T) composite method, include the dispersion corrected global hybrid B3LYP-D3(BJ),25−27 two range-separated functionals ωB97X-V28 and ωB97M-V29 related to the best performing meta-GGA, the original double-hybrid with dispersion correction B2PLYP-D3(BJ),9,27 and the range-separated functional with nonlocal correlation LC-VV10.8 Somewhat surprisingly, the dispersionless double hybrid PBE0-210 also performs very well. A more detailed look into the top performers in the four categories of complexes comprising the NCDT16 data set is shown in Figure 2. Detailed results for all methods are listed in the Supporting Information. The ranking of the methods varies significantly, with only B3LYP-D3(BJ)25−27 achieving a “top 10” result in three out of the four categoriesits performance for the H-bonded data set is average. The lowest MAE for the five Ne-containing complexes is significantly higher than in the

methods used in the current work is shown in Table 1. For full references of all listed methods, see the Supporting Information.

3. RESULTS AND DISCUSSION The mean absolute errors (MAEs) of selected noncovalent bond lengths (|ΔreSE|) averaged over the 16 complexes comprising the NCDT16 data set (|ΔreSE|) are shown in Figure 1 along with their associated standard errors (σ(|ΔreSE|)/ 16 ). The best performer is the meta-GGA functional B97M-V,23 with |ΔreSE| = 48 mÅ. This corresponds to an accuracy better than several double hybrid functionals and even the considerably more expensive CCSD(T) method. The only other method with a MAE below 50 mÅ is the double-hybrid DSD-PBEPBE-D3(BJ).24 Other methods achieving lower MAE 4896


Article


significantly overpredicting the H2CO−HCN angle compared to both MP2 + δ(T) and the semiexperimental uncertainty. 3.1. Interaction Energies of the NCDT16 Data Set. Counterpoise-corrected interaction energies of the 16 weakly bound clusters have been calculated at the semiexperimental geometries, as well as at the optimized geometries, with the ten best-performing methods. The full results are listed in the Supporting Information. The MAE values with respect to the MP2 + δ(T) results are shown in Figure 4. The average

other three categories, with meta-GGA’s performing comparably poorly: the best meta-GGA hybrid M06-2X-D330,31 is only rank 17. On the other hand, for H-bonded complexes, the more recent meta-GGA’s of Truhlar and co-workers32,33 top the fieldthe overall MAE of the range-separated meta-GGA MN12-SX,34 which is the top method in this category, is hampered only by the large MAE for the Ne−Ne−N2O trimer. On the other hand, the highest ranked double hybrid for the Hbonded complexes is the PBE0-2 functional10 at number 13. The other two classifications show more mixed results. The excellent performance of dispersion-corrected GGA’s for π−π bonded systems is somewhat surprising, especially in the case of the fourth best N12-D3(BJ)35,36 functional which has a 10× larger overall |ΔreSE| at ∼270 mÅ. Also notable are the rather average results of MP2 + δ(T) for π−π interactions, at rank 30. The absolute deviations of the calculated bond lengths from the semiexperimental equilibrium values (|ΔrSE e |) are shown for the top density functional B97M-V and the MP2 + δ(T) methods in Figure 3. The reported uncertainties in the

Figure 4. MAE values of counterpoise-corrected interaction energies obtained with the listed methods with respect to the MP2 + δ(T) data at rSE e geometries (blue) and at optimized geometries (orange).

difference between Eint values at rSE e and at the MP2 + δ(T)optimized geometries is just below 1.1 kJ/mol, with the highest difference (10.8 kJ/mol or 43%) obtained for the HF−H2CO complex. This increase in interaction energies and associated errors should be considered when further components of binding energies are required, as zero point vibrational corrections have to be obtained at the minima. The MN15-D3(BJ) functional performs rather well for both sets of interaction energies and it is the best out of the 20 methods studied. Among the best methods for noncovalent geometries (top part of Figure 4), the error in the interaction energies at the minima increases compared to the error at the rSE e structures, mainly due to a further overprediction of the interaction energies in H-bonded complexes. Methods without empirical or nonlocal dispersion terms generally underpredict the interactions for dispersion-dominated complexes, even when MP2 correlation is included (cf. PBE0-2). The only exception among the functionals shown is M05-2X, which on average overpredicts Eint of Ne-containing complexes by about 0.5 kJ/mol. The large absolute error in CCSD(T) data compared to the MP2 + δ(T) results is a systematic underprediction of the interaction energies due to the smaller basis set used, possibly further exacerbated by the application of counterpoise corrections.37 3.2. Correlations with Other Data Sets. The performance of the methods listed in Table 1 for the accuracy of noncovalent bond lengths has been correlated with the performance observed for five other benchmark data sets, with the overall results shown in Table 2. The six data sets can

Figure 3. Absolute deviations of the noncovalent bond lengths calculated with B97M-V (gray) and MP2 + δ(T) (red) compared to the semiexperimental uncertainties (error bars).

semiexperimental bond lengths11 are also included. For Necontaining complexes, both methods struggle with the Ne− OCS complex and with the position of Ne in the Ne−Ar−N2O and Ne−Ar−HCl trimers. In other dispersion dominated complexes, the rather poor performance of the composite method for the Kr−OCS dimer could be due to the frozen core approximation or neglection of relativistic effects in the basis sets used. The results of the two methods for π−π interactions is comparable, with the large deviation in one of the CS2−OCS parameters due to a different angle between the two monomers. In H-bonded complexes, the composite method generally performs better than B97M-V, with the latter method 4897


Article


poor description of dispersion in the three-parameter functionals in all rare gas complexes (see the Supporting Information for details). However, with dispersion corrections, the B3LYP-D3(BJ) method25−27 is one of the top methods. The latter two correlations are shown in graphical form in

Table 2. Pearson’s Correlation Coefficient between Studied Benchmark Data Sets NCDT16 A21G A24 NCCE31 CCse22

A21G

A24

NCCE31

CCse22

CS20

0.94

0.84 0.87

0.56 0.62 0.74

0.21 0.21 0.31 0.54

0.24 0.19 0.38 0.48 0.81

Figure 5. Only methods with |ΔreSE| < 200 mÅ in the NCDT16

be split into three pairs: noncovalent semiexperimental and theoretical equilibrium bond lengths (NCDT1611 and A21G,3 respectively), counterpoise-corrected interaction energies (A243 and NCCE3114), and covalent semiexperimental and theoretical equilibrium bond lengths (CCse2213 and CS20,14 respectively). Despite using different basis sets, our results calculated for the five data sets are generally comparable to literature data, obtaining qualitative agreement with the A21G data of Mardirossian and Head-Gordon,28 an agreement within 0.64 kJ/mol with the results of Peverati and Truhlar14 and Roch38 for the A24 and NCCE31 data sets, or within 1 mÅ for the CS20 results of Peverati and Truhlar.14 The average errors of the methods for the two data sets of noncovalent bond lengths are strongly correlated, with ρNCDT16, A21G = 0.94. The MAEs obtained for the A21G data set are much lower than for the NCDT16 data set, as in the A21G benchmark the whole bond matrix is assessed. The only bond length present in both data setsthe r(C−C) of the Tshaped HCCH−HCCH dimerdiffers by 70 mÅ, with CCSD(T)/CBS underpredicting the rSE e . This difference is consistent with the |ΔreSE| ∼ 65 mÅ from our CCSD(T) and MP2 + δ(T) results. The correlation between the two data sets of counterpoisecorrected interaction energies is also strong (ρA24, NCCE31 = 0.74). Similarly to the NCDT16 data above, our counterpoisecorrected CCSD(T)/may-cc-pVTZ calculation consistently overpredicts the reference CP-CCSD(T)/CBS binding energies with |ΔE int| = 1.5 kJ/mol. As our extrapolated MP2 + δ(T) results have a much lower |ΔE int| = 0.1 kJ/mol, the difference is attributed to the basis set incompleteness error. The correlation between the covalent bond length data sets (ρCCse22, CS20 = 0.81) is weaker than for the noncovalent bond data sets. The SOGGA11-X39 method has been previously highlighted as an excellent functional for covalent bond lengths.40 However, for noncovalent bond lenghts, its performance is poor

Figure 5. Correlation diagrams for performance in prediction of noncovalent bonds (|ΔreSE|, NCDT16 data set): (A) correlation with MAE of counterpoise-corrected interaction energies (|ΔE int|, A24); SE |, CCse22). The (B) correlation with MAE of covalent bonds (|Δrcov dashed line shows average performance of the data set. For clarity, only selected methods are shown. Colors are as in Figure 1.

data set are listed, and only the best-performing dispersion correction is shown for each functional, if it is better than the parent functional. While the B97M-V23 functional is the best studied method for noncovalent bond lengths, the related range-separated ωB97M-V29 offers a more balanced performance for accurate interaction energies (Figure 5A, |ΔE int| = 0.3

(|ΔreSE| > 750 mÅ, 300 mÅ with D3(BJ)) even when compared to nondispersion corrected single hybrids, such as PBE0 (|ΔreSE| = 182 mÅ). The lack of transferability between errors in covalent and noncovalent bonds is also confirmed by the low correlation coefficient ρNCDT16,CCse22 = 0.21. On the other hand, the comparably high correlation coefficient between interaction energies and noncovalent bonds (ρNCDT16,A24 = 0.84) shows that accurate description of energetics generally corresponds to accurate geometriesthe notable exceptions with poorly predicted geometries despite |ΔE int| < 1 kJ/mol are ωB97XD,41 B97-D3,6,42 TPSS-D3(BJ),27,43 and M08-HX(-D3).44,45 Curiously, the two half-and-half hybrid functionals25 (BH and H with 50% Slater exchange and BH and HLYP with 50% B88 exchange) perform much better than both three-parameter functionals (B3PW9125 and B3LYP25,26). This is due to the

SE |=4 kJ/mol) and covalent bond lengths (Figure 5B, |Δrcov mÅ). The most balanced double hybrid functional remains DSD-PBEPBE-D3(BJ),24 as its somewhat large |ΔE int| of 0.5 kJ/mol in the A24 data set is offset by its remarkably good performance for the NCCE31 data set. For nonmeta-GGA single hybrid functionals, good overall functionals are the rangeseparated ωB97X-V28 and global hybrid B3LYP-D3(BJ).25−27 With the exception of B97M-V, the meta-GGA functionals achieve worse results than dispersion-corrected GGA functionals. The most balanced GGA functional is HCTH/120-

4898


Article


describe dispersion interactions well, the inclusion of D2

SE D3(BJ)27,46 with |ΔE int| = 1.3 kJ/mol and |Δrcov | = 9 mÅ, as the dispersion-corrected forms of PBE and BLYP have comparably large deviation in covalent bond lengths

correction5 already reduces |ΔreSE| at least by a factor of 7. However, the D3 and D3(BJ) corrections6,27 perform much better, with further reduction in MAE’s by a factor of 2. The the revised version (D3m(BJ))48 shows no systematic improvement compared to the original Becke−Johnson-damped version (D3(BJ)). However, the most accurate method for improving noncovalent bond lengths is the nonlocal correlation in VV10,8 with 4 out of 5 methods including this component performing better than CCSD(T). In light of this, the recently developed spin-component scaled double hybrids with nonlocal correlation in place of the empirical dispersion term (ie. DSDPBEPBE-NL) are likely to achieve good results;49 however, the availability of such methods is limited. On the other hand, the nonlocal correlation approach has a higher computational cost than empirical dispersion corrections, especially as the implementation of analytical gradients for these functionals is not widespread. It would be interesting to compare the current results to methods including with the recently developed D4 correction;7 unfortunately, their availability is also currently limited. The importance of developing density functionals with systematic consideration of dispersion corrections can also be highlighted. The dispersion-corrected spin-component-scaled double hybrids of Kozuch and Martin24,50 as well as the recent functionals of Mardirossian and Head-Gordon23,28,29 were systematically parametrized with dispersion corrections, and as a result show good performance for dispersion dominated complexes. However, some dispersionless functionals trained on data sets including noncovalent interactions, such as the functionals from the Truhlar group,30,33,34,51 are excellent for H-bonded systems even without empirical dispersion terms. While applying such further corrections to these functionals is possible,36,45 it does not lead to a systematic improvement of noncovalent geometries.

SE (|Δrcov | > 16 mÅ). 3.3. Climbing Jacob’s Ladder. The large amount of collected data allows us to investigate the effects of systematic improvement of DFT methods by climbing the so-called “Jacob’s ladder”. The resulting average errors in noncovalent bond lengths and interaction energies are shown in Figure 6. In

Figure 6. Performance of related functionals in noncovalent interactions: (A) bond lengths; (B) interaction energies. Color denotes functional family (green, BLYP; red, PBE; blue, B97; cyan, meta-B97) and dispersion corrections (yellow, empirical; black, VV10). Vertical lines delineate GGA’s from global single hybrids, range-separated hybrids, double hybrids, and spin-component scaled double hybrids.

the BLYP-based series (BLYP-D3(BJ), B3LYP-D3(BJ), B2PLYP-D3(BJ), and DSD-BLYP-D3(BJ); green), a significant improvement in prediction of the structure is achieved by inclusion of HF exchange, while the inclusion of MP2 correlation shows mixed results. In the PBE-based series (PBE-D3(BJ), PBE0-D3(BJ), ωPBE-D3(BJ), LC-VV10, PBE02, and DSD-PBEPBE-D3(BJ); red), inclusion of HF-exchange leads to a better prediction of interaction energies in addition to the structural parameters, especially in the range-separated form and with inclusion of the VV10 nonlocal correlation. The dispersionless PBE0-2 double hybrid functional performs well for structures, but interaction energies are not improved from the related single hybrids, as a consequence of the lack of dispersion corrections. In the B97-based data (blue), the results are mixed. Without dispersion corrections, the range-separated hybrid ωB97X performs quite poorly for interaction energies worse than the original B97-D2 form. The empirical dispersion correction included in ωB97X-D improves the agreement in interaction energies but significantly degrades the structural performance; inclusion of the VV10 nonlocal correlation in ωB97X-V remedies this problem. 3.4. Role of Dispersion Corrections. It is widely accepted47 that correctly scaled dispersion corrections can greatly improve the agreement in noncovalent interaction energies. For noncovalent bond lengths, near-MP2 quality can be achieved with dispersion corrected GGAs.27 This statement is further confirmed by our overall results (cf. Figure 1) and is especially pronounced for π−π interactions (cf. Figure 2). For the “worst case” of BLYP and B3LYP functionals, which do not

4. CONCLUSIONS In the current work, we have benchmarked 89 ab initio methods using the NCDT16 data set11 to establish their performance for prediction of noncovalent bond lengths, and five other reference data sets to analyze correlations with accuracy for interaction energies and covalent bonds. The best overall method is the ωB97M-V range-separated meta-GGA hybrid functional of Mardirossian and Head-Gordon,29 thanks to its balanced performance for noncovalent structures, interaction energies and covalent bond lengths. The best results in noncovalent structures were achieved by the related meta-GGA B97M-V.23 Other recommended methods include the single hybrid B3LYP-D3(BJ)25−27 with excellent results for complexes with binding energies below 10 kJ/mol and its widespread availability, and the double hybrid DSD-PBEPBED3(BJ)24 which is considerably more accurate in interaction energies. For π−π interactions, the B97-D3 method6,42 performs rather well, while for H-bonded complexes the newer Minnesota functionals from the MN12 family32,34 and MN15-D3(BJ)33,45 are excellent performers. The results also highlight that the ranking of a method for covalent bonds is not necessarily transferable to noncovalent bonds, while a low average error in interaction energies generally correlates with the quality of prediction of noncovalent bond lengths. Additionally, at the highest level of theory used in the current work (MP2/CBS + δCCSD(T)/34899


Article

The Journal of Physical Chemistry A ζ), the differences in interaction energies between the semiexperimental and optimized geometries can be as large as 40% (HF−H2CO complex). The MAE between the two sets of structures of 62 mÅ confirms the importance of validation of computational methods against experimental observables. Furthermore, systematic improvements of the density functional approximations by range-separation, MP2 correlation and empirical dispersion generally improve the agreement with reference interaction energies. However, the most consistent functionals for noncovalent bond lengths contain the nonlocal correlation from the VV108 functional, and are systematically parametrized with the dispersion terms taken into account.

■

correction (DFT-D) for the 94 elements H-Pu. J. Chem. Phys. 2010, 132, 154104. (7) Caldeweyher, E.; Bannwarth, C.; Grimme, S. Extension of the D3 dispersion coefficient model. J. Chem. Phys. 2017, 147, 034112. (8) Vydrov, O. A.; Van Voorhis, T. Nonlocal van der Waals density functional: The simpler the better. J. Chem. Phys. 2010, 133, 244103. (9) Grimme, S. Semiempirical hybrid density functional with perturbative second-order correlation. J. Chem. Phys. 2006, 124, 034108. (10) Chai, J. D.; Mao, S. P. Seeking for reliable double-hybrid density functionals without fitting parameters: The PBE0-2 functional. Chem. Phys. Lett. 2012, 538, 121−125. (11) Kraus, P.; Obenchain, D. A. D.; Frank, I. Benchmark-quality semiexperimental structural parameters of van der Waals complexes. J. Phys. Chem. A 2018, 122, 1077−1087. (12) Grimme, S.; Brandenburg, J. G.; Bannwarth, C.; Hansen, A. Consistent structures and interactions by density functional theory with small atomic orbital basis sets. J. Chem. Phys. 2015, 143, 054107. (13) Piccardo, M.; Penocchio, E.; Puzzarini, C.; Biczysko, M.; Barone, V. Semi-experimental equilibrium structure determinations by employing B3LYP/SNSD anharmonic force fields: Validation and application to semirigid organic molecules. J. Phys. Chem. A 2015, 119, 2058−2082. (14) Peverati, R.; Truhlar, D. G. Quest for a universal density functional: the accuracy of density functionals across a broad spectrum of databases in chemistry and physics. Philos. Trans. R. Soc., A 2014, 372, 20120476. (15) Parrish, R. M.; et al. Psi4 1.1: An open-source electronic structure program emphasizing automation, advanced libraries, and interoperability. J. Chem. Theory Comput. 2017, 13, 3185−3197. (16) Lehtola, S.; Steigemann, C.; Oliveira, M. J.; Marques, M. A. Recent developments in LIBXC − A comprehensive library of functionals for density functional theory. SoftwareX 2018, 7, 1−5. (17) Rappoport, D.; Furche, F. Property-optimized Gaussian basis sets for molecular response calculations. J. Chem. Phys. 2010, 133, 134105. (18) Hellweg, A.; Rappoport, D. Development of new auxiliary basis functions of the Karlsruhe segmented contracted basis sets including diffuse basis functions (def2-SVPD, def2-TZVPPD, and def2-QVPPD) for RI-MP2 and RI-CC calculations. Phys. Chem. Chem. Phys. 2015, 17, 1010−1017. (19) Papajak, E.; Zheng, J.; Xu, X.; Leverentz, H. R.; Truhlar, D. G. Perspectives on basis sets beautiful: Seasonal plantings of diffuse basis functions. J. Chem. Theory Comput. 2011, 7, 3027−3034. (20) Halkier, A.; Helgaker, T.; Jørgensen, P.; Klopper, W.; Olsen, J. Basis-set convergence of the energy in molecular Hartree−Fock calculations. Chem. Phys. Lett. 1999, 302, 437−446. (21) Kendall, R. A.; Dunning, T. H.; Harrison, R. J. Electron affinities of the first-row atoms revisited. Systematic basis sets and wave functions. J. Chem. Phys. 1992, 96, 6796−6806. (22) Woon, D. E.; Dunning, T. H. Gaussian basis sets for use in correlated molecular calculations. X. The atoms aluminum through argon revisited. J. Chem. Phys. 1993, 98, 1358−1371. (23) Mardirossian, N.; Head-Gordon, M. Mapping the genome of meta-generalized gradient approximation density functionals: The search for B97M-V. J. Chem. Phys. 2015, 142, 074111. (24) Kozuch, S.; Martin, J. M. L. Spin-component-scaled double hybrids: An extensive search for the best fifth-rung functionals blending DFT and perturbation theory. J. Comput. Chem. 2013, 34, 2327−2344. (25) Becke, A. D. A new mixing of Hartree−Fock and local densityfunctional theories. J. Chem. Phys. 1993, 98, 1372−1377. (26) Stephens, P. J.; Devlin, F. J.; Chabalowski, C. F.; Frisch, M. J. Ab initio calculation of vibrational absorption and circular dichroism spectra using density functional force fields. J. Phys. Chem. 1994, 98, 11623−11627. (27) Grimme, S.; Ehrlich, S.; Goerigk, L. Effect of the damping function in dispersion corrected density functional theory. J. Comput. Chem. 2011, 32, 1456−1465.

ASSOCIATED CONTENT

S Supporting Information *

The Supporting Information is available free of charge on the ACS Publications website at DOI: 10.1021/acs.jpca.8b03345. Detailed references for all functionals (PDF) Spreadsheet with the MAE values for all six data sets and all complexes in the NCDT16 data set (XLSX) Archive of PYTHON driver routines for performing the benchmarks (ZIP) Archive of optimized structures with selected methods (ZIP)

■

AUTHOR INFORMATION

Corresponding Author

*(P.K.) E-mail: [email protected]. ORCID

P. Kraus: 0000-0002-4359-5003 Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS P.K. would like to thank Dr. S. Lehtola for his help with LIBXC, Dr. D. G. A. Smith and Dr. L. A. Burns for their help with ψ4, and Dr. D. A. Obenchain for his support and helpful comments. This work was partially carried out on the Leibniz Universität Hannover computer cluster, which is funded by the Leibniz Universität Hannover, the Lower Saxony Ministry of Science and Culture (MWK), and the German Research Association (DFG).

■

REFERENCES

(1) Zhao, Y.; Truhlar, D. G. Benchmark databases for nonbonded interactions and their use to test density functional theory. J. Chem. Theory Comput. 2005, 1, 415−432. (2) Hohenstein, E. G.; Sherrill, C. D. Wavefunction methods for noncovalent interactions. Wiley Interdiscip. Rev. Comput. Mol. Sci. 2012, 2, 304−326. (3) Ř ezác,̌ J.; Hobza, P. Describing noncovalent interactions beyond the common approximations: How accurate is the “gold standard,” CCSD(T) at the complete basis set limit? J. Chem. Theory Comput. 2013, 9, 2151−2155. (4) Ř ezác,̌ J.; Hobza, P. Benchmark calculations of interaction energies in noncovalent complexes and their applications. Chem. Rev. 2016, 116, 5038−5071. (5) Grimme, S. Semiempirical GGA-type density functional constructed with a long-range dispersion correction. J. Comput. Chem. 2006, 27, 1787−1799. (6) Grimme, S.; Antony, J.; Ehrlich, S.; Krieg, H. A consistent and accurate ab initio parametrization of density functional dispersion 4900


Article

The Journal of Physical Chemistry A (28) Mardirossian, N.; Head-Gordon, M. ωB97X-V: A 10-parameter, range-separated hybrid, generalized gradient approximation density functional with nonlocal correlation, designed by a survival-of-thefittest strategy. Phys. Chem. Chem. Phys. 2014, 16, 9904. (29) Mardirossian, N.; Head-Gordon, M. ωB97M-V: A combinatorially optimized, range-separated hybrid, meta-GGA density functional with VV10 nonlocal correlation. J. Chem. Phys. 2016, 144, 214110. (30) Zhao, Y.; Truhlar, D. G. The M06 suite of density functionals for main group thermochemistry, thermochemical kinetics, noncovalent interactions, excited states, and transition elements: Two new functionals and systematic testing of four M06-class functionals and 12 other function. Theor. Chem. Acc. 2008, 120, 215−241. (31) Goerigk, L.; Grimme, S. A thorough benchmark of density functional methods for general main group thermochemistry, kinetics, and noncovalent interactions. Phys. Chem. Chem. Phys. 2011, 13, 6670−6688. (32) Peverati, R.; Truhlar, D. G. An improved and broadly accurate local approximation to the exchange-correlation density functional: The MN12-L functional for electronic structure calculations in chemistry and physics. Phys. Chem. Chem. Phys. 2012, 14, 13171− 13174. (33) Yu, H. S.; He, X.; Li, S. L.; Truhlar, D. G. MN15: A Kohn-Sham global-hybrid exchange-correlation density functional with broad accuracy for multi-reference and single-reference systems and noncovalent interactions. Chem. Sci. 2016, 7, 5032−5051. (34) Peverati, R.; Truhlar, D. G. Screened-exchange density functionals with broad accuracy for chemistry and solid-state physics. Phys. Chem. Chem. Phys. 2012, 14, 16187−16191. (35) Peverati, R.; Truhlar, D. G. Exchange-correlation functional with good accuracy for both structural and energetic properties while depending only on the density and its gradient. J. Chem. Theory Comput. 2012, 8, 2310−2319. (36) Goerigk, L. Treating London-dispersion effects with the latest Minnesota density functionals: Problems and possible solutions. J. Phys. Chem. Lett. 2015, 6, 3891−3896. (37) Mentel, M.; Baerends, E. J. Can the counterpoise correction for basis set superposition effect be justified? J. Chem. Theory Comput. 2014, 10, 252−267. (38) Roch, L. M.; Baldridge, K. K. Dispersion-corrected spincomponent-scaled double-hybrid density functional theory: Implementation and performance for non-covalent interactions. J. Chem. Theory Comput. 2017, 13, 2650−2666. (39) Peverati, R.; Truhlar, D. G. Communication: A global hybrid generalized gradient approximation to the exchange-correlation functional that satisfies the second-order density-gradient constraint and has broad applicability in chemistry. J. Chem. Phys. 2011, 135, 191102. (40) Brémond, É.; Savarese, M.; Su, N. Q.; Pérez-Jiménez, Á . J.; Xu, X.; Sancho-García, J. C.; Adamo, C. Benchmarking density functionals on structural parameters of small-/medium-sized organic molecules. J. Chem. Theory Comput. 2016, 12, 459−465. (41) Chai, J.-D.; Head-Gordon, M. Long-range corrected hybrid density functionals with damped atom-atom dispersion corrections. Phys. Chem. Chem. Phys. 2008, 10, 6615−6620. (42) Becke, A. D. Density-functional thermochemistry. V. Systematic optimization of exchange-correlation functionals. J. Chem. Phys. 1997, 107, 8554−8560. (43) Tao, J.; Perdew, J. P.; Staroverov, V. N.; Scuseria, G. E. Climbing the density functional ladder: Non-Empirical meta-generalized gradient approximation designed for molecules and solids. Phys. Rev. Lett. 2003, 91, 146401. (44) Zhao, Y.; Truhlar, D. G. Exploring the limit of accuracy of the global hybrid meta density functional for main-group thermochemistry, kinetics, and noncovalent interactions. J. Chem. Theory Comput. 2008, 4, 1849−1868. (45) Goerigk, L.; Hansen, A.; Bauer, C. A.; Ehrlich, S.; Najibi, A.; Grimme, S. A look at the density functional theory zoo with the advanced GMTKN55 database for general main group thermochem-

istry, kinetics and noncovalent interactions. Phys. Chem. Chem. Phys. 2017, 19, 32184−32215. (46) Boese, A. D.; Doltsinis, N. L.; Handy, N. C.; Sprik, M. New generalized gradient approximation functionals. J. Chem. Phys. 2000, 112, 1670−1678. (47) Kruse, H.; Goerigk, L.; Grimme, S. Why the standard B3LYP/631G* model chemistry should not be used in DFT calculations of molecular thermochemistry: Understanding and correcting the problem. J. Org. Chem. 2012, 77, 10824−10834. (48) Smith, D. G. A.; Burns, L. A.; Patkowski, K.; Sherrill, C. D. Revised damping parameters for the D3 dispersion correction to density functional theory. J. Phys. Chem. Lett. 2016, 7, 2197−2203. (49) Kesharwani, M. K.; Karton, A.; Martin, J. M. Benchmark ab initio conformational energies for the proteinogenic amino acids through explicitly correlated methods. Assessment of density functional methods. J. Chem. Theory Comput. 2016, 12, 444−454. (50) Kozuch, S.; Martin, J. M. L. DSD-PBEP86: in search of the best double-hybrid DFT with spin-component scaled MP2 and dispersion corrections. Phys. Chem. Chem. Phys. 2011, 13, 20104−20107. (51) Peverati, R.; Truhlar, D. G. Improving the accuracy of hybrid meta-GGA density functionals by range separation. J. Phys. Chem. Lett. 2011, 2, 2810−2817.

4901


Density Functional Theory for Microwave Spectroscopy of Non

Recommend Documents