Systematic Model Chemistries Based on Density Functional Theory: Comparison with traditional Models and with Experiment

J.M. Seminario (Editor) Recent Developments and Applications of Modern Density Functional Theory

Theoretical and Computational Chemistry, Vol. 4 9 1996 Elsevier Science B.V. All rights reserved.

679

S y s t e m a t i c M o d e l C h e m i s t r i e s B a s e d on D e n s i t y F u n c t i o n a l T h e o r y : Comparison with traditional Models and with Experiment Michael J. Frisch, Gary W. Trucks and James R. Cheeseman Lorentzian, Inc., 140 Washington Avenue, North Haven, CT 06472, USA Model chemistries based on DFT methods are compared to one another and to those based on traditional ab initio methods such as Hartree-Fock and MP2 for a variety of job types including energies, structures, vibrational frequencies and related properties and NMR chemical shifts. In general, model chemistries based on hybrid functionals perform best of all DFT models, many times achieving equal or better accuracy to the more computationally expensive MP2 theory.

1. INTRODUCTION Kohn-Sham density functional theory (DFT) [1, 2] has been used increasingly frequently in recent years to model molecular systems and chemical reactions. Very good results have been reported by many different researchers on a wide variety of compounds and problems [3-25]. In such cases, various DFT methods have proven to provide an effective tradeoff between computational cost and accuracy. However, as with any computational method, model chemistries based on DFT methods require systematic calibration in order to determine their range of applicability and expected accuracy for previously-unstudied molecular systems. This chapter will examine such results. Current implementations of DFT methods replace the exchange-correlation energy by a one-electron integral of the electron spin densities and possibly their gradients; functionals involving only the density itself are termed "local," while those involving both the density and its gradient are termed "gradient-corrected." Many different functionals have been proposed, and separate functionals have been devised for the exchange and correlation portions of this term. These pure functionals are the ones which will be considered: SVWN5, comprised of the standard local exchange functional [1, 2, 26] and the local correlation functional of Vosko, Wilk and Nusair (VWN5, denoted functional V in the paper) [27].

680 9 SVWN3, comprised of the standard local exchange functional [1, 2, 26] and an alternate local correlation functional of Vosko, Wilk and Nusair (VWN3, denoted functional III in the paper) [27]. 9 BLYP, which combines Becke's 1988 gradient-corrected exchange functional (B88) [28] with the gradient-corrected correlation functional of Lee, Yang and Parr (LYP) [29, 30]. In addition to these "pure" functionals, Becke has introduced functionals comprised of a mixture of Hartree-Fock and DFT exchange along with a gradientcorrected DFT correlation functional. Becke has proposed several of hybrid functionals [31-33], one of which has the following form [31]:

Exc

=

]~LSDA

--'xc + ao(Exexact

. -. , ~ L S D A

--.x

_ A~BSS ) + axLU~x + acL~EcPWgl

(~)

