Minimally augmented Karlsruhe basis sets

Author's personal copy Theor Chem Acc (2011) 128:295–305 DOI 10.1007/s00214-010-0846-z REGULAR ARTICLE Minimally augmented Karlsruhe basis sets Jingjing Zheng • Xuefei Xu • Donald G. Truhlar Received: 1 October 2010 / Accepted: 26 October 2010 / Published online: 1 December 2010 ! Springer-Verlag 2010 Abstract We propose an extension of the basis sets proposed by Ahlrichs and coworkers at Karlsruhe (these basis sets are designated as the second-generation default or ‘‘def2’’ basis sets in the Turbomole program). The Karlsruhe basis sets are very appealing because they constitute balanced and economical basis sets of graded quality from partially polarized double zeta to heavily polarized quadruple zeta for all elements up to radon (Z = 86). The extension consists of adding a minimal set of diffuse functions to a subset of the elements. This yields basis sets labeled minimally augmented or with ‘‘ma’’ as a prefix. We find that diffuse functions are not quite as important for the def2 basis sets as they are for Pople basis sets, but they are still necessary for good results on barrier heights and electron affinities. We provide assessments and validations of this extension for a variety of data sets and representative cases. We recommend the new ma-TZVP basis set for general-purpose applications of density functional theory. Keywords Electronic structure ! Basis sets ! Density functional theory ! Bond dissociation energies ! Barrier heights ! Electron affinities ! Ionization potentials ! Noncovalent interactions ! Diffuse functions ! Minimally augmented basis set ! Double zeta ! Triple zeta ! Quadruple zeta ! ma-TZVP ! DBH24/08 database ! S22A database Electronic supplementary material The online version of this article (doi:10.1007/s00214-010-0846-z) contains supplementary material, which is available to authorized users. J. Zheng ! X. Xu ! D. G. Truhlar (&) Department of Chemistry and Supercomputing Institute, University of Minnesota, Minneapolis, MN 55455-0431, USA e-mail: truhlar@umn.edu 1 Introduction A large number of one-electron basis sets have been proposed for use in electronic structure calculations, and although some practitioners gravitate to the most familiar basis sets for all their work, many other quantum chemists, including our group, routinely struggle with the question of which basis to use whenever they start a project or decide to make their calculations more accurate (by selecting a larger or more appropriate basis set) or more affordable for larger systems (by selecting a smaller basis set). Popular collections of basis sets are those developed by the groups of Pople [1], Dunning [2, 3], and Ahlrichs [4, 5]. The present article is concerned with the most recent collection of basis sets from the Ahlrichs group [5]. These basis sets were developed as a second generation of default basis sets for the popular TURBOMOLE program [6] and are called the def2 basis sets, in particular def2-SV(P), def2-SVP, def2-TZVP, def2-TZVPP, def2-QZVP, and def2-QZVPP. In these abbreviations, SV denotes split valence (another name for valence double zeta), TZV denotes valence triple zeta, QZV denotes valence quadruple zeta, (P) denotes partially polarized, P denotes polarized, and PP denotes heavily polarized. At the risk of oversimplifying their recommendation, we summarize the recommendation of Weigend and Ahlrichs by saying that they recommend P-type basis sets for density functional theory (DFT) and PP-type basis sets for correlated wave function calculations. The present article is primarily concerned with DFT, although the basis sets proposed here can also be used with wave function theory (WFT). We note that all the basis sets described here can be used in any electronic structure package, not just TURBOMOLE; for example, our group uses them routinely in Gaussian and MOLPRO. 123 Author's personal copy 296 Theor Chem Acc (2011) 128:295–305 Table 1 Smallest s and p exponential parameters in some common basis sets Zeta Basis Group Diffuse? Carbon s Oxygen p s p Double 6-31G* Pople No 0.168714 0.168714 0.270006 0.270006 Double 6-31?G* Pople Yes 0.0438 0.0438 0.0845 0.0845 Double cc-pVDZ Dunning No 0.1596 0.1517 0.3023 0.2753 Double aug-cc-pVDZ Dunning Yes 0.0469 0.04041 0.07896 0.06856 Double def2-SV(P) Ahlrichs No 0.130902 0.152686 0.255308 0.276415 Double def2-SVP Ahlrichs No 0.130902 0.152686 0.255308 0.276415 Triple Triple 6-311G(2d,p) 6-311?G(2d,p) Pople Pople No Yes 0.145585 0.0438 0.145585 0.0438 0.255611 0.0845 0.255611 0.0845 Triple cc-pVTZ Dunning No 0.1285 0.1209 0.2384 0.214 Triple aug-cc-pVTZ Dunning Yes 0.04402 0.03569 0.07376 0.05974 Triple def2-TZVP Ahlrichs No 0.095164 0.100568 0.185045 0.174784 Triple def2-TZVPP Ahlrichs No 0.095164 0.100568 0.185045 0.174784 Quadruple cc-pVQZ Dunning No 0.1111 0.1007 0.2067 0.175 Quadruple aug-cc-pVQZ Dunning Yes 0.04145 0.03218 0.06959 0.05348 Quadruple def2-QZVP Ahlrichs No 0.107399 0.075984 0.197727 0.128640 Quadruple def2-QZVPP Ahlrichs No 0.107399 0.075984 0.197727 0.128640 Parameters with digits beyond the millionths place are round to the nearest 0.000001 for this table The def2 basis sets are very appealing because they are designed to provide consistent accuracy across the whole periodic table and they are available for all elements up to radon (Z = 86) [5], whereas this is not true for the more commonly used basis sets from the Pople and Dunning groups. We have begun using the def2 basis sets for these reasons and have obtained generally satisfactory (often excellent) results, but in our opinion, two aspects of the design of these basis sets required further study, namely that (1) they were designed without provision of diffuse functions and (2) they do not allow for scalar relativistic effective core potentials in the fourth period of the periodic table (especially the 3d and 4p blocks), where relativistic effects are beginning to become chemically important. (They do include a relativistic effective core potential starting with Rb, at the beginning of the fifth row). The present article is concerned with issue (1). The Turbomole manual [7] states that ‘‘Diffuse functions should only be added if really necessary. E.g., for small anions or treatment of excited states use: TZVP instead of SVP ? diffuse….’’. The manual recommends adding diffuse functions only for excited states of small molecules or excited states with (a partial) Rydberg character, and for such calculations, the manual recommends switching to Dunning’s augmented basis sets. Using the augmented Dunning basis sets usually constitutes costly overkill for DFT [8, 9]; furthermore, this does not correspond to our experience that adding diffuse functions to double zeta basis sets is more important than going to a 123 triple zeta basis set for anions and often for other properties [10], and we decided to make further investigations; these investigations are reported here. The first step is to examine the basis sets themselves. The ‘‘diffuseness’’ of a basis set is measured by the value of the smallest exponential parameter for each angular momentum, with smaller exponential parameters leading to a basis capable of expanding a more diffuse charge distribution. As an example, Table 1 shows the smallest exponential parameters for s and p subshells of several popular basis sets for carbon and oxygen. The table shows that the def2 subshells are usually but not always slightly more diffuse than standard nondiffuse basis sets, but not nearly as diffuse as standard diffuse ones. The observation that the def2 basis sets are slightly more diffuse than the unaugmented sets of Pople or Dunning but less diffuse than their augmented ones raises the question of whether they are diffuse enough that the recommendation in the manual should be followed. This is the motivation for the systematic tests presented here. 2 Augmented Karlsruhe basis sets In order to make diffuse versions of the def2 basis sets, diffuse s and p Gaussian functions were added to the polarized def2 Turbomole basis sets series [5], i.e., split valence (SV) bases, triple zeta valence (TZV) bases, and quadruple zeta valence (QZV) bases on all atoms except H. After some experimentation (see Appendix), we settled on Author's personal copy Theor Chem Acc (2011) 128:295–305 a universal prescription for the exponential parameters of the added functions, in particular we took them to be equal to the exponential parameter of the most diffuse s or p function in the corresponding def2 basis divided by a factor of 3. A prefix ‘‘ma’’ is used to replace the ‘‘def2’’ prefix to denote this augmentation. For example, maTZVPP is the basis def2-TZVPP plus a set of diffuse s and p functions on all atoms except H. Adding only diffuse s and p functions follows the strategy of constructing ‘‘?’’ and ‘‘maug’’ basis sets that was recommended for DFT calculations in previous studies [8, 9], and we call this minimal augmentation; thus ‘‘ma’’ may be read as minimally augmented. Note that def2-QZVP and def2-QZVPP are identical for elements studied in Sect. 3 (and in fact for all elements up to Ar [Z = 18]), and hence their minimally augmented counterparts are also identical, but they are not identical for many other atoms, e.g., Cs. Dividing the most diffuse exponential parameter by three sometimes produces a repeating decimal. In such cases, we round the parameter to the nearest 0.00000000001. 3 Diffuse functions for nonmetals H–Cl 3.1 Computational details The minimally augmented def2 basis sets (the ‘‘ma’’ … basis sets) together with standard def2 basis sets were tested against four databases: the 24 diverse barrier heights database [11, 12] (DBH24/08), the IP13/3, and EA13/3 databases with respectively 13 ionization potentials and 13 electron affinities [13], and the S22A database [14–16] (S22A) containing the noncovalent interaction energies of 22 small diverse complexes. The DBH24/08 database is a collection of four sub-databases: HATBH6 containing six barrier heights for heavy-atom transfer reactions with neutral reagents (e.g., H ? N2O ? OH ? N2), NSBH6 containing six barrier heights for nucleophilic substitution reactions of anions (e.g., F- ? CH3Cl ? FCH3 ? Cl-), UABH6 containing six barrier heights for unimolecular and association reactions (e.g., H ? C2H4 ? C2H5), and HTBH6 containing six barrier heights for hydrogen-atomtransfer reactions (e.g., H ? H2S ? H2 ? HS). The IP13/ 3 and EA13/3 databases contain atoms, diatoms, and triatoms (e.g., O, Cl2, and PH2). The S22A database is composed of three sub-databases: HB7A with seven hydrogenbonded complexes (e.g., formamide dimer), D8A with eight complexes bound predominantly by dispersion-like interactions, including p–p stacking (e.g., methane dimer and adenine–thymine stack), and M7A with seven complexes whose interaction involves a mixture of hydrogen bonding and dispersion-like interactions (e.g., phenol 297 dimer). All data in all databases considered here are relative Born–Oppenheimer energies (electronic energy plus nuclear repulsion) without vibrational energy; that is we are calculating relative energies of points on potential energy surfaces. The heaviest atom in any of the databases in Sect. 3 is Cl (Z = 17). The calculations of noncovalent interaction energies are reported both with and without the Boys-Bernardi [17] counterpoise correction (CpC); other calculations are reported only uncorrected, as usual. The M06-2X [18] density functional was chosen for most of the calculations in Sect. 3 because of its efficiency and accuracy; the use of a higher accuracy density functional helps to keep the basis set errors from being obscured by the density functional approximation. In order to check the generalization and validity of our conclusions obtained by using the M06-2X functional, some additional calculations with other density functionals, i.e., M06, [18] M06-L [19], xB97 [20], xB97X [20], and xB97X-D [21] were also performed for comparison. The Minnesota functionals used here have percentages of Hartree–Fock exchange from 0 in M06-L to 27 in M06 and to 54 in M06-2X, and the xB97, xB97X, and xB97X-D have 100% Hartree–Fock exchange at large interelectronic separation, and from 0 (xB97) to 15.8 (xB97X) and 22.2 (xB97X-D) at small interelectronic separation. In addition to these density functional calculations, some calculations employing WFT, in particular the Møller–Plesset second-order perturbation approximation (MP2) [22], were also carried out. The results of the MP2 calculations are given only in supporting information. We note that correlated WFT calculations generally require larger basis sets, especially higher-angular-momentum functions, than DFT calculations. In all cases, the same standard geometries were used for calculations with all density functionals and all basis sets and for the MP2 calculations as well. The DBH24, IP13/3, and EA13/3 databases use QCISD/MG3 geometries as explained in the original presentations of these databases [11–13]. In the S22A database, the best estimates of the energies in the original S22 database [13] were improved by the extrapolating all CCSD(T)/complete basis set interaction energies using larger basis sets for the CCSD(T) component of the computation in a recent Takatani et al. [16] work. The S22A geometries are optimized by MP2/ccpVTZ, with or without counterpoise correction (CpC) or by CCSD(T)/cc-pVTZ or CCSD(T)/cc-pVQZ without CpC. An Internet site specifies the optimization method for each individual complex [23]. For comparison with the results obtained with the unaugmented and augmented def2 basis sets, we also calculated some of the quantities using basis sets [10, 24–30] based on the work of Pople and coworkers, in particular MG3S and MG3T. The MG3S basis set is equivalent to a 123 Author's personal copy 298 Theor Chem Acc (2011) 128:295–305 6-311?G(3d2f,2df,2p) basis for H–Si (Z = 1–14) and is an improved version of the 6-311?G(3d2f) basis for P–Ar (Z = 15–18). The MG3T basis is the same as MG3S except that all diffuse functions are removed. Note that S denotes ‘‘semidiffuse’’ (same number of diffuse functions as plus [1], maug [8, 9], or ma), and T denotes ‘‘tight’’ (no diffuse functions). All the electronic structure calculations were performed by using the Gaussian09 package [31] except that some of the MP2 calculations for the S22A database were carried out by the Molpro 2008.1 package [32]. 3.2 Results In all cases, we give two kinds of errors: mean signed error (MSE, which is useful for showing systematic errors) and mean unsigned error (MUE, which is the same as mean absolute error). Tables 2 and 3 give errors for DBH24/08 and its sub-databases. Tables 4 and 5 give errors for S22A and its sub-databases. Table 6 shows errors for EA13/3 and IP13/3. 3.3 Discussion Table 2 shows that the errors in barrier heights for the M06-2X functional are almost always decreased when diffuse functions are added. The errors without diffuse functions are especially large for def2-SV(P) and for the nucleophilic reactions of anions. At the same time, it is noteworthy that the relative errors in omitting diffuse functions with def2-TZVP and def2-TZVPP relative to maTZVP and ma-TZVPP are much smaller than the differences between MG3T and MG3S. The latter observation may be related to the overtightness of the Pople-group unaugmented basis sets, as has been remarked previously [27, 33, 34]. Furthermore, the manual recommendation mentioned in the introduction has some truth to it, in that the errors are sometimes lower with def2-TZVP than with ma-SVP, with a trend that is generally opposite of the analogous trend [35] with Pople-group basis sets. A relevant consideration though is the relative cost of these two options, that is, increasing zeta versus adding s and p diffuse functions to nonhydrogenic atoms. Table 7 shows that the def2-TZVP basis set contains many more contracted basis functions than the ma-SVP basis set. Similarly, Tables 2 and 7 show that the ma-TZVP basis set gives about the same accuracy as the def2-QZVPP one, but with about half as many contracted functions. Note that the number of diffuse functions in the fully augmented Dunning-group basis sets (denoted ‘‘aug’’) is larger than the number in the present minimal augmentation strategy (denoted in various articles as plus, semiduffuse, maug, or ma, depending on the exponential parameters and the context). It is consistent with the recommendations of Weigend and Ahlrichs [5] that def2-TZVPP is not much more accurate than def2TZVP for density functional calculations, and this is also true with the minimally augmented versions. Table 3 shows larger errors with other density functionals, but similar trends in most respects. One interesting observation in both Tables 2 and 3 is that ma-TZVP is more accurate on average than def2-TZVPP, but Table 7 shows that it involves a smaller number of contracted functions. Both Tables 2 and 3 indicate that there is no advantage in going beyond ma-TZVP for density functional calculations. The four larger basis sets than this have about the same errors or, in the one case of def2-TZVPP, even large errors. Table 2 Mean signed and unsigned errors in barrier heights with the M06-2X density functional (kcal/mol) Basis set HATBH6 NSBH6 UABH6 HTBH6 DBH24/08 MSE MUE MSE MUE MSE MUE MSE MUE MUE MG3T -8.74 9.15 -7.10 17.50 0.57 1.69 -0.74 1.50 7.46 MG3S -0.02 0.73 0.60 0.86 0.37 1.09 -0.49 1.24 0.98 def2-SV(P) -0.99 4.05 -4.71 7.97 0.89 1.38 -1.51 2.44 3.96 ma-SV(P) -0.90 3.59 1.07 1.96 1.20 1.25 -1.10 1.69 2.12 def2-SVP -0.69 2.95 -4.79 7.98 0.47 1.14 -1.41 1.60 3.42 ma-SVP -0.42 2.60 1.24 1.94 0.74 1.51 -0.95 1.15 1.80 0.09 0.98 -1.44 2.60 0.17 0.90 -0.36 1.24 1.43 def2-TZVP ma-TZVP -0.21 0.70 1.04 1.04 0.18 0.92 -0.34 1.24 0.98 def2-TZVPP 0.16 0.81 -1.41 2.55 0.20 1.05 -0.46 1.21 1.41 ma-TZVPP def2-QZVPP -0.11 -0.19 0.52 0.45 1.08 -0.09 1.08 1.31 0.21 0.16 1.07 1.03 -0.45 -0.51 1.21 1.13 0.97 0.98 ma-QZVPP -0.21 0.38 1.09 1.09 0.15 1.02 -0.48 1.10 0.90 123 Author's personal copy Theor Chem Acc (2011) 128:295–305 299 Table 3 Mean signed and unsigned errors in barrier heights obtained using five density functionals (kcal/mol) Basis set HATBH6 NSBH6 UABH6 HTBH6 MSE DBH24/08 MSE MUE MSE MUE MSE MUE MUE MUE MG3T -12.22 13.19 -10.60 16.71 0.86 2.42 -4.22 4.34 9.16 MG3S -5.85 6.87 -3.35 3.35 0.52 1.77 -4.11 4.21 4.05 M06-L def2-SV(P) -5.99 9.35 -8.95 9.92 0.74 1.84 -4.21 4.35 6.36 ma-SV(P) -6.34 8.81 -3.27 3.27 0.84 2.09 -4.16 4.24 4.60 def2-SVP -5.71 8.34 -9.58 10.35 0.53 2.35 -4.68 4.81 6.46 ma-SVP def2-TZVP -5.94 -5.62 7.80 7.10 -3.44 -5.18 3.44 5.18 0.60 0.42 2.58 1.57 -4.57 -3.96 4.65 4.12 4.62 4.49 ma-TZVP -5.98 6.91 -2.54 2.54 0.40 1.57 -3.99 4.14 3.79 def2-TZVPP -5.60 6.82 -5.26 5.26 0.44 1.59 -4.01 4.05 4.43 ma-TZVPP -5.95 6.77 -2.61 2.61 0.43 1.59 -4.03 4.06 3.76 def2-QZVPP -5.89 6.64 -3.26 3.26 0.44 1.57 -3.72 3.75 3.80 ma-QZVPP -5.91 6.62 -2.42 2.42 0.42 1.57 -3.69 3.72 3.58 8.00 M06 MG3T -11.33 11.70 -9.52 15.86 0.79 2.52 -1.81 1.95 MG3S -3.62 4.06 -1.61 1.64 0.54 1.91 -1.53 1.66 2.32 def2-SV(P) -4.96 6.04 -6.86 8.06 0.54 1.73 -2.09 2.66 4.62 ma-SV(P) -4.91 5.58 -1.45 2.37 0.86 2.12 -1.70 2.17 3.06 def2-SVP -4.50 5.15 -7.27 8.33 0.35 2.44 -2.34 2.36 4.57 ma-SVP -4.29 4.90 -1.42 2.25 0.62 2.81 -1.93 1.93 2.97 def2-TZVP -3.10 4.02 -3.41 3.49 0.41 1.62 -1.29 1.53 2.66 ma-TZVP def2-TZVPP -3.44 -3.03 4.08 3.78 -0.58 -3.54 1.11 3.54 0.41 0.46 1.64 1.80 -1.25 -1.42 1.48 1.52 2.08 2.66 ma-TZVPP -3.33 3.88 -0.68 1.13 0.47 1.82 -1.38 1.47 2.07 def2-QZVPP -3.32 3.80 -1.71 1.71 0.52 1.70 -1.14 1.39 2.15 ma-QZVPP -3.32 3.76 -0.76 1.27 0.50 1.72 -1.10 1.34 2.02 MG3S 2.42 3.12 1.49 2.09 2.50 3.08 -0.35 2.23 2.63 def2-TZVP 2.52 3.26 -0.86 3.32 2.38 2.93 -0.26 2.27 2.95 ma-TZVP 2.26 2.98 2.00 2.31 2.39 2.94 -0.24 2.25 2.62 def2-TZVPP 2.56 3.04 -0.86 3.38 2.43 3.05 -0.33 2.14 2.90 ma-TZVPP 2.32 2.81 2.02 2.28 2.43 3.06 -0.32 2.15 2.57 MG3S 0.80 2.32 0.83 1.32 1.86 2.34 -1.22 2.18 2.04 def2-TZVP 0.87 2.46 -1.52 2.90 1.70 2.14 -1.14 2.22 2.43 ma-TZVP 0.63 2.19 1.44 1.52 1.70 2.15 -1.13 2.22 2.02 def2-TZVPP 0.96 2.22 -1.51 2.90 1.75 2.31 -1.18 2.10 2.38 0.73 1.96 1.46 1.47 1.75 2.33 -1.17 2.09 1.96 -1.31 2.19 0.10 0.67 1.03 2.04 -2.04 2.24 1.78 xB97 xB97X ma-TZVPP xB97X-D MG3S def2-TZVP -1.13 2.23 -2.11 2.59 0.87 1.85 -1.92 2.21 2.22 ma-TZVP -1.37 1.96 0.82 0.94 0.87 1.86 -1.89 2.19 1.74 def2-TZVPP -1.04 2.00 -2.10 2.55 0.91 2.01 -1.96 2.13 2.18 ma-TZVPP -1.27 1.74 0.84 0.92 0.91 2.03 -1.93 2.11 1.70 123 Author's personal copy 300 Theor Chem Acc (2011) 128:295–305 Table 4 Mean signed and unsigned errors (kcal/mol) in the S22A noncovalent interaction database and its sub-databases with the M06-2X functional Basis HB7A MSE D8A MUE MSE M7A MUE MSE S22A MUE MSE MUE Not counterpoise corrected MG3T -1.67 1.67 -1.45 1.46 -0.69 0.69 -1.28 1.28 MG3S def2-SV(P) 0.35 -2.96 0.68 2.96 -0.68 -1.30 0.70 1.31 -0.17 -0.71 0.26 0.73 -0.19 -1.64 0.55 1.65 ma-SV(P) -0.94 0.94 -1.56 1.60 -0.74 0.74 -1.10 1.11 def2-SVP -2.88 2.88 -1.29 1.31 -0.63 0.63 -1.58 1.59 ma-SVP -0.43 0.48 -1.32 1.38 -0.42 0.43 -0.75 0.79 def2-TZVP 0.24 0.81 -0.12 0.23 -0.08 0.36 0.01 0.46 ma-TZVP 0.55 0.89 -0.03 0.16 0.02 0.35 0.17 0.45 def2-TZVPP 0.15 0.68 -0.04 0.25 0.04 0.34 0.05 0.42 ma-TZVPP 0.46 0.69 0.07 0.18 0.16 0.34 0.22 0.40 MG3T 0.76 0.76 0.02 0.15 0.21 0.32 0.32 0.40 MG3S 0.80 0.85 -0.13 0.17 0.19 0.31 0.27 0.43 def2-SV(P) 0.51 0.96 0.38 0.38 0.42 0.42 0.43 0.58 ma-SV(P) 0.67 1.17 -0.08 0.32 0.31 0.37 0.28 0.61 def2-SVP 0.52 0.71 0.39 0.39 0.41 0.41 0.44 0.50 ma-SVP def2-TZVP 0.53 0.68 0.77 0.88 -0.08 0.21 0.27 0.25 0.27 0.24 0.32 0.38 0.22 0.37 0.44 0.49 Counterpoise corrected ma-TZVP 0.76 0.90 0.21 0.25 0.26 0.39 0.40 0.50 def2-TZVPP 0.55 0.72 0.24 0.27 0.28 0.39 0.35 0.45 ma-TZVPP 0.59 0.72 0.23 0.26 0.28 0.39 0.36 0.45 Although it is not the main purpose of this article, it is also interesting to compare ma-TZVP to MG3S, with which we were very satisfied in the past. Table 7 shows that ma-TZVP is smaller, but Table 2 shows it is equally accurate for M06-2X, and Table 3 shows it is more accurate for M06-L and M06. Since it is also defined for a larger number of atoms than MG3S, it is a reasonable choice for general-purpose density functional calculations of barrier heights. Table 4 shows that diffuse functions decrease the counterpoise correction (that is, the results without CpC are closer to the corrected results when diffuse functions are included), especially for def2-SV(P) and def2-SVP. This is a very significant advantage of these minimally augmented basis sets since counterpoise corrections are known to be problematic [36]. The effect of diffuse functions on the uncorrected interactions energies is large for the triple zeta MG3S basis but small for the ma-TZVP basis, again indicating that diffuse functions are less important with the Ahlrichs-group basis sets than with the Pople-type basis sets. Table 5 shows similar trends to Table 4. Comparison of the results in the tables to those in Table 7 shows that the most efficient basis set for density 123 functional calculations of noncovalent interaction energies is def2-TZVP. A comment is in order on the good performance of the xB97X-D density functional in Table 5. This density functional includes an empirical molecular mechanics term for damped dispersion interactions, so the good performance comes partly from the empirical term and only partly from the density functional itself. The table shows that this strategy, originally proposed by Toennies and coworkers [37, 38], is very successful for noncovalent interaction energies, but caution is advised for any method containing empirical functional forms. Table 6 shows that electron affinities are much worse without diffuse functions, and ionization potentials are slightly better. This is consistent with our previous work with other kinds of basis sets. As for barrier heights, maTZVP appears to be a good choice for electron affinities, and as for noncovalent interactions, def2-TZVP appears to be a good choice for ionization potentials. It is noteworthy that ma-TZVP is more accurate than MG3S for both ionization potentials and electron affinities. Since boron atom is not included in the databases discussed above, we made a test of the adequacy of ma-TZVP Author's personal copy Theor Chem Acc (2011) 128:295–305 301 Table 5 Mean signed and unsigned errors (kcal/mol) in the S22A noncovalent interaction database and its sub-databases with the xB97X-D and xB97X density functionals Basis HB7A D8A MSE MUE M7A MSE MUE S22A MSE MUE MSE MUE xB97X Not counterpoise corrected MG3S def2-TZVPP -0.45 -0.58 0.64 0.76 0.72 1.03 0.85 1.10 0.03 0.09 0.34 0.39 0.13 0.22 0.62 0.76 ma-TZVPP -0.25 0.67 1.12 1.18 0.19 0.38 0.39 0.76 MG3S -0.08 0.64 1.12 1.20 0.30 0.37 0.48 0.76 def2-TZVPP -0.14 0.68 1.31 1.35 0.34 0.40 0.54 0.83 ma-TZVPP -0.10 0.66 1.29 1.33 0.33 0.40 0.54 0.82 Counterpoise corrected xB97X-D Not counterpoise corrected MG3S -0.42 0.42 -0.85 0.85 -0.47 0.47 -0.59 0.59 def2-TZVPP -0.45 0.45 -0.34 0.34 -0.32 0.32 -0.37 0.37 ma-TZVPP -0.10 0.20 -0.22 0.22 -0.18 0.19 -0.17 0.20 0.03 0.13 -0.32 0.32 -0.12 0.14 -0.15 0.20 -0.01 0.18 -0.06 0.17 -0.07 0.14 -0.04 0.16 0.05 0.18 -0.07 0.17 -0.06 0.14 -0.03 0.16 Counterpoise corrected MG3S def2-TZVPP ma-TZVPP Table 6 Mean signed and unsigned errors (in kcal/mol) for electron affinities and ionization potentials in the IP13/3 and EA13/3 databases with the M06-2X functional Table 7 Number of contracted basis functions/primitive Gaussian functions of two largest systems studied in the DBH24/08 database and the S22A database, respectively Basis set Basis set ClCH3!!!Cl- Adenine thymine complex MG3T EA IP MSE MUE MSE MUE MG3T 157/237 669/1033 10.67 10.74 0.59 2.09 MG3S 169/249 745/1109 56/111 68/123 288/519 464/595 MG3S 1.26 2.09 1.08 2.52 cc-pVTZ 9.85 9.85 0.70 2.39 def2-SV(P) ma-SV(P) aug-cc-pVTZ 0.53 1.55 1.12 2.70 def2-SVP 65/120 321/552 77/132 397/628 def2-SV(P) 23.18 23.18 0.90 2.97 ma-SVP ma-SV(P) -0.10 2.99 3.38 4.31 def2-TZVP 123/213 655/1057 3.09 ma-TZVP 135/225 731/1133 4.21 def2-TZVPP 157/240 743/1156 169/252 819/1232 def2-SVP ma-SVP 23.19 -0.18 23.19 2.90 0.75 3.28 def2-TZVP 7.23 7.23 0.41 2.04 ma-TZVPP ma-TZVP 0.88 1.64 0.59 2.17 def2-QZVPP 287/448 1413/2166 299/460 1489/2242 def2-TZVPP 7.18 7.18 0.45 2.05 ma-QZVPP ma-TZVPP 0.81 1.56 0.64 2.19 def2-QZVPP 3.54 3.60 0.71 2.33 We use spherical harmonic polarization functions, not Cartesian functions, for all calculations in this article ma-QZVPP 0.92 1.57 0.79 2.37 for boron. In particular, the test is the electron affinity of boron atom, which is a very severe test since electron affinity is the property most sensitive to diffuse functions and since an unbonded atom, unlike an atom in a molecule, has no neighboring atoms whose basis functions overlap it, thereby diminishing the need for diffuse functions. We used the M06-2X functional to calculate the boron electron affinity with the ma-TZVP, aug-cc-pVTZ, and aug-ccpV5Z basis sets. The ma-TZVP, aug-cc-pVTZ, and aug-ccPV5Z basis sets give the boron electron affinity as 5.60, 5.55, and 5.88 kcal/mol, respectively. These values may be compared with the experimental value of 6.45 kcal/mol 123 Author's personal copy 302 [39] and the def2-TZVP value of 1.42 kcal/mol. We conclude that the ma-TZVP basis is also reasonable for boron. The MP2 calculations in supporting information show results for DBH24/08, S22A, EA13/3, and IP13/3 for the interested reader. MP2 is less accurate for barrier heights and noncovalent interactions than the modern density functionals included in the present paper, and so the MP2 results are only of interest for specialized purposes. Theor Chem Acc (2011) 128:295–305 For testing purposes, we specify two other basis sets by dividing the elements into three categories: (i) metals and semimetals, (ii) nonhydrogenic nonmetals, and (iii) hydrogen. A prefix ‘‘nma’’ (denoting ‘‘nonmetals minimally augmented’’) indicates that diffuse s and p functions are added only to elements in category (ii), e.g., N, O, F, and Cl. A prefix ‘‘mma’’ (denoting ‘‘metals minimally augmented’’) indicates that diffuse s and p functions are added only to elements in category (i), e.g., the metals Li and Cs and the semimetals Sb and Bi. 4 Diffuse functions for other elements 4.2 Results and discussion Section 3 addressed elements with Z = 1, 6–9, and 15–17, none of which are metals or semimetals, plus boron (Z = 5). Several questions arise next: should we add diffuse s and p subshells to all nonhydrogenic nonmetals? Should we also add them to semimetals and/or metals, and—if not—what is the distinction between elements that need diffuse functions and those that do not? To answer these questions, calculations were performed for selected molecules including metal or semimetal elements with minimally augmented def2 basis sets, to compare with the results from standard def2 basis sets; in particular, we examined the dissociation energy (De) of strongly ionic LiCl and LiO and less ionic SbN and BiN and the electron affinities (EA) of strongly ionic LiCl and CsF and less ionic SbN and CsH. No molecules involving the elements in period four are included here because the conclusions based on comparison to experiments for molecules containing those elements could be affected by omission of scalar relativistic effects. 4.1 Computational details The broadly applicable M06 density functional has been used for most of the calculations of this section. For comparison, the M06-2X, B3LYP [40], B1LYP [41–43], and BLYP [44, 45] functionals were employed for some of the calculations. The geometries of LiCl (2.04365733 Å) and LiO (1.6069033 Å) come from Database/3 [12], where they were obtained by QCISD/MG3 calculations. Geometries optimized at the CCSD(T)/cc-pV5Z-PP level by Peterson [47] were used in calculations of SbN and BiN. For CsF (2.3756 Å) and CsH (2.5114 Å), the bond length was optimized by the M06 density functional with the maQZVPP basis set. The ma-QZVPP basis set is very large and is used as a standard for comparison to smaller basis sets, except for CsH where we use QZVPP& defined as ma-QZVPP further augmented by diffuse s and p functions on H (again with exponents determined by the factor-of-3 prescription). For the convenience of the reader, the largest basis in each case of Tables 8 and 9 is labeled NCBS for ‘‘nearly complete basis set.’’ 123 Table 8 shows the dissociation energy for LiO, LiCl, SbN, and BiN calculated by the M06 and other density functionals and accurate [46, 47] results. The corresponding NCBS limit is also indicated as a separate row for convenience of the reader. Because the def2-QZVP and def2QZVPP basis sets are the same for all elements in these four molecules, the NCBS limit is same as ma-QZVP. M06 and M06-2X give similar results for LiO and LiCl. Inspection of Table 8 shows a large deviation of B1LYP from BLYP for SbN and BiN, which we have previously used [46] as a diagnostic for significant multireference character, which is why M06-2X should not be applied to these systems. In contrast, both M06 and M06-2X appear accurate for LiO and LiCl. However, the applicability of the various density functionals is a side point. The main point is that all results in Table 8 indicate that, except for SVP results for LiO and LiCl, adding diffuse s and p functions on both elements forming the bond improves the accuracy compared with the NCBS limit, but much less than going to the bigger def2 basis set. Thus, we suggest that one does not need diffuse functions on either element for De calculations. The bonds investigated here vary from strongly ionic to weakly ionic, and furthermore, the conclusion agrees with that drawn for bond energies in two previous studies [10, 11], and they answer the question of whether diffuse functions on nonmetals might be more important when they are bonded to metals (because diffuse functions might be more important for atoms that have negative partial charge in a highly ionic bond). Therefore, we can recommend as a general conclusion that no diffuse functions are needed for calculating bond dissociation energies. Table 9 shows the calculated electron affinities of LiCl, SbN, CsF, and CsH. The results show that adding diffuse functions on both metal and nonmetal elements is indispensable in EA calculations. For the strongly ionic bonds, Li–Cl and Cs-F, diffuse functions are more effective on metal elements than on nonmetal elements. Adding diffuse functions only on metal elements does almost as much good as adding diffuse functions on both elements. This is Author's personal copy Theor Chem Acc (2011) 128:295–305 Table 8 The calculated dissociation energies (De, kcal/ mol) for LiO, LiCl, SbN, and BiN by the M06 and other density functionals and accurate results a ma-QZVPP b From the MLBE21105 database [46] 303 LiO M06 LiCl M06 SbN M06 BiN M06 def2-SVP 83.79 112.74 79.64 70.01 ma-SVP 88.61 115.13 80.70 71.17 def2-TZVP 81.85 112.11 86.89 78.03 ma-TZVP 82.31 112.30 86.85 78.11 def2-QZVP 84.16 113.31 88.35 79.05 ma-QZVP 84.26 113.28 88.46 79.10 NCBSa 84.26 113.28 88.46 79.10 M06-2X M06-2X B3LYP, B1LYP, BLYP B3LYP, B1LYP, BLYP def2-SVP 80.93 113.23 75.19, 68.73, 93.59 59.62, 52.51, 79.36 ma-SVP 84.76 115.88 76.22, 69.90, 94.05 60.81, 54.14, 80.66 def2-TZVP 82.81 111.97 81.58, 75.33, 98.97 66.46, 59.88, 85.82 ma-TZVP 83.26 112.11 82.01, 75.79, 99.42 67.10, 60.54, 86.48 def2-QZVP 84.67 113.41 82.95, 76.72, 100.36 67.80, 61.23, 87.15 ma-QZVP 84.90 113.48 82.98, 76.76, 100.38 67.86, 61.30, 87.20 NCBSa 84.90 113.48 82.98, 76.76, 100.38 67.86, 61.30, 87.20 Experiment Experiment WFT WFT c CCSD(T)/cc-pV5Z-PP from Peterson [47] Accurate 82.0 b Table 9 The calculated electron affinities (EA, kcal/mol) of LiCl, SbN, CsF, and CsH LiCl SbN M06-2X B3LYP M06 CsF CsH B3LYP B3LYP B3LYP def2-SVP -3.44 -8.22 -5.04 -5.92 nma-SVP -7.78 -12.45 -12.45 -16.36 -6.04 -5.92 mma-SVP ma-SVP -9.66 -10.98 -13.51 -11.39 -16.47 -6.72 -15.69 -17.60 -22.22 -7.58 -7.41 -7.41 def2-TZVP -6.48 -10.56 -13.46 -18.05 -4.97 -5.83 nma-TZVP -7.58 -11.88 -16.12 -21.02 -5.46 -5.83 mma-TZVP -11.26 -14.95 -15.59 -20.50 -6.74 -7.36 ma-TZVP -11.69 -15.73 -17.61 -22.77 -7.09 -7.36 def2-QZVP -9.91 -14.01 -15.60 -20.90 -6.11 -6.73 nma-QZVP -10.16 -14.34 -16.19 -21.73 -6.29 -6.73 mma-QZVP -11.30 -15.33 -16.38 -21.91 -6.95 -7.47 ma-QZVP -11.45 -15.56 -17.39 -23.01 -7.07 -7.47 -9.91 -14.01 -15.60 -20.90 -6.12 -6.78 nma-QZVPP -10.16 -14.34 -16.19 -21.73 -6.27 -6.78 mma-QZVP -11.30 -15.33 -16.38 -21.91 -6.95 -7.50 ma-QZVPP -11.45 -15.56 -17.39 -23.01 -7.05 -7.50 NCBSa -11.45 -15.56 -17.39 -23.01 -7.05 -7.55 def2-QZVPP -7.39 -3.57 a For LiCl, SbN, and CsF, NCBS denotes ma-QZVPP; for CsH, NCBS denotes the QZVPP & basis set, where a set of s and p diffuse functions has been added to both Cs and H a result of the almost full valence shells of Cl and F in these strongly ionic bonds. For the less ionic bond Sb-N, the effects of diffuse functions on the semi-metal element and b 113.9 82.33 c 68.10c the nonmetal element are similar, and both are important and have to be included for reliable calculations of EA values. From Table 9, we can notice that the ma-TZVP basis set already approaches the NCBS limit, and it gives a remarkably better result than the more expensive def2QZVP basis set. The same conclusion is obtained as in Sect. 3: ma-TZVP is an excellent choice for EA calculations. In Table 9, for CsH, the NCBS results include a set of s and p diffuse functions not only on Cs but also on H. The added diffuse functions on H do not change the calculated EA very much; hence, we confirm that no diffuse functions are needed for H, even for the electron affinity of the metal hydrides where H is bonded to the most electropositive metal (considering only Z B 86, which is the scope of this paper). 5 Conclusions We conclude that adding diffuse s and p functions on nonhydrogenic atoms is a more efficient way to increase the accuracy of density functional calculations of barrier heights and electron affinities with def2 basis sets than either increasing the valence space (increasing zeta) or adding extra polarization functions (going from def2TZVP to def2-TZVPP). The new ma-TZVP basis set presented here is highly recommended for barrier height 123 Author's personal copy 304 Theor Chem Acc (2011) 128:295–305 Table 10 Mean signed unsigned errors in barrier heights (kcal/mol) using two different schemes for augmentation Basis set HATBH6 MSE NSBH6 MUE MSE UABH6 MUE MSE HTBH6 MUE MSE DBH24/08 MUE MUE maGE-TZVPPa M06-2X -0.13 0.50 1.09 1.09 0.21 1.07 -0.45 1.21 0.97 M06 -3.36 3.90 -0.65 1.15 0.46 1.82 -1.38 1.47 2.08 M06-L -5.96 6.78 -2.59 2.59 0.42 1.58 -4.03 4.06 3.75 ma-TZVPPb M06-2X -0.11 0.52 1.08 1.08 0.21 1.07 -0.45 1.21 0.97 M06 M06-L -3.33 -5.95 3.88 6.77 -0.68 -2.61 1.13 2.61 0.47 0.43 1.82 1.59 -1.38 -4.03 1.47 4.06 2.07 3.76 a Diffuse functions obtained by geometric extension (GE): the exponential parameter of the added s or p diffuse function is equal to the exponential parameter of the most diffuse s or p function in def2-TZVPP divided by the ratio of the exponents of the most two most diffuse s or p functions in the def2-TZVPP b The exponential parameter of the added s or p diffuse function is equal to the exponential parameter of the most diffuse s or p function in def2TZVPP divided by 3 and electron affinity calculations with density functional theory. Diffuse functions are less important for noncovalent interactions and ionization potentials than for barrier heights and electron affinities, and def2-TZVP appears adequate for bond energies, noncovalent interactions, and ionization potentials. The def2-TZVP basis is also adequate for some of the barrier heights; for example, it has good average performance for hydrogen-atom-transfer reactions. For large systems or systems involving anions, it is recommended to try ma-SVP before increasing the size of the valence space to valence triple zeta. There is little advantage in going to quadruple zeta for density functional theory, at least for all of the cases in this paper. The Turbomole manual recommends taking diffuse functions, when needed, from the aug-cc-pV…Z basis sets, but these are only available up to Kr (Z = 36), not including K and Ca (Z = 19 and 20). Furthermore, the augcc-pV…Z sets include higher-angular-momentum diffuse functions that we have found, both here and in previous work [8, 9], to be unnecessary for most DFT calculations, and they include diffuse functions on H, which are not recommended. Therefore, we have defined minimally augmented def2 basis sets (ma basis sets) for all elements up to radon (Z = 86), and the present tests show that maTZVP appears to be an excellent choice for general-purpose applications of density functional theory, even better than the MG3S basis with which we previously had very good success. We recommend ma-TZVP for general use and ma-SVP or 6-31?G(d, p) when ma-TZVP is unaffordable, and the tests presented provide a validation for omitting diffuse functions for properties for which they are less important. 123 Acknowledgments This work was supported in part by the U. S. Department of Energy, Office of Basic Energy Sciences, under grant no. DE-FG02-86ER13579 and by the Air Force Office of Scientific Research under grant no. FA9550-08-1-0183. Appendix Table 10 shows errors in the DBH24/08 barrier heights for three density functionals with two different schemes for extending the basis sets; the errors are nearly the same when the basis is extended by a geometric series in the exponential parameters or by the simpler scheme of a factor of 3. These and other less systematic considerations led us to conclude that the simple scheme of dividing by 3 is adequate for most purposes, and we therefore adopted that scheme for our standard definition of the augmented def2 basis sets. References 1. Hehre WJ, Radom L, Schleyer P, Pople JA (1986) Ab Initio Molecular Orbital Theory. New York, John Wiley & Sons 2. Dunning T H Jr (1989) J Chem Phys 90:1007 3. Wilson AK, Woon DE, Peterson KA, Dunning T H Jr (1999) J Chem Phys 110:7667 4. Weigend F, Furche F, Ahlrichs R (2003) J Chem Phys 19:12753 5. Weigend F, Ahlrichs R (2005) Phys Chem Chem Phys 7:3297 6. Ahlrichs R et al (2010) Turbomole—program package for ab initio electronic structure calculations. http://www.turbomole. com. Accessed 21 April, 2010 7. Ahlrichs R et al (2010) Turbomole user’s manual, Version 6.0, Feburary 3, 2009. http://www.cosmologic.de/data/DOK.pdf. Accessed 30, April, 2010 8. Papajak E, Leverentz H, Zheng J, Truhlar DG (2009) J Chem Theory Comput 5:1197 Author's personal copy Theor Chem Acc (2011) 128:295–305 9. Papajak E, Truhlar DG (2010) J Chem Theory Comput 6:597 10. Lynch BJ, Zhao Y, Truhlar DG (2003) J Phys Chem A 107:1384 11. Zheng J, Zhao Y, Truhlar DG (2007) J Chem Theory Comput 3:569 12. Zheng J, Zhao Y, Truhlar DG (2009) J Chem Theory Comput 5:808 13. Lynch BJ, Truhlar DG (2003) J Phys Chem A 107:3898 14. Jurecka P, Sponer J, Cerny J, Hobza P (2006) Phys Chem Chem Phys 8:1985 15. Marchetti O, Werner H-J (2009) J Phys Chem A 113:11580 16. Takatani T, Hohenstein EG, Malagoli M, Marshall MS, Sherrill CD (2010) J Chem Phys 132:144104 17. Boys SF, Bernardi F (1970) Mol Phys 19:553 18. Zhao Y, Truhlar DG (2008) Theor Chem Acc 120:215 19. Zhao Y, Truhlar DG (2006) J Chem Phys 125:194101 20. Chai J-D, Head-Gordon M (2008) J Chem Phys 128:084106 21. Chai J-D, Head-Gordon M (2008) Phys Chem Chem Phys 10:6615 22. Møller C, Plesset MS (1934) Phys Rev 46:618 23. Hobza P (2010) http://www.begdb.com. Accessed 30 April, 2010 24. Krishnan R, Binkley JS, Seeger R, Pople JA (1980) J Chem Phys 72:650 25. McLean AD, Chandler GS (1980) J Chem Phys 72:5639 26. Clark T, Chandrasekhar J, Spitznagel GW, Schleyer PV (1983) J Comput Chem 4:294 27. Frisch MJ, Pople JA, Binkley JS (1984) J Chem Phys 80:3265 28. Curtiss LA, Raghavachari K, Redfern PC, Rassolov V, Pople JA (1998) J Chem Phys 109:7764 29. Curtiss LA, Redfern PC, Raghavachari K, Rassolov V, Pople JA (1999) J Chem Phy 110:4703 30. Fast PL, Sanchez ML, Truhlar DG (1999) Chem Phys Lett 306:407 31. Frisch MJ, Trucks GW, Schlegel HB, Scuseria GE, Robb MA, Cheeseman JR, Scalmani G, Barone V, Mennucci B, Petersson GA, Nakatsuji H, Caricato M, Li X, Hratchian HP, Izmaylov AF, Bloino J, Zheng G, Sonnenberg JL, Hada M, Ehara M, Toyota K, Fukuda R, Hasegawa J, Ishida M, Nakajima T, Honda Y, Kitao O, Nakai H, Vreven T, Montgomery Jr. JA, Peralta JE, Ogliaro F, 305 32. 33. 34. 35. 36. 37. 38. 39. 40. 41. 42. 43. 44. 45. 46. 47. Bearpark M, Heyd JJ, Brothers E, Kudin KN, Staroverov VN, Kobayashi R, Normand J, Raghavachari K, Rendell A, Burant JC, Iyengar SS, Tomasi J, Cossi M, Rega Millam NJ, Klene M, Knox JE, Cross JB, Bakken V, Adamo C, Jaramillo J, Gomperts RE, Stratmann O, Yazyev AJ, Austin R, Cammi C, Pomelli, JW, Ochterski R, Martin RL, Morokuma K, Zakrzewski VG, Voth GA, Salvador P, Dannenberg JJ, Dapprich S, Daniels AD, Farkas O, Foresman JB, Ortiz, JV, Cioslowski J, Fox DJ (2009) Gaussian 09. Revision A 02, Gaussian Inc., Wallingford, CT Werner H-J, Knowles PJ, Lindh R, Manby FR, Schutz M, Celani P, Korona T, Mitrushenkov A, Rauhut G, Adler TB, Amos RD, Bernhardsson A, Berning A, Cooper DL, Deegan MJO, Dobbyn AJ, Eckert F, Goll E, Hampel C, Hetzer G, Hrenar T, Knizia,G, Koeppl C, Liu Y, Lloyd AW, Mata RA, May AJ, McNicholas SJ, Meyer W, Mura ME, Nicklass A, Palmieri P, Pfueger K, Pitzer R, Reiher M, Schumann U, Stoll H, Stone AJ, Tarroni R, Thorsteinsson T, Wang M, Wolf A (2008) MOLPRO, version 2008.1,Universität Stuttgart, Stuttgart, Germany Del Bene JE, Aue DH, Shavitt I (1992) J Am Chem Soc 114:1631 Lynch BJ, Zhao Y, Truhlar DG (2005) J Phys Chem A 109:1643 Lynch BJ, Zhao Y, Truhlar DG (2003) J Phys Chem A 107:1384 Alvarez-Idaboy JR, Galano A (2010) Theo Chem Acc 126:75 Tang KT, Toennies JP (1984) J Chem Phys 80:3726 Paesani F, Gianturco F, Lewerenz M, Toennies JP (1999) J Chem Phys 111:6897 Scheer M, Bilodeau RC, Haugen HK (1998) Phys Rev Lett 80:2562 Stephens PJ, Devlin FJ, Chabalowski CF, Frisch MJ (1994) J Phys Chem 98:11623 Lee C, Yang W, Parr RG (1988) Phys Rev B 37:785 Becke AD (1988) Phys Rev A 38:3098 Adamo C, Barone V (1999) J Chem Phys 110:6158 Adamo C, Cossi M, Barone V (1999) Theochem 493:147 Ernzerhof M, Scuseria GE (1999) J Chem Phys 110:5029 Schultz NE, Zhao Y, Truhlar DG (2005) J Phys Chem A 109:11127 Peterson KA (2003) J Chem Phys 119:11099 123