Accurate DFT-D3 Calculations in a Small Basis Set - Journal of

Article pubs.acs.org/JCTC

Accurate DFT-D3 Calculations in a Small Basis Set Jiří Hostaš†,‡ and Jan Ř ezác*̌ ,† †

Institute of Organic Chemistry and Biochemistry, Czech Academy of Sciences, 166 10 Prague, Czech Republic Department of Chemistry, Institute for Quantum Science and Technology, and Centre for Molecular Simulation, University of Calgary, 2500 University Drive NW, Calgary, Alberta T2N 1N4, Canada

‡

S Supporting Information *

ABSTRACT: Calculations of interaction energies of noncovalent interactions in small basis sets are affected by the basis set superposition error and dispersion-corrected DFT-D methods and are thus usually parametrized only for triple-ζ and larger basis sets. Nevertheless, some smaller basis sets could also perform well. Among many combinations tested, we obtained excellent results with the DZVP-DFT basis and newly parametrized D3 dispersion correction. The accuracy of interaction energies and geometries is close to significantly more expensive calculations.

1. INTRODUCTION In the last decades, the density functional theory (DFT), a method dominating solid-state physics, was brought into the spotlight in the fields of chemistry.1 Substantial increase in the accuracy of both its energetic and geometric descriptions has led to the development of its more sophisticated applications in, e.g., energy decomposition, QM/MM, and vibrational analysis.2−4 Affordable scaling has enabled its use for assemblies up to several hundred atoms, opening space in the areas of catalysis, supramolecular chemistry, molecular biology, or biotechnology.5 With the increasing number of applications to large systems, it has become apparent that an accurate description of London dispersion is of crucial importance.6,7 Nowadays, there have been continuous efforts to improve the description of attractive long-range van der Waals interactions which is well reflected in the number of new DFT functionals addressing this issue published each year.8,9 We have chosen to follow an alternative approach to compensate for the lack of dispersion energy, namely, a posteriori calculated empirical correction term.10 The most successful version of this approach, Grimme’s D3 correction, has been parametrized for use with a large def2-QZVP basis set in order to prevent any problems caused by the basis set superposition error (BSSE).11 Some parameters for this correction exist also for a triple-ζ basis set. It has been assumed that in a smaller basis set the error of the DFT calculation due to the BSSE is too large to yield acceptably accurate interaction energies. This problem is manifested most strongly in hydrogen bonds whose interaction energy can be overestimated by a factor of more than two. Other dispersion corrections have also been parametrized for smaller basis sets, but the results are not satisfactory.12 There is no systematic study of dispersion correction for double-ζ basis © 2017 American Chemical Society

sets, and only several more or less empirical approaches have been developed over the past years.13,14 This includes the 3c corrected HF and DFT methods in a minimal or double-ζ basis which combine the D3 dispersion with an additional correction for the BSSE.14,15 Nevertheless, semiempirical quantum mechanical (SQM) methods with corrections for noncovalent interactions16−19 usually offer comparable or better accuracy than DFT in a small basis set, and their computational cost is lower. However, there are cases where the SQM methods fail and full DFT calculations are needed. The SQM methods are often parametrized for a limited set of elements only, and their high accuracy holds only for well-behaved molecules close to the ones used in the parametrization. The same applies to the DFT-D and HF-D calculations in a minimal basis, notwithstanding additional corrections including a correction for BSSE (as shown below in Section 3.3). It is clear that there is a room for a method, preferably based on DFT, that would be more accurate and more robust while being significantly more efficient than DFT-D in large basis sets. Our experience suggests that such method has to use at least double-ζ basis set. On the other hand, using basis sets larger than that reduces the computational efficiency considerably. To find the best setup, we have searched a library of basis sets ranging from a minimal basis set to split valence (SV) and double-ζ (DZ) basis sets, selected those with the smallest BSSE, and parametrized the D3 correction for them in order to identify the best setup. The best results were achieved with the DZVP-DFT basis which outperforms other basis sets of similar size and more Received: April 5, 2017 Published: July 17, 2017 3575

DOI: 10.1021/acs.jctc.7b00365 J. Chem. Theory Comput. 2017, 13, 3575−3585

Article

Journal of Chemical Theory and Computation approximate methods by a large margin. This basis set was constructed specifically for the use with DFT, and it exhibits significantly smaller BSSE than other basis sets of similar size. For this basis set, we have reparametrized multiple variants of the D3 correction on a larger training set. In addition to the original damping functions used by Grimme (zero damping and Becke−Johnson damping), we also test two recently introduced ones (the optimized power damping20 and C-six-only damping21). In all cases, we parametrized the dispersion correction both with and without the optional three-body term. The resulting methods have then been tested on multiple data sets of benchmark interaction energies as well as on geometries of noncovalent complexes. We have found that the best results were achieved with BLYP and B3LYP functionals, which is in perfect agreement with the much broader survey of functionals in DFT-D3 with large basis set.22 Among the damping functions, the new optimized power (OP) damping works best, but it is closely followed by the Becke−Johnson (BJ) damping which is readily available in any software that implements DFT-D3. Because of the computational efficiency of GGA functionals, we do recommend the BLYP-D3(OP)/DZVP-DFT or BLYP-D3(BJ)/DZVP-DFT setups for general use. They are as accurate as analogous calculations in the triple-ζ basis set, which would be notably more expensive.

R 0AB =

∑

AB n = 6,8,10

Cn ,AB n rAB

fd , n (rAB)

+ (a1R 0AB + a 2)n

(5)

(6)

2.3. Three-Body Contribution. The dispersion is highly additive and a pairwise formulation of the correction can yield very accurate results. However, it can be extended with manybody terms of which the leading three-body (3-body) term is by far the most important. It can be introduced as an extension of the pairwise dispersion correction using the Axilrod−TellerMuto approximation.26,27 It has been included in the D3 correction as an optional extension10 which takes the form of E3‐body = − ∑ fd , n (rAB , rBC , rCA)E ABC ABC

E ABC =

(2)

(7)

C9ABC(3cos θa cos θb cos θc + 1) (rABrBCrCA)3

(8)

where θa, θb, and θc describe the mutual angles between three respective atoms in a triangle ABC with interatomic distances rAB, rBC, and rCA. The fd,n(rAB,rBC,rCA) refers to the zero damping function (eq 2) without any adjustable parameters and the CABC 9 dispersion coefficient is estimated from the C6 coefficients as

n snrAB n rAB

β

r ABn + (a1R 0AB + a 2)βn

6 ⎡ ⎤ RAB a1 ⎥ 6 fdD,63(CSO) = ⎢s6 + AB [1 + exp(RAB − 2.5R 0 )] ⎦ RAB + (2.52 )6 ⎣

(1)

Later, revised rational damping to finite values for small interatomic distances developed by Becke and Johnson (BJ) was adopted.23−25 It enters the calculations as11 fd , n (rAB) =

r ABn