This functional, denoted B3PW91, is a mixture of exact (Hartree-Fock) and DFT exchange (the latter is Becke's gradient-corrected exchange functional [28]) and the 1991 gradient-corrected correlation functional of Perdew and Wang (PW91) [34]. The three constants were determined by fitting to 116 calculations drawn from the G1 molecule set [35] and its components (atomization energies, ionization potentials, proton affinities and first-row atomic energies). Becke arrived at values of %=0.20, ax=0.72, and ac=0.81 for the three parameters. Two hybrid functionals are considered in this work: B3PW91, as defined in Equation 1. 9 B3LYP, a similar hybrid functional correlation functional as follows:

constructed

EB3LYP= I~LSDA (E~F _ I~LSD^ _ ,,~BSS EcVWN3 xc ~'x +a0 -~x ) + a x ~ x + +ac(Ec

LYP

using

--

the

EcVWN3 )

LYP

(2)

The parameters are again those suggested by Becke. Practical model chemistries based on DFT methods will combine a pure or hybrid functional with a specific basis set, creating a computational p r o c e d u r e which may be generally applied to chemical problems and whose results may be compared for related compounds. However, most systematic studies of density functional methods were completed prior to the widespread use of hybrid functionals. Our goal here is to consider calibration studies which include the whole range of DFT functionals using a variety of basis sets in an effort to determine which are most applicable to the chemical problems of greatest current interest.

681 2. GEOMETRIES AND ENERGIES

The choice of a model for the optimized molecular structure is a natural starting point for any chemical investigation undertaken with electronic structure theory. Thus, this section begins by considering the performance of several model chemistries involving DFT methods for geometry optimizations. It then goes on to discuss models for predicting high-accuracy energies and thermochemical properties of molecular systems. 2.1 Optimized Geometries Large scale, systematic studies of optimized geometries using the full range of DFT functionals are still in preparation by several researchers as of this writing. Nevertheless, the performance of a range of functionals on several compounds known to cause difficulty for traditional ab initio methods can be considered. Murray and coworkers reported results from calculations for this set of several small molecules known to be troublesome at the Hartree-Fock level for a variety of DFT methods, using the 6-31G(d), DZ2P and TZ2P basis sets [14]. Their optimized geometries agree very well with the corresponding calculations in Table I (the HF, MP2, SVWN5 and BLYP optimizations were run in order to use identical basis sets for all calculations). In addition, the table provides results using the B3LYP and B3PW91 hybrid functionals. In general, using the larger, triple-zeta basis set decreases the predicted bond lengths (the sole exception is the bond length in F2 for the BLYP functional). For these systems, Hartree-Fock bond lengths are too short, sometimes quite significantly so. MP2 theory improves on the Hartree-Fock values, often producing bond lengths in excellent agreement with experiment for the larger basis set. Like MP2, the models employing pure DFT functionals tend to predict bond lengths that are longer than the observed values. The models using hybrid functionals produce errors whose sign is less predictable but generally are smaller in magnitude than those for the pure functionals; their values are comparable to or more accurate than the MP2 predictions for all of the bond lengths not involving fluorine. In virtually all cases, the models including the hybrid functionals perform significantly better than the models using SVWN5 and BLYP (the one exception is hydrogen molecule where BLYP with the larger basis set performs about as well as the hybrid functionals). For these molecules, the accuracies of B3LYP and B3PW91 are comparable.

Both Hartree-Fock theory and MP2 provide very good agreement with experiment for the bond angle in formaldehyde. The DFT models using the smaller basis set all significantly underestimate the bond angle in this structure. The bond angle computed by the DFT-based models with the larger basis set tends to be slightly too small (and this is especially true for the hybrid functionals).

682 Table 1 Optimized geometries for several small molecules using traditional ab initio and DFT-based model chemistries (bond lengths in A; angles in degrees).

~-R~ HF MP2 SVWN5 BLYP B3LYP B3PW91 Experiment"

6-31G(d,p) 0.732 0.734 0.765 0.747 0.743 0.743

6-311+G(2d,2p) 0.734 0.737 0.766 0.747 0.743 0.744 0.741

6-31G(d,p) 0.901 0.921 0.933 0.937 0.925 0.922

F2: R~v

N~:R~ HF MP2 SVWN5 BLYP B3LYP B3PW91 Experiment"

HF MP2 SVWN5

6-311+G(2d,2p) 0.897 0.918 0.932 0.933 0.923 0.920 0.917

6-31G(d,p) 1.078 1.131 1.111 1.118 1.106 1.105

6-311+G(2d,2p) 1.067 1.114 1.096 1.104 1.092 1.091 1.098 CO: Rco 6-31G(d,p) 6-311+G(2d,2p) 1.114 1.104 1.151 1.138 1.142 1.128

6-31G (d,p) 1.345 1.421 1.389 1.434 1.403 1.393

6-311 +G (2d,2p) 1.330 1.410 1.388 1.437 1.400 1.388 1.412

Formaldehyde: Rco 6-31G(d,p) 6-311+G(2d,2p) 1.184 1.179 1.220 1.213 1.207 1.199

BLYP

1.150

1.138

1.218

1.213

B3LYP B3PW91 Experiment"

1.138 1.137

1.126 1.125

1.207 1.205

1.201 1.199

HF MP2 SVWN5 BLYP B3LYP B3PW91 Experiment" "

Reference [14].

1.128 Formaldehyde: Rc, 6-31G (d,p) 6-311 +G (2d,2p) 1.093 1.091 1.100 1.098 1.126 1.121 1.120 1.114 1.110 1.105 1.111 1.107 1.099

1.203 Formaldehyde: Ac, c 6-31G (d,p) 6-311 +G (2d,2p) 115.7 116.4 115.5 116.7 115.1 116.3 114.9 116.1 115.2 116.1 115.3 116.1 116.5

683 Table 2 gives the absolute average deviation from experiment for the bond lengths in these compounds for the models under consideration as well as the standard deviation of the absolute errors. MP2 with the larger, triple-zeta basis set and the models including the hybrid functionals have the smallest overall error, one third to a quarter that of Hartree-Fock. The models using the pure functionals have overall errors comparable to Hartree-Fock with the smaller basis sets and less than half that of Hartree-Fock for the larger basis set. However, when the standard deviations are examined, which are generally comparable in size to the average errors themselves, it is evident that these differences are not statistically significant. FOOF is another molecule which has historically proven to be a challenge for successive electronic structure methods. This molecule has an unusually long FO bond length, indicating a fairly weak interactions. FOOF is a case where MP2 theory performs rather poorly, and CCSD(T) or QCISD(T) with a large basis set is necessary for good agreement with experiment. Amos and coworkers [3] have published geometries for FOOF optimized with the LDA and BLYP functionals using the TZ2P and TZ2Pf basis sets. These optimizations were run with the 6-311+G(2d) and 6-311+G(2df) basis sets and obtained results in good agreement with Amos and coworkers. Table 2 Average absolute bond length differences with experiment (/~) and their standard deviations for several small molecules using traditional ab initio and DFT-based model chemistries. J

,,

6-31G(d,p) I ARe~Vl Std. Dev. HF MP2 SVWN5 BLYP B3LYP B3PW91 ,,

,

0.022 0.013 0.017 0.018 0.007 0.008

0.021 0.011 0.008 0.006 0.003 0.006 ,

6-311+G(2di2p) I AReXVl Std. Dev. 0.028 0.006 0.013 0.013 0.005 0.007

0.025 0.006 0.011 0.007 0.004 0.008

,,

Geometry optimizations were also performed with the B3LYP and B3PW91 hybrid functionals. Table 3 compares our results to Scuseria's CCSD(T) results [36] for the optimization of FOOF; Hartree-Fock and MP2 results with a very large basis set are also given for reference.

684 Table 3 ptimized structure of FOOF using various model chemistries (bond lengths in ; bond angles in degrees).

~

SVWN5 BLYP B3LYP B3PW91 CCSD(T)" Experiment"

6-311+G(2d) a Roo Rot AFoo DFooF 1.189 1.568 110.7 88.6 1.208 1.643 111.6 89.7 1.227 1.530 109.4 88.3 1.221 1.515 109.4 88.3 1.216 1.614 108.9 87.8 1.217 1.575 109.5 87.5

HF/6-311+G(3df) MP2/6-311+G(3df)

6-311 +G (2df) b Roo 1.187 1.205 1.226 1.221 1.218 1.217

RoF 1.559 1.634 1.522 1.506 1.589 1.575

AFoo 110.7 111.5 109.3 109.4 108.9d 109.5

DFooF 88.4 89.6 88.1 88.1 87.8d 87.5

1.296 1.143

1.352 1.672

106.3 111.4

85.1 89.3

" TZ2P for the CCSD(T) calculation. b TZ2Pf for the CCSD(T) calculation. ' Reference [36]. d These parameters were held fixed during the geometry optimization [36].

All the DFT-based models perform somewhat less accurately than CCSD(T), especially with respect to the O-F bond length, but their results are nevertheless substantially superior to Hartree-Fock and MP2 for this difficult case. Note that the addition of f functions to the basis set does not generally improve the predicted structure of this molecule for these model chemistries. The larger basis set results in a slight lengthening of the O-O bond and a shortening of the O-F bond, so it is of substantial help only when the latter bond length is overestimated at the TZ2P level. Some additional geometry optimization results are discussed as part of o u r consideration of accurate energies predictions using DFT-based model chemistries, the subject of the next section. 2.2 Thermochemical Results This section considers the performance of density functional theory for predicting thermochemical quantities: atomization energies, electron affinities, ionization potentials, and proton affinities. A standard series of reference calculations for which very accurate experimental values are known is used for this purpose [35, 37, 38]. This molecule set is commonly referred to as the G2 set, and it consists of 55 atomization energies, 38 ionization potentials, 25 electron affinities, and 7 proton affinities; it and its component subsets are frequently used to test new methods. This test set of molecules has certain inherent limitations. For example, all of the component molecules are very small systems containing only first and second row atoms. Only a limited number of bond types are represented (there are no ring systems). Nevertheless, it is useful for comparing the relative accuracies of

685 different model chemistries. The computations in the molecule set are difficult to model very accurately, and thus a model's performance on the set is also very likely to be a lower bound on its true accuracy. Foresman and coworkers have performed all of the calculations in the molecule set for over 30 different model chemistries [39]. Table 4 presents some of their overall results for various model chemistries. (Seminario performed a related study comparing LDA atomization energies to those computed by G2 theory [18].) Table 4 Overall results on the 125 calculations in the G2 molecule set for various model chemistries (all values are in kcal/mol)." ,,

B3LYP/6-311+G(2d,p) / / B3LYP / 6-311 +G (2d,p) BLYP/6-311+G(2d,p) / / BLYP/6-311+G(2d,p) MP2/6-311+G(2d,p) / / MP2 / 6-311 +G (2d,p) SVWN5 / 6-311 +G(2d,p) / / SVWN5/6-311 +G (2d,p) SVWN3 / 6-311 +G(2d,p) / / SVWN3 / 6-311 +G (2d,p) BLYP/6-31+G(d,p) / / BLYP/6-31+G(d,p) B3LYP/6-31+G(d,p) / / B3LYP/6-31+G(d,p) B3LYP/6-31G(d) / / B3LYP/6-31G(d) MP2/6-31+G(d,p) / / MP2/6-31+G(d,p) HF / 6-31+G(d,p) / / HF/6-31+G(d,p) HF/6-31G(d) / / HF/6-31G(d)

Mean Absolute Standard Deviation Deviation 3.1 3.0

Maximum Errors Positive 13.6

Negative -19.7

3.9

3.2

14.3

-15.9

8.9

7.8

18.3

-39.2

18.1

19.8

81.0

-10.1

24.9

19.2

89.3

-10.4

3.9 ......

3.2

15.2

-i5.2

3.9

4.2

17.6

-33.8

7.9

9.5

12.2

-54.2

11.4

8.1

15.6

-44.0

46.7

40.6

10.1

-179.8

51.0

41.2

11.5

-184.2 ,,

"

Reference [39].

In each case, the geometry of each system w a s first optimized at the given level, and vibrational frequencies were computed in order to predict the zeropoint energy (ZPE), which was scaled using the scaling factors of Wong [25] and added to the predicted total energy. Stability calculations were performed for

686 each system for every model chemistry to ensure that the lowest energy wavefunction of the proper electronic state had been located. The columns in the table give the mean absolute deviation (MAD) from experiment (the average difference between the calculated and experimental values ignoring sign), the standard deviation of the absolute errors, and the largest positive and negative errors for any calculation in the set. The MAD gives a measure of the average magnitude of error with respect to experiment, the standard deviation gives a sense of how widely the model's performance varied across the various calculations, and the maximum errors indicate each model's worst-case accuracy for this molecule set. The five models in the top portion of the table all use the large 6-311+G(2d,p) basis set, and accordingly may be used to compare the overall results for these methods somewhat independently of basis set effects. B3LYP performs best of all of them, with a MAD of only 3.1 kcal/mol. It's maximum error is around 20 kcal/mol, however, indicating that it still requires calibration for new classes of systems. The performance of the BLYP functional is worse than that of the hybrid functional with a MAD of 5.0 kcal/mol, and its largest error also greater. The results for MP2 theory at the same basis set are significantly inferior to those for B3LYP and BLYP. The two local DFT functionals perform very poorly on this set of calculations, having errors at least twice as large as MP2. The models in the lower section of Table 4 use smaller basis sets, and these results are provided for comparison purposes. The B3LYP method using the 631+G(d,p) basis set performs almost as well as B3LYP with the larger basis set, although the maximum error increases substantially. BLYP using the 6-31+G(d,p) basis set has a lower MAD than BLYP with the larger basis set although the other statistics are essentially the same for the two BLYP-based model chemistries. B3LYP using the modest 6-31G(d) basis set is still more accurate than MP2 with the larger basis set. The two Hartree-Fock model chemistries indicate the unsuitability of Hartree-Fock theory for calculations at this desired level of accuracy. The 6-311+G(2d,p) basis set is rather large, and will be out of reach for geometry optimizations and frequency calculations for many large systems of chemical interest. For such cases, it is often desirable to perform the geometry optimization and frequency calculation at a lower level of theory than the final energy calculation. Overall G2 set results for several models based u p o n B3LYP/6-31G(d), HF/3-21G(d) and AM1 geometries and frequencies are given in Table 5, along with the result for Gaussian-2 theory for comparison purposes [38].

687 Table 5 Overall results on the 125 calculations in the G2 molecule set for various model chemistries at several specified geometry optimization and frequency calculation levels (all values are in kcal/mol), a

G2 Theory b B3LYP/6-311+G(3df,2df,2p) / / B3LYP/6-31G(d) B3LYP/6-311+G(2d,p) / / B3LYP/6-31G(d) B3LYP/6-31+G(d,p) / / B3LYP/6-31G(d) B3LYP/6-31G(d) / / B3LYP/6-31G(d) MP2/6-31 l+G(2d,p) / / B3LYP/6-31G(d) B3LYP/6-311+G(2d,p) / / HF/3-21G(d) B3LYP/6-31G(d) / / HF/3-21G(d) B3LYP/6-31G(d)//AM1 c AM1 / / A M 1 ~ HF/6-31G+(d,p)//AM1 ~ HF/6-31G(d)//AM1 ~

Maximum Errors Mean Absolute Standard Deviation Deviation Positive Negative 5.1 -3.5 1.2 0.9 12.5 -9.3 2.7 2.6 3.2

3.0

13.6

-20.1

4.0

4.2

17.6

-33.9

7.9

9.5

12.2

-54.2

8.9

7.8

29.7

-39.2

3.2

3.0

13.8

-21.2

8.0

9.4

9.4

-54.2

10.5 18.8 49.4 54.2

11.3 16.9 43.1 43.1

14.7 47.8 8.0 8.6

-54.2 -95.5 -206.1 -207.2

Reference [39]. b Derived from data in Reference [38]. c AM1 zero-point energies were not scaled. The calculations involving Na, Li and Mg were excluded from models using AM1. The non-zero energy for H § predicted by AM1 was computed and taken into account in the calculation of proton affinities. "

Using B3LYP/6-311+G(3df,2df,2p) energies computed at the B3LYP/6-31G(d) geometry and corrected by the (scaled) ZPE at the optimized geometry produces very good overall results across the computations in the G2 set. In fact, this model chemistry had the lowest MAD of all of those considered in this study. (The basis set notation specifies different sets of polarization functions for second-row and third-row and higher heavy atoms.) This basis set is similar to one recommended by Bauschlicher and Partridge [40]. The B3LYP/6-311+G(2d,p) model chemistry at the B3LYP/6-31G(d) geometry also produces good results, although the maximum error is almost twice that

688 using the larger basis set. Using smaller basis sets for the final energy calculation correspondingly increases the mean absolute deviation and maximum errors. The MP2/6-311+G(2d,p) model performs less accurately than the hybrid functionals for all of the models using the B3LYP/6-31G(d) geometries. For systems where B3LYP/6-31G(d) is too expensive for a geometry optimization and frequency calculation, Hartree-Fock geometries and frequencies using the 6-31G(d) or 3-21G(d) basis sets may be a viable alternative. Two examples of the latter are given in the middle section of the table. The final section of the table gives results for three models using AM1 geometries and frequencies. Such models are not very accurate, but they may be all that is practical for some very large molecules. For such systems, using B3LYP for a final single point energy calculation substantially improves the overall accuracy of the thermochemistry. These data suggest some interesting observations about methods for geometry optimizations and frequency calculations when thermochemical predictions is the final calculation goal. Table 6 summarizes the relevant results. Table 6 B3LYP-based model chemistries with the best overall results on the 125 calculations in the G2 molecule set (all values are in kcal/mol)."

B3LYP/6-311+G(3df,2df,2p) / / B3LYP/6-31G(d) B3LYP/6-311+G(2d,p) / / B3LYP / 6-311 +G (2d,p) B3LYP/6-311+G(2d,p) / / B3LYP/6-31G(d) B3LYP/6-311+G(2d,p) / / HF/3-21G(d) B3LYP/6-31+G(d,p) / / B3LYP/6-31+G(d,p) B3LYP/6-31+G(d,p) / / B3LYP/6-31G(d) B3LYP/6-31G(d) / / B3LYP/6-31G(d) ,

"

Maximum Errors Mean Absolute Standard Deviation Deviation Positive Negative 2.7 2.6 12.5 -9.3 3.1

3.0

13.6

-19.7

3.2

3.0

13.6

-20.1

3.2

3.0

13.8

-21.2

3.9

4.2

17.6

-33.8

4.0

4.2

17.6

-33.9

7.9

9.5

12.2

-54.2

.,

Reference [39].

Note that the overall results for B3LYP/6-311+G(2d,p) final energies are not very sensitive to the basis set used for the optimized geometry. Similarly, the two models using B3LYP/6-31+G(d,p) total energies are essentially unchanged by the

689 addition of hydrogen polarization and heavy-atom diffuse functions in the geometry optimization and frequency calculation. However, these additional functions do make a significant difference in the final energy calculation (compare the final two models). Table 7 lists the results for various model chemistries on the atomization energy calculations in the G2 molecule set. Table 7 Results on the 55 atomization energy calculations in the G2 molecule set for various model chemistries (all values are in kcal/mol)."

Mean Absolute B3LYP/6-311+G(3df,2df,2p) / / B3LYP/6-31G(d) B3LYP/6-311+G(2d,p) / / B3LYP / 6-311G + (2d,p) B3LYP/6-311+G(2d,p) / / B3LYP/6-31G(d) B3LYP/6-311+G(2d,p) / / HF/3-21G(d) BLYP/6-31+G(d,p) / / BLYP/6-31+G(d,p) BLYP/6-311+G(2d,p) / / BLYP / 6-311 +G (2d,p) B3LYP/6-31+G(d,p) / / B3LYP/6-31+G(d,p) B3LYP/6-31+G(d,p) / / B3LYP/6-31G(d) B3LYP/6-31G(d) / / B3LYP/6-31G(d) B3LYP/6-31G(d)// HF/3-21G(d) B3LYP/6-31G(d)//AM1 b SVWN5/6-311+G(2d,p) / / SVWN5 / 6-311 +G (2d,p) SVWN3/6-311+G(2d,p) / / SVWN3/6-311+G(2d,p)

Maximum Errors

Standard

Deviation Deviation 2.3 2.1

Positive 8.1

Negative -9.3

3.0

3.0

8.0

-19.7

3.2

3.1

8.0

-20.1

3.3

3.1

8.1

-21.2

4.0

3.4

13.5

-15.2

4.5

3.5

14.3

-15.9

4.4

5.0

8.1

-33.8

4.6

5.1

8.3

-33.9

5.4

5.0

8.1

-31.5

5.7

5.4

8.3

-35.0

8.6 34.0

8.2 20.8

8.0 81.0

-44.2 -0.6

37.6

22.8

89.3

-0.4

,,

" b

Reference [39]. AM1 zero-point energies were not scaled. The calculations involving Na, Li and Mg were excluded from the model using AM1.

690 The relative performance of the various model chemistries is quite similar to that of the overall results. The major exception occurs with the two LSDA functionals which perform much more poorly for atomization energies than they do for other calculation types (and hence on the overall set). As was true for the overall results, the 6-311+G(3df,2df,2p)//B3LYP/6-31G(d)is the most accurate. Table 8 similarly lists the results for the ionization potential calculations. The ordering of the various model chemistries is very similar to that for the atomization energies. Table 8 Results on the 38 ionization potential calculations in the G2 molecule set for various model chemistries (all values are in kcal/mol), a

B3LYP/6-311+G(3df,2df,2p) / / B3LYP/6-31G(d) B3LYP/6-311+G(2d,p) / / B3LYP/6-31G(d) B3LYP/6-311+G(2d,p) / / HF/3-21G(d) B3LYP/6-311+G(2d,p) / / B3LYP / 6-311 + G (2d,p) B3LYP/6-31+G(d,p) / / B3LYP/6-31G(d) B3LYP/6-31+G(d,p) / / B3LYP/6-31+G(d,p) B3LYP/6-31G(d) / / B3LYP/6-31G(d) B3LYP/6-31G(d) //HF/3-21G(d) BLYP/6-311+G(d2,p) / / BLYP/6-311+G(1d,p) BLYP/6-31+G(d,p) / / BLYP/6-31+G(d,p) SVWN5/6-311+G(2d,p) / / SVWN 5 / 6-311 +G (2d,p) B3LYP/6-31G(d)//AM1 b SVWN3/6-311+G(2d,p) / / SVWN3/6-311+G(2d,p)

Mean Absolute Deviation 3.6

Standard Deviation 3.2

Maximum

Errors

Positive 12.5

Negative -4.8

3.6

3.2

13.2

-4.8

3.6

3.1

12.0

-5.6

3.7

3.2

13.1

-5.2

4.0

3.3

15.7

-6.0

4.1

3.3

15.7

-6.0

4.1

3.5

12.2

-13.8

4.2 4.4

3.3 3.7

8.5 12.2

-13.8 -9.0

4.7

3.1

12.4

-9.4

4.9

3.4

13.5

-6.9

5.2 14.9

6.0 5.1

14.7 26.5

-31.2 -10.4

" Reference [39].

b AM1 zero-point energies were not scaled. The calculations involving Na, Li and Mg were excluded from the model using AM1.

691 Once again the model chemistries based on local functionals perform significantly less well than the gradient-corrected functionals. The hybrid functionals continue to perform better than the pure ones. Table 9 shows similar trends for electron affinities and proton affinities, although BLYP performs slightly better than B3LYP for electron affinities. For proton affinities, BLYP and B3LYP perform equally well when using the 631+G(d,p) basis set, but the hybrid functional has a significant advantage when the larger 6-311+G(2d,p) basis set is used. 2.3 Using DFT Geometries within Compound Model Chemistries Several groups are investigating incorporating DFT methods into composite methods for computing high accuracy total energies such as Gaussian-2 theory [38]. For example, Bauschlicher and Partridge have examined a modification of the G2(MP2) method which uses B3LYP structures and vibrational frequencies [41]. Morokuma and coworkers have proposed a modified G2 theory incorporating coupled cluster energies and DFT geometries and frequencies [42]. Bauschlicher and Partridge denote their method G2(B3LYP/MP2/CC). Figure 1 summarizes the steps involved in the original G2(MP2) method [43], which simplifies G2 by computing the basis-set-extension corrections at the MP2 level rather than the MP4 level (in the interests of reducing computation time and disk storage requirements), as well as the variations proposed by Bauschlicher and Partridge. Figure 1 Summary of the G2(MP2) methods and the G2(B3LYP/MP2/CC) modifications to the procedure. . . . . . . . . . . . . . . Step G2(MP2) Method 1 HF/6-31G(d) Opt./Freq. for ZPE 2 MP2/6-31G(d) Opt. for final structure 3 QCISD(T)/6-311G(d,p) energy 4 MP2/6-311+G(3df,2p) energy 5 Higher level correction (HLC) computed from # of valence (x and ~ electrons ECaa~2) = E~SDm + (EMP2(step4)_EMP2(Step3)) + 0.8929(ZPE HE) + HLC ....

G2(B3LYP/MP2/CC) Variations B3LYP/6-31G(d) Opt./Freq. for ZPE Use B3LYP/6-31G(d) geometry CCSD(T)/6-311G(d,p) energy Same Same

EG2(B3LYP/MP2/Cc)= ECCSDC r) + (EMP2~S,~p4)_EM~,~pS)) + 0.98(ZPE B3LYP)+ HLC .

