5274
J. Phys. Chem. 1996, 100, 5274-5280
Electron Densities of Homonuclear Diatomic Molecules As Calculated from Density Functional Theory Jian Wang,† Leif A. Eriksson,‡ Benny G. Johnson,§ and Russell J. Boyd*,† Department of Chemistry, Dalhousie UniVersity, Halifax, NoVa Scotia, Canada B3H 4J3, Department of Physics, UniVersity of Stockholm, Box 6730, S-113 85 Stockholm, Sweden, and Q-Chem Inc., 317 Whipple Street, Pittsburgh, PennsylVania 15218 ReceiVed: April 11, 1995; In Final Form: NoVember 29, 1995X
The electron densities of H2, Li2, N2, O2, and F2 have been calculated by density functional theory (DFT) methods with various exchange and correlation functionals such as SVWN, BLYP, B3LYP, BP86, B3P86, and B3PW91 and compared with the results from quadratic configuration interaction calculations including all single and double substitutions (QCISD) using Dunning’s correlation consistent polarized valence doubleand triple-ζ basis sets augmented with the corresponding diffuse functions (aug-CC-PVDZ and aug-CCPVTZ). The DFT methods yield electron densities and Laplacians of the densities in good agreement with the QCISD results. The gradient-corrected functionals improve upon the densities calculated by use of the local spin density approximation (SVWN) relative to the corresponding QCISD densities. All gradientcorrected functionals are generally insensitive to the grids used in the study. The basis set, however, has a significant effect on the electron density and shows a strong dependence on the choice of the DFT functionals.
Introduction (DFT)1-3
is an attractive alternative Density functional theory method for the study of the electronic structures of molecules. Gradient-corrected density functional methods have been successfully applied to many systems, including intramolecular or intermolecular hydrogen bonded systems, in which the electron correlation effects have to be taken into account with great care.4,5 However, there is no straightforward way to improve the approximations and corrections which form the basis of practical density functional theory. Many approximate functionals have been developed, directly or indirectly, by intuition. We have suggested that the electron density can, in principle, be used to optimize and develop functionals for use in density functional theory.6-8 The idea behind this approach is that the electron density is a physically observable quantity and any theoretical method should converge to the same electron density in its exact limit. If there exists a high-quality electron density, the DFT functionals can be reparametrized or even reformalized so that they yield an electron density in better agreement with the reference density. Since all ground state molecular properties in DFT are unique functions of the corresponding electron density, functionals which produce an accurate electron density are expected to yield good results for other molecular properties, including energies and related quantities. Ideally, the reference density should be taken from experimental results. In reality, however, it is very difficult to obtain a reliable electron density from a direct experimental measurement. As a second option, the electron density calculated from a high-level ab initio method can be used as a reference. Since a knowledge of the ability of existing DFT functionals to yield an accurate electron density is essential for the reoptimization or reformalization of DFT functionals, the purpose of this paper is to investigate electron densities calculated by various DFT methods. The DFT electron densities of five homonuclear diatomic molecules (H2, Li2, N2, O2, and †
Dalhousie University. University of Stockholm. § Q-Chem Inc. X Abstract published in AdVance ACS Abstracts, February 15, 1996. ‡
0022-3654/96/20100-5274$12.00/0
F2) are compared with the results from quadratic configuration interaction calculations including all single and double substitutions (QCISD). For a homonuclear diatomic molecule, there is no charge transfer between the atoms, and accordingly the choice of the basis set is not as critical as in a heteronuclear polyatomic molecule. Therefore, we use Dunning’s correlation consistent polarized valence double- and triple-ζ basis sets augmented with the corresponding diffuse functions9 (aug-CCPVDZ and aug-CC-PVTZ), and we focus on the effects of the exchange and correlation functionals and the numerical grid accuracy on the total electron density. Computational Details The density functional calculations reported in this work were performed using the Gaussian 92/DFT program,10 which employs Gaussian-type orbitals (GTO) as basis sets to solve the spin-unrestricted Kohn-Sham (KS) equations. The form of the local spin density approximation used is that derived by Vosko, Wilk, and Nusair11 together with the Slater exchange term (SVWN). Becke’s 198812 and his three-parameter13 hybrid exact (Hartree-Fock) exchange with local and gradient-corrected exchange functionals together with Perdew’s 1986,14 Perdew and Wang’s 1991,15 and Lee, Yang, and Parr’s16 gradient-corrected correlation functionals were used in the gradient-corrected DFT calculations. Hereafter, BP86 and BLYP refer to Becke’s 1988 exchange functional together with Perdew’s 1986 and Lee, Yang, and Parr’s correlation functionals, respectively, while B3P86, B3PW91, and B3LYP refer to Becke’s three-parameter exchange functional with Perdew’s 1986, Perdew and Wang’s 1991, and Lee, Yang, and Parr’s correlation functionals, respectively. The DFT electron densities were calculated from the SCF converged Kohn-Sham orbitals. The fine grid (75 radial shells and 302 angular points per shell) and tight SCF convergence criterion were used in most of the calculations. All electron densities in this paper were calculated at the experimental equilibrium structures (Table 1) and are in atomic units (au). The basis sets are Dunning’s aug-CC-PVDZ and aug-CC-PVTZ9 except for Li2, for which the 6-311G(d,p) basis set was used. The Dunning basis sets are denoted by PVDZ and PVTZ in this paper. © 1996 American Chemical Society
Electron Densities of Homonuclear Diatomic Molecules
J. Phys. Chem., Vol. 100, No. 13, 1996 5275
TABLE 1: Equilibrium Bond Lengths (Å) and Vibrational Frequencies (cm-1) Obtained Using the aug-CC-PVDZ Basis Set HF
MP2
QCISD
SVWN
BLYP
BP86
B3P86
B3LYP
B3PW91
expt
0.745 4555
0.755 4460
0.761 4349
0.779 4166
0.766 4270
0.767 4264
0.760 4374
0.761 4360
0.741 4372
0.741 4401
2.784 337
2.749 342
2.693 348
2.701 340
2.715 333
2.739 330
2.714 359
2.705 355
2.670 377
2.670 351
1.078 2736
1.132 2157
1.115 2378
1.109 2401
1.117 2329
1.116 2342
1.103 2463
1.104 2447
1.098 2457
1.098 2360
1.163 1974
1.235 1417
1.213 1620
1.207 1633
1.234 1498
1.225 1546
1.201 1686
1.209 1643
1.207 1685
1.207 1580
1.338 1216
1.427 934
1.427 919
1.388 1045
1.438 936
1.423 966
1.390 1050
1.403 1020
1.392 1054
1.417 891
H2 re υ Li2a re υ N2 re υ O2 re υ F2 re υ a
Using the 6-311G(d,p) basis set.
Figure 1. Electron density differences in H2 along the internuclear axis.
Figure 3. QCISD-B3PW91 electron density differences in Li2: (a) electron density difference contours; (b) electron density differences along the internuclear axis.
Figure 2. QCISD-SVWN electron density differences in Li2: (a) electron density difference contours; (b) electron density differences along the internuclear axis. Solid lines in the contours correspond to positive values, while dashed lines relate to negative values. All densities are plotted in a plane that includes the nuclei. The outer contour of the electron 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.
In the MP2 and QCISD calculations, the generalized electron densities, which are obtained by means of the Z-vector method, were used. Previous studies17,18 have demonstrated that the Z-vector method represents the wave function response to a perturbation and is thus more closely tied to experimental observations. All occupied molecular orbitals were included in the QCISD calculations. The equilibrium bond lengths and corresponding vibrational frequencies are listed in Table 1. In general the Hartree-Fock
(HF) method yields shorter internuclear distances and higher vibrational frequencies than the second-order Møller-Plesset perturbation theory (MP2) and QCISD methods. Density functional results agree favorably with experimental and other ab initio results for these molecules. Equilibrium geometries calculated using Becke’s three-parameter exchange functional are a clear improvement over those based on Becke’s 1988 exchange functional. Some improvement in the vibrational frequencies is also observed. However, the various DFT vibrational frequencies are not as accurate as the QCISD values. This indicates that the DFT potential energy surfaces are still not quite as good as the QCISD ones. Electron Densities The QCISD-SVWN and QCISD-B3PW91 electron density differences are shown in Figures 1-9 for H2, Li2, N2, O2, and F2. The electron density differences are illustrated by means of the difference contours in a plane containing the nuclei (except for H2) and by plots of the electron density differences along the internuclear axis, the latter providing additional information on the magnitude of the differences. Relative to the QCISD density, the SVWN method in general corresponds to a more delocalized density distribution, which is characterized by less electron density near the nuclei (solid lines in the difference contours) and slightly more electron density farther from the nuclei and in the internuclear bonding region (dashed lines). The gradient-corrected B3PW91 func-
5276 J. Phys. Chem., Vol. 100, No. 13, 1996
Wang et al.
Figure 4. QCISD-SVWN electron density differences in N2: (a) electron density difference contours; (b) electron density differences along the internuclear axis.
Figure 7. QCISD-B3PW91 electron density differences in O2: (a) electron density difference contours; (b) electron density differences along the internuclear axis.
Figure 5. QCISD-B3PW91 electron density differences in N2: (a) electron density difference contours; (b) electron density differences along the internuclear axis.
Figure 8. QCISD-SVWN electron density differences in F2: (a) electron density difference contours; (b) electron density differences along the internuclear axis.
Figure 6. QCISD-SVWN electron density differences in O2: (a) electron density difference contours; (b) electron density differences along the internuclear axis.
Figure 9. QCISD-B3PW91 electron density differences in F2: (a) electron density difference contours; (b) electron density differences along the internuclear axis.
tionals improve the SVWN electron density toward the QCISD density. This is evident from the small sizes of the QCISDB3PW91 difference contours and the small magnitudes of the differences. The SVWN density is smaller than the QCISD density at all points between the nuclei in H2 (Figure 1). Also it is interesting to note that the B3PW91 functionals predict a larger density at the nuclei than the QCISD density for Li2
(Figure 3) and N2 (Figure 5). The latter results are qualitatively opposed to the results for O2 and F2 (Figures 7 and 9). That the SVWN density is more delocalized than the QCISD results in all five molecules can be attributed to the local spin density approximation. This approximation is based on the parametrization of a homogeneous electron gas model and therefore allows for a large delocalization of the electrons. The
Electron Densities of Homonuclear Diatomic Molecules
J. Phys. Chem., Vol. 100, No. 13, 1996 5277 do so to a lesser extent. In general, however, all five gradientcorrected functionals lead to qualitatively similar features in Figure 1, in which there is a maximum at the midpoint of the bond and minima near the nuclei. On the other hand, they all differ from the QCISD-SVWN curve, where there is a minimum at the midpoint and maxima at the nuclei. On the basis of the electron density, B3PW91 improves slightly the results for the H2 molecule relative to other gradient-corrected functionals. Relative to the other functionals, B3P86 and B3PW91 yield electron densities for Li2 closer to the QCISD results (Figure 10). BLYP largely overestimates the electron density at the nuclei. One interesting feature is that all five gradient-corrected functionals, including BP86, overestimate the electron density in the internuclear region with respect to the QCISD density (Figure 10). A similar observation holds for N2, as seen in Figure 11. All gradient-corrected functionals underestimate the electron density at the nuclei in O2 (Figure 12) relative to the QCISD results. Becke’s three-parameter exchange functional does not improve the results significantly over his 1988 functional. A similar observation holds for F2 (Figure 13). The Effect of Numerical Accuracy
Figure 10. QCISD-DFT electron density differences along the internuclear axis in Li2.
gradient-corrected DFT methods with the Becke (1993) and Perdew and Wang (1991) functionals (also including the BLYP, B3LYP, BP86, and B3P86 functionals, as shown below) improve upon the local approximation and result in a significant increase of the electron density at the nuclei. The Effects of Exchange and Correlation Functionals on the Electron Density As seen in the previous discussion, the electron density is somewhat sensitive to the functionals used in the DFT calculations. The effects of the exchange and correlation functionals (XC) on the electron density are investigated now in more detail with the focus on the gradient corrections mentioned above. All discussions are based on calculations using the PVTZ basis set. For the sake of convenience, the QCISD density is chosen as the reference density. Figure 1 shows the electron density difference of H2 along the internuclear axis with respect to various functionals. Relative to the QCISD results, BLYP pulls electrons toward the nuclei and the bonding region. The B3PW91, B3P86 (not shown, but nearly identical to B3PW91), and B3LYP functionals
In the Gaussian 92/DFT program, the Coulomb integrals are calculated directly, rather than using an auxiliary basis set to fit the electron density. Thus, only the exchange and correlation parts need to be fitted numerically. The accuracy of the numerical fitting will directly affect the resulting density. Three different sets of grid points were used to compare the resulting densities. The first set has 50 radial shells and 194 angular points per shell and is referred to as the normal grid. The fine grid is defined as having 75 radial shells and 302 angular points per shell, while the extrafine grid has 99 radial shells and 302 angular points per shell, which results in 30 000 grid points per atom. The discussion in this section is restricted to the gradientcorrected functionals only. Overall, the electron densities are not particularly sensitive to the grids beyond the normal one. The maximum density difference obtained with different sets of grids is less than 0.001 for the five molecules. This is consistent with the observations of Daul et al.19 In N2 (Figure 14), for example, the greatest change in the electron density due to the different grids is about 0.000 12. The extrafine grid does not improve the BLYP and B3LYP densities (Figure 14a,b) substantially. With the B3P86 functionals, the fine and extrafine grids have qualitatively different effects. Relative to the normal grid, the fine grid tends to
Figure 11. QCISD-DFT electron density differences along the internuclear axis in N2.
5278 J. Phys. Chem., Vol. 100, No. 13, 1996
Wang et al.
Figure 12. QCISD-DFT electron density differences along the internuclear axis in O2.
Figure 13. QCISD-DFT electron density differences along the internuclear axis in F2.
increase the B3P86 electron density near the nuclei, but the extrafine grid has the opposite effect (Figure 14c). The same feature was observed with the BP86 and B3PW91 functionals (not shown in Figure 14). In terms of the magnitude of the difference, the grid size has less effect on the density with B3P86, BP86, and B3PW91 than with B3LYP and BLYP. The Effect of the Basis Set on the Electron Densities Despite the fact that there is no charge transfer in a homonuclear diatomic molecule, an investigation of basis set effects can still provide some insight into the balance between various factors in the DFT calculations. All discussions in this section are based on calculations using the PVDZ and PVTZ basis sets, with the SVWN, BLYP, B3LYP, BP86, B3P86, and B3PW91 functionals. Figures 15-18 illustrate the PVTZPVDZ density differences. The fine grid was used in all calculations. Li2 is not included because Dunning basis sets are not available for Li. Figure 15 shows the PVTZ-PVDZ density differences in H2 as calculated with three different functionals. The PVTZ basis set leads to more electron density near the nuclei and in the bonding region. All six functionals, including SVWN, exhibit nearly the same basis set effects. In contrast to H2, the effect of the basis set on the N2 density depends on the choice of the functionals. B3PW91 and B3P86 have the least effect (Figure 16), while the SVWN functionals show the greatest effect. Whereas the PVTZ basis set increases the electron density in the bonding region of N2 with all functionals, it tends to decrease the electron density near the nuclei with SVWN, which is in contrast to the results from the
BP86 or B3P86 calculations. Our results suggest that the exchange functional is more sensitive to the quality of the basis set: BP86 and BLYP show a similar basis set dependence, whereas there are differences, especially near the nuclei, between the B3P86 and BP86 results. The same conclusion is obtained by comparing the B3P86 and B3LYP results with the B3LYP and BLYP results. As in the case of N2, the gradient-corrected functionals are less sensitive to the basis set than the SVWN functionals in O2 and F2 (Figures 17 and 18). Furthermore, all four gradientcorrected functionals show essentially the same basis set dependence in these two molecules. It should be pointed out that the quality of the basis set has a much more pronounced effect on the electron density than the grid accuracy. In O2, for example, the greatest density change due to varying the grid size in the B3LYP calculations is less than 0.000 12 au, while it is more than 0.75 au due to the basis set. The Laplacian of the Electron Density Another way to view the electron density is by means of the Laplacian of the density.20 In the region of space where the Laplacian has negative values, the electronic charges are concentrated, while charges are depleted in those regions in space where the Laplacian is positive. The Laplacian of the electron density in F2 calculated with the PVDZ basis set is shown in Figure 19. The Laplacian of the B3LYP electron density (Figure 19a) is basically the same as the QCISD (Figure 19b) results. This further illustrates that the electron densities calculated from these two methods are very close to each other.
Electron Densities of Homonuclear Diatomic Molecules
J. Phys. Chem., Vol. 100, No. 13, 1996 5279
Figure 16. FPVTZ-FPVDZ electron density differences along the internuclear axis in N2. The basis set effect for BLYP is essentially the same as that illustrated for BP86. The basis set effect for B3PW91 and B3LYP is very similar to that shown for B3P86.
Figure 17. FPVTZ-FPVDZ electron density differences along the internuclear axis in O2. The basis set effect for all other gradient-corrected functionals is similar to that shown for B3PW91.
Figure 14. Electron density differences along the internuclear axis in N2 as calculated by use of different grids: (a) Fextrafine-Fnormal with BLYP functionals; (b) Fextrafine-Fnormal with B3LYP functionals; (c) FextrafineFnormal with B3P86 functionals (solid line) and Ffine-Fnormal with B3P86 functionals (dashed line).
Figure 18. FPVTZ-FPVDZ electron density differences along the internuclear axis in F2. The basis set effect for all other gradient-corrected functionals is similar to that shown for B3PW91.
Figure 15. FPVTZ-FPVDZ electron density differences along the internuclear axis in H2.
Similar results are observed with the SVWN, BLYP, BP86, B3P86, and B3PW91 functionals for all five molecules. As pointed out before, since there is no electron transfer between the atoms in a homonuclear diatomic molecule, the Laplacian of the electron density, which describes the concentration and depletion of the electron density, does not provide a sufficiently sensitive measure of the difference in the density calculated by various methods. Conclusions DFT electron densities calculated with various exchange and correlation functionals have been studied in detail and compared with the QCISD densities, using the correlation consistent polarized valence double- and triple-ζ basis sets augmented with
diffuse functions. For the homonuclear diatomic molecules, H2, Li2, N2, O2, and F2, the DFT methods yield electron densities and their Laplacians in very close agreement with high-level ab initio results away from the nuclei. The major differences between the DFT and QCISD densities occur near the nuclei, where SVWN overestimates the densities, while the nonlocal functionals such as BLYP, B3LYP, BP86, B3P86, and B3PW91 result in significant improvement to the QCISD results. It should be noted that while the SVWN method overestimates the electron density near the nuclei, it underestimates the electron density at the nuclei. All functionals are generally insensitive to the choice of the numerical grid size, and there appears to be no need to go beyond the normal grid in most of the calculations. On the other hand, the basis set has a much more important effect on the electron density with the local DFT functionals, showing a greater dependence on the basis set than the gradient-corrected functionals. Acknowledgment. We thank the Natural Sciences and Engineering Research Council of Canada (R.J.B.) and the
5280 J. Phys. Chem., Vol. 100, No. 13, 1996
Figure 19. Contours of the Laplacian of the electron density in F2 as calculated by use of the PVDZ basis set: (a) B3LYP; (b) QCISD.
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) Parr, R. G.; Yang, W. Density-Functional Theory of Atoms and Molecules; Oxford University Press: New York, 1989.
Wang et al. (2) Labanowski, J. K.; Andzelm, J. W. Density Functional Methods in Chemistry; Springer-Verlag: New York, 1991. (3) Ziegler, T. Chem. ReV. 1991, 91, 651. (4) Oie, T.; Topol, I. A.; Burt, S. K. J. Phys. Chem. 1994, 98, 1121. (5) Brooks, C. L., III; Case, D. A. Chem. ReV. 1993, 93, 2487. (6) Wang, J.; Eriksson, L. A; Boyd, R. J.; Shi, Z.; Johnson, B. G. J. Phys. Chem. 1994, 98, 1844. (7) Wang, J.; Shi, Z.; Boyd, R. J.; Gonzalez, C. A. J. Phys. Chem. 1994, 98, 6988. (8) Boyd, R. J.; Wang, J.; Eriksson, L. A. In Recent AdVances in Density Functional Methods; Chong, D. P., Ed.; World Scientific Publishing: Singapore, 1995; pp 369-401. (9) Dunning, T. H., Jr. J. Chem. Phys. 1989, 90, 1007. (10) 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. (11) Vosko, S. H.; Wilk, L.; Nusair, M. Can. J. Phys. 1980, 58, 1200. (12) Becke, A. D. Phys. ReV. 1988, A38, 3098. (13) (a) Becke, A. D. J. Chem. Phys. 1993, 98, 1372. (b) Becke, A. D. J. Chem. Phys. 1993, 98, 5648. (c) Stephens, P. J.; Devlin, F. J.; Chabalowski, C. F.; Frisch M. J. J. Phys. Chem. 1994, 98, 11623. (14) Perdew, J. P.; Wang, Y. Phys. ReV. 1986, B33, 8800. (15) (a) Perdew, J. P. In Electronic Structure of Solids; Ziesche, P., Eschrig, H., Eds.; Akademie Verlag: Berlin, 1991. (b) Perdew, J. P.; Chevary, J. A.; Vosko, S. H.; Jackson, K. A.; Pederson, M. A.; Singh, D. J.; Fiolhais, C. Phys. ReV. B 1992, 46, 6671. (c) Perdew, J. P.; Wang, Y. Unpublished results. (16) Lee, C.; Yang, W.; Parr, R. G. Phys. ReV. 1988, B37, 785. (17) Handy, N. C.; Schaefer, H. F., III. J. Chem. Phys. 1984, 81, 5031. (18) Wiberg, K. B.; Hadad, C. M.; LePage, T. J.; Breneman, C. M.; Frisch, M. J. J. Phys. Chem. 1992, 96, 671. (19) Daul, C. A.; Goursot, A.; Salahub, D. R. In NATO ARW Proceedings on Numerical Grid Methods and Their Applications to Schro¨ dinger’s Equation; Cerjan, C., Ed.; Kluwer Academic Publishers: Dordrecht, 1993. (20) Bader, R. F. W. Atoms in Molecules-A Quantum Theory; Oxford University Press: New York, 1990.
JP951023G