where β6 is the additional adjustable parameter and β8 = β6 + 2. When β6 is set to 6, OP damping is equivalent to the BJ damping, and for higher values of β6, it gradually approaches a function similar to the zero damping. This additional flexibility brings slightly higher accuracy when compared to its predecessor, the BJ damping. On the contrary, the dispersion correction scheme introduced by Schröder et al. is motivated by reduction of the number of empirical parameters in hand with truncation of expansion at n = 6; therefore, it is referred to as C-six-only (CSO) damping. It employs sigmoidal interpolation function of the form 1 + 1/(R2AB + const) in an attempt to approximate higher order terms while keeping the damping to finite dispersion energy at zero distance (eq 6). It delivers results similar to the BJ damping21 while reducing the number of method-dependent parameters from three to one (a1). There were three additional parameters considered in the original work, but they were found to be almost independent of the DFT functional used so that they enter the final formula as constants:

sn 1 + 6(rAB/(sR , nR 0AB))−αn

(4)

β

fdD, n3(OP) (rAB) =

where the expansion is terminated at n = 8, A and B are indexes of the atoms in the system, Cn,AB are dispersion coefficients derived from time-dependent polarizabilities for each atom pair AB, and rAB is the respective interatomic distance. For the description of dispersion between atoms at shorter distances where part of the dispersion is covered by the correlation term in the DFT functional, the damping function fd,n(rAB) has to be used. Grimme had introduced two versions of the damping function. In the so-called zero damping approach, it is defined as a function of tabulated cutoff radii RAB 0 for an AB pair, the radius scaling parameter sR,n (where sR,6 is used as an adjustable parameter and sR,8 = 1), a global scaling factor s6 = 1, the steepness parameters α6 = 14 and α8 = 16, and is given by fd , n (rAB) =

C6,AB

Here, a1 and a2 are free fit parameters. Note that the number of adjustable parameters is two (s8 and sR,6) and three (s8, a1, a2) for zero and BJ damping, respectively. 2.2. Alternative Damping Functions. Recently, other authors have introduced modifications of the damping function in the D3 methodology. Head-Gordon et al. generalized the BJ damping by optimizing the exponent n in eq 3.20 The resulting damping function, named optimized power (OP) damping, takes the form of

2. COMPUTATIONAL METHODS 2.1. D3 Dispersion Correction. This work is based on the D3 dispersion correction developed by Grimme which is described in detail in the literature.10,12 We use the formulation of the correction and the atomic parameters from Grimme’s work, fitting only the few method-specific global parameters in the damping function. The dispersion takes the form of a damped pairwise interatomic potential Edisp = −∑

C8,AB

(3)

C9ABC = − C6,ABC6,ACC6,BC

where the cutoff radii are calculated as 3576

(9) DOI: 10.1021/acs.jctc.7b00365 J. Chem. Theory Comput. 2017, 13, 3575−3585

Article

Journal of Chemical Theory and Computation