.

.

.

.

.

.

.

.

.

692 Table 9 Results on the 25 electron affinity calculations (top section) and 7 proton affinity calculations in the G2 molecule set for various model chemistries (all values are in kcal/mol)." ,,,

,

BLYP/6-311+G(2d,p) / / BLYP/6-311+G(2d,p) B3LYP/6-311+G(2d,p) / / B3LYP/6-311 +G (2d,p) B3LYP/6-311+G(3df,2df,2p) / / B3LYP/6-31G(d) BLYP/6-31+G(d,p) / / BLYP/6-31+G(d,p) B3LYP/6-311+G(2d,p) / / B3LYP/6-31G(d) B3LYP/6-311+G(2d,p) / / HF/3-21G(d) B3LYP/ 6-31+G(d,p) / / B3LYP/6-31+G(d,p) SVWN5/6-311+G(2d,p) / / SVWN5/6-311+G(2d,p) SVWN3/6-311+G(2d,p) / / SVWN3/6-311+G(2d,p) B3LYP/6-31i+G(3df,2df,2p) / / B3LYP/6-31G(d) B3LYP/6-311+G(2d,p) / / B3LYP/6-311 +G (2d,p) B3LYP/6-311+G(2d,p) / / B3LYP/6-31G(d) B3LYP/6-311+G(2d,p) / / HF/3-21G(d) BLYP/6-311+G(2d,p) / / BLYP/ 6-311 +G (2d,p) BLYP/6-31+G(d,p) / / BLYP/6-31+G(d,p) B3LYP/6-31+G(d,p) / / B3LYP/6-31+G(d,p) SVWN3/6-311 +G(2d,p) / / SVWN3/6-311+G(2d,p) SVWN5/6-311+G(2d,p) / / SVWN5/6-311+G(2d,p) " R e f e r e n c e [39].

