Article Cite This: J. Phys. Chem. A XXXX, XXX, XXX−XXX
pubs.acs.org/JPCA
Density Functional Theory and the Basis Set Truncation Problem with Correlation Consistent Basis Sets: Elephant in the Room or Mouse in the Closet? David Feller*,† and David A. Dixon‡ †
Department of Chemistry, Washington State University, Pullman, Washington 99164-4630, United States Department of Chemistry, The University of Alabama, Shelby Hall, Box 870336, Tuscaloosa, Alabama 35487-0336, United States
‡
S Supporting Information *
ABSTRACT: Two recent papers in this journal called into question the suitability of the correlation consistent basis sets for density functional theory (DFT) calculations, because the sets were designed for correlated methods such as configuration interaction, perturbation theory, and coupled cluster theory. These papers focused on the ability of the correlation consistent and other basis sets to reproduce total energies, atomization energies, and dipole moments obtained from “quasi-exact” multiwavelet results. Undesirably large errors were observed for the correlation consistent basis sets. One of the papers argued that basis sets specifically optimized for DFT methods were “essential” for obtaining high accuracy. In this work we reexamined the performance of the correlation consistent basis sets by resolving problems with the previous calculations and by making more appropriate basis set choices for the alkali and alkaline-earth metals and second-row elements. When this is done, the statistical errors with respect to the benchmark values and with respect to DFT optimized basis sets are greatly reduced, especially in light of the relatively large intrinsic error of the underlying DFT method. When judged with respect to high-quality Feller-Peterson-Dixon coupled cluster theory atomization energies, the PBE0 DFT method used in the previous studies exhibits a mean absolute deviation more than a factor of 50 larger than the quintuple zeta basis set truncation error.
I. INTRODUCTION S. R. Jensen et al.1 recently compared the performance of Gaussian-type orbitals (GTOs), numeric atom-centered centered orbitals (NAOs), and augmented plane wave methods (APW) against a collection of benchmark quality, “quasi-exact” density functional theory (DFT) multiwavelet (MW) results for 211 first- and second-row molecules. The MW total energies were reported to be accurate to ∼1 μEh, although we know of no independent verification of this claim. Comparisons were limited to total energies, atomization energies, and dipole moments. Three DFT variants were selected from among the hundreds in the literature: the local density SVWN method,2,3 the generalized gradient approximation Perdew−Burke− Ernzerhof (PBE),4 and hybrid PBE05 methods. The first of the three comparison criteria (total energies) is seldom of interest to chemists. While in principle the total energy of an atom can be deduced from the sum of all ionization potentials or a molecule’s total energy can be based on the sum of all bond dissociation energies, the complete set of experimental data is seldom available with sufficient accuracy to yield an accurate value. Nor is high precision in total energies a prerequisite for high accuracy in thermochemical or spectroscopic properties due to exploitable cancellation of error. Two molecules (CH3CH2O and CCH) were excluded from the © XXXX American Chemical Society
PBE0 test set results due to convergence problems. The title of the S. R. Jensen et al.1 paper includes the phrase “The Elephant in the Room of Density Functional Theory”. While this metaphor was never explicitly explained, it appears to refer to the role of GTO basis set truncation error in DFT calculations. Although most of the S. R. Jensen et al.1 comparisons were performed using up through the aug-cc-pV5Z basis set, for Li, Be, Na, and Mg, the authors substituted the aug-cc-pVQZ basis set, because they were unaware of the availability of the appropriate quintuple zeta sets. S. R. Jensen et al.1 conclude that “moderate sized” GTO basis sets display errors in total energies that, on average, exceed 10 kcal/mol and that even “very large” basis sets show significant outliers. The GTO calculations were performed with NWChem.6 A second, even more recent “elephant” paper based on the same test set of molecules compared the performance of the DFT-optimized polarization consistent (pc) basis sets with the correlation consistent sets using the same three criteria employed in the earlier paper. In the second paper, entitled “How Large is the Elephant in the Density Functional Room?”, Received: January 12, 2018 Revised: February 19, 2018 Published: February 20, 2018 A
DOI: 10.1021/acs.jpca.8b00392 J. Phys. Chem. A XXXX, XXX, XXX−XXX
Article
The Journal of Physical Chemistry A F. Jensen7 concludes that “standard” Gaussian basis sets, meaning not optimized for DFT, are “unsuitable for achieving high accuracy”. In a study of noncovalent intermolecular interactions, Witte et al.8 remarked that the correlation consistent basis set sequences “converge remarkably slowly to the basis set limit for DFT.” However, it should be pointed out that the study of Witte suffers from the same poor choice of correlation consistent basis sets for second-row elements as the two “elephant” papers. It is not clear from the publication whether F. Jensen7 adopted any of the correlation consistent results from the work of S. R. Jensen et al.1 or if they were recomputed. The reported mean absolute deviations (MADs) in the two papers are similar to one another but not identical. For example, the F. Jensen7 and S. R. Jensen et al.1 PBE0/ aVQZ MAD values are 0.709 and 0.582 kcal/mol, respectively. Although F. Jensen7 does not report aug-cc-pV5Z statistics, at the cc-pV5Z level the pcseg-4 basis set outperforms the correlation consistent set by almost a factor of 3 in terms of the MAD, 0.167 (V5Z) versus 0.058 (pcseg4). In much earlier work, Bauschlicher and Partridge studied the sensitivity of DFT calculations to the choice of basis set.9 Their investigation compared the ability of a variety of Pople-style and correlation consistent basis sets (aug-cc-pVTZ/aug-ccpV(D+d)Z) to reproduce experimental, zero-point-inclusive atomization energies for 55 first- and second-row compounds. The authors concluded that, for the empirically adjusted, hybrid B3LYP10,11 functional, the 6-311+G(3df,2p) basis set slightly outperformed the aug-cc-pVTZ (2.20 vs 2.59 kcal/mol MAD), but it was not obvious if atomic spin−orbit corrections were included when computing the atomization energies. The experimental uncertainty for some of the molecules exceeded 3 kcal/mol, thus limiting the significance of the statistical analysis. Several years later, Raymond and Wheeler12 reexamined using B3LYP with a wider range of correlation consistent basis sets (aug-cc-pVnZ, n = D, T,Q). They focused on the ability of B3LYP to reproduce experimental structures, atomization energies, dipole moments, electron affinities, and harmonic frequencies for a small collection of di- and triatomics. The MAD for atomization energies relative to experiment obtained with the aug-cc-pVQZ basis set was 4.6 kcal/mol, well outside the nominal chemical accuracy limit of ±1 kcal/mol. An attempt was made to extrapolate the properties to the complete basis set (CBS) limit, but for atomization energies the convergence failed to exhibit uniform, monotonic convergence. Martell et al. studied six gradientcorrected functionals with two Pople-style basis sets and two correlation consistent basis sets (cc-pVDZ and cc-pVTZ).13 MAD values for atomization energies across a collection of 44 molecules fell into the 3.1 to 14.2 kcal/mol range. Boese et al.14 assessed the impact of basis set selection on a modified Hamprecht-Cohen-Tozer-Handy (HCTH)15 functional, as well as a long list of other functionals, using a combination of Poplestyle and correlation consistent basis sets up through cc-pVQZ. They recommended the use of Pople-style basis sets over the correlation consistent sets. The older Dunning TZ2P basis set16 provided the lowest overall error when using generalized gradient approximation (GGA) functionals. Experimental atomization energies, electron affinities, proton affinities, and ionization potentials were used to evaluate the basis sets and functionals, with root mean squared (RMS) errors falling into the 4.6−46.9 kcal/mol range. In a series of seven papers, Wilson and co-workers extensively examined the correlation consistent basis set family
in the context of DFT calculations.17−23 In the first of these, they focused on the atomization energies of three small secondrow molecules (SO2, ClO2, and CCl) with the B3LYP and B3PW9124 functionals and concluded, in agreement with Bauschlicher and Partridge,9 that tight d functions are important for accelerating basis set convergence. Smooth convergence to the basis set limit was found. The next paper in the series examined the performance of the correlation consistent basis sets up through aug-cc-pV5Z on a collection of 17 atomization energies for closed-shell molecules using six functionals. As with many of the earlier studies, experimental values were used as the reference point. Smooth convergence to the basis set limit was again observed. The later papers in this series stressed: (1) the importance of tight d functions for second-row molecules, (2) that the cc-pVnZ correlation consistent basis sets perform about as well as the pc-n basis sets for B3LYP atomization energies using experimental reference values, (3) that correcting for basis set superposition error can help in cases where irregular basis set convergence is observed, and (4) that recontraction of the correlation consistent basis sets can also help ensure smooth convergence in DFT calculations. The goal of the present work is to re-examine the performance of the correlation consistent basis sets for DFT atomization energies using one of the functions employed in the two “elephant” papers. To isolate the basis set truncation factor from multiple other issues involved with comparisons to experiment, we will adopt the test set of molecules and the MW results of S. R. Jensen et al.1 as our reference. The MW results were obtained with the MRChem program.
II. APPROACH PBE0 calculations were performed with Gaussian 0925 using the default “fine” grid. Tests of the “ultrafine” grid produced total energies that differed by ∼1 μEh. For very large basis sets, it is necessary to disable the default Davidson-style26 basis transformation in G09 intended to improve the speed of the integrals evaluation (NoBasisTransform), as well as to request higher accuracy in the two-electron integrals (Acc2E = 14) to achieve consistent results. If Gaussian senses near linear dependency in the basis set, it performs a rectangular transformation to remove eigenfunctions of the overlap matrix with very small eigenvalues. We decreased the threshold to 1 × 10−7 (IOp(3/59 = 7)). Failure to change the program defaults can result in energies that are wrong by as much as 1 mEh. It is not entirely clear whether the definition of the PBE0 functional in Gaussian corresponds exactly to the definition used in the MRChem27 program used for obtaining the MW results. To limit the scope of this project, our analysis will be limited to the PBE0 functional, unlike the two “elephant” papers,1,7 which also included the PBE functional. It is unlikely that the conclusions reached here will be limited only to this functional. With aV(7+d)Z quality contracted basis sets, the Gaussian 09 total energies still differed by ∼1 mEh from the MW energies, an amount approximately half as large as the difference observed with the aug-pc-4 basis sets. NWChem and G09 energies obtained with the same basis set differed by ∼0.1 mEh, indicating the limit of reproducibility between these codes. Uncontracting the 7 zeta quality basis sets brought the G09 energies to within 0.1 mEh (0.06 kcal/mol) of the MW values. These uncontracted basis sets required the use of the Gaussian SuperFineGrid. B
DOI: 10.1021/acs.jpca.8b00392 J. Phys. Chem. A XXXX, XXX, XXX−XXX
Article
The Journal of Physical Chemistry A
Table 1. Statistical PBE0 Atomization Energy Error Metricsa
Statistical analyses were limited to the diffuse function augmented correlation consistent basis sets of Dunning, Peterson, and co-workers.28−48 These basis sets are conventionally denoted aug-cc-pVnZ, n = D, T, Q, ...10, but for the sake of brevity we will shorten the labels to aVnZ. As is the convention with the correlation consistent basis sets, only the spherical forms of the basis functions were used. Note, in particular, that alkali and alkaline-earth metal basis sets for Li, Be, Na, and Mg up through aV5Z are available and were taken from Pascher et al.49 Second-row basis sets with an additional tight d function are denoted aV(n+d)Z. On the basis of the work cited above, it is expected that tight d functions will be important for second-row molecules. Calculations were performed through at least the aV5Z/aV(5+d)Z basis sets. In some cases basis sets up through 7 zeta were used. Comparisons between the correlation consistent aVnZ/aV(n +d)Z basis sets and the polarization consistent pc-n basis sets are complicated by the use of different numbers of functions at each basis set index (n) level. For example, the aug-pcseg-4 first-row basis sets consist of a (21s 13p 7d 4f 3g 2h) primitive set contracted to [7s 6p 7d 4f 3g 2h], while the aV5Z basis set is (15s 8p 5d 4f 3g 2h) contracted to [7s 6p 5d 4f 3g 2h]. Thus, the aug-pcseg-4 basis set would be expected to yield results slightly closer to the basis set limit than the aV5Z basis set due to the inclusion of more d functions in the former. The aV(n +d)Z sets only add a single tight d function. All geometries were taken from the work of S. R. Jensen et al.1 While the two “elephant” papers also included an examination of dipole moments, in this study we restrict ourselves to considering atomization energies only, since they are among the most difficult properties for electronic structure methods to reproduce accurately. To obtain a better estimate of the performance of correlation consistent basis sets, the test set molecules with the largest 50 errors relative to the MW results were re-examined. A thorough review of the largest aV5Z atomization energy errors reveals that all were associated with molecules containing second-row elements or alkali and alkaline-earth metals. First-row compounds usually differed by less than 0.1 kcal/mol. PBE0 total energies and atomization energies for basis sets VnZ, aVnZ, and aV(n+d)Z, n = D, T, Q, and 5 are provided in Supporting Information tables. A limited number of aV(6+d)Z and aV(7+d)Z results are also listed.
basis set
MAD
RMS
abs max
aV5Zb
0.211
0.548
5.223
aV5Zb aV(T+d)Ze aV(Q+d)Ze aV(5+d)Ze CBS(aVQ5Z)
0.151 1.153 0.249 0.082 0.079 (0.065)f
0.303 1.375 0.437 0.116 0.137 (0.084)
aug-pcseg-4
0.053 (0.063)
aV(6+d)Zh
0.029
2.281 3.841 2.639 0.509 0.984 (0.218) 0.23 (0.23) 0.104
0.046
source S. R. Jensen et al.c this workd this work this work this work this work F. Jenseng this work
a
In kilocalories per mole, on the basis of the multiwavelet reference set of Jensen et al.1 bNumber of comparisons = 210. cS. R. Jensen et al.1 d Replacing the results for Al2 and Si2 and using the aV5Z basis sets for alkali and alkaline-earth elements. eNumber of aV(n+d)Z values = 47. f Values in parentheses correspond to ignoring Mg compounds, where the CBS extrapolation greatly overestimates the atomization energies. g F. Jensen.7 The values in parentheses correspond to complete basis set extrapolations. hNumber of aV(6+d)Z comparisons = 7.
Figure 1. Atomization energy mean absolute deviations (kcal/mol).
basis set level, replacing the erroneous S. R. Jensen et al.1 values for Al2 and Si2, as well as using the correct aV5Z alkali metal and alkaline-earth basis sets, reduced the MAD from 0.211 to 0.151 kcal/mol. The impact of tight d functions on the atomization energies is significant, in part due to the presence of 81 main-group second-row compounds in the test set. As mentioned earlier, tight d functions were found to be important by Bauschlicher and Partridge9 and by Wang and Wilson.17 For the cases in which the aVnZ basis set was replaced by aV(n+d)Z, the MAD is reduced by almost a factor of 2. Note that the aV(n+d)Z MAD and RMS values listed in Table 1 almost certainly represent upper bounds on the values corresponding to a complete replacement of all second-row data. Applying a simple two-point CBS extrapolation formula50−54
III. RESULTS AND DISCUSSION Several of the PBE0/aV5Z and aVQZ atomization errors appearing in the Supporting Information (SI) file attached to the S. R. Jensen et al.1 paper appeared unexpectedly large. For example, the absolute errors for Al2 (3Πu), Si2 (3Σg−), CH2NHO-zwitter (1A), and CH2NHO-s (1A′) were reported to be 3.583 (aV5Z), 5.223 (aV5Z), 7.074 (aVQZ), and 3.584 (aVQZ) kcal/mol, respectively. In our own calculations we find the absolute errors to be 0.007, 0.056, 0.032, and 0.051 kcal/ mol, respectively. Other randomly selected aV5Z atomization energies agreed with the present results within the uncertainty of the potentially slightly different definitions of the PBE0 functional in NWChem and Gaussian 09. It was also noticed that the geometry for CH3CN (acetonitrile) contained an error related to the geometrical parameters of the methyl group. Table 1 lists the PBE0 atomization energy statistical error metrics (MAD, root-mean-square (RMS), and absolute maximum deviations) obtained at various levels of theory. These results are displayed graphically in Figure 1. At the aV5Z
E(Smax) = ECBS + A /(Smax + 1/2)4
(1)
can reduce the MAD even further if it is judiciously applied. Our experience indicates that the formula works well for basis sets in the aVTZ to aV6Z range with coupled cluster theory or configuration interaction energies. However, it has not been calibrated for DFT, and its sphere of applicability may be C
DOI: 10.1021/acs.jpca.8b00392 J. Phys. Chem. A XXXX, XXX, XXX−XXX
Article
The Journal of Physical Chemistry A
additional convergence beyond aV(5+d)Z is possible. While other aV(6+d)Z and aV(7+d)Z calculations would have been possible, it was felt that additional cases were not needed for making this point. For example, in the case of S2 (3Σg−) the aV(7+d)Z basis set yielded a dissociation energy that differed by only 0.001 kcal/mol from the MW value. For Al2 (3Πu) the difference was 0.007 kcal/mol. Altering the Gaussian 09 program defaults as described previously should enable calculations at the aV(6+d)Z level for most of the systems in the test set, although that was not attempted for this study. Finally, the impact of core/valence functions was tested in a series of calculations on S2 (3Σg−) that added the tight weighted core valence (CV) functions to the aV(n+d)Z basis sets. Although the improvements in the dissociation energies were uniform across the series of basis sets, they were relatively modest at ∼0.1 kcal/mol. Given the large number of CV functions it is more economical to increase the basis set index, n, if higher accuracy is desired, rather than add CV functions. Compared to orbital-based methods such as MP2 or CCSD(T), the addition of tight functions in DFT calculation has only a minor impact.
limited. In the present study, the formula proved effective at estimating the PBE0/CBS limit in most cases. An exception was observed for Mg compounds. The aVnZ sequence of atomization energies failed to exhibit monotonic convergence to the CBS limit, as represented by the MW results. Since the raw aV5Z atomization energies overestimated the MW values, the CBS extrapolation for Mg compounds produced atomization energies that were significantly too large. Eliminating the Mg compounds from the reference set produced the smaller error metrics given in parentheses in Table 1. It is not known if similar behavior would be found with other DFT methods. Using basis sets without the added diffuse functions avoids this problem. Overall, if aVnZ basis sets are used for alkali metals and alkaline-earth elements and tight d function aV(n+d)Z basis sets are used for main-group second-row elements (Al − Ar), we find that the correlation consistent sets are only slightly poorer at reproducing MW-DFT results than the diffuse function augmented, segmented polarization consistent (augpcseg-n) basis sets. The extent to which the larger number of functions contained in the aug-pcseg-n basis sets, compared to the aVnZ sets, contributes to these findings, as opposed to the explicit optimization at the DFT level of theory, is unknown. The small difference between correlation consistent and pc basis sets found in this study is in sharp contrast to the difference indicated by F. Jensen7 (see his Table 2). To give one specific example, the atomization error for SF6 at the aVQZ level is given as 10.18 kcal/mol in both elephant papers. But substituting the aV(Q+d)Z basis set reduces the error to 1.32 kcal/mol. At the aV(5+d)Z level the error is reduced even further, to 0.50 kcal/mol. To assess the intrinsic accuracy of the PBE0 functional, a statistical analysis was performed using high-quality FellerPeterson-Dixon55−57 atomization energies, excluding atomic spin−orbit and zero-point vibrational effects as the best available data. The FPD composite approach involves the use of CCSD(T) with large basis sets, CBS extrapolation, corrections for core/valence and scalar relativistic effects, and inclusion of higher-order correlation recovery beyond CCSD(T). It has been shown capable of reproducing well-established experimental thermochemical properties to an accuracy of less than 0.3 kcal/mol (296 comparisons).58 When the comparison set is restricted to experimental values with uncertainties below 0.2 kcal/mol, the MAD for the FPD results falls to 0.09 kcal/ mol (RMS = 0.13 kcal/mol). The error metrics for PBE0 versus FPD were 3.88 (MAD), 5.62 (RMS), and 18.25 (Max) kcal/ mol. Consequently, the accuracy of the PBE0 functional, which was selected for use in the first elephant paper because it was found to be “relatively accurate for atomization energies”, is found to provide considerably worse than so-called “chemical accuracy” (±1 kcal/mol). For this test set, even the small aV(T +d)Z basis set was found to be capable of sufficient accuracy to match the capabilities of the PBE0 functional relative to FPD results. Pursuing highly precise (from a basis set perspective) results that are ultimately inaccurate (relative to much higher level theory) makes little sense. The presence of systematic errors, even in cases with high precision, can render the nominally precise result of limited value.59 The wide breadth of the correlation consistent sets (up through aV10Z for some elements) allows researchers to reduce the basis set truncation error to very small values in cases where that is desirable. A limited number of aV(6+d)Z and aV(7+d)Z calculations were performed to demonstrate that
IV. CONCLUSIONS By employing the correct aV5Z basis sets for the alkali and alkaline-earth metal compounds and by using basis sets with tight d functions for second-row elements, it has been shown that correlation consistent basis sets perform nearly as well as the larger pc DFT-optimized basis sets for the atomization energies associated with the collection of molecules used in the two earlier studies. Multiple errors were discovered in the results of the S. R. Jensen et al. study.1 On the basis of the current results, the “elephant in the room” can more appropriately be characterized as a mouse. The comparison of the current results with those of the two elephant papers illustrates the importance of using the best available finite basis set values when numerical limit methods, such as the multi wavelet approach, are involved. It is helpful to the computational chemistry community when key input data such as basis sets are in carefully maintained databases, for example, the initial versions of the EMSL database library.60,61 This issue is clearly important for the precision of the calculation as well as for benchmark comparisons whenever basis sets are involved. Dipole moments were not considered in this study. This property is much less sensitive to the quality of the basis set as long as diffuse functions are present. We expect the addition of tight d functions to have a small (