proteins for parametrization of semiempirical methods from which we here use only the equilibrium geometries (IHB15).17 Finally, we pass from interaction energies to conformation energies of a tripeptide which are determined by an interplay of noncovalent interactions and constraints due to the covalent structure, especially torsional potentials. The set comes from our earlier work34 and was also included in the GMTKN database35,36 under the name PCONF. However, the original benchmark interaction energies were calculated in basis sets that are not up to the current standards. To make this data set compatible with the data sets of interaction energies, we have recalculated it with a setup identical to the one used for the S66x8 data set28 (i.e., a composite CCSD(T)/CBS scheme built from MP2/CBS energy extrapolated from aug-cc-pVTZ and aug-cc-pVQZ basis sets and from a CCSD(T) correction calculated in the aug-cc-pVDZ basis). The new benchmark conformation energies are reported in the Supporting Information in Table S1. In the validation, we work solely with root-mean-square error (RMSE) which is a more robust error measure than the mean absolute error. To compare the error between different data sets, we calculate a corresponding relative error as RMSE divided by the average magnitude of the interaction energy in the respective data set (the coefficient of variance of RMSE, for simplicity referred to as relative RMSE). In the PCONF data set, the conformation energies are made relative to the average conformer energy in the set, and the relative RMSE is calculated as RMSE divided by the range of these relative values. 2.5. Parametrization Protocol. The parametrization was carried out by the means of minimization of an error function using its gradient with respect to the parameters computed numerically. In the case of a training set consisting of a single, consistent data set, we choose the root-mean-square error (RMSE) as the error function. When the training set is built from multiple data sets with different properties, the error function is built as a sum of relative RMSE errors calculated separately in each data set. As there is a large number of combinations of basis sets, DFT functionals, and variants of the D3 correction, we divide the parametrization procedure into two steps. The first step is a preliminary parametrization where we consider multiple basis sets but only the two canonical forms of the D3 correction, i.e. using only zero and BJ damping. In this step, the S66x8 data set with excluded hydrogen bonds is used as the training set. In the second step, we narrow the focus down to the single best performing basis set identified earlier. In this step, we explore all the combinations of the tested DFT functionals, the four damping functions described above, and we repeat the fitting with and without the 3-body term. Here, we use larger training set consisting of the S66x8, X40, and L7 data sets. 2.6. Software Used. All of the DFT calculations presented in this article were performed using the Turbomole 7.0 package employing the resolution of identity approximation.37 HF-3c and PBEh-3c calculations were performed in ORCA package.38 The PM6-D3H4X17 calculations were performed using the MOPAC2016 program.39 The basis sets not available in the package were downloaded from the EMSL Basis Set Exchange Web site.40 The Cuby framework (http://cuby4.molecular.cz) was used to automate the calculations of the data sets.41 The parameters developed in this work are available as a patch applicable to the dftd3 software42 (version 3.2, BJ damping only) and are going to be included in the new version

The importance of the 3-body term in the dispersion correction remains a subject of active debate. As the D3 correction was parametrized without it, using the 3-body dispersion energy leads to slightly worse results in small molecules. The authors did not recommend using the 3-body dispersion term in general applications.10 It was also found that the contribution of dispersion is smaller than the intrinsic errors in the description of 3-body effects in DFT.6 On the other hand, it is obvious that it improves the description of dispersion in larger systems (where the number of 3-body contributions grows faster than the number of pairs of atoms). On the other hand, this contribution is only very short-ranged, and it could thus be effectively accounted for in the parametrization of the pairwise dispersion correction. In this work, we test following three variants: First, we neglect the 3-body dispersion completely. Second, we parametrize the correction without it, but we apply it in the validation. Third, it is already used in the parametrization. The last approach requires reliable benchmark data for large systems to be used in parametrization and validation which was made possible only recently when more accurate calculations became available for sufficiently large model systems. 2.4. Data Sets Used. The parametrization and subsequent testing of the dispersion correction relies on accurate benchmark interaction energies. Here, we use multiple established data sets of benchmark data that cover various aspects of the method. In most of them, the benchmark results were obtained using composite CCSD(T)/CBS calculations with estimated accuracy of 1−2%.2 The basis for the parametrization is the S66x8 data set, a set of 66 dissociation curves of noncovalent complexes featuring hydrogen bonds, dispersion, and mixed-type interactions of organic molecules.28 The use of dissociation curves is crucial for proper description of the whole potential energy surface and not only the equilibrium geometries. The S66 set of equilibrium geometries29 is then used in a subsequent test of the accuracy because it features more accurate benchmark energies. In the extended parametrization, we use also the X40 data set30 which was built as an extension of the S66 data set covering halogenated compounds. When used in the parametrization, it ensures better transferability of the resulting methods to different elements. When working with the 3-body term in the correction, it is necessary to include larger systems in the parametrization and in testing. For the former, we use the L7 data set which features the largest systems accessible to the CCSD(T) or equivalent QCISD(T) calculations.31 For validation, we use the S12L data set32 of interaction energies in supramolecular complexes. However, we have found that the original interaction energies are not accurate enough for detailed testing of the DFT-D methods. Instead, we use the latest recalculation of the data set by Hesselmann and Korona who used DFT-SAPT to obtain more reliable and consistent benchmark interaction energies.33 For the validation, we use several additional data sets. For general testing in smaller model systems, we use the set of 207 dimers of organic molecules in geometries taken from crystals (3B69 dimers).6 This is not only a large and diverse test set well suitable for testing of methods developed using the S66 data set but it also covers different binding motifs and nonequilibrium geometries. To investigate the effects of BSSE in ionic interactions, we also include a set of hydrogen bonds of charged residues. We have previously developed 15 dissociation curves of hydrogen bonds covering charged moieties found in 3577


Article


Figure 1. Root-mean-square error (RMSE) of the preliminary parametrization of DFT-D3 (for def2-QZVP, the original parameters are used). The results are plotted separately for each data set (X40, S66, and L7) and damping function (zero = a standard “zero-damping” formula and BJ = rational damping to finite values according to Becke and Johnson). Average of the individual errors is plotted in black.

sets of comparable size. It is clear that this basis set is ideally suited for the intended purpose. 3.2. Preliminary Parametrization of D3 Dispersion Correction. In the second step, we performed a basic parametrization of the D3 dispersion correction for the DZVP-DFT basis as well as for several other basis sets selected for comparison (def2-TZVP, def2-SVP, and 6-31G*). The parametrization of dispersion was performed on the subset of the S66x8 data set, containing only dispersion-dominated complexes (i.e., the geometries and interaction energies of dispersion-bound complexes at points along dissociation curves). This is a training set sufficient for parametrization of the dispersion correction for well-behaved methods, and it is close to the original parametrization of D3 for the def2-QZVP basis which was performed on the complete S66x8 data set. Note that we intentionally excluded systems containing hydrogen bonds from the parametrization where the BSSE is the largest. Such parametrization would lead to underestimated dispersion in all the other systems and collapsing two terms with diverse asymptotic decay, dispersion, and BSSE, into one correction. In this preliminary parametrization, we used only the BJ and zero damping. In some combinations of functionals and basis sets, the parameters of the BJ damping function reach unphysical values, which results in positive dispersion energy, balancing the overstabilization caused by the large BSSE. These combinations have been omitted from the discussion. The zero damping is more robust and generally more suitable for methods where the interaction energy may be overestimated even in the underlying method. The transferability and quality of parametrization have been tested on the data sets X40 and L7.30,31 The former data set plays a prominent role in the validation because it contains halogen elements not present in the training set. The L7 data set mostly comprises dispersion-bound complexes of larger size which makes it possible to assess the transferability of the correction to large systems. We have utilized the RMSE as the most important statistic criterion for ranking the methods. The errors for all the tested setups are summarized as a plot of the RMSE (Figure 1). The combinations of basis sets and functionals are ordered by their complexity, starting with the

of the Cuby framework which also implements the additional damping functions.41,43

3. RESULTS AND DISCUSSION 3.1. Search for Optimal Basis Set. In the first step, we searched a large library of basis sets in order to identify those that exhibit acceptably small basis set superposition error. The basis sets considered were STO-3G, MINI, MINI3, MINIS, MINI3S (minimal basis), 6-31G, 6-31G*, 6-31G**, SV, SV(P), SVP, def2-SV(P), def2-SVP, def2-SVPD, MIDI, MIDIX, DZ, DZP, DZVP, DZVP-DFT, cc-pVDZ, and aug-cc-pVDZ (splitvalence and double-ζ).44−51 These basis sets were combined with multiple common DFT functionals: BLYP,52 PBE53 (GGA functionals), B3LYP,54 PBE055 (hybrid functionals), and TPSS56 (meta-GGA). The list of DFT functionals is not exhaustive, but the accuracy of the calculations here is determined mostly by the quality of the basis set while the difference between the functionals is less important. Additionally, our choice of the functionals is based on previous experience with DFT-D3 in larger basis sets. To maximize the efficiency, it is preferred to use a pure GGA functional. In the detailed survey of DFT-D3 methods, Grimme reported that very small error (less than 0.3 kcal/mol in the S22 data set) can be achieved with the BLYP functional.22 Additionally, Sherill with co-workers studied the accuracy of D3 correction for large basis sets where they showed that BLYP and B3LYP are superior to other tested functionals when applied to noncovalent interactions.57 To these, we add several more functionals for comparison. Our results confirm their qualities and show that the relative performance of different functionals in the small basis is very similar to the calculations in large basis so that previous evidence supporting their choice can be used. In this initial screening, the magnitude of the BSSE was estimated from the error of interaction energies of hydrogenbonded systems from the S66 data set. All minimal basis sets yielded poor results therefore basis set of at least double-ζ size had to be used in our method which does not include specific BSSE correction. Among the remaining basis sets, the DZVPDFT basis exhibited results significantly better than other basis 3578


Article


well the gap between more approximate but less robust semiempirical methods or methods relying on large empirical corrections and the more accurate but also significantly more expensive DFT-D in large basis. Even though BLYP-D3/DZVP-DFT is somewhat less accurate than DFT-D in large basis, it performs consistently in different kinds of systems. The less expensive approaches work well for some of them but fail in others, which indicates that they are not robust enough for general applications without previous validation. Finally, we compare DFT-D3 in DZVP-DFT with another DFT-D approach using double-ζ basis, PBEh-3c.14 This method works reasonably well in small model complexes, but the error grows significantly when we pass to larger systems. The basis set used here, modified def2SVP basis (def2-mSVP), has very large BSSE. That is corrected by a separate empirical correction of the same form as used in HF-3c. These empirical corrections for BSSE and other errors due to insufficient size of the basis set can work well is some cases, but their transferability is limited because they do not depend on the actual electronic structure of the interacting molecules. It is thus better to avoid this issue by using a basis set that is less susceptible to BSSE rather than relying on an additional correction. We demonstrate this on the dissociation curve of the water dimer which features a hydrogen bond, an interaction motif prone to large BSSE. Here, the BSSE of several combinations of method and basis set was evaluated more rigorously using the counterpoise scheme of Boys and Bernardi.59 The results are plotted in Figure 3. In the DFT calculations combining the BLYP functional with different basis sets (left panel of the plot), the BSSE generally behaves as expected, growing with the decreasing size of the basis set. However, there are large differences among the double-ζ basis sets. The BSSE of the DZVP-DFT basis is only negligibly larger than of the def2-TZVP basis. On the other hand, the def2-SVP basis yields BSSE comparable to the minimal basis MINI. The right panel of Figure 3 illustrates the behavior of the HF3c and PBEh-3c methods. Here, we plot both the actual BSSE in the QM calculation (solid line) and its value after adding the corrections for BSSE and basis set size that are part of the 3c corrections. The actual BSSE is rather large, although smaller than in the BLYP calculations in the same basis. The empirical corrections work well in the case of HF-3c but overcorrect in PBEh-3c. The resulting error in the latter case remains very large. Moreover, the distance dependence of PBEh-3c becomes unphysical which is the most probable reason for the failure of the method in larger systems discussed above. On the basis of all these results, it is clear that the DZVPDFT combined with BLYP yields exceptionally small BSSE for a basis of this size. That translates into very good description of noncovalent interactions. We thus prefer to use this basis set and couple it with unmodified dispersion correction, rather than depending on additional corrections for basis set size effects whose robustness is questionable. 3.4. Final Parametrization for DZVP-DFT Basis Set. The small training set used in the preliminary parametrization is not sufficient for a robust parametrization of damping functions with larger number of parameters. This is already indicated by the problems with parametrization of the BJ damping described above, and we found it even more difficult in the case of the OP damping. We thus performed the final round of parametrization with larger and more diverse training set. Now, we focus only

large def2-QZVP coupled with the original parametrization of the correction. The def2-TZVP basis stands for a triple-ζ basis sets from the same series. Finally, the DZVP-DFT basis is compared to two other basis sets of similar size that are commonly used, 6-31G* and def2-SVP. Generally, the BJ damping has proved to be superior over zero damping, although this behavior is strongly system-, basis set-, and functional-dependent. There is no simple increase in accuracy with the complexity of the functional (decreasing from left to right in the plot). The BLYP and B3LYP functionals yielded on average more accurate results than PBE, TPSS, and PBE0, whereas the BLYP has lower computational demands when compared to the hybrid B3LYP functional. The final results support the special qualities of the DZVP-DFT basis set coupled with the BLYP functional having the average RMSE 0.142 kcal/mol which is comparable to the triple-ζ basis set. This exceptional behavior is most likely caused by the more diffuse character of the basis set (although it outperforms even the DZ basis set with additional diffuse functions) combined with parametrization specifically focused on its use in DFT calculations. 3.3. Analysis of Basis Set Superposition Error in DZVP-DFT Basis Set. At this stage, we compare the best resulting method, BLYP-D3/DZVP-DFT, with both more efficient methods (PM6-D3H4X17,58 and HF-3c15) as well as to more accurate DFT-D (BLYP-D3/def2-QZVP). The PM6D3H4X is a semiempirical method that yields the best description of noncovalent interactions in its class because of the use of advanced corrections for dispersion and hydrogen and halogen bonding. The HF-3c method does not rely on the semiempirical approximations, but the use of the minimal basis set (a modified MINI basis) has to be compensated by empirical corrections (in addition to a D3 dispersion correction). For this validation, we use the X40 and L7 data sets that test the transferability of the methods to other elements and to larger systems. Note that these two data sets are included in the final parametrization. The errors in these data sets are plotted in Figure 2. It is obvious that the DFT-D3 in small basis fills

Figure 2. Comparison of the preliminary parametrization of the BLYPD3/DZVP-DFT method with both more approximate approaches and with DFT-D in large basis. The root-mean-square error (RMSE) of the selected methods in the S66, X40, and L7 data sets is plotted in green, red, and blue, respectively. Since the PM6-D3H4X method does not contain any specific correction for strong hydrogen bonds in halogenated compounds, we omit these complexes in the X40 data set. 3579


Article


Figure 3. Basis Set Superposition Error (BSSE) for geometries along dissociation curve of water dimer calculated for DFT in different basis sets (left) and for the HF-3c and PBEh-3c methods (right, without and with the empirical corrections). Vertical line marks the equilibrium distance.

on the DZVP-DFT basis, which was identified as the best candidate, and we explore more variants of the D3 correction. The training was constructed from three data sets. First, we use the complete S66x8 data set which covers most common interactions at different distances. Second, we add the L7 data set which ensures that the parametrization is transferable to larger systems. Third, we also include the X40 set containing halogens in order to avoid a parametrization specific only to organic elements. The error measure minimized in the parametrization is constructed as a sum of the relative errors in these data sets to which all the three data sets contribute with equal weight. More robust training set allows us not only to parametrize the correction with the BJ damping but also the OP damping which has one more parameter. To make the study complete, we have included also the CSO version of the D3 correction which uses only one global parameter. The inclusion of larger noncovalent complexes (the L7 data set) in the training set also allows us performing the parametrization of the D3 correction with the 3-body term. This approach should produce a consistent set of parameters that balances well the description of both small and large systems. This is an important difference compared to the original Grimme’s parametrization in which the parametrization was performed without the 3-body term. That resulted into a small decrease in accuracy in small systems when the 3-body term is used. Combining the five functionals, four damping functions, and the optional use of the 3-body term resulted in 40 parametrization runs. The resulting parameters for the BLYP functional are listed in Table 1; parameters for all the functionals are provided in the Supporting Information in Table S2. All the parametrizations successfully converged, and there was only one case where the parameters became negative. Surprisingly, this did not occur in the BJ damping that was problematic in the preliminary parametrization but in the zero damping. Here, the s8 parameter became slightly negative in the case of PBE0 functional when the 3-body term was not used.

Table 1. Final Parameters for BLYP-D3/DZVP-DFT with Different Variants of Damping Functionsa Damping BJ BJ/3b CSO CSO/3b OP OP/3b zero zero/3b a

sr,6

1.0838 1.0937

s8

a1

a2

3.3122 3.9502

0.4866 0.5083 0.8977 0.9873 0.5137 0.5162

4.686 4.7204

0.9237 1.1769 1.1805 1.3347

3.3564 3.4447

β

19.6287 18.1201

Values are dimensionless or in au.

The quality of the fit (which is indicated by the optimized value of the error) provides a valuable insight into the data. First, we average the quality of the fit over all the variants of the D3 correction for each functional. The best fit was obtained for BLYP closely followed by B3LYP. Following that, the fit was slightly worse for TPSS and considerably worse for PBE0 and PBE. The values of the overall error (which is a relative number) are, in this order, 0.52, 0.55, 0.60, 0.68, and 0.79. The same analysis can be performed for the different damping functions, averaging over all the functionals. Here, the best fit was obtained with the OP damping which is most flexible, closely followed by the BJ damping. As follows is the CSO damping, which is somewhat surprising because it uses a single parameter only. The worst fit was obtained with the zero damping. It is interesting that for each damping function slightly better fit was obtained without the 3-body term. Although this difference was very small, it suggests that the 3body term is rather problematic (possibly because we followed the Grimme’s approach with no adjustable parameters in it). Likewise, it is interesting to compare new parameters for the DZVP-DFT basis set with the original parametrization of D3(BJ) for the large basis set. We do this for the BLYP functional. For the sake of comparison, we reoptimized the parameters for the def2-QZVP basis set using the same protocols we applied to the smaller basis sets, and we also list the parameters for several other basis sets obtained in the 3580


Article

Journal of Chemical Theory and Computation Table 2. Comparison of Parameters in Becke−Johnson Damping Function for BLYP Functional and Different Basis Sets As Obtained by Parametrization on Different Training Setsa Basis set

Parametrization

a1

a2

s8

def2-QZVP

Grimme Preliminary Final Final, with 3b Preliminary Preliminary Final Final, with 3b Preliminary Preliminary

0.430 0.470 0.467 0.489 0.539 0.385 0.487 0.508 0.421 0.145

4.230 3.740 3.976 4.033 3.335 4.133 4.686 4.720 3.833 4.860

2.700 2.097 2.392 2.840 1.865 1.773 3.312 3.950 0.914 0.656

def2-TZVP DZVP-DFT

def2-SVP 6-31G* a

Table 3. Error of 15 Best Performing Variants of DFT-D3/ DZVP-DFT Method Obtained as a Relative Error Averaged over Four Validation Data Sets

Values are dimensionless or in au.

preliminary parametrization. The parameters are listed in Table 2. First, comparison of the original parametrization for def2QZVP and the results of our parametrizations for the same basis sets shows that they are all reasonably similar. The results of our final protocol are close to the original parametrization. In a large basis set where BSSE is negligible, the method seems to be very robust, and the parameters do not depend much on the training set used. Second, we pass to smaller double-ζ basis sets ordered in the table by increasing error. Here, the BSSE is being compensated mainly by the reduction of the s8 scaling factor which becomes substantially smaller in the case of def2SVP and 6-31G* basis sets where the BSSE becomes very large. Third, when the 3-body term is used in the parametrization, the value s8 increases in order to compensate the 3-body effects which are on average repulsive. Finally, the parameters for the DZVP-DFT basis are significantly more similar to the ones for the larger basis sets than to those for other double-ζ basis sets. This confirms our hypothesis based on the previous analysis of the magnitude of the BSSE: The DZVP-DFT basis set works so well because of the exceptionally small BSSE, rather than because of the compensation of the BSSE by unphysical parametrization of the dispersion correction. 3.5. Validation of DFT-D3/DZVP-DFT Energies. The parametrization of the D3 correction is validated on four data sets not used in its development. The accuracy of interaction energies in smaller systems is tested in the 3B69 dimers set which provides a diverse test set covering also the nonequilibrium geometries. The S12L data set is used to test the applicability of the method to larger molecular systems. As a test sensitive specifically to BSSE, we use the IHB15 data set which features ionic H-bonds, strong polar interactions where BSSE could be very large. Finally, to extend the validation beyond interaction energies, we add the PCONF data set of conformation energies of a FGG tripeptide. As a measure of the overall accuracy of all the variants of the DFT-D3/DZVP-DFT method, we use the average of the relative RMSE in these data sets. To discern the general trends, all the variants of the method can be ordered according to the average error in the tests sets, analyzing the factors that affect this ranking. It is obvious that among the DFT functionals BLYP works best. Among the top 15 variants out of 60 (top 25%, listed in Table 3), BLYP is present seven times (and also occupies the first two places), followed by TPSS (four times), B3LYP (three times), and

a

Functional

Damping

3-bodya

Average error (%)

BLYP BLYP B3LYP TPSS BLYP B3LYP TPSS BLYP TPSS BLYP BLYP PBE0 TPSS BLYP B3LYP

BJ CSO BJ BJ OP CSO CSO OP OP OP BJ OP OP CSO OP

−/Yes −/Yes −/Yes −/Yes −/Yes −/Yes −/Yes Yes/Yes −/Yes −/− Yes/Yes −/Yes Yes/Yes Yes/Yes −/Yes

11.6 11.6 11.8 11.9 12.1 12.2 12.5 12.7 12.8 13.2 13.2 13.2 13.3 13.3 13.4

Used in parametrization/in validation.

PBE0 (once). As a pure GGA functional, BLYP is also the fastest in the group, so it is an obvious choice for practical use. Among the damping functions, the OP damping works best (present 7 times in the top 15 variants), followed by the BJ and CSO damping which are present four times. Zero damping scored considerably worse (in the second half of the rankings). The success of the OP damping can be easily explained by its greater flexibility. More surprising is the relatively good performance of the CSO approach which uses only one adjustable parameter. It is apparent here that the other parameters which were set to the fixed values in the original work are truly method independent. Even more interesting is the role of the 3-body term. The best results were obtained when it was not considered in the fitting but used in the testing (10 out of 15 cases, including the seven best variants). Above, we have shown that the use of the 3-body term makes the fit worse, yet it is important in the large systems included in the test set. The best solution is thus to follow the common practice in DFT-D3; it is to use only the pairwise parametrization and add the 3-body term in large systems without changing the parameters. This approach is followed by the parameter sets obtained with the 3-body term. Among the top 15 methods, there is only one which does not use the 3-body term at all. It uses the OP damping where the flexibility of the damping function allows a fit that better balances the description of small and large systems. In summary, we recommend using the BLYP-D3/DZVPDFT setup with BJ damping parametrized without the 3-body term and applying this term when the size of the studied system justifies it. This set of parameters was added to the dftd3 software42 where it can be readily used; we provide a patch for the source code in the Supporting Information. In the next paragraphs, we look more closely at the individual results obtained by the BLYP functional (listed in Table 4). The table lists not only the final BLYP-D3/DZVP-DFT results but also the results obtained with the preliminary parametrization and with the parameters for the def2-QZVP basis. We also include the results calculated in the def2-QZVP basis for comparison. The complete tables of results for the remaining functionals are available in the Supporting Information in Tables S3−S6. 3581


Article

Journal of Chemical Theory and Computation Table 4. Errors (RMSE, kcal/mol) of Variants of BLYP-D3 Calculations in def2-QZVP and DZVP-DFT Basis Setsa 3-bodyb

S66x8

S66

BJ BJ

−/− −/Yes

0.22 0.22

0.24 0.22

BJ BJ

Yes/Yes −/−

0.23 0.23

0.21 0.21

BJ BJ

−/− −/Yes

0.80 0.78

0.93 0.89

BJ BJ

−/− −/Yes

0.68 0.68

0.73 0.72

−/− Yes/Yes −/Yes −/− Yes/Yes −/Yes −/− Yes/Yes −/Yes −/− Yes/Yes −/Yes

0.71 0.70 0.75 0.66 0.64 0.69 0.64 0.64 0.66 0.74 0.76 0.75

0.60 0.58 0.65 0.60 0.59 0.64 0.53 0.53 0.54 0.74 0.79 0.73

Damping

BJ BJ BJ CSO CSO CSO OP OP OP zero zero zero

X40

L7

S12L

PCONF

IHB15

Averagec

0.71 0.71

0.74 0.69

0.63 0.63

15.7% 13.1%

0.71 0.72

0.69 0.72

0.63 0.63

13.0% 13.7%

1.03 1.02

0.59 0.51

2.51 2.51

24.6% 21.7%

0.95 0.95

0.46 0.40

2.45 2.44

17.5% 15.0%

0.86 0.86 0.86 0.86 0.87 0.87 0.82 0.82 0.82 0.90 0.93 0.90

0.36 0.35 0.35 0.35 0.33 0.32 0.36 0.34 0.31 0.39 0.39 0.35

1.95 1.98 1.95 2.03 2.10 2.03 1.87 1.88 1.87 2.30 2.38 2.30

13.5% 13.2% 11.6% 13.6% 13.3% 11.6% 13.2% 12.7% 12.1% 15.2% 15.5% 15.1%

3B69 dimers

def2-QZVP, original parameters 0.39 1.63 4.76 0.38 0.90 2.52 def2-QZVP, reparametrization using our protocol 0.32 0.97 2.36 0.33 0.96 2.79 DZVP-DFT, parameters for def2-QZVP 1.52 4.35 10.66 1.51 3.13 8.35 DZVP-DFT, preliminary parametrization 1.39 1.38 5.05 1.38 1.23 3.05 DZVP-DFT, final version 1.01 0.88 3.00 1.01 0.95 2.74 1.01 1.50 1.13 1.13 0.84 3.04 1.19 0.90 2.79 1.13 1.56 1.20 0.91 1.01 2.98 0.92 0.99 2.69 0.90 1.64 2.29 1.07 1.96 3.71 1.13 2.18 3.77 1.07 2.73 3.94

a Errors are listed for the training and validation data sets. The overall average relative error is calculated from the validation data sets only. bUsed in parametrization/in validation. cAverage of relative errors in the four validation sets.

in the S12L data set are surprisingly good. This is an important result as this methodology is likely to be used for large systems. The only case where the DZVP-DFT data set works considerably worse than the def2-QZVP basis is the data set of ionic H-bonds. In this case, the BSSE is considerably larger than in the neutral systems and becomes the limiting factor. However, the interaction energies in this data set are correspondingly larger when compared to neutral systems so that the relative error is still acceptably small (around 10%). 3.6. Testing DFT-D3/DZVP-DFT in Geometry Optimizations. In the second validation step, we analyze the performance of this method in terms of geometry optimization. The use of a smaller basis set size for geometry optimization prior to single-point calculations in a larger basis set is a common practice in many studies. It is generally assumed that the geometries are less sensitive to the basis set size. Since we are only interested in the effect of the basis set size in DFT calculations, we take geometries optimized with BLYP-D3/ def2-QZVP as a benchmark (in the S66 data set, rescaled MP2 geometries are used). They are used for the validation of smaller basis sets both combined with the BLYP functional and the reparametrized D3 dispersion. Table 5 summarizes a statistical analysis of geometry optimizations in the S66 data set performed with def2-TZVP, def2-SVP, and DZVP-DFT basis sets with the newly developed parameters (for the former two basis sets, the parameters from the preliminary parametrization are used). It is evident that the results strongly depend on the quality of the basis set. For example, when passing within the same basis set family from a triple-ζ (def2-TZVP) to double-ζ (def2-SVP) size basis set, the RMSD increases seven times. On the contrary, the DZVP-DFT basis set yielded satisfactory

We start the discussion with the calculations in the def2QZVP basis. When the original Grimme’s parametrization is used, the overall error in the validation sets amounts to 15.7% and drops to 13.1% when the 3-body term is used. It improves the results in all the data sets, most notably in the large systems in the L7 and S12L sets. The parameters for this basis set obtained by the protocol described in this paper improve the results only slightly which indicates that the original parametrization is sufficient despite being based on smaller training set. Interestingly, it is possible to obtain very similar results both with and without the 3-body term when the validation is consistent with the parametrization. This suggests that the 3body term can be effectively accounted for at the pairwise level as long as large systems are included in the training set. When the original parameters for the def2-QZVP are used for the DZVP-DFT basis, the results are not very satisfactory (with the average errors above 20%). Compared to the final parametrization, the worst cases are the large systems in the S12L set. Although the preliminary parametrization improves the results substantially, the errors in the large systems remain relatively high. The final parametrization which included the large systems (the L7 data set) in the parametrization improves the description of the large systems consistently. The differences among the damping functions are very small, with the exception of the zero damping which performs considerably worse. The overall error in the test sets is similar to the calculations in the def2-QZVP basis, but this is caused mainly by the surprising improvement of the peptide conformation energies which is likely caused by an error compensation. In the S12L and 3B69 dimers, the results are similar to the ones in the larger basis but slightly worse. On the other hand, in absolute terms, the errors 3582


Article


considered. However, the BJ damping yields comparably good results and is readily available in many software packages.

Table 5. Average Root-Mean-Square Deviation (average RMSD) and Statistical Analysis of Closest Intermolecular Distance in Complexes (both in Å) def2-TZVP

DZVP-DFT

Overall difference to benchmark geometries 0.010 0.023 Errors in the closest intermolecular distance Mean signed error −0.011 0.012 Mean unsigned error 0.011 0.028

■

def2-SVP

RMSD

The Supporting Information is available free of charge on the ACS Publications website at DOI: 10.1021/acs.jctc.7b00365. Table of updated benchmark energies for the PCONF data set, table of parameters for all final variants of the D3 correction described in this paper, tables of the results obtained using these parameters and a patch with the new parameters for the dftd3 program. (PDF)

0.070 −0.044 0.054

results, delivering three times smaller RMSD than another basis set of similar size (def2-SVP). Besides the overall root-mean-square deviation, it is also interesting to look at the closest intermolecular distances that are most affected by the BSSE. As expected, these closest contacts are underestimated in calculations where BSSE is large, such as in the def2-SVP basis set. The DZVP-DFT basis performs significantly better with systematic error (calculated as mean unsigned error) of only 0.012 Å. 3.7. Computational Times. Finally, we discuss the computational demands of the calculations. Table 6 shows

■

CB[5] = 90 atoms

CB[8] = 144 atoms

6-31G* def2-SVP DZVP-DFT def2-TZVP def2-QZVP

5.8 6.4 7.7 25.0 108.5

6.8 8.0 9.3 31.1 161.2

AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]. ORCID