,,

MAD 2.5

,

,,,,

,,,,,,

,

,

Standard Maximum Errors Deviation Positive Negative 2.5 11.4 -4.5

2.6

2.9

13.6

-1.7

2.6

2.5

10.5

-7.0

2.8

3.1

15.2

-4.3

2.9

3.0

13.6

-7.0

2.9

3.0

13.8

-2.6

3.0

3.7

17.6

-1.9

6.6

3.3

15.1

0.0

17.3

3.5

26.9

0.0

" 1'.0

0.6

0.6

-2.2

1.2

0.7

1.2

-2.4

1.3

0.6

0.7

-2.3

1.4

0.9

1.0

-2.9

1.7

1.5

0.7

-4.3

1.9

0.9

2.7

-2.7

2.0

1.0

3.0

-3.5

6.0

2.0

0.0

-9.4

6.5

2.0

0.0

--10.1

693 Bauschlicher and Partridge's variations are designed to better address systems which require an electron correlation treatment to obtain accurate vibrational frequencies and zero-point energy and to reduce the computational requirements of the overall method by replacing the MP2 optimization with the less expensive but comparably accurate B3LYP calculation. The CCSD(T) coupled cluster theory method also replaces QCISD(T). These researchers tested their approach on the 55 atomization energies in the G2 molecule set. Table 10 summarizes their results for the B3LYP/6311+G(3df,2p) model chemistry and for the G2(MP2/B3LYP/CC) method, as well as for the Complete Basis Set method known as CBS-Q, another compound model chemistry designed to produce highly accurate thermochemistry results [44, 45]. Table 10 Results on the 55 atomization energy calculations in the G2 molecule set for variations of the G2(MP2) model chemistry (all values are in kcal/mol). Mean Maximum Absolute Standard Absolute Deviation Deviation Error CBS-Q" 0.8 0.6b 2.3 G2(MP2)" 1.3 0.9b 4.2 G2(B3LYP/MP2)C 1.3 1.0 3.7 G2(B3LYP / MP2 / CC) ~ 1.3 0.9 3.1 B3LYP / 6-311+G(3df,2p) * 2.2 1.9 8.1 Reference [45]. b Standard deviation calculated from the reported RMS value. c Reference [41]; ZPEs scaled by 0.989. "

