J. Phys. Chem. 1996, 100, 6317-6324
6317
Electron Densities of Several Small Molecules As Calculated from Density Functional Theory Jian Wang,† Benny G. Johnson,‡ Russell J. Boyd,*,† and Leif A. Eriksson§ Department of Chemistry, Dalhousie UniVersity, Halifax, NoVa Scotia, Canada B3H 4J3, Q-Chem Inc., 317 Whipple Street, Pittsburgh, PennsylVania 15218, and Department of Physics, UniVersity of Stockholm, Box 6730, S-113 85 Stockholm, Sweden ReceiVed: October 4, 1995; In Final Form: December 20, 1995X
The electron densities of CH3F, CO, CO2, H2CO, H2O, HCN, HF, and NH3 have been calculated by density functional theory (DFT) methods with various exchange and correlation functionals and are compared with QCISD (quadratic configuration interaction method including single and double substitutions) results. The DFT methods yield electron densities and Laplacians of the densities in good agreement with the QCISD results. As expected, the gradient-corrected functionals including B3LYP, B3P86, BP86, and BLYP improve the densities calculated by the local spin density approximation relative to the corresponding QCISD densities. Relative to his 1988 functional, Becke’s three-parameter exchange functional clearly improves the results for equilibrium geometries, vibrational frequencies, dipole moments, and electron densities. For all gradientcorrected functionals, the results are generally insensitive to the numerical accuracy of the grid. The basis set, however, has a significant effect on the electron density, which also shows a large dependence on the particular choice of functionals used in the DFT calculation.
Introduction Although density functional theory (DFT)1-3 is one of the important methods for the study of the electronic structures of molecules,4-6 it is difficult to improve the approximations and corrections embedded therein. We have previously suggested that an accurate representation of the electron density can serve as the basis for the development of functionals for use in DFT.7-10 The advantage of this approach is that the electron density is a scalar quantity covering all space and can thus provide much more information than integrated quantities such as the total energy. Also, the electron density is more sensitive to the errors introduced by use of an approximate wave function than the total energy. An accurate density will in general result in good values for the energy as well as other related quantities, whereas a low energy does not have to be associated with an accurate total density. Since it is not always feasible to obtain an accurate experimental electron density distribution of a molecule, we propose to use a high-level ab initio density as an alternative. To do so, the similarities and differences between the DFT and high-level ab initio densities need to be systematically investigated. In a previous study9 on homonuclear diatomic molecules, we have shown that the electron density and its Laplacian as calculated from DFT methods are in very good agreement with QCISD results. The gradient-corrected functionals improve the local spin density approximation (LSDA) densities relative to the corresponding QCISD densities. In this paper, the study of the DFT electron densities is extended to several heteronuclear molecules: CH3F, CO, CO2, H2CO, H2O, HCN, HF, and NH3. This set of compounds is taken from the G2 data set,11 including neutral molecules with only first-row nonmetal atoms. The differences between the DFT and QCISD densities and between their Laplacians are plotted, and the effects on the electron density of the exchange and correlation functionals, numerical grid accuracy, and basis set are examined. In addition, the properties of the DFT electron densities at the bond critical points12 are compared with the QCISD results. †
Dalhousie University. Q-Chem Inc. § University of Stockholm. X Abstract published in AdVance ACS Abstracts, March 15, 1996. ‡
0022-3654/96/20100-6317$12.00/0
Computational Details The Gaussian 92/DFT program13 was used in this study. The form of the local spin density approximation used is that derived by Vosko et al.14 for the correlation together with the Slater exchange term (SVWN). Becke’s 198815 or his threeparameter16 gradient-corrected exchange functionals were used together with either Perdew’s 198617 or Lee et al.’s18 gradientcorrected correlation functionals. Hereafter, BP86 and BLYP refer to Becke’s 1988 exchange functional together with Perdew’s 1986 or Lee et al.’s correlation functionals, respectively, while B3P86 and B3LYP refer to Becke’s threeparameter hybrid exchange functional with Perdew’s 1986 or Lee et al.’s correlation functionals. The actual form of the B3 hybrid exchange functional is that implemented in the Gaussian 92/DFT program and differs slightly from the original formulation.19 The DFT electron densities were directly calculated from the SCF-converged Kohn-Sham orbitals. The fine grid (75 radial shells and 302 angular points per shell) and tight SCF convergence criteria were used in most of the calculations. All electron densities in this paper were calculated at the experimental equilibrium structures and are in atomic units (au). The generalized electron densities, which are obtained by means of the Z-vector method,20,21 were used in the MP2 and QCISD calculations. Since the electron density calculated from the DFT and ab initio methods is generally very sensitive to the basis set used in the calculations, Dunning’s correlation consistent polarized valence double- and triple-ζ basis sets augmented with the corresponding diffuse functions22 (aug-CC-PVDZ and aug-CCPVTZ) were used throughout. Electron Densities and Their Laplacians The DFT densities are generally very similar to the QCISD results for the molecules studied here. The main differences between the DFT and QCISD densities occur near the nuclei, where the DFT methods in general yield a more diffuse electron density distribution relative to the QCISD method. Among the five DFT methods (SVWN, BP86, BLYP, B3P86, and B3LYP) examined here, SVWN, as expected, gives the worst performance in comparison with the QCISD results, while the © 1996 American Chemical Society
6318 J. Phys. Chem., Vol. 100, No. 15, 1996
Wang et al.
TABLE 1: Theoretical and Experimental Dipole Moments (D)a CH3F CO H2CO H2O HCN HF NH3 a
SVWN
BLYP
BP86
B3P86
B3LYP
HF
MP2
QCISD
exptl27
1.744 0.229 2.291 1.854 3.009 1.780 1.533
1.729 0.188 2.258 1.796 2.957 1.734 1.477
1.699 0.220 2.239 1.798 2.954 1.734 1.482
1.802 0.112 2.400 1.859 3.037 1.791 1.528
1.817 0.093 2.406 1.854 3.037 1.788 1.523
2.116 0.257 2.856 2.000 3.295 1.930 1.625
1.948 0.300 2.347 1.868 3.029 1.800 1.534
1.955 0.110 2.401 1.862 3.030 1.795 1.528
1.85 0.112 2.33 1.85 2.98 1.82 1.47
Theoretical values were calculated at experimental geometries with the aug-CC-PVDZ basis set.
Figure 1. Electron density difference contours in HCN. (a) QCISDSVWN, (b) QCISD-BP86, and (c) QCISD-B3P86. Solid lines correspond to positive values (excess QCISD density), while dashed lines relate to negative values (excess DFT density). The outermost contour of the charge density difference equals 0.001 au. The succeeding contours increase in value in the order 2 × 10n, 4 × 10n, and 8 × 10n with n beginning at -3 and increasing by unity. The same set of contours is used throughout the paper.
gradient-corrected functionals improve upon the SVWN results. Becke’s three-parameter hybrid exchange functional leads to better results than his 1988 functional. Taking the linear molecule HCN as an example, the SVWN functionals underestimate the density near the nuclei of hydrogen, carbon, and nitrogen (solid lines in Figure 1a) relative to the QCISD density. In addition, SVWN overestimates the density in the outer regions of the heavy atoms (carbon and nitrogen) and in the C-N bond region. The gradient-corrected functionals yield densities in better agreement with the QCISD density than the SVWN results. This is illustrated by the QCISD-BP86 and QCISD-B3P86 density difference contours in Figure 1b and 1c. The maximum differences of the densities are about 2.0 au for QCISD-SVWN, 0.40 au for QCISD-BP86, and 0.35 au for QCISD-B3P86, respectively. Of the four gradient-corrected functionals considered, B3P86 leads to the best agreement with the QCISD density (Figure 1c), while B3LYP is the second best (not shown in Figure 1). Becke’s three-parameter hybrid exchange functional yields more accurate densities than his 1988 functional, as shown in Figure 1b and 1c. The observation that the four gradient-corrected functionals yield densities in better agreement with the QCISD densities than the SVWN results also holds for nonlinear molecules such as H2CO (Figure 2). The B3P86 density agrees best with the QCISD density (Figure 2c) in terms of the difference contours and the magnitude of the differences. B3LYP yields slightly worse results than B3P86 in terms of the magnitude of the density differences (Figure 2b).
Figure 2. Electron density difference contours in H2CO. (a) QCISDSVWN, (b) QCISD-B3LYP, and (c) QCISD-B3P86. Solid lines correspond to positive values, while dashed lines relate to negative values. The contour plane is the molecular plane.
The differences between the DFT and high-level ab initio densities are closely related to the differences of other quantities. Table 1 compares the experimental dipole moments of the chosen molecules with the theoretical values calculated at the experimental geometries with the aug-CC-PVDZ basis set. Overall, the DFT dipole moments agree well with the experimental and QCISD results, which indicates that the DFT densities are, in general, very close to the high-level ab initio results. The improvement gained by use of Becke’s threeparameter exchange functional is also demonstrated in Table 1. For hydrogen fluoride, for example, the dipole moments calculated by means of the five DFT methods agree with the QCISD result to within 0.06 D and with the experimental values to within 0.09 D. For this diatomic molecule, the more delocalized the charge distribution, the smaller the molecular dipole moment. Relative to the QCISD density, the BP86 density distribution (Figure 3b) is more delocalized over the H-F bond region and in the outer part of fluorine; consequently, it results in a smaller dipole moment than the QCISD one. The SVWN density distribution (Figure 3a), on the other hand, is less delocalized than the BP86 results (compare the solid lines in Figure 3a and 3b). The SVWN dipole moment is thus larger than the BP86 value. However, the SVWN density is more delocalized in the outer part of fluorine relative to the QCISD density, which causes its dipole moment to be smaller than the QCISD results. Similarly, the B3P86 functionals yield a more delocalized density (Figure 3c) than the QCISD density distribution, which gives rise to a relatively small dipole moment in the B3P86 calculation. However, since the density difference between the QCISD and B3P86 results is smaller than QCISD-
Electron Densities of Several Small Molecules
J. Phys. Chem., Vol. 100, No. 15, 1996 6319
Figure 3. Electron density difference contours in HF. (a) QCISDSVWN, (b) QCISD-BP86, and (c) QCISD-B3P86. Solid lines correspond to positive values (excess QCISD density), while dashed lines relate to negative values (excess DFT density).
SVWN and QCISD-BP86, the B3P86 dipole moment actually agrees better with the QCISD and experimental dipole moments than the other four DFT methods. Another way to view the density is by means of its Laplacian, which indicates where the electron density is locally concentrated or depleted.12 Negative values of the Laplacian correspond to a local concentration of electron density. Electron density is depleted in those regions where the Laplacian is positive. The Laplacian distribution recovers the valence-shell structure of a molecule and can be used to study the active site of a chemical reaction.12,23 It is surprising, however, that the Laplacian distributions calculated from the five DFT methods used in this study are nearly the same as that calculated from the QCISD results. In CH3F, for example, there are almost no differences between the SVWN and QCISD Laplacians (Figure 4a and 4b). The Laplacian of the Hartree-Fock density shown in Figure 4c is obviously different from the QCISD results in the region near fluorine. This indicates that even the SVWN method yields a better description of the electron density than the Hartree-Fock method. As in the case of CH3F, the DFT Laplacian distributions are in very good agreement with the QCISD results for CO, CO2, H2O, HCN, HF, and NH3. A minor exception is observed for H2CO. In this molecule, the SVWN functionals lead to a negative Laplacian which is nearly contiguous over the valence region of carbon and oxygen (Figure 5a), while the QCISD density shows a clear separation in the region, as shown by the solid lines in Figure 5b. The B3P86 and B3LYP functionals do not show any marked improvements over the SVWN results. Effect of the Grid Size on the Electron Density Three sets of grids were used in order to study the dependence of the electron density on the choice of the grid size. The first set has 50 radial shells and 194 angular points per shell and is referred to as a normal grid. The fine grid is defined as having 75 radial shells and 302 angular points per shell, while the extrafine grid corresponds to 99 radial shells and 302 angular points per shell. The discussion in this section is restricted to the gradient-corrected functionals only. The electron densities are in general not sensitive to the choice of the grid size. The effect on the electron density due to changing the grid is less than 0.2 au in all cases and less than
Figure 4. Contours of the Laplacian of the electron density in CH3F. (a) SVWN, (b) QCISD, and (c) Hartree-Fock. Solid lines correspond to negative values, while dashed lines relate to positive values. The contour scale is the same as in Figure 1. The densities are plotted in a CHF plane. The F atom is on the left, and the H atom is in the upper right.
Figure 5. Contours of the Laplacian of the electron density in H2CO. (a) SVWN and (b) QCISD. Solid lines correspond to negative values, while dashed lines relate to positive values. The contour scale is the same as in Figure 1. The contour plane is the molecular plane. The oxygen atom is on the right.
0.001 au in many cases. Also, the use of a larger grid does not necessarily lead to any significant change in the calculated dipole moments, equilibrium geometries, and vibrational frequencies. The explicit details of the grid effect are illustrated here by choosing HCN and CO as model systems. In HCN, the maximum change in electron density due to changing the grid is about 0.00012 au (Figure 6). The fine and extrafine grids have the same effects on the density with all five functionals. However, the absolute grid effects on the densities depend on the functional used in the calculations, as illustrated by Figure 6a and 6b. Relative to the normal grid results, the fine (or extrafine) grid decreases the B3P86 density at the nitrogen atom by 1 order of magnitude more than it does in the case of the B3LYP functionals. The same observation holds for the BLYP vs BP86 results. In CO, on the other hand,
6320 J. Phys. Chem., Vol. 100, No. 15, 1996
Wang et al.
Figure 8. B3LYP density difference along the internuclear axis in CO2 as calculated by using the extrafine and normal grids.
Figure 6. Electron density differences along the internuclear axis in HCN as calculated by using different grids. (a) Fextrafine(B3LYP) Fnormal(B3LYP) and (b) Fextrafine(B3P86) - Fnormal(B3P86).
Figure 9. Electron density differences along the internuclear axis in HF as calculated by using the aug-CC-PVTZ and aug-CC-PVDZ basis sets. (a) FPVTZ(B3P86) - FPVDZ(B3P86) and (b) FPVTZ(SVWN) - FPVDZ(SVWN).
BLYP, BP86, and B3P86, the greatest change in the density due to the variation of the grid size is less than 0.000 12 au.
Figure 7. Electron density differences along the internuclear axis in CO as calculated by using different grids. (a) Fextrafine(B3LYP) - Fnormal(B3LYP); (b) Ffine(B3P86) - Fnormal(B3P86) and (c) Fextrafine(B3P86) Fnormal(B3P86).
the fine and extrafine grids do not have the same effect on the B3P86 density (Figure 7b and 7c), although they still show the same effect on the density when the B3LYP functionals are used (Figure 7a). It is interesting to note that in the B3LYP calculations for the CO2 molecule, the fine and extrafine grids lead to an increase of the density of up to 0.18 au at the oxygen nuclei (Figure 8), while in the calculations using other functionals including
Effects of Basis Sets on the Electron Densities The basis set has been previously shown to have a considerable effect on the DFT electron densities of homonuclear diatomic molecules.9 In this section, electron densities obtained with the aug-CC-PVDZ and aug-CC-PVTZ basis sets are studied in detail for the five functionals SVWN, BP86, BLYP, B3P86, and B3LYP. The fine grid is used in all calculations. Since the density differences of linear molecules are much easier to show pictorially, HF, HCN, and CO will be discussed in detail below. Figure 9 shows the density differences calculated by using the aug-CC-PVDZ and aug-CC-PVTZ basis sets in the HF molecule. The better quality basis set (aug-CC-PVTZ) in the B3P86 calculations leads to a change in the density of up to 1.0 au at fluorine (Figure 9a). The density at hydrogen remains almost unchanged on this scale. The same observation holds for the other three gradient-corrected functionals. The results with the SVWN functionals, however, are more sensitive to the basis set than the gradient-corrected functionals in terms of the magnitude of the density difference (Figure 9b): its maximum density difference due to the basis set is about 2.0 au. In HCN, on the other hand, the gradient-corrected functionals are as sensitive to the basis set as the SVWN functionals (Figure 10). Relative to the aug-CC-PVDZ basis set, the aug-CC-PVTZ basis set leads to a decrease in the densities at the heavy atoms when the gradient-corrected functionals are used, whereas the opposite effect is observed in the SVWN calculations. Fur-
Electron Densities of Several Small Molecules
J. Phys. Chem., Vol. 100, No. 15, 1996 6321
TABLE 2: Position of Bond Critical Points, rc (Å), and the Electron Density, Gc, and Its Laplacian,a ∇2Gc (au), at the Bond Critical Points As Calculated by Various Theoretical Methods SVWN
BLYP
BP86
B3P86
B3LYP
HF
MP2
QCISD
CH3F C-F bond rc(C) Fc ∇2Fc C-H bond rc(H) Fc ∇2Fc
0.481 0.233 -0.084
0.481 0.232 -0.109
0.478 0.232 -0.092
0.469 0.230 -0.023
0.471 0.230 -0.042
0.443 0.219 0.355
0.463 0.226 0.028
0.462 0.226 0.034
0.375 0.272 -1.023
0.394 0.278 -1.029
0.393 0.276 -1.019
0.393 0.279 -1.043
0.394 0.280 -1.049
0.407 0.286 -1.106
0.396 0.278 -1.043
0.400 0.278 -1.028
rc(C) Fc ∇2Fc
0.374 0.489 1.313
0.373 0.487 1.379
0.372 0.485 1.424
CO 0.372 0.486 1.451
0.372 0.487 1.417
0.369 0.486 1.584
0.371 0.476 1.630
0.371 0.479 1.576
rc(C) Fc ∇2Fc
0.392 0.454 0.424
0.390 0.455 0.429
0.390 0.453 0.495
CO2 0.389 0.456 0.467
0.390 0.453 0.495
0.385 0.464 0.403
0.387 0.446 0.641
0.388 0.452 0.496
HCN C-H bond rc(H) Fc ∇2Fc C-N bond rc(C) Fc ∇2Fc rc(H) Fc ∇2Fc
0.331 0.278 -1.157
0.349 0.284 -1.169
0.348 0.283 -1.164
0.343 0.285 -1.195
0.345 0.286 -1.197
0.336 0.291 -1.285
0.343 0.283 -1.179
0.349 0.283 -1.169
0.396 0.476 0.075
0.394 0.476 0.162
0.393 0.474 0.205
0.392 0.475 0.259
0.392 0.477 0.225
0.387 0.475 0.511
0.392 0.469 0.267
0.391 0.469 0.338
0.187 0.357 -2.178
0.195 0.366 -2.088
0.194 0.365 -2.098
H2O 0.191 0.366 -2.161
0.192 0.366 -2.150
0.185 0.368 -2.344
0.185 0.369 -2.354
0.192 0.364 -2.103
H2CO C-O bond rc(C) Fc ∇2Fc C-H bond rc(H) Fc ∇2Fc
0.400 0.411 0.264
0.398 0.411 0.291
0.397 0.409 0.345
0.396 0.410 0.394
0.396 0.411 0.353
0.390 0.409 0.611
0.395 0.403 0.460
0.394 0.404 0.477
0.383 0.264 -0.969
0.400 0.269 -0.968
0.399 0.268 -0.961
0.399 0.271 -0.991
0.399 0.271 -0.994
0.409 0.279 -1.071
0.403 0.270 -0.986
0.406 0.270 -0.975
rc(H) Fc ∇2Fc
0.152 0.363 -3.304
0.160 0.372 -3.062
0.159 0.372 -3.090
HF 0.155 0.371 -3.244
0.156 0.371 -3.215
0.144 0.370 -3.774
0.153 0.365 -3.239
0.156 0.368 -3.150
0.251 0.336 -1.515
0.246 0.341 -1.592
0.246 0.332 -1.489
0.251 0.334 -1.485
NH3 N-H bond rc(H) Fc ∇2Fc a
0.242 0.327 -1.465
0.253 0.335 -1.496
0.251 0.334 -1.485
0.249 0.335 -1.509
Negative values of the Laplacian correspond to a local concentration of electron density; positive values correspond to a local depletion.
Figure 10. Electron density differences along the internuclear axis in HCN as calculated by using the aug-CC-PVTZ and aug-CC-PVDZ basis sets.
Figure 11. Electron density differences along the internuclear axis in CO as calculated by using the aug-CC-PVTZ and aug-CC-PVDZ basis sets.
thermore, the gradient-corrected methods lead to a greater change of the density at carbon than at nitrogen, while SVWN leads to a greater change at nitrogen than at carbon. Hence, the effect of the basis set on the density is dependent on the
choice of functionals. This is further illustrated in the CO molecule, as shown in Figure 11. Relative to the aug-CC-PVDZ basis set, the aug-CC-PVTZ basis set mainly increases the densities at the C and O nuclei in the SVWN calculations, while
6322 J. Phys. Chem., Vol. 100, No. 15, 1996
Wang et al.
TABLE 3: Theoretical and Experimental Geometriesa SVWN
BLYP
BP86
B3P86
1.383 1.109 109.0
1.422 1.105 108.2
1.411 1.107 108.4
1.392 1.098 108.6
1.137
1.146
1.146
1.133
B3LYP
HF
MP2
QCISD
exptl27
1.402 1.098 108.4
1.376 1.088 108.5
1.409 1.098 108.3
1.407 1.101 108.4
1.383 1.100 108.4
1.134
1.111
1.150
1.144
1.128
1.179
1.142
1.180
1.173
1.162
1.185 1.099 121.8
1.224 1.112 121.6
1.218 1.114 121.7
1.208 1.116 121.8
0.965 104.8
0.944 105.9
0.966 103.8
0.965 104.1
0.958 104.5
1.075 1.157
1.065 1.134
1.078 1.183
1.080 1.170
1.065 1.153
0.926
0.900
0.925
0.923
0.917
1.019 106.8
1.004 107.5
1.020 106.3
1.022 106.1
1.012 106.7
CH3F r(CF) r(CH) ∠(HCF)
CO r(CO)
CO2 r(CO)
1.170
1.181
1.179
1.165
r(CO) r(CH) ∠(HCO)
1.206 1.129 121.8
1.219 1.123 121.9
1.217 1.125 121.9
H2CO 1.205 1.207 1.114 1.114 121.8 121.8
r(OH) ∠(HOH)
0.975 104.5
0.975 104.0
0.974 103.9
0.963 104.6
r(CH) r(CN)
1.086 1.162
1.081 1.168
1.083 1.169
1.075 1.156
r(HF)
0.935
0.936
0.935
0.923
1.028 106.6
1.028 106.0
1.028 105.0
1.017 106.5
H2O
HCN
HF NH3 r(NH) ∠(HNH) a
Bond distances in angstroms, bond and dihedral angles in degrees. Theoretical values were obtained with the aug-CC-PVDZ basis set.
in the B3P86 calculations, the aug-CC-PVTZ set decreases the density at C and increases the density at O. The electron density is more sensitive to the choice of the basis set than to the grid, and therefore, the quality of the basis set has a much more pronounced effect than the grid on dipole moments, equilibrium geometries, vibrational frequencies, and other quantities. Since DFT methods require much less CPU time and disk space than conventional post-Hartree-Fock methods, choosing a better quality basis set in the DFT calculation is thus feasible and should be preferable for systems in which electron correlation has to be taken into account with great care. Note that a careful selection of a basis set is also very important in conventional ab initio methods.24 The differences in unpaired spin densities (FR-β(r)) at various nuclei in radicals have recently been computed using these same basis sets.25 It was concluded that the PVDZ basis set is quite inadequate in describing radical hyperfine properties, whereas the aug-CC-PVTZ basis set represents a significant improvement, especially when decontracting the s shell. Electron Density at the Bond Critical Points In the atoms-in-molecules theory,12 the bond critical point (BCP) is defined as a point where ∇F ) 0 and the associated Hessian matrix of the density has two negative eigenvalues and one positive eigenvalue. The properties of the electron density at a BCP, such as its position and the values of the density and its Laplacian at the BCP, can be used as a measure of interatomic interaction in a molecule.12 No attempt has been made to discuss these interactions in detail in this work. Instead, we have chosen to compare the results calculated from various ab initio (HF, MP2, and QCISD) and DFT (SVWN, BP86, BLYP, B3P86, and B3LYP) methods (Table 2). Relative to the QCISD results for the C-N, C-O, and C-F bonds in HCN, CO, CO2, and CH3F, all currently employed DFT functionals, shift the bond critical points closer to carbon. The greatest shifts of the BCP in the DFT calculations are 0.005 Å for the C-N bond, 0.003 Å for the C-O bonds, and 0.019 Å for the C-F bond, respectively. The DFT densities at the BCP agree with the corresponding QCISD values to within 0.007 au for the C-N bond, 0.010 for the C-O bonds, and
0.007 au for the C-F bond, respectively. The consistently larger density at the BCP in the DFT calculations than at the QCISD level indicates that the DFT methods lead to a greater density than the QCISD method in the bonding regions. In other words, the DFT densities are more delocalized than the QCISD density. This can be further illustrated by the Laplacian of the density at the BCP. For example, all five DFT functionals correspond to a negative ∇2F for the C-F bond in the CH3F molecule, which indicates that electrons are locally concentrated at the BCP. The QCISD density (and the HF and MP2 densities) has a positive ∇2F, which corresponds to a depletion of electrons at this point. Of the five DFT functionals, SVWN leads to the worst agreement with the QCISD results for the properties at the BCP. The better overall performance of the B3LYP and B3P86 calculations is obvious from Table 2. It is interesting to note that the properties of the electron density at the BCP obtained at the MP2 level are not always comparable to the results calculated using the gradient-corrected DFT methods. Similarly, as for the bonds between heavy atoms, for C-H, N-H, and O-H bonds, the SVWN method predicts the BCP to be closer to the hydrogen than the QCISD and gradientcorrected DFT methods. Also, the SVWN density at the BCP deviates more from the QCISD results than the gradientcorrected DFT methods. The improvement of the gradientcorrected methods, and in particular the B3LYP and B3P86 functionals, over the SVWN functionals is clearly shown in Table 2. In CH3F, for example, SVWN predicts the BCP of the C-H bond to be closer to the hydrogen by 0.025 Å than that in QCISD, while the B3LYP and B3P86 functionals predict the BCP to be only 0.006 and 0.007 Å closer to the hydrogen, respectively. The densities at the BCP of the C-H bond in CH3F are 0.272, 0.280, 0.279, and 0.278 au for SVWN, B3LYP, B3P86, and QCISD, respectively. Equilibrium Geometries and Vibrational Frequencies Tables 3 and 4 list the equilibrium geometries and vibrational frequencies calculated by various DFT and ab initio methods with the aug-CC-PVDZ basis set. In general, the HF method yields shorter internuclear distances and higher vibrational
Electron Densities of Several Small Molecules
J. Phys. Chem., Vol. 100, No. 15, 1996 6323
TABLE 4: Theoretical and Experimental Vibrational Frequencies (cm-1)a exptl HF
MP2
QCISD
obsd28,29
CH3F 1025 1446 3039 1169 1458 3132
1132 1582 3199 1280 1587 3293
1036 1470 3091 1186 1490 3204
1040 1471 3060 1186 1484 3159
1049 1464 2964 1182 1467 3006
2201
CO 2187
2403
2072
2130
2143
636 1302 2332
674 1370 2425
CO2 636 1302 2332
772 1497 2538
656 1305 2379
675 1339 2353
667 1333 2349
1459 1719 2796 1130 1197 2831
1466 1742 2792 1146 1204 2849
1512 1823 2908 1194 1246 2978
H2CO 1514 1801 2892 1192 1246 2960
1630 1979 3106 1326 1348 3185
1527 1727 2978 1188 1252 3062
1527 1769 2953 1187 1255 3028
1500 1746 2783 1167 1249 2843
1563 1764 2944 1191 1287 3009
1538 3704 3815
1580 3647 3748
1582 3677 3789
1620 3834 3946
H2O 1617 3799 3909
1744 4130 4238
1622 3804 3938
1647 3816 3932
1595 3657 3756
1648 3832 3943
Π Σ Σ
714 2147 3363
666 2085 3285
693 2103 3363
735 2200 3461
HCN 727 2189 3450
852 2401 3622
702 1991 3454
716 2125 3438
712 2097 3311
727 2129 3442
Σ
3984
3926
3946
4113
HF 4070
4467
4082
4103
3962
4139
1026 3477 1638 3611
NH3 1018 3455 1638 3583
1103 3685 1767 3818
1045 3481 1649 3635
1071 3461 1661 3596
950 3337 1627 3444
1022 3506 1691 3577
SVWN
BLYP
BP86
B3P86
A1 A1 A1 E E E
1065 1376 2949 1125 1386 3047
949 1378 2955 1114 1404 3006
997 1396 2956 1130 1409 3053
1058 1446 3048 1173 1457 3149
Σ
2158
2089
2099
Πu Σg Σu
645 1342 2401
619 1290 2292
A1 A1 A1 B1 B2 B2
1447 1792 2778 1134 1194 2841
A1 A1 B2
A1 A1 E E a
953 3377 1557 3516
1015 3344 1602 3459
1021 3348 1595 3478
B3LYP
harm.30
2170
Theoretical values were obtained with the aug-CC-PVDZ basis set.
frequencies than MP2, which tends to overcorrect the HartreeFock results. Note that the experimental harmonic vibrational frequencies should be used for comparison whenever available, since the theoretical values are based on the harmonic approximation. Density functional methods perform fairly well relative to experimental and ab initio results. The gradient-corrected DFT calculations using Becke’s three-parameter exchange functional are a clear improvement over Becke’s 1988 exchange functional for equilibrium geometries and vibrational frequencies. In fact, the B3P86 and B3LYP functionals yield equilibrium geometries and vibrational frequencies comparable with the QCISD results. In general, the present results obtained with the aug-CC-PVDZ basis set are very similar to the SVWN, BLYP, HF, MP2, and QCISD equilibrium geometries and vibrational frequencies obtained with the smaller 6-31G(d) basis set.26 (The effect of using the larger basis set is usually less than 0.005 Å, 0.5°, and 50 cm-1 for bond lengths, bond angles, and frequencies, respectively, at all levels of theory.) Conclusions The electron densities calculated by density functional methods, employing various exchange and correlation functionals, have been studied in detail and compared with the QCISD electron densities. The correlation consistent polarized valence double- and triple-ζ basis sets augmented with the corresponding diffuse functions have been used throughout. For heteronuclear diatomic and polyatomic molecules like CH3F,
CO, CO2, H2CO, H2O, HCN, HF, and NH3, the DFT methods yield electron densities and Laplacians of the densities in very close agreement with the high level ab initio results. The differences between the DFT and QCISD densities mainly occur at the nuclei, where SVWN underestimates the densities, while the gradient-corrected functionals such as BP86, BLYP, B3P86, and B3LYP result in improvements toward the QCISD results. Of the four gradient-corrected functionals, the B3P86 and B3LYP functionals yield dipole moments, equilibrium geometries, and vibrational frequencies in the best agreement with the experimental and QCISD results. All four functionals are in general insensitive to the numerical grid accuracy, and it appears that there is no need to go beyond the normal grid in most DFT calculations. On the other hand, the basis set has a much more pronounced effect on the electron density. As a consequence of the smaller requirements for CPU time and disk space in DFT calculations than in post-Hartree-Fock methods, it is possible within the DFT formalism to choose a better quality basis set. The significant improvement in the results for equilibrium geometries, vibrational frequencies, dipole moments, and electron densities gained by using Becke’s three-parameter hybrid exchange functional (B3) requires further discussion. B3 is a semiempirical hybrid of Hartree-Fock and DFT functionals. The better performance of the B3LYP and B3P86 functionals is to be expected since B3 mixes two methods which give densities that are too diffuse (DFT) and too compact (HF). Therefore, B3 gets better answers often by cancellation of errors.
6324 J. Phys. Chem., Vol. 100, No. 15, 1996 Another way to look at this is in terms of the spurious Coulomb self-interaction introduced by the approximate DFT functionals. This is a large contributing factor as to why the Kohn-Sham densities are too diffuse. The B3 method partially remedies the problem by mixing in HF exchange, which is selfinteraction-free, and reducing the approximate exchange-correlation component. An alternative approach would be to correct for the self-interaction directly either by attempting to develop functionals which are self-interaction-free, or by incorporating a self-interaction correction into the SCF calculation directly. Acknowledgment. We thank the Natural Sciences and Engineering Research Council of Canada (R.J.B.) and the Swedish Natural Science Research Council (L.A.E.) for financial support. J. W. thanks the Killam Trust for a Predoctoral Scholarship. References and Notes (1) Hohenberg, P.; Kohn, W. Phys. ReV. B 1964, 136, 864. (2) Kohn, W.; Sham, L. J. Phys. ReV. A, 1965, 140, 1133. (3) Parr, R. G.; Yang, W. Density-Functional Theory of Atoms and Molecules; Oxford University: New York, 1989. (4) Labanowski, J. K., Andzelm, J. W., Eds. Density Functional Methods in Chemistry; Springer-Verlag: New York, 1991. (5) Ziegler, T. Chem. ReV. 1991, 91, 651. (6) Politzer, P., Seminario, J. M., Eds. Theoretical and Computational Chemistry; Elsevier Scientific: Amsterdam, 1995; Vol. 2. (7) Wang, J.; Eriksson, L. A.; Boyd, R. J.; Shi, Z.; Johnson, B. G. J. Phys. Chem. 1994, 98, 1844. (8) Wang, J.; Shi, Z.; Boyd, R. J.; Gonzalez, C. A. J. Phys. Chem. 1994, 98, 6988. (9) Wang, J.; Eriksson, L. A.; Johnson, B. G.; Boyd, R. J. J. Phys. Chem., in press. (10) Boyd, R. J.; Wang, J.; Eriksson, L. A. In Recent AdVances in Density Functional Methods; Chong, D. P., Ed.; World Scientific: Singapore, 1995; Vol. 1, pp 369-401.
Wang et al. (11) Curtiss, L. A.; Raghavachari, K.; Trucks, G. W.; Pople, J. A. J. Chem. Phys. 1991, 94, 7221. (12) Bader, R. F. W. Atoms in MoleculessA Quantum Theory; Oxford University: New York, 1990. (13) Frisch, M. J.; Trucks, G. W.; Head-Gordon, M.; Gill, P. M. W.; Wong, M. W.; Foresman, J. B.; Johnson, B. G.; Schlegel, H. B.; Robb, M. A.; Replogle, E. S.; Gomperts, R.; Andres, J. L.; Raghavachari, K.; Binkley, J. S.; Gonzalez, C.; Martin, R. L.; Fox, D. J.; DeFrees, D. J.; Baker, J.; Stewart, J. J. P.; Pople, J. A. Gaussian 92, Gaussian, Inc.: Pittsburgh, PA, 1992. (14) Vosko, S. H.; Wilk, L.; Nusair, M. Can. J. Phys. 1980, 58, 1200. (15) Becke, A. D. Phys. ReV. 1988, A23, 3098. (16) Becke, A. D. J. Chem. Phys. 1993, 98, 5648. (17) Perdew, J. P.; Wang, Y. Phys. ReV. 1986, B33, 8800. (18) Lee, C.; Yang, W.; Parr, R. G. Phys. ReV. 1988, B37, 785. (19) The VWN3 parameterization of the uniform electron gas correlation is used in the Gaussian 92/DFT program. The definition of the functional and the parameters are described in: Stephens, P. J.; Devlin, F. J.; Chabalowski, C. F.; Frisch, M. J. J. Phys. Chem. 1994, 98, 11623. (20) Handy, N. C.; Schaefer, H. F., III. J. Chem. Phys. 1984, 81, 5031. (21) Wiberg, K. B.; Hadad, C. M.; LePage, T. J.; Breneman, C. M.; Frisch, M. J. J. Phys. Chem. 1992, 96, 671. (22) Dunning, T. H., Jr. J. Chem. Phys. 1989, 90, 1007. (23) (a) Shi, Z.; Boyd, R. J. J. Phys. Chem. 1991, 95, 4698. (b) Shi, Z.; Boyd, R. J. J. Chem. Phys. 1988, 88, 4375. (24) See, for example: Cramer, C. J.; Truhlar, D. G. J. Am. Chem. Soc. 1994, 116, 3892. (25) Gauld, J. W.; Eriksson, L. A.; Radom, L. To be published. (26) Johnson, B. G.; Gill, P. M. W.; Pople, J. A. J. Chem. Phys. 1993, 98, 5612. (27) Lide, D. R. Handbook of Chemistry and Physics, 72nd ed.; CRC: Boca Raton, FL, 1991. (28) Herzberg, G. Molecular Spectra and Molecular Structure I. Spectra of Diatomic Molecules, 2nd ed.; Van Nostrand, New York, 1950. (29) Pople, J. A.; Schlegel, H. B.; Krishnan, R.; Defrees, D. J.; Binkley, J. S.; Frisch; M. J.; Whiteside, R. A.; Houk, K. N.; Hehre, W. J. Int. J. Quantum Chem. Symp. 1981, 15, 269. (30) Hehre, W. J.; Radom, L.; Schleyer, P. v. R.; Pople, J. A. Ab initio Molecular Orbital Theory; Wiley, New York, 1986.
JP952944U