Jiří Hostaš: 0000-0001-6652-602X Jan Ř ezác:̌ 0000-0001-6849-7314 Notes

The authors declare no competing financial interest.

■

Table 6. Computational Times for Cucurbit[n]uril System, where n = 5, 8 with 90 and 144 Atoms, Respectivelya Basis set/system size

ASSOCIATED CONTENT

S Supporting Information *

ACKNOWLEDGMENTS This work has been supported by the Czech Science Foundation by Grant No. P208/16-11321Y and is a part of the Research Project RVO 61388963 of the Institute of Organic Chemistry and Biochemistry, Czech Academy of Sciences.

■

REFERENCES

(1) Burke, K. Perspective on density functional theory. J. Chem. Phys. 2012, 136, 150901. (2) Rezac, J.; Hobza, P. Benchmark Calculations of Interaction Energies in Noncovalent Complexes and Their Applications. Chem. Rev. 2016, 116, 5038−5071. (3) Misquitta, A. J.; Jeziorski, B.; Szalewicz, K. Dispersion energy from density-functional theory description of monomers. Phys. Rev. Lett. 2003, 91, 4. (4) Brunk, E.; Rothlisberger, U. Mixed Quantum Mechanical/ Molecular Mechanical Molecular Dynamics Simulations of Biological Systems in Ground and Electronically Excited States. Chem. Rev. 2015, 115, 6217−6263. (5) Polly, R.; Werner, H. J.; Manby, F. R.; Knowles, P. J. Fast Hartree-Fock theory using local density fitting approximations. Mol. Phys. 2004, 102, 2311−2321. (6) Rezac, J.; Huang, Y.; Hobza, P.; Beran, G. J. O. Benchmark Calculations of Three-Body Intermolecular Interactions and the Performance of Low-Cost Electronic Structure Methods. J. Chem. Theory Comput. 2015, 11, 3065−3079. (7) Grimme, S. Density functional theory with London dispersion corrections. Wiley Interdiscip. Rev.: Comput. Mol. Sci. 2011, 1, 211−228. (8) 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, 23. (9) 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, 52. (10) Grimme, S.; Antony, J.; Ehrlich, S.; Krieg, H. A consistent and accurate ab initio parametrization of density functional dispersion correction (DFT-D) for the 94 elements H-Pu. J. Chem. Phys. 2010, 132, 19. (11) Grimme, S.; Ehrlich, S.; Goerigk, L. Effect of the Damping Function in Dispersion Corrected Density Functional Theory. J. Comput. Chem. 2011, 32, 1456−1465.

a

BLYP functional. Intel Xeon E5630 2.53 GHz, 8 cores, 5.8 GB RAM per core. Values in minutes.

the computational times obtained for systems increasing in size and for a given basis set. The steep step between double- and triple-ζ basis sets clearly rationalizes our focus on DZ-size basis sets. The calculations in DZVP-DFT were the slowest among the double-ζ basis sets; however, the differences were practically negligible.

4. CONCLUSIONS To conclude, it has been shown that the DFT method combined with a double-ζ size basis set can yield surprisingly accurate results if the basis is chosen well. Among the tested combinations of functionals and basis sets, the DZVP-DFT basis and BLYP functional are recommended for both geometry optimizations and energy calculations in large systems. The accuracy of this setup is as good as DFT-D3 calculations in triple-ζ basis sets, but the computational cost is significantly lower. The method was tested specifically for transferability to large systems, and the results obtained in the S12L data set are very promising. The superior performance of the DZVP-DFT basis can be attributed mainly to the fact that this basis set exhibits much smaller BSSE than other basis sets of comparable size. The resulting DFT-D3 method is more reliable than using a basis set that exhibits larger BSSE and adding another empirical correction for that such as in the PBEh-3c method which uses modified def2-SVP basis. Among the variants of the D3 correction, the OP damping worked best because it is the most flexible damping function 3583


Article


(35) Goerigk, L.; Grimme, S. Efficient and Accurate Double-HybridMeta-GGA Density Functionals-Evaluation with the Extended GMTKN30 Database for General Main Group Thermochemistry, Kinetics, and Noncovalent Interactions. J. Chem. Theory Comput. 2011, 7, 291−309. (36) Goerigk, L.; Grimme, S. A General Database for Main Group Thermochemistry, Kinetics, and Noncovalent Interactions - Assessment of Common and Reparameterized (meta-)GGA Density Functionals. J. Chem. Theory Comput. 2010, 6, 107−126. (37) TURBOMOLE V7.0 2015, a development of University of Karlsruhe and Forschungszentrum Karlsruhe GmbH, 1989−2007, TURBOMOLE GmbH, since 2007. http://www.turbomole.com (accessed July 2017). (38) Neese, F. The ORCA program system. Wiley Interdiscip. Rev.: Comput. Mol. Sci. 2012, 2, 73−78. (39) Stewart, J. J. P. MOPAC2012, 2012. Stewart Computational Chemistry, Colorado Springs, CO. http://OpenMOPAC.net (accessed July 2017). (40) Basis Set Exchange. https://bse.pnl.gov/bse/portal (accessed February 7, 2017). (41) Rezac, J. Cuby: An Integrative Framework for Computational Chemistry. J. Comput. Chem. 2016, 37, 1230−1237. (42) DFT-D3 Website. https://www.chemie.uni-bonn.de/pctc/ mulliken-center/software/dft-d3/ (accessed July 2017). (43) Rezac, J. Cuby 4, software framework for computational chemistry. http://cuby4.molecular.cz/ (accessed July 2017). (44) Weigend, F.; Ahlrichs, R. Balanced basis sets of split valence, triple zeta valence and quadruple zeta valence quality for H to Rn: Design and assessment of accuracy. Phys. Chem. Chem. Phys. 2005, 7, 3297−3305. (45) Godbout, N.; Salahub, D. R.; Andzelm, J.; Wimmer, E. Optimization of Gaussian-type basis sets for local spin density functional calculations. Part I. Boron through neon, optimization technique and validation. Can. J. Chem. 1992, 70, 560−571. (46) Frisch, M. J.; Pople, J. A.; Binkley, J. S. Self-consistent molecular orbital methods 25. Supplementary functions for Gaussian basis sets. J. Chem. Phys. 1984, 80, 3265−3269. (47) Andzelm, J.; Huzinaga, S.; Klobukowski, M.; Radzio, E. Model potential study of the interactions in Ar2, Kr2 and Xe2 dimers. Mol. Phys. 1984, 52, 1495−1513. (48) Hehre, W. J.; Stewart, R. F.; Pople, J. A. Self-Consistent Molecular-Orbital Methods. I. Use of Gaussian Expansions of SlaterType Atomic Orbitals. J. Chem. Phys. 1969, 51, 2657−2664. (49) Hehre, W. J.; Ditchfield, R.; Pople, J. A. Sef-Consistent Molecular Orbital Methods. XII. Further Extensions of Gaussian-Type Basis Sets for Use in Molecular Orbital Studies of Organic Molecules. J. Chem. Phys. 1972, 56, 2257. (50) Ahlrichs, R.; Bar, M.; Haser, M.; Horn, H.; Kolmel, C. Electronic structure calculations on workstation computers: The program system turbomole. Chem. Phys. Lett. 1989, 162, 165−169. (51) Hariharan, Pc; Pople, J. A. The influence of polarization functions on molecular orbital hydrogenation energies. Theor. Chim. Acta 1973, 28, 213−222. (52) Lee, C. T.; Yang, W. T.; Parr, R. G. Development of the ColleSalvetti correlation-energy formula into a functional of the electron density. Phys. Rev. B: Condens. Matter Mater. Phys. 1988, 37, 785−789. (53) Perdew, J.; Burke, K.; Ernzerhof, M. Generalized gradient approximation made simple. Phys. Rev. Lett. 1996, 77, 3865−3868. (54) Becke, A. D. Densityfunctional thermochemistry. III. The role of exact exchange. J. Chem. Phys. 1993, 98, 5648−5652. (55) Adamo, C.; Barone, V. Toward reliable density functional methods without adjustable parameters: The PBE0 model. J. Chem. Phys. 1999, 110, 6158−6170. (56) Tao, J. M.; Perdew, J. P.; Staroverov, V. N.; Scuseria, G. E. Climbing the density functional ladder: Nonempirical meta-generalized gradient approximation designed for molecules and solids. Phys. Rev. Lett. 2003, 91, 4.

(12) Grimme, S.; Hansen, A.; Brandenburg, J. G.; Bannwarth, C. Dispersion-Corrected Mean-Field Electronic Structure Methods. Chem. Rev. 2016, 116, 5105−5154. (13) DiLabio, G. A.; Koleini, M.; Torres, E. Extension of the B3LYPdispersion-correcting potential approach to the accurate treatment of both inter- and intra-molecular interactions. Theor. Chem. Acc. 2013, 132, 13. (14) 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, 19. (15) Sure, R.; Grimme, S. Corrected small basis set Hartree-Fock method for large systems. J. Comput. Chem. 2013, 34, 1672−1685. (16) Korth, M.; Pitonak, M.; Rezac, J.; Hobza, P. A Transferable HBonding Correction for Semiempirical Quantum-Chemical Methods. J. Chem. Theory Comput. 2010, 6, 344−352. (17) Rezac, J.; Hobza, P. Advanced Corrections of Hydrogen Bonding and Dispersion for Semiempirical Quantum Mechanical Methods. J. Chem. Theory Comput. 2012, 8, 141−151. (18) Miriyala, V. M.; Rezac, J. Description of Non-Covalent Interactions in SCC-DFTB Methods. J. Comput. Chem. 2017, 38, 688−697. (19) Christensen, A. S.; Kubar, T.; Cui, Q.; Elstner, M. Semiempirical Quantum Mechanical Methods for Noncovalent Interactions for Chemical and Biochemical Applications. Chem. Rev. 2016, 116, 5301− 5337. (20) Witte, J.; Mardirossian, N.; Neaton, J. B.; Head-Gordon, M. Assessing DFT-D3 Damping Functions Across Widely Used Density Functionals: Can We do Better? J. Chem. Theory Comput. 2017, 13, 2043−2052. (21) Schroeder, H.; Creon, A.; Schwabe, T. Reformulation of the D3(Becke-Johnson) Dispersion Correction without Resorting to Higher than C-6 Dispersion Coefficients. J. Chem. Theory Comput. 2015, 11, 3163−3170. (22) Goerigk, L.; Kruse, H.; Grimme, S. Benchmarking Density Functional Methods against the S66 and S66x8 Datasets for NonCovalent Interactions. ChemPhysChem 2011, 12, 3421−3433. (23) Becke, A. D.; Johnson, E. R. A density-functional model of the dispersion interaction. J. Chem. Phys. 2005, 123, 154101. (24) Johnson, E. R.; Becke, A. D. A post-Hartree-Fock model of intermolecular interactions. J. Chem. Phys. 2005, 123, 024101. (25) Becke, A. D.; Johnson, E. R. Exchange-hole dipole moment and the dispersion interaction. J. Chem. Phys. 2005, 122, 154104. (26) Axilrod, P. M.; Teller, E. Interaction of the van der Waals type between three atoms. J. Chem. Phys. 1943, 11, 299−300. (27) Muto, Y. Force between nonpolar molecules. Proc. Phys. Math. Soc. Japan 1943, 17, 629−631. (28) Rezac, J.; Riley, K. E.; Hobza, P. S66: A Well-balanced Database of Benchmark Interaction Energies Relevant to Biomolecular Structures. J. Chem. Theory Comput. 2011, 7, 2427−2438. (29) Rezac, J.; Riley, K. E.; Hobza, P. Extensions of the S66 Data Set: More Accurate Interaction Energies and Angular-Displaced Nonequilibrium Geometries. J. Chem. Theory Comput. 2011, 7, 3466−3470. (30) Rezac, J.; Riley, K. E.; Hobza, P. Benchmark Calculations of Noncovalent Interactions of Halogenated Molecules. J. Chem. Theory Comput. 2012, 8, 4285−4292. (31) Sedlak, R.; Janowski, T.; Pitonak, M.; Rezac, J.; Pulay, P.; Hobza, P. Accuracy of Quantum Chemical Methods for Large Noncovalent Complexes. J. Chem. Theory Comput. 2013, 9, 3364−3374. (32) Grimme, S. Supramolecular Binding Thermodynamics by Dispersion-Corrected Density Functional Theory. Chem. - Eur. J. 2012, 18, 9955−9964. (33) Hesselmann, A.; Korona, T. Intermolecular symmetry-adapted perturbation theory study of large organic complexes. J. Chem. Phys. 2014, 141, 094107. (34) Valdes, H.; Pluhackova, K.; Pitonak, M.; Rezac, J.; Hobza, P. Benchmark database on isolated small peptides containing an aromatic side chain: comparison between wave function and density functional theory methods and empirical force field. Phys. Chem. Chem. Phys. 2008, 10, 2747−2757. 3584


Article

Journal of Chemical Theory and Computation (57) 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. (58) Rezac, J.; Hobza, P. A halogen-bonding correction for the semiempirical PM6 method. Chem. Phys. Lett. 2011, 506, 286−289. (59) Boys, S. F.; Bernardi, F. The calculation of small molecular interactions by the differences of separate total energies. Some procedures with reduced errors. Mol. Phys. 1970, 19, 553−566.

3585


Accurate DFT-D3 Calculations in a Small Basis Set - Journal of

Recommend Documents