Bauschlicher and Partridge considered their two changes to the G2(MP2) procedure separately. Using the B3LYP/6-31G(d) geometry and ZPE--denoted G2(B3LYP/MP2)--reduced the maximum absolute error slightly. Replacing QCISD(T) with CCSD(T) further reduced this maximum error, although neither substitution had any significant effect on the MAD or its standard deviation for this series of calculations. Thus, this approach does reduce the computational intensity of G2(MP2) without a decrease in accuracy. None of the G2 variations is as accurate as CBS-Q, and CBS-Q is significantly less computationally expensive than those based on G2 theory. The atomization energies computed via the B3LYP/6-311+G(3df,2p) model chemistry are not as accurate and the maximum error is about twice as large as those of the compound models incorporating high-level and exl~ensive traditional electronic structure methods (the cost of which scales as N'). This result is also in line with those for other B3LYP-based model chemistries examined in the previous section of this chapter.

694 3. MOLECULAR AND VIBRATIONAL PROPERTIES Molecular properties may also be predicted with DFT-based methods, including both ones derived from the electron density and the vibrational frequencies of the molecule and its related properties. 3.1 M o l e c u l a r P r o p e r t i e s

De Proft and coworkers have computed atomic populations according to several schemes using several different ab initio and DFT methods [46]. Their preliminary examination at the QCISD/cc-pVDZ level indicated that charges computed with the atomic polar tensor (APT) method [47] were the most sensitive to the effects of electron correlation. Accordingly, they computed charges with this scheme for a variety of methods. Their results are summarized in Table 11. This table compare the results for several theoretical methods (all using the ccpVDZ basis set) to the QCISD/cc-pVDZ charges. All of the methods involving electron correlation achieve reasonable agreement with the high-level results, and the differences between them are much smaller than the difference between the Hartree-Fock and QCISD results. For this set of molecules, the smallest overall difference from the QCISD results was exhibited by the hybrid functionals. Table 11 Mean absolute deviation of APT charges from the QCISD results for 15 first-row molecules using a variety of methods and the cc-pVDZ basis set (units of atomic charge). HE 0.104

MP2 0.046

~vIAD(AQCISD)". LDA BLYP 0.034 0.032

. . . . B3LYP B3PW91 0.013 0.019

" Reference [46].

Table 12 summarizes De Proft and coworkers' comparison of the dipole moment predicted by B3LYP and QCISD, using the 6-31G(d) and cc-pVDZ basis sets, for a set of 12 compounds containing first-row atoms to the observed values for this quantity. With the smaller basis set, B3LYP and QCISD perform comparably for this set of compounds. Agreement with experiment improves with the larger basis set for both methods, but the improvement is greater for QCISD, and the B3LYP average absolute deviation is twice as large as that for QCISD.

695 Table 12 Mean absolute deviation from experiment of the predicted dipole moment and its standard deviation for 12 first-row compounds using the B3LYP and QCISD methods and the cc-pVDZ basis set (values in debye).

B3LYP/6-31G(d) QCISD/6-31G(d) B3LYP / cc-pVDZ QCISD / cc-pVDZ b

MAD(AExp) ~ 0.197 0.196 0.116 0.059

Standard Deviation" 0.129 0.166 0.112 0.049

" Computed from data in Reference [46]. b Only 10 of 12 QCISD values are reported at this basis set. Wiberg and coworkers have investigated spin polarization in heterosubstituted allyl radicals using a range of ab initio and DFT methods [23]. Goodexperimental values exist only for acetyl radical; Wiberg and coworkers' results for this compound as well as some smaller calculations of our own are given in Table 13. Table 13 Predicted spin polarizations for acetyl radical using several different model chemistries. O

C 1 (away from C2 substituent) (near substituent) 6-31G(d) 6-311G(d,p)" 6-31G(d) 6-311G(d,p) a 6-31G(d) 6-311G(d,p) a 0.651 1.091 -0.615 HF 0.023 0.039 0.994 0.980 0.127 0.117 MP2 0.386 0.292 0.864 0.911 -0.154 -0.122 B3LYP 0.334 0.275 0.964 1.002 -0.187 -0.168 QCISD Experiment" 0.2 0.7 " Reference [23].

Based on these calculations, spin densities appear to be another property for which the hybrid B3LYP functional produces excellent agreement with the highlevel ab initio method, in this case QCISD. The spin densities in acetyl radical is a problem where MP2 theory performs poorly, predicting a value for the oxygen atom that is much too small and a value for the second carbon atom that has the wrong sign. Spin contamination is the major factor reducing the accuracy of the Hartree-Fock and MP2 results.

696 3.2 Vibrational Frequencies and Zero-Point Energies Vibrational frequencies computed with ab initio methods are known to suffer from systematic errors. Accordingly, computed frequencies are typically scaled by standard factors before they are compared with experimental observations [48, 49]. Frequencies used to compute zero-point energies and other thermal corrections to the total energy are also scaled, and the optimal scaling factors for each of these uses differ from one another. Wong [25] and Scott and Radom [17] have investigated the performance of several DFT models for predicting vibrational frequencies and zero-point energies, and they have independently computed optimal scaling factors for each model chemistry. Where their investigations overlapped, there is very good agreement between the two studies. Their findings are also in agreement with the scale factors for zero-point energies computed by Bauschlicher and Partridge [41]. Table 14 summarizes the computed scaling factors for a variety of ab initio and DFT model chemistries based on the 6-31G(d) basis set, as well as the errors that remain in the computed frequencies and zero-point energies after scaling.

Table 14 Optimal scaling factors and overall RMS errors after scaling for fundamental frequencies and zero-point energies for various model chemistries (RMS errors in cm 1 for frequencies and kJ/mol for ZPEs). Zero-point Energies Fundamental Frequencies Scale RMS Error Scale RMS Error Factor Factor HF/6-31G(d) b 0.8953 50 0.9182 0.46 0.9833 43 1.0079 0.57 SVWN3/6-31G(d)" B-LYP/6-31G(d) d 0.9940 45 1.0057 0.44 0.9613 34 0.9826 0.30 B3LYP/6-31G(d) d 0.9559 34 0.9760 0.28 B3P86/6-31G(d) d 0.9573 34 0.9774 0.28 B3PW91/6-31G(d) ~ MP2(FulI)/6-31G(d) c 0.9427 61 0.9643 0.50 0.9434 63 0.9676 0.50 MP2 (FC) / 6-31G(d) b QCISD/6-31G(d) b 0.9537 37 0.9805 0.38 ,, ,,

,

..,

.

.

.

.

.

" Reference [25]. b Reference [17]. Reference [50]. d Frequency data from Reference [25] and ZPE data from Reference [17].

697 Both studies used a test set of 1066 frequencies of 122 molecules (known as the Z1 set) to determine frequency scaling factors. Table 14 reports the ZPE scale factors of Scott and Radom, who used a set of 89 frequencies from 39 molecules to compute them. All of the computed scaling factors are near unity, indicating a correlation of modest remaining errors in the harmonic force field with the (neglected) anharmonic contributions. The DFT-based model chemistries perform very well. Those using pure functionals produce overall RMS errors of about the same magnitude as traditional ab initio methods such as Hartree-Fock and MP2. The hybrid functionals produce significantly smaller overall RMS errors, with values on the same order as QCISD. These results demonstrate that DFT methods, and especially those based on hybrid functionals, are likely to be reliable for predicting vibrational frequencies and properties derived from them. Indeed, the B3LYP and B3PW91 models produce the best results of all of those considered. 3.3 Vibrational Circular Dichroism

Vibrational circular dichroism (VCD)is exhibited by chiral molecules [51], and can accordingly be used as a direct method for determining the absolute configuration of a chiral molecule. To do so in practice requires the accurate prediction of VCD intensities; these are determined by vibrational rotational strengths, which are themselves in turn determined by vibrational electric and magnetic dipole transition moments. Computation of VCD spectra within the harmonic approximation thus requires the calculation of the Harmonic Force Field (HFF), the Atomic Polar Tensors (APTs), and the Atomic Axial Tensors (AATs) [52, 53]. The development of a method for predicting VCD spectra which employs density functional theory to compute each of the required components of the vibrational rotational strengths has recently been reported [10, 54, 55]. This implementation includes a new, efficient DFT methodology for the calculation of AATs using direct methods and gauge-including atomic orbital (GIAO) basis sets. Systematic calibration studies on a large set of molecules are still underway. However, preliminary results suggest that hybrid functionals produce excellent agreement with experimental VCD spectra. For example, Table 15 compares the rotational strengths predicted by Hartree-Fock theory and three DFT-based methods to available observed values for 6,8-dioxabicyclo[3.2.1]octane.

698 Table 15 Calculated and experimental rotational strengths for the (+)-(1R,5S) enantiomer of 6,8-dioxabicyclo[3.2.1]octane for various methods using the TZ2P basis set (values in 10-~ e s u 2 cm2). iCundamental 38 33

HF" 6.0 -16.5

LSDA" 10.6 0.5

BLYP" 9.3 -9.0

B3LYP" 8.5 -11.9

Exp? 19.4 -18.8

3z

o.s

-5.9

6.s

8.s

24.9

31 30 29 28 26 23 21 20 19 18 17 16 15 14 13 12

13.5 7.4 -8.0 2.8 -13.0 -13.1 0.2 73.7 -44.4 19.1 -37.4 -2.7 51.2 -12.0 12.8 12.0

-7.4 32.5 -13.6 12.0 -22.7 4.8 -12.7 12.6 -19.4 -13.2 -26.9 38.7 28.1 21.2 -15.6 3.1

4.9 12.3 2.7 -1.7 -9.9 5.5 36.6 6.8 7.5 -82.1 -63.3 9.0 57.2 66.1 -9.2 23.1

-13.8 25.7 -0.7 -1.4 -11.9 7.5 6.1 12.6 6.8 -44.1 -45.4 -2.2 74.8 47.8 -44.0 30.6

-29.8 52.3 -15.2 11.6 -16.4 13.6 -9.6 19.3 20.3 -62.3 -80.0 6.4 88.0 63.3 -74.0 32.5

" Reference [55]. b Reference [56]. The B3LYP spectra best reproduces the experimental spectra for this molecule. In fact, this calculation yields a nearly unambiguous assignment of fundamentals 10-38. The BLYP spectra is significantly less accurate than the B3LYP results, and the differences are even greater for the LSDA spectra. For this problem, the functionals can be rated LSDA < BLYP < B3LYP in terms of relative accuracy. Finally, the poor Hartree-Fock results illustrate the importance of electron correlation for this property. 3.4 N M R Chemical Shifts Electronic structure methods have been applied to the prediction of NMR properties for several years [57, 58], and recently DFT methods using GIAOs have also proven effective for predicting nuclear magnetic resonance shielding tensors and calculating NMR chemical shifts [7, 59, 60]. Various methods for c o m p u t i n g NMR chemical shifts are compared for a set of 16 molecules in Table 16.

699 The MP2 results are the most accurate. However, the gradient-corrected DFTbased models represent a significant improvement over the Hartree-Fock results. Interestingly, the raw differences between the predicted values and experiment indicate that the DFT-based methods generally predict shifts which are too deshielded, while Hartree-Fock theory makes errors in both directions. These shifts were computed using a very large basis set. For general use, we recommend the B3LYP/6-311+G(2d,p) model for computing NMR properties of chemical systems, using a geometry optimized at the B3LYP/6-31G(d) level. For systems where it is computationally practical, this level represents a reasonable trade-off between accuracy and cost, since B3LYP is significantly more costeffective than MP2. This model includes the effects of electron correlation in the calculation and provides quantitatively good agreement with experiment. Table 17 summarizes results computed with these models for a set of 32 ~C and 5 ~N chemical shifts from 22 molecules. Even the Hartree-Fock NMR predictions are adequate for assignment of most experimental spectra, and the B3LYP results are quite accurate.

Table 16 Mean absolute deviation, standard deviation of the absolute error, and maximum absolute deviation of ~3C,~SN and ~70 chemical shifts a from experiment for a set of 16 molecules b using a variety of methods and the QZ2P basis set (values in ppm). Overall (25 shifts)

~C only (17 shifts)

~N only (5 shifts)

MAD Std. Dev. Max. Error MAD Std. Dev. Max. Error MAD Std. Dev. Max. Error

I-t.~ 19.0 17.9 65.2

LSDA c 20.2 11.6 44.8

BI~YP~ 12.1 8.0 30.2

B3LYP~ B3PW91 r 14.1 13.2 10.6 10.6 39.6 38.8

MP2 'a''' 4.0 5.2 23.9

8.9 7.4 30.8

14.7 7.3 26.3

8.1 4.7 18.6

8.0 5.5 17.3

7.1 5.2 17.3

1.6 1.8 6.5

47.8 12.1 65.2

32.8 10.8 44.8

20.6 6.8 30.2

28.4 6.5 39.6

26.6 7.0 38.8

11.8 6.4 23.9

~3C shifts relative to CH 4, ~SN shifts relative to NH3. and ~70 shifts relative to H20. b Geometries optimized with MP2/TZ2P. c Reference [71. a MP2 statistics based on data from Reference [58]. "

700 Table 17 Mean absolute deviation, standard deviation of the absolute errors, and maximum absolute deviation of ~3C and ~SN chemical shifts" from experiment for a set of 22 molecules b using the model chemistries recommended in the text (values in ppm). -

,

BC only (32 shifts)

MAD Std. Dev. Max. Error

~N only (5 shifts)

MAD Std. Dev. Max. Error

,,

HF/6-31G(d) ~ 9.3 5.9 22.4

l~3LYP/6_311 +G(2d,p),

12.5 7.2 19.7

16.9 3.0 21.6

3.6 2.2 8.4

" ~C shifts relative to TMS, and 15N shifts relative to NH 3. b Geometries optimized with B3LYP/6-31G(d). " Reference [7].

Unfortunately, many interesting molecules are substantially too large for their NMR properties to be practically studied with these DFT-based methods. For such systems, Hartree-Fock theory, using the modest 6-31G(d) basis set, may provide a reasonable alternative approach to computing ~3C chemical shifts (HF/6-31G(d) BC shifts are only slightly less accurate than ones computed with the HF/6-311+G(2d,p) model chemistry [7]). Taxol provides an appropriate test case (this molecule's structure is illustrated in Figure 2). The structure of this molecule was optimized using the HF/STO-3G model, and its NMR properties were computed at the HF/6-31G(d) level. The RMS error of the computed values for the 41 13C chemical shifts for which experimental data exists was 6.4 ppm, and the mean absolute deviation was 5.0 ppm. The maximum error for carbon shifts is 18.7 ppm. These results appear to be sufficiently accurate to be of aid in experimental peak assignments. The DFT functionals used here to predict NMR properties do not currently include a specific magnetic-dependent term. It is possible that the inclusion of such a term would improve upon these results.

701 Figure 2 Structure of the taxol molecule. O

H H3C~

~

f

H

i

0

0

OH

o H~c

,'#

4. IMPORTANCE AND CALIBRATION OF INTEGRATION GRIDS In contrast to conventional ab initio methods, which rely on analytic algorithms for their principal calculation types, DFT calculations currently proceed by employing three-dimensional numerical quadrature techniques, which involves numerical integration over a three dimensional grid of points in space [61-67]. Such grids are defined as a number of radial shells around each atom, each of which contains a set number of integration points. There is a strong correlation between the number and choice of grid points and both the accuracy and cost of the quantum chemical calculation. Grids are usually denoted as a parenthesized pair of integers. For example, the designation (75,302) refers to a grid having 75 radial shells around each atom each of which contains 302 points. Pruning a uniform gridma grid containing the same number of angular points at each radial distance--to eliminate some of the points in specific regions of space provides a means of reducing the computational cost for a calculation with negligible loss in accuracy. Pruned grids are reduced from the full, uniform grid

702 form so that fewer points are used on the shells near the core, where the density is nearly spherically symmetric, and far from the nucleus, where the density is small. Thus, pruned grids are densest in the region of the atom where properties are changing most rapidly. Available pruned grids differ in three main ways: their computational accuracy, the degree to which they exhibit rotational invariance and their range of applicability by element. In this section, two widely-used pruned grids are compared [68]: the SG1 grid suggested by Gill and coworkers [69], which is a pruned (50,194) grid, and the (75,302)p grid, a pruned version of (75,302). 4.1

Accuracy

The choice of an integration grid can influence the accuracy of the final DFT calculation. Table 18 gives the overall accuracy differences between several commonly-used grids with respect to the energies computed with the very large (96,32,64) grid of Murray and coworkers [70] (a grid of 96 radial shells, each composed of a 32x64 uniform angular grid). Table 18 Mean absolute deviation, standard deviation of the absolute error, and maximum absolute errors for the differences in predicted total energy with respect to the (96,32,64) grid for a set of 36 molecules a containing first and second-row atoms, computed with the BLYP/6-31G(d) model chemistry using various integration grids (energy differences in kcal/mol), b

MAD Standard Deviation Maximum Error

(99,302) 0.009 0.017 0.079

(75,302) (75,302)p 0.009 0.008 0.017 0.016 0.078

0.078

SG1 (50,194) [(50,194)p] 0.043 0.044 0.059 0.062 0.265

0.280

" Tight-convergence geometry optimizations were performed at the MP2/6-31G(d) level for systems with up to one heavy atom and at the HF/6-31G(d) level for those with two or more heavy atoms. b Reference [68]. "

Although all of these energy differences are relatively small for these small test systems, there is a significant difference in overall accuracy between the SG1 and (75,302)p grids. Larger differences are to be expected with large molecules. SG1 has a MAD about 5 times as large as (75,302)p, and their maximum errors also differ by a factor of 4. Both grids perform just as well as the larger uniform grids from which they were pruned.

703 4.2 R o t a t i o n a l

Invariance

All finite grids suffer from some degree of rotational invariance: changing the orientation of the molecular system under investigation can alter the predicted energy. (75,302)p exhibits relatively little rotational invariance while for SG1 this effect can be substantial. Table 19 gives the results for an example case which illustrates this difference, A14P4. Table 19 Difference in computed total energy between two orientations of A14P4 for two integration grids, computed at the BLYP/6-31G(d) level (values in kcal/mol)." _

. ,

,,

AEorientation

SG1 [(50,192)p] (75,302)p .

"

.

.

0.354 0.079

.

.

.

.

.

.

Reference [68].

The change in energy arising solely from the differing orientation of the molecule is over 4 times as large for SG1 as for the (75,302)p grid. Once again, the error produced by SG1 is small but significant. 4.3 E x t e n s i b i l i t y

Ideally, the same integration grid should be used for all DFT calculations undertaken, regardless of the atom types within the molecular system under investigation. Initially, both SG1 and (75,302)p were defined for the first two rows of the periodic table. (75,302)p has been extended to the third row. However, SG 1 cannot be so extended in any straightforward way since its parent grid, the uniform (50,194)grid, performs very poorly on third-row atoms. Table 20 illustrates this effect. Table 20 Difference in predicted total energy with respect to the (96,32,64) grid for two third-row atoms, using the BLYP/ 6-31G(d) model chemistry (values in kcal/mol)." .

.

Mn Zn

.

.

.

(99,302) 0.000 0.000 ,

"

.

.

.

(75,302ip 0.008 0.009 .

Reference [68].

.

.

.

.

.

.

.

.

.

.

.

i50,i94) -0.727 -0.660

704 The energies predicted by the (75,302) grid differ from those of the large grid only by microhartrees. In contrast, those predicted by the full, uniform (50,194) grid differ by over 0.6 kcal/mol, indicating that this grid is unsuitable for compounds involving third-row and higher atoms even before pruning.

5. CONCLUSION Model chemistries based on gradient-corrected DFT methods generally provide good reliability for a wide range of electronic structure calculations. The performance of models based on local DFT methods is both less accurate and more erratic. While some previous studies have claimed that LSDA is better than the gradient-corrected functionals for some properties and the reverse is true for others, we find that DFT methods using hybrid functionals are in general as good as or better than either type of pure functional. Hybrid functionals perform very well in predicting the energies, geometries, vibrational frequencies of chemical systems, as well as many molecular properties including atomic populations, dipole moments, spin densities, VCD spectra and NMR chemical shifts. They also predict thermochemical properties of systems with very good accuracy when combined with an accurate geometry and a large basis set. In these cases, such model chemistries perform as well as or better than traditional, more expensive electron correlation methods such as MP2. In the case of NMR chemical shifts, DFT methods represent an improvement over HartreeFock theory, but they are not as accurate MP2 for this property. The various hybrid functionals do not differ substantially in accuracy. In particular, B3LYP and B3PW91 appear to be comparable in accuracy for the calculation types considered here. The same parameters, chosen on the basis of thermochemical calculations at fixed geometries for one functional can be used with other functionals and give improved accuracy over pure functionals for a wide range of calculations. This reflects the fact that inclusion of exact exchange is based on sound underlying physics [31, 32] and is not simply an ad-hoc fitting procedure.

Acknowledgments The authors thank ,ZEleen Frisch for her help in preparing this chapter. All results reported for the first time in this work were computed with the Gaussian 94 program [71].

705 References

P. Hohenberg and W. Kohn, Phys. Rev., 136, B864 (1964). W. Kohn and L. J. Sham, Phys. Rev., 140, Al133 (1965). R.D. Amos, C. W. Murray and N. C. Handy, Chem. Phys. Lett., 202, 489 (1993). J. Andzelm and E. Wimmer, J. Chem. Phys., 96, 1280 (1992). J. Andzelm, C. Sosa and R. A. Eades, J. Phys. Chem., 97, 4664 (1993). J.E. Carpenter and C. P. Sosa, J. Mol. Struct., (1993). J.R. Cheeseman, M. J. Frisch, G. W. Trucks and T. A. Keith, J. Chem. Phys., 104, 5497 (1996) 8. S.J. Collins and P. J. O'Malley, Chem. Phys. Lett., 228, 246 (1994). 9. S.M. Colwell, C. W. Murray, N. C. Handy and R. D. Amos, Chem. Phys. Lett., 210, 261 (1993). 10. F.J. Devlin, P. J. Stephens, J. R. Cheeseman and M. J. Frisch, J. Am. Chem. Soc., in press (1996). 11. N. C. Handy, P. E. Maslen, R. D. Amos, J. S. Andrews, C. W. Murray and G. J. Laming, Chem. Phys. Lett., 197, 506 (1992). 12. K. Laasonen, M. Parrinello, R. Car, C. Lee and D. Vanderbilt, Chem. Phys. Lett., 207, 208 (1993 ). 13. R. Merkle, A. Savin and H. Preuss, Chem. Phys. Lett., 194, 32 (1992). 14. C. W. Murray, G. J. Laming, N. C. Handy and R. D. Amos, Chem. Phys. Lett., 199, 551 (1992). 15. T. V. Russo, R. L. Martin and P. J. Hay, J. Chem. Phys., (1993). 16. B. G. Johnson, P. M. W. Gill and J. A. Pople, J. Chem. Phys., 98, 5612 (1993). 17. A. P. Scott and L. Radom, J. Phys. Chem., submitted (1996). 18. J. M. Seminario, Chem. Phys. Lett., 206, 547 (1993). 19. C. Sosa, J. Andzelm, B. C. Elkin, E. Wimmer, K. D. Dobbs and D. A. Dixon, J. Phys. Chem., 96, 6630 (1992). 20. C. Sosa, C. Lee, G. Fitzgerald and R. A. Eades, Chem. Phys. Lett., 211, 265 (1993). 21. P.J. Stephens, F. J. Devlin, C. F. Chabalowski and M. J. Frisch, J. Phys. Chem., 98, 11623 (1994). 22. D. J. Tozer and C. P. Sosa, J. Chem. Phys., submitted (1995). 23. K. B. Wiberg, J. R. Cheeseman, J. W. Ochterski and M. J. Frisch, J. Am. Chem. S0c., 117, 6535 (1995). 24. K. B. Wiberg and J. W. Ochterski, J. Comp. Chem., accepted (1996). 25. M.W. Wong, Chem. Phys. Lett., in press (1996). 26. J. C. Slater, Quantum Theory of Molecular and Solids. Vol. 4: The Self-Consistent Field for Molecular and Solids, McGraw-Hill: New York, 1974. 27. S. H. Vosko, L. Wilk, and M. Nusair, Can. J. Phys., 58, 1200 (1980).

1. 2. 3. 4. 5. 6. 7.

706 28. 29. 30. 31. 32. 33. 34. 35. 36. 37. 38. 39. 40. 41. 42. 43. 44. 45. 46. 47. 48. 49.

50. 51. 52. 53. 54. 55. 56.

A. D. Becke, Phys. Rev. A, 38, 3098 (1988). C. Lee, W. Yang and R. G. Parr, Physical Review B, 37, 785 (1988). B. Miehlich, A. Savin, H. Stoll and H. Preuss, Chem. Phys. Lett., 157, 200 (1989). A. D. Becke, J. Chem. Phys., 98, 5648 (1993). A. D. Becke, J. Chem. Phys., 98, 1372 (1993). A. Becke, J. Chem. Phys., 104, 1040 (1996). J. P. Perdew and Y. Wang, Phys. Rev. B, 45, 13244 (1992). J.A. Pople, M. Head-Gordon, D. J. Fox, K. Raghavachari and L. A. Curtiss, J. Chem. Phys., 90, 5622 (1989). G. E. Scuseria, ]. Chem. Phys., 94, 442 (1991). L. A. Curtiss, C. Jones, G. W. Trucks, K. Raghavachari and J. A. Pople, J. Chem. Phys., 93, 2537 (1990). L. A. Curtiss, K. Raghavachari, G. W. Trucks and J. A. Pople, J. Chem. Phys., 94, 7221 (1991). J. B. Foresman, tE. Frisch, J. W. Ochterski and M. J. Frisch, in preparation (1996). C. W. Bauschlicher Jr. and H. Partridge, Chem. Phys. Lett., 240, 533 (1995). C.W. Bauschlicher Jr. and H. Partridge, J. Chem. Phys., 103, 1788 (1995). A. M. Mebel, K. Morokuma, and M.C. Lin, J. Chem. Phys., 103, 7414 (1995). L. A. Curtiss, K. Raghavachari and J. A. Pople, J. Chem. Phys., 98, 1293 (1993). G. A. Petersson, T. G. Tensfeldt and J. A. Montgomery Jr., J. Chem. Phys., 94, 6091 (1991). J. W. Ochterski, G.A. Petersson, and J.A. Montgomery, J. Chem. Phys., 104, 2598 (1996). F. De Proft, J. M. L. Martin and P. Geerlings, Chem. Phys. Lett., 250, 393 (1996). J. Cioslowski, J. Am. Chem. Soc., 111, 8333 (1989). J. A. Pople, R. Krishnan, H. B. Schlegel and J. S. Binkley, Int. J. Quant. Chem., Quant. Chem. Symp., 13, 325 (1979). J. A. Pople, IL Krishnan, H. B. Schlegel, D. DeFrees, J. S. Binkley, M. J. Frisch, R. F. Whiteside, R. F. Hour and W. J. Hehre, Int. J. Quant. Chem., Quant. Chem. Symp., 15, 269 (1981). J. A. Pople, A. P. Scott, M. W. Wong and L. Radom, Isr. J. Chem., 33, 345 (1993). P.J. Stephens and M. A. Lowe, Ann. Rev. Phys. Chem., 36, 213 (1985). P.J. Stephens, J. Phys. Chem., 89, 748 (1985). P.J. Stephens, J. Phys. Chem., 91, 1712 (1987). J. IL Cheeseman, C. S. Ashvar, M. J. Frisch, F. J. Devlin and P. J. Stephens, Chem. Phys. Lett., 252, 211 (1996). P. J. Stephens, C. S. Ashvar, F. J. Devlin, J. R. Cheeseman and M. J. Frisch, Mol. Phys., 88, in press (1996). T. Eggimann, R. A. Shaw and H. Weiser, J. Phys. Chem., 95, 591 (1991).

707 57. 58. 59. 60. 61. 62. 63. 64. 65. 66. 67. 68. 69. 70. 71.

K. Wolinski, J. F. Hilton and P. Pulay, J. Am. Chem. Soc., 112, 8251 (1990). J. Gauss, J. Chem. Phys., 99, 3629 (1993). G. Schreckenbach and T. Ziegler, J. Phys. Chem., 99, 606 (1995). A. M. Lee, N. C. Handy and S. M. Colwell, J. Chem. Phys., 103, 10095 (1995). A. D. Becke, J. Chem. Phys., 88, 2547 (1988). V. I. Lebedev, Zh. Vychisl. Mat. Mat. Fiz., 15, 48 (1975). V. I. Lebedev, Zh. Vychisl. Mat. Mat. Fiz., 16, 293 (1976). V. I. Lebedev, in Proc. Conf. Diff. Eqn. Numer. Math., Novosibivsk, 1978, Ed. S. L. Sobolev, Nauka: Novosibivsk, 1980, 110. V. I. Lebedev, Russian Acad. Sci. Dokl. Math., 45, 587 (1992). N.C. Handy and S. F. Boys, Theo. Chim. Acta, 31, 195 (1973). C. W. Murray, N. C. Handy and G. J. Laming, Mol. Phys., 78, 997 (1993). G. W. Trucks and M. J. Frisch, in preparation (1996). P. M. W. Gill, B. G. Johnson and J. A. Pople, Chem. Phys. Lett., 209, 506 (1993). C. W. Murray, N. C. Handy and G. J. Laming, Molecular Physics, 78, 997 (1993). M. J. Frisch, G. W. Trucks, H. B. Schlegel, P. M. W. Gill, B. G. Johnson, M. A. Robb, J. R. Cheeseman, T. A. Keith, G. A. Petersson, J. A. Montgomery, K. Raghavachari, M. A. A1-Laham, V. G. Zakrzewski, J. V. Ortiz, J. B. Foresman, J. Cioslowski, B. B. Stefanov, A. Nanayakkara, M. ChaUacombe, C. Y. Peng, P. Y. Ayala, W. Chen, M. W. Wong, J. L. Andres, E. S. Replogle, R. Gomperts, R. L. Martin, D. J. Fox, J. S. Binkley, D. J. Defrees, J. Baker, J. P. Stewart, M. HeadGordon, C. Gonzalez and J. A. Pople, Gaussian 94, Revision D.3, Gaussian, Inc.: Pittsburgh, PA, 1995